Merge pull request #7 from chamberm/magnet-tracking

Merge early october

src/dataset/Maximas.cpp

@@ -6,6 +6,7 @@
 #include "Maximas.h"
 
 #include "DatasetManager.h"
+#include "RTTrackingHelper.h"
 #include "../Logger.h"
 #include "../gfx/ShaderHelper.h"
 #include "../gui/MyListCtrl.h"
@@ -158,11 +159,13 @@ bool Maximas::load( nifti_image *pHeader, nifti_image *pBody )
     if( m_originalAxialOrientation == 0 )
     { 
         flipAxis( X_AXIS, true );
+        RTTrackingHelper::getInstance()->setMaximaFlip(Vector(-1,1,1));
         flipAnat( X_AXIS );     
     }
     if( m_originalSagOrientation == 0 )
     {
         flipAxis( Y_AXIS, true );
+        RTTrackingHelper::getInstance()->setMaximaFlip(Vector(1,-1,1));
         flipAnat( Y_AXIS );
     }
 
@@ -534,3 +537,13 @@ void Maximas::updatePropertiesSizer()
     m_pRadMapOnSphere->Hide();
     m_pRadMainAxis->Hide();
 }
+
+bool Maximas::save( wxXmlNode *pNode, const wxString &rootPath ) const
+{
+    assert( pNode != NULL );
+
+    pNode->SetName( wxT( "dataset" ) );
+    DatasetInfo::save( pNode, rootPath );
+
+    return true;
+}

src/dataset/Maximas.h

@@ -18,6 +18,7 @@ class Maximas : public Glyph
     //Maximas();
     Maximas( const wxString &filename );
     virtual ~Maximas();
+    bool save( wxXmlNode *pNode, const wxString &rootPath ) const;
 
     std::vector< std::vector<float> > *getMainDirData()                       { return &m_mainDirections;           };
 

src/dataset/RTTFibers.cpp

@@ -35,7 +35,7 @@ RTTFibers::RTTFibers()
     m_FAThreshold( 0.20f ),
     m_angleThreshold( 35.0f ),
     m_step( 1.0f ),
-    m_GMstep( 15 ),
+    m_GMstep( 5 ),
     m_nbSeed ( 10.0f ),
     m_nbMeshPt ( 0 ),
     m_puncture( 0.2f ),
@@ -47,11 +47,14 @@ RTTFibers::RTTFibers()
     m_isHARDI( false ),
     m_countGMstep( 0 ),
 	m_stop( false ),
+    m_and( true ),
     m_pExcludeInfo( NULL ),
+    m_pIncludeInfo( NULL ),
     m_pSeedMapInfo( NULL ),
     m_pGMInfo( NULL ),
     m_render( true ),
     m_steppedOnceInsideChildBox( false ),
+    m_steppedOnceIntoAND( false ),
     m_prune(true)
 {
     m_bufferObjectsRTT = new GLuint[2];
@@ -219,7 +222,7 @@ void RTTFibers::seed()
                             }
 						    //Track both sides
 							performHARDIRTT( seed,  1, pointsF, colorF); //First pass
-                            draw = m_render;
+                            draw = m_render && m_and;
 						    performHARDIRTT( seed, -1, pointsB, colorB); //Second pass
                         }
                         else
@@ -229,7 +232,7 @@ void RTTFibers::seed()
 						    performDTIRTT( Vector(x,y,z), -1, pointsB, colorB); //Second pass
                         }
 
-						if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && !m_stop && (m_render || draw))
+						if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && (m_render || draw )&& (draw || m_and))
 						{
                             bool keepRight = false;
                             bool keepLeft = false;
@@ -273,6 +276,7 @@ void RTTFibers::seed()
                                 insertPointsForTractoDriven(pointsF, pointsB);
                             }
                             m_steppedOnceInsideChildBox = false;
+                            m_steppedOnceIntoAND = false;
 						}
 					}
 				}
@@ -320,7 +324,7 @@ void RTTFibers::seed()
                             }
 						    //Track both sides
 							performHARDIRTT( seed,  1, pointsF, colorF); //First pass
-                            draw = m_render;
+                            draw = m_render && m_and;
 						    performHARDIRTT( seed, -1, pointsB, colorB); //Second pass
                         }
                         else
@@ -330,7 +334,7 @@ void RTTFibers::seed()
 						    performDTIRTT( Vector(x,y,z), -1, pointsB, colorB); //Second pass
                         }
 
-						if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && !m_stop && (m_render || draw))
+						if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && (m_render || draw )&& (draw || m_and))
 						{
                             bool keepRight = false;
                             bool keepLeft = false;
@@ -375,6 +379,7 @@ void RTTFibers::seed()
                                 insertPointsForTractoDriven(pointsF, pointsB);
                             }
                             m_steppedOnceInsideChildBox = false;
