Annotation of /trunk/test/tracto.diderot

Parent Directory | Revision Log

---

Revision 56 - (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 :			vec3 guide = seedSign * principleEvec(F@seedPoint);
22 :			real aniso = 0.0;
23 :
24 :			update
25 :			{
26 :			if (inside (pos,F)) {
27 :			tensor[3,3] ten = F@pos;
28 :			aniso = CL(ten); // "CL" is a linear anisotropy measure
29 :			if (aniso < clmin) {
30 :			// terminate because anisotropy went too low
31 :			stabilize;
32 :			}
33 :			// else
34 :			vec3 evec = principleEvec(ten);
35 :			if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign
36 :			evec = principleEvec(F @ (pos + 0.5*stepSz*evec));
37 :			if (dot(evec, guide) < 0) { evec = -evec; } // fix eigenvector sign
38 :			guide = evec;
39 :			pos = pos + stepSz*evec;
40 :			stepNum = 1 + stepNum;
41 :			if (stepNumMax > 0 && stepNum > stepNumMax) {
42 :			// terminate because we took too many steps
43 :			stabilize;
44 :			}
45 :			} else {
46 :			// terminate because we went out of bounds
47 :			stabilize;
48 :			}
49 :			}
50 :
51 :
52 :			}

---

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