Home My Page Projects Code Snippets Project Openings diderot
 Summary Activity Tracker Tasks SCM

# SCM Repository

[diderot] Annotation of /trunk/test/tracto.diderot
 [diderot] / trunk / test / tracto.diderot

# Annotation of /trunk/test/tracto.diderot

Revision 248 - (view) (download)

 1 : glk 56 // tracto.diderot 2 : // 3 : // simple tractography; this only traces half the path, depending on 4 : // seedSign; caller will half to piece these together into single path 5 : // uses midpoint method (RK2) for integration 6 : // 7 : 8 : input string dataFile; // name of dataset 9 : input real stepSz; // size of steps 10 : input real CLmin; // minimum value of CL anisotropy measure 11 : input int stepNumMax; // max number of steps allowed, or zero if no limit 12 : 13 : jhr 104 image(3)[3,3] img = load (dataFile); 14 : glk 56 15 : jhr 248 //field#1(3)[3,3] F = img ⊛ bspln3; 16 : field#2(3)[3,3] F = img ⊛ bspln3; 17 : glk 56 18 : actor Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1 19 : { 20 : int stepNum = 0; 21 : vec3 pos = seedPoint; 22 : glk 58 //// principleEvec : tensor[3,3] -> vec3 23 : //// * : real x vec3 -> vec3 24 : jhr 92 vec3 guide = real(seedSign) * principleEvec(F@seedPoint); 25 : glk 56 real aniso = 0.0; 26 : 27 : update 28 : { 29 : if (inside (pos,F)) { 30 : glk 58 //// @ : field x vec -> tensor[3,3] 31 : glk 56 tensor[3,3] ten = F@pos; 32 : glk 58 //// CL : tensor[3,3] -> real 33 : glk 56 aniso = CL(ten); // "CL" is a linear anisotropy measure 34 : glk 58 //// < : real x real -> bool 35 : jhr 92 if (aniso < CLmin) { 36 : glk 56 // terminate because anisotropy went too low 37 : stabilize; 38 : } 39 : // else 40 : vec3 evec = principleEvec(ten); 41 : glk 58 //// dot : vec3 x vec3 -> real 42 : //// - : vec3 -> vec3 43 : jhr 92 if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign 44 : glk 58 //// * : real x real -> real 45 : //// * : real x vec3 -> vec3 46 : //// + : vec3 x vec3 -> vec3 47 : glk 56 evec = principleEvec(F @ (pos + 0.5*stepSz*evec)); 48 : jhr 92 if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign 49 : glk 56 guide = evec; 50 : pos = pos + stepSz*evec; 51 : glk 58 //// + : int x int -> int 52 : glk 56 stepNum = 1 + stepNum; 53 : glk 58 //// > : real x real -> bool 54 : //// && : bool x bool -> bool 55 : glk 56 if (stepNumMax > 0 && stepNum > stepNumMax) { 56 : // terminate because we took too many steps 57 : stabilize; 58 : } 59 : } else { 60 : // terminate because we went out of bounds 61 : stabilize; 62 : } 63 : } 64 : 65 : 66 : }

 root@smlnj-gforge.cs.uchicago.edu ViewVC Help Powered by ViewVC 1.0.0