// 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 string dataFile; // name of dataset
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
image(3)[] img = load (dataFile);
field#1(3)[3,3] F = convolve (bspln3, img);
actor 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;
}
}
}