+                            m_steppedOnceIntoAND = false;
 						}
 					}
 				}
@@ -421,7 +426,7 @@ void RTTFibers::seed()
                             }
 						    //Track both sides
 							performHARDIRTT( seed,  1, pointsF, colorF); //First pass
-                            draw = m_render;
+                            draw = m_render && m_and;
 						    performHARDIRTT( seed, -1, pointsB, colorB); //Second pass
 						}
 						else
@@ -431,7 +436,7 @@ void RTTFibers::seed()
 							performDTIRTT( Vector(x,y,z), -1, pointsB, colorB); //Second pass
 						}
 
-						if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && !m_stop && (m_render || draw))
+						if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && (m_render || draw ) && (draw || m_and))
 						{
                             bool keepRight = false;
                             bool keepLeft = false;
@@ -476,6 +481,7 @@ void RTTFibers::seed()
                                 insertPointsForTractoDriven(pointsF, pointsB);
                             }
                             m_steppedOnceInsideChildBox = false;
+                            m_steppedOnceIntoAND = false;
 						}
 					}
 				}
@@ -506,7 +512,7 @@ void RTTFibers::seed()
                 {
 					//Track both sides
 					performHARDIRTT( Vector(positions[k].x,positions[k].y,positions[k].z),  1, pointsF, colorF); //First pass
-                    draw = m_render;
+                    draw = m_render && m_and;
 					performHARDIRTT( Vector(positions[k].x,positions[k].y,positions[k].z), -1, pointsB, colorB); //Second pass
                 }
                 else
@@ -516,7 +522,7 @@ void RTTFibers::seed()
 					performDTIRTT( Vector(positions[k].x,positions[k].y,positions[k].z), -1, pointsB, colorB); //Second pass
                 }
 
-				if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && !m_stop && (m_render || draw))
+				if( (pointsF.size() + pointsB.size())/3 * getStep() > getMinFiberLength() && (pointsF.size() + pointsB.size())/3 * getStep() < getMaxFiberLength() && (m_render || draw) && (draw || m_and))
 				{
                     bool keepRight = false;
                     bool keepLeft = false;
@@ -560,6 +566,7 @@ void RTTFibers::seed()
                         insertPointsForTractoDriven(pointsF, pointsB);
                     }
                     m_steppedOnceInsideChildBox = false;
+                    m_steppedOnceIntoAND = false;
 				}
             }
         }
