Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] Diff of /trunk/test/tracto.diderot
 [diderot] / trunk / test / tracto.diderot

# Diff of /trunk/test/tracto.diderot

revision 56, Mon Apr 26 00:57:58 2010 UTC revision 104, Tue Jun 1 16:04:19 2010 UTC
# Line 10  Line 10
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
14
15  field#1(3)[3,3] F = convolve (bspln3, img);  //field#1(3)[3,3] F = convolve (bspln3, img);
16    field#2(3)[3,3] F = convolve (bspln3, img);
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;

Legend:
 Removed from v.56 changed lines Added in v.104