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 |
{ |
{ |
20 |
int stepNum = 0; |
int stepNum = 0; |
21 |
vec3 pos = seedPoint; |
vec3 pos = seedPoint; |
22 |
vec3 guide = seedSign * principleEvec(F@seedPoint); |
//// principleEvec : tensor[3,3] -> vec3 |
23 |
|
//// * : real x vec3 -> vec3 |
24 |
|
vec3 guide = real(seedSign) * principleEvec(F@seedPoint); |
25 |
real aniso = 0.0; |
real aniso = 0.0; |
26 |
|
|
27 |
update |
update |
28 |
{ |
{ |
29 |
if (inside (pos,F)) { |
if (inside (pos,F)) { |
30 |
|
//// @ : field x vec -> tensor[3,3] |
31 |
tensor[3,3] ten = F@pos; |
tensor[3,3] ten = F@pos; |
32 |
|
//// 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 |
if (aniso < clmin) { |
//// < : real x real -> bool |
35 |
|
if (aniso < CLmin) { |
36 |
// terminate because anisotropy went too low |
// terminate because anisotropy went too low |
37 |
stabilize; |
stabilize; |
38 |
} |
} |
39 |
// else |
// else |
40 |
vec3 evec = principleEvec(ten); |
vec3 evec = principleEvec(ten); |
41 |
if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign |
//// dot : vec3 x vec3 -> real |
42 |
|
//// - : vec3 -> vec3 |
43 |
|
if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign |
44 |
|
//// * : real x real -> real |
45 |
|
//// * : real x vec3 -> vec3 |
46 |
|
//// + : 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 |
52 |
stepNum = 1 + stepNum; |
stepNum = 1 + stepNum; |
53 |
|
//// > : real x real -> bool |
54 |
|
//// && : bool x bool -> bool |
55 |
if (stepNumMax > 0 && stepNum > stepNumMax) { |
if (stepNumMax > 0 && stepNum > stepNumMax) { |
56 |
// terminate because we took too many steps |
// terminate because we took too many steps |
57 |
stabilize; |
stabilize; |