Annotation of /trunk/test/tracto.diderot

Parent Directory | Revision Log

---

Revision 58 - (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 :			image(3)[] img = load (dataFile);
14 :
15 :			field#1(3)[3,3] F = convolve (bspln3, img);
16 :
17 :			actor Tracto (vec3 seedPoint, int seedSign) // seedSign == +1 or -1
18 :			{
19 :			int stepNum = 0;
20 :			vec3 pos = seedPoint;
21 :	glk	58	//// principleEvec : tensor[3,3] -> vec3
22 :			//// * : real x vec3 -> vec3
23 :	glk	56	vec3 guide = seedSign * principleEvec(F@seedPoint);
24 :			real aniso = 0.0;
25 :
26 :			update
27 :			{
28 :			if (inside (pos,F)) {
29 :	glk	58	//// @ : field x vec -> tensor[3,3]
30 :	glk	56	tensor[3,3] ten = F@pos;
31 :	glk	58	//// CL : tensor[3,3] -> real
32 :	glk	56	aniso = CL(ten); // "CL" is a linear anisotropy measure
33 :	glk	58	//// < : real x real -> bool
34 :	glk	56	if (aniso < clmin) {
35 :			// terminate because anisotropy went too low
36 :			stabilize;
37 :			}
38 :			// else
39 :			vec3 evec = principleEvec(ten);
40 :	glk	58	//// dot : vec3 x vec3 -> real
41 :			//// - : vec3 -> vec3
42 :	glk	56	if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign
43 :	glk	58	//// * : real x real -> real
44 :			//// * : real x vec3 -> vec3
45 :			//// + : vec3 x vec3 -> vec3
46 :	glk	56	evec = principleEvec(F @ (pos + 0.5*stepSz*evec));
47 :			if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign
48 :			guide = evec;
49 :			pos = pos + stepSz*evec;
50 :	glk	58	//// + : int x int -> int
51 :	glk	56	stepNum = 1 + stepNum;
52 :	glk	58	//// > : real x real -> bool
53 :			//// && : bool x bool -> bool
54 :	glk	56	if (stepNumMax > 0 && stepNum > stepNumMax) {
55 :			// terminate because we took too many steps
56 :			stabilize;
57 :			}
58 :			} else {
59 :			// terminate because we went out of bounds
60 :			stabilize;
61 :			}
62 :			}
63 :
64 :
65 :			}

---

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