// tracto.diderot // // simple tractography; this only traces half the path, depending on // seedSign; caller will half to piece these together into single path // uses midpoint method (RK2) for integration // input real stepSz; // size of steps input real CLmin; // minimum value of CL anisotropy measure input int stepNumMax; // max number of steps allowed, or zero if no limit input image(3)[3,3] img; //field#1(3)[3,3] F = img ⊛ bspln3; field#2(3)[3,3] F = img ⊛ bspln3; strand Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1 { int stepNum = 0; vec3 pos = seedPoint; //// principleEvec : tensor[3,3] -> vec3 //// * : real x vec3 -> vec3 vec3 guide = real(seedSign) * principleEvec(F@seedPoint); real aniso = 0.0; update { if (inside (pos,F)) { //// @ : field x vec -> tensor[3,3] tensor[3,3] ten = F@pos; //// CL : tensor[3,3] -> real aniso = CL(ten); // "CL" is a linear anisotropy measure //// < : real x real -> bool if (aniso < CLmin) { // terminate because anisotropy went too low stabilize; } // else vec3 evec = principleEvec(ten); //// dot : vec3 x vec3 -> real //// - : vec3 -> vec3 if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign //// * : real x real -> real //// * : real x vec3 -> vec3 //// + : vec3 x vec3 -> vec3 evec = principleEvec(F @ (pos + 0.5*stepSz*evec)); if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign guide = evec; pos = pos + stepSz*evec; //// + : int x int -> int stepNum = 1 + stepNum; //// > : real x real -> bool //// && : bool x bool -> bool if (stepNumMax > 0 && stepNum > stepNumMax) { // terminate because we took too many steps stabilize; } } else { // terminate because we went out of bounds stabilize; } } }
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: { // terminate because we went out of bounds stabilize; } } }