// 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)[3,3] img = image (dataFile);

//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;
        }
    }

}

initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];