Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/pure-cfg/test/tracto.diderot
ViewVC logotype

Annotation of /branches/pure-cfg/test/tracto.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 477 - (view) (download)

1 : glk 56 // tracto.diderot
2 :     //
3 :     // simple tractography; this only traces half the path, depending on
4 :     // seedSign; caller will half to piece these together into single path
5 :     // uses midpoint method (RK2) for integration
6 :     //
7 :    
8 :     input string dataFile; // name of dataset
9 :     input real stepSz; // size of steps
10 :     input real CLmin; // minimum value of CL anisotropy measure
11 :     input int stepNumMax; // max number of steps allowed, or zero if no limit
12 :    
13 : jhr 104 image(3)[3,3] img = load (dataFile);
14 : glk 56
15 : jhr 248 //field#1(3)[3,3] F = img ⊛ bspln3;
16 :     field#2(3)[3,3] F = img ⊛ bspln3;
17 : glk 56
18 :     actor Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1
19 :     {
20 :     int stepNum = 0;
21 :     vec3 pos = seedPoint;
22 : glk 58 //// principleEvec : tensor[3,3] -> vec3
23 :     //// * : real x vec3 -> vec3
24 : jhr 92 vec3 guide = real(seedSign) * principleEvec(F@seedPoint);
25 : glk 56 real aniso = 0.0;
26 :    
27 :     update
28 :     {
29 :     if (inside (pos,F)) {
30 : glk 58 //// @ : field x vec -> tensor[3,3]
31 : glk 56 tensor[3,3] ten = F@pos;
32 : glk 58 //// CL : tensor[3,3] -> real
33 : glk 56 aniso = CL(ten); // "CL" is a linear anisotropy measure
34 : glk 58 //// < : real x real -> bool
35 : jhr 92 if (aniso < CLmin) {
36 : glk 56 // terminate because anisotropy went too low
37 :     stabilize;
38 :     }
39 :     // else
40 :     vec3 evec = principleEvec(ten);
41 : glk 58 //// dot : vec3 x vec3 -> real
42 :     //// - : vec3 -> vec3
43 : jhr 92 if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign
44 : glk 58 //// * : real x real -> real
45 :     //// * : real x vec3 -> vec3
46 :     //// + : vec3 x vec3 -> vec3
47 : glk 56 evec = principleEvec(F @ (pos + 0.5*stepSz*evec));
48 : jhr 92 if (dot(evec, guide) < 0.0) { evec = -evec; } // fix eigenvector sign
49 : glk 56 guide = evec;
50 :     pos = pos + stepSz*evec;
51 : glk 58 //// + : int x int -> int
52 : glk 56 stepNum = 1 + stepNum;
53 : glk 58 //// > : real x real -> bool
54 :     //// && : bool x bool -> bool
55 : glk 56 if (stepNumMax > 0 && stepNum > stepNumMax) {
56 :     // terminate because we took too many steps
57 :     stabilize;
58 :     }
59 :     } else {
60 :     // terminate because we went out of bounds
61 :     stabilize;
62 :     }
63 :     }
64 :    
65 :    
66 :     }

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0