@@ -797,7 +804,7 @@ Vector RTTFibers::advecIntegrateHARDI( Vector vin, const std::vector<float> &sti
     Vector vMagnet(0,0,0);
     float angleMin = 360.0f;
     float angle = 0.0f;
-    float puncture = m_vinvout;
+    float g = m_vinvout;
     float wm = m_pMaskInfo->at(s_number);
     float gm = 0;
     float F = 0;
@@ -813,7 +820,7 @@ Vector RTTFibers::advecIntegrateHARDI( Vector vin, const std::vector<float> &sti
     bool isMagnetOn = RTTrackingHelper::getInstance()->isMagnetOn();
     if(isMagnetOn)
     {  
-        vMagnet = magneticField(vin, sticks, s_number, pos, vOut, F);
+        vMagnet = magneticField(vin, sticks, s_number, pos, vOut, F, g);
     }
     else
     {
@@ -853,12 +860,12 @@ Vector RTTFibers::advecIntegrateHARDI( Vector vin, const std::vector<float> &sti
     //}
 
     //Weight between in and out directions. Magnet will also be weighted by distance.
-    Vector res = (1-F)*((1.0 - puncture ) * vin + puncture * vOut) + F * vMagnet;
+    Vector res = (1-F)*((1.0 - g) * vin + g * vOut) + F * vMagnet;
 
     return res;
 }
 
-Vector RTTFibers::magneticField(Vector vin, const std::vector<float> &sticks, float s_number, Vector pos, Vector& vOut, float& F) 
+Vector RTTFibers::magneticField(Vector vin, const std::vector<float> &sticks, float s_number, Vector pos, Vector& vOut, float& F, float& G) 
 {
     Vector final = vin;
     bool alreadyAffected = false;
@@ -914,11 +921,12 @@ Vector RTTFibers::magneticField(Vector vin, const std::vector<float> &sticks, fl
                         if( angle < angleMin )
                         {
                             angleMin = angle;
-                            final = v1;
+                            vOut = v1;
                         }      
                     }
                 }
-                F = selObjs[b]->getStrength();
+                //F = selObjs[b]->getStrength();
+                //G = 0;
                 alreadyAffected = true;
             } 
             else
@@ -1234,14 +1242,14 @@ std::vector<float> RTTFibers::pickDirection(std::vector<float> initialPeaks, boo
     std::vector<float> draftedPeak;
     if(!initWithDir)
     {
-
-	    float norms[3];
+	    unsigned int nbPeaks = initialPeaks.size()/3;
+	    std::vector<float> norms;
 	    float sum = 0.0f;
 
-	    for(unsigned int i=0; i < initialPeaks.size()/3; i++)
+	    for(unsigned int i=0; i < nbPeaks; i++)
         {
             Vector v1(initialPeaks[i*3],initialPeaks[i*3+1], initialPeaks[i*3+2]);
-		    norms[i] = v1.getLength();
+            norms.push_back(v1.getLength());
 		    sum += norms[i];
 	    }
 
@@ -1319,7 +1327,8 @@ std::vector<float> RTTFibers::pickDirection(std::vector<float> initialPeaks, boo
 bool RTTFibers::checkExclude( unsigned int sticksNumber)
 {
 	bool res = true;
-	if(m_pExcludeInfo != NULL)
+
+    if(m_pExcludeInfo != NULL && RTTrackingHelper::getInstance()->isNotMapOn())
 	{
 		if(m_pExcludeInfo->at(sticksNumber) != 0)
 		{
@@ -1328,6 +1337,17 @@ bool RTTFibers::checkExclude( unsigned int sticksNumber)
 		}
 	}
 
+    if(m_pIncludeInfo != NULL && RTTrackingHelper::getInstance()->isAndMapOn())
+    {
+        if(!m_steppedOnceIntoAND)
+            m_and = false;
+        if(m_pIncludeInfo->at(sticksNumber) != 0)
+        {
+            m_and = true;
+            m_steppedOnceIntoAND = true;
+        }
+    }
+
 	return res;	
 }
 

src/dataset/RTTFibers.h

@@ -35,7 +35,7 @@ class RTTFibers
     FMatrix trilinearInterp( float fx, float fy, float fz );
     Vector advecIntegrate( Vector vin, const FMatrix &tensor, Vector e1, Vector e2, Vector e3, float tensorNumber );
     Vector advecIntegrateHARDI( Vector vin, const std::vector<float> &sticks, float peaksNumber, Vector pos );
-    Vector magneticField( Vector vin, const std::vector<float> &sticks, float peaksNumber, Vector pos, Vector& vOut, float& F);
+    Vector magneticField( Vector vin, const std::vector<float> &sticks, float peaksNumber, Vector pos, Vector& vOut, float& F, float& G);
 
     void clearFibersRTT();
 
@@ -66,7 +66,9 @@ class RTTFibers
     void insertPointsForTractoDriven( std::vector<float> pointsF, std::vector<float> pointsB);
 
 	void setExcludeInfo( Anatomy* info )                              { m_pExcludeInfo = info; }
+    void setIncludeInfo( Anatomy* info )                              { m_pIncludeInfo = info; }
 	bool checkExclude(unsigned int sticksNumber);						  
+    void setAnd(bool val) { m_and = val;}
 
     float getFAThreshold()                       { return m_FAThreshold; }
     float getAngleThreshold()                    { return m_angleThreshold; }
@@ -118,6 +120,7 @@ class RTTFibers
 	DatasetInfo *m_pShellInfo;
     Anatomy     *m_pMaskInfo;
 	Anatomy     *m_pExcludeInfo;
+    Anatomy     *m_pIncludeInfo;
 	Anatomy     *m_pSeedMapInfo;
     Anatomy     *m_pGMInfo;
     Vector       m_initVec;
@@ -129,8 +132,10 @@ class RTTFibers
 	float m_alpha;
     int m_currentSeedBoxID;
 	bool m_stop;
+    bool m_and;
     bool m_render;
     bool m_steppedOnceInsideChildBox;
+    bool m_steppedOnceIntoAND;
     bool m_prune;
     GLuint*     m_bufferObjectsRTT;
 

src/dataset/RTTrackingHelper.cpp

@@ -19,9 +19,12 @@ RTTrackingHelper::RTTrackingHelper()
 	m_isSeedFromfMRI ( false ),
     m_isTractoDrivenRSN( false ),
     m_isGMallowed( false ),
+    m_isNotMapOn( false ),
+    m_isAndMapOn( false ),
 	m_id( 0 ),
     m_isSrcAlpha( true ),
     m_isRandomInit( true ),
+    flippedAxes( 1, 1, 1),
     m_bufferID( 0 )
 {
 }

src/dataset/RTTrackingHelper.h

@@ -29,6 +29,8 @@ class RTTrackingHelper
 	bool isSeedFromfMRI() const { return m_isSeedFromfMRI; }
     bool isTractoDrivenRSN() const { return m_isTractoDrivenRSN; }
     bool isGMAllowed() const { return m_isGMallowed; }
+    bool isNotMapOn() const {return m_isNotMapOn;}
+    bool isAndMapOn() const {return m_isAndMapOn;}
     bool isSrcAlpha() const { return m_isSrcAlpha; }
     bool isRandomInit() const { return m_isRandomInit; }
 
@@ -38,7 +40,7 @@ class RTTrackingHelper
     void setRTTDirty( bool dirty )            { m_isRTTDirty = dirty; }
     void setRTTActive( bool active )          { m_isRTTActive = active; }
 	void setSeedFromfMRI( bool seedFromfMRI ) { m_isSeedFromfMRI = seedFromfMRI; }
-    void setMaximaFlip( Vector flip)          { flippedAxes = flip; }
+    void setMaximaFlip( Vector flip)          { flippedAxes.x = flip.x; flippedAxes.y = flip.y; flippedAxes.z = flip.z; }
     Vector getMaximaFlip()          { return flippedAxes; }
 
     bool toggleInterpolateTensors() { return m_interpolateTensors = !m_interpolateTensors; }
@@ -50,6 +52,8 @@ class RTTrackingHelper
     bool toggleMagnet()             { return m_isMagnetOn = !m_isMagnetOn; }
     bool toggleSrcAlpha()           { return m_isSrcAlpha = !m_isSrcAlpha; }
     bool toggleRandomInit()         { return m_isRandomInit = !m_isRandomInit; }
+    bool toggleNotMap()             { return m_isNotMapOn = !m_isNotMapOn;}
+    bool toggleAndMap()             { return m_isAndMapOn = !m_isAndMapOn;}
 
     bool togglePlayStop()             { return m_isPaused = !m_isPaused; }
     bool isTrackActionPlaying() const { return m_isTrackActionPlaying; }
@@ -59,6 +63,8 @@ class RTTrackingHelper
 	void setShellSeed(bool enabled) {m_isShellSeeds = enabled;}
 	void setSeedMap(bool enabled) {m_isSeedMap = enabled;}
     void setGMmap(bool enabled) {m_isGMallowed = enabled;}
+    void setNotmap(bool enabled) {m_isNotMapOn = enabled;}
+    void setAndmap(bool enabled) {m_isAndMapOn = enabled;}
 
 	int generateId() {return (++m_id);}
     bool toogleTractoDrivenRSN()          { return m_isTractoDrivenRSN = !m_isTractoDrivenRSN; }
@@ -97,6 +103,8 @@ class RTTrackingHelper
 	bool m_isSeedFromfMRI;
     bool m_isTractoDrivenRSN;
     bool m_isGMallowed;
+    bool m_isNotMapOn;
+    bool m_isAndMapOn;
 	int m_id;
     bool m_isSrcAlpha;
     bool m_isRandomInit;

src/gui/PropertiesWindow.cpp

@@ -405,9 +405,7 @@ void PropertiesWindow::OnSliderIntensityThresholdMoved( wxCommandEvent& WXUNUSED
 
             std::vector< Vector > positions = s->m_tMesh->getVerts();
             float shellSeedNb = positions.size();
-            RTTrackingHelper::getInstance()->m_pTxtTotalSeedNbBox->SetValue(wxString::Format( wxT( "%.1f"), shellSeedNb) );
-
-
+            RTTrackingHelper::getInstance()->m_pTxtTotalSeedNbBox->SetValue(wxString::Format( wxT( "%.1f"), shellSeedNb) );      
         }
 
         // This slider will set the Brightness level. Currently only the glyphs uses this value.
@@ -501,6 +499,7 @@ void PropertiesWindow::OnFitToAnat( wxCommandEvent& WXUNUSED(event) )
         if( pFibers != NULL )
         {
             pFibers->fitToAnat(false);
+            pFibers->updateLinesShown();
         }
     }
     else

src/gui/SceneManager.cpp

@@ -602,7 +602,7 @@ bool SceneManager::loadOldVersion( wxXmlNode * pRoot, const wxString &rootPath
                 }
                 else
                 {
-                    wxFileName fullDatasetPath( rootPath + wxFileName::GetPathSeparator() + path );
+                    wxFileName fullDatasetPath( path );
                     index = DatasetManager::getInstance()->load( fullDatasetPath.GetFullPath(), extension );
                 }