Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

# Annotation of /trunk/test/tracto.diderot

 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 : image(3)[] img = load (dataFile); 14 : 15 : field#1(3)[3,3] F = convolve (bspln3, img); 16 : 17 : actor Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1 18 : { 19 : int stepNum = 0; 20 : vec3 pos = seedPoint; 21 : glk 58 //// principleEvec : tensor[3,3] -> vec3 22 : //// * : real x vec3 -> vec3 23 : glk 56 vec3 guide = seedSign * principleEvec(F@seedPoint); 24 : real aniso = 0.0; 25 : 26 : update 27 : { 28 : if (inside (pos,F)) { 29 : glk 58 //// @ : field x vec -> tensor[3,3] 30 : glk 56 tensor[3,3] ten = F@pos; 31 : glk 58 //// CL : tensor[3,3] -> real 32 : glk 56 aniso = CL(ten); // "CL" is a linear anisotropy measure 33 : glk 58 //// < : real x real -> bool 34 : glk 56 if (aniso < clmin) { 35 : // terminate because anisotropy went too low 36 : stabilize; 37 : } 38 : // else 39 : vec3 evec = principleEvec(ten); 40 : glk 58 //// dot : vec3 x vec3 -> real 41 : //// - : vec3 -> vec3 42 : glk 56 if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign 43 : glk 58 //// * : real x real -> real 44 : //// * : real x vec3 -> vec3 45 : //// + : vec3 x vec3 -> vec3 46 : glk 56 evec = principleEvec(F @ (pos + 0.5*stepSz*evec)); 47 : if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign 48 : guide = evec; 49 : pos = pos + stepSz*evec; 50 : glk 58 //// + : int x int -> int 51 : glk 56 stepNum = 1 + stepNum; 52 : glk 58 //// > : real x real -> bool 53 : //// && : bool x bool -> bool 54 : glk 56 if (stepNumMax > 0 && stepNum > stepNumMax) { 55 : // terminate because we took too many steps 56 : stabilize; 57 : } 58 : } else { 59 : // terminate because we went out of bounds 60 : stabilize; 61 : } 62 : } 63 : 64 : 65 : }