10 |
input real CLmin; // minimum value of CL anisotropy measure |
input real CLmin; // minimum value of CL anisotropy measure |
11 |
input int stepNumMax; // max number of steps allowed, or zero if no limit |
input int stepNumMax; // max number of steps allowed, or zero if no limit |
12 |
|
|
13 |
image(3)[] img = load (dataFile); |
image(3)[3,3] img = load (dataFile); |
14 |
|
|
15 |
field#1(3)[3,3] F = convolve (bspln3, img); |
//field#1(3)[3,3] F = img ⊛ bspln3; |
16 |
|
field#2(3)[3,3] F = img ⊛ bspln3; |
17 |
|
|
18 |
actor Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1 |
actor Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1 |
19 |
{ |
{ |
21 |
vec3 pos = seedPoint; |
vec3 pos = seedPoint; |
22 |
//// principleEvec : tensor[3,3] -> vec3 |
//// principleEvec : tensor[3,3] -> vec3 |
23 |
//// * : real x vec3 -> vec3 |
//// * : real x vec3 -> vec3 |
24 |
vec3 guide = seedSign * principleEvec(F@seedPoint); |
vec3 guide = real(seedSign) * principleEvec(F@seedPoint); |
25 |
real aniso = 0.0; |
real aniso = 0.0; |
26 |
|
|
27 |
update |
update |
32 |
//// CL : tensor[3,3] -> real |
//// CL : tensor[3,3] -> real |
33 |
aniso = CL(ten); // "CL" is a linear anisotropy measure |
aniso = CL(ten); // "CL" is a linear anisotropy measure |
34 |
//// < : real x real -> bool |
//// < : real x real -> bool |
35 |
if (aniso < clmin) { |
if (aniso < CLmin) { |
36 |
// terminate because anisotropy went too low |
// terminate because anisotropy went too low |
37 |
stabilize; |
stabilize; |
38 |
} |
} |
40 |
vec3 evec = principleEvec(ten); |
vec3 evec = principleEvec(ten); |
41 |
//// dot : vec3 x vec3 -> real |
//// dot : vec3 x vec3 -> real |
42 |
//// - : vec3 -> vec3 |
//// - : vec3 -> vec3 |
43 |
if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign |
if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign |
44 |
//// * : real x real -> real |
//// * : real x real -> real |
45 |
//// * : real x vec3 -> vec3 |
//// * : real x vec3 -> vec3 |
46 |
//// + : vec3 x vec3 -> vec3 |
//// + : vec3 x vec3 -> vec3 |
47 |
evec = principleEvec(F @ (pos + 0.5*stepSz*evec)); |
evec = principleEvec(F @ (pos + 0.5*stepSz*evec)); |
48 |
if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign |
if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign |
49 |
guide = evec; |
guide = evec; |
50 |
pos = pos + stepSz*evec; |
pos = pos + stepSz*evec; |
51 |
//// + : int x int -> int |
//// + : int x int -> int |