Tue, 18 Mar 2014 01:27:15 +0000
Update PTdecode to handle output from other Ptouch drivers
1 /*
2 #
3 # File : dtmri_view.cpp
4 # ( C++ source file )
5 #
6 # Description : A viewer of Diffusion-Tensor MRI volumes (medical imaging).
7 # This file is a part of the CImg Library project.
8 # ( http://cimg.sourceforge.net )
9 #
10 # Copyright : David Tschumperle
11 # ( http://www.greyc.ensicaen.fr/~dtschump/ )
12 #
13 # License : CeCILL v2.0
14 # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
15 #
16 # This software is governed by the CeCILL license under French law and
17 # abiding by the rules of distribution of free software. You can use,
18 # modify and/ or redistribute the software under the terms of the CeCILL
19 # license as circulated by CEA, CNRS and INRIA at the following URL
20 # "http://www.cecill.info".
21 #
22 # As a counterpart to the access to the source code and rights to copy,
23 # modify and redistribute granted by the license, users are provided only
24 # with a limited warranty and the software's author, the holder of the
25 # economic rights, and the successive licensors have only limited
26 # liability.
27 #
28 # In this respect, the user's attention is drawn to the risks associated
29 # with loading, using, modifying and/or developing or reproducing the
30 # software by the user in light of its specific status of free software,
31 # that may mean that it is complicated to manipulate, and that also
32 # therefore means that it is reserved for developers and experienced
33 # professionals having in-depth computer knowledge. Users are therefore
34 # encouraged to load and test the software's suitability as regards their
35 # requirements in conditions enabling the security of their systems and/or
36 # data to be ensured and, more generally, to use and operate it in the
37 # same conditions as regards security.
38 #
39 # The fact that you are presently reading this means that you have had
40 # knowledge of the CeCILL license and that you accept its terms.
41 #
42 */
44 #include "CImg.h"
45 using namespace cimg_library;
47 // The lines below are necessary when using a non-standard compiler as visualcpp6.
48 #ifdef cimg_use_visualcpp6
49 #define std
50 #endif
51 #ifdef min
52 #undef min
53 #undef max
54 #endif
56 // Compute fractional anisotropy (FA) of a tensor
57 //-------------------------------------------
58 template<typename T> float get_FA(const T& val1, const T& val2, const T& val3) {
59 const float
60 l1 = val1>0?val1:0, l2 = val2>0?val2:0, l3 = val3>0?val3:0,
61 lm = (l1+l2+l3)/3,
62 tr2 = 2*( l1*l1 + l2*l2 + l3*l3 ),
63 ll1 = l1-lm,
64 ll2 = l2-lm,
65 ll3 = l3-lm;
66 if (tr2>0) return (float)std::sqrt( 3*(ll1*ll1 + ll2*ll2 + ll3*ll3)/tr2 );
67 return 0;
68 }
70 // Insert an ellipsoid in a CImg 3D scene
71 //----------------------------------------
72 template<typename t,typename tp,typename tf,typename tc>
73 void insert_ellipsoid(const CImg<t>& tensor,const float X,const float Y,const float Z,const float tfact,
74 const float vx, const float vy, const float vz,
75 CImgList<tp>& points, CImgList<tf>& faces, CImgList<tc>& colors,
76 const unsigned int res1 = 20, const unsigned int res2 = 20) {
78 // Compute eigen elements
79 const float l1 = tensor[0], l2 = tensor[1], l3 = tensor[2], fa = get_FA(l1,l2,l3);
81 CImg<> vec = CImg<>::matrix(tensor[3],tensor[6],tensor[9],
82 tensor[4],tensor[7],tensor[10],
83 tensor[5],tensor[8],tensor[11]);
84 const int
85 r = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[3]),255.0f),
86 g = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[4]),255.0f),
87 b = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[5]),255.0f);
89 // Define mesh points
90 const unsigned int N0 = points.size;
91 for (unsigned int v=1; v<res2; v++)
92 for (unsigned int u=0; u<res1; u++) {
93 const float
94 alpha = (float)(u*2*cimg::valuePI/res1),
95 beta = (float)(-cimg::valuePI/2 + v*cimg::valuePI/res2),
96 x = (float)(tfact*l1*std::cos(beta)*std::cos(alpha)),
97 y = (float)(tfact*l2*std::cos(beta)*std::sin(alpha)),
98 z = (float)(tfact*l3*std::sin(beta));
99 points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(x,y,z)).mul(CImg<tp>::vector(vx,vy,vz)));
100 }
101 const unsigned int N1 = points.size;
102 points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(0,0,-l3*tfact)));
103 points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(0,0,l3*tfact)));
104 points[points.size-2](0)*=vx; points[points.size-2](1)*=vy; points[points.size-2](2)*=vz;
105 points[points.size-1](0)*=vx; points[points.size-1](1)*=vy; points[points.size-1](2)*=vz;
107 // Define mesh triangles
108 for (unsigned int vv=0; vv<res2-2; vv++)
109 for (unsigned int uu=0; uu<res1; uu++) {
110 const int nv = (vv+1)%(res2-1), nu = (uu+1)%res1;
111 faces.insert(CImg<tf>::vector(N0+res1*vv+nu,N0+res1*nv+uu,N0+res1*vv+uu));
112 faces.insert(CImg<tf>::vector(N0+res1*vv+nu,N0+res1*nv+nu,N0+res1*nv+uu));
113 colors.insert(CImg<tc>::vector(r,g,b));
114 colors.insert(CImg<tc>::vector(r,g,b));
115 }
116 for (unsigned int uu=0; uu<res1; uu++) {
117 const int nu = (uu+1)%res1;
118 faces.insert(CImg<tf>::vector(N0+nu,N0+uu,N1));
119 faces.insert(CImg<tf>::vector(N0+res1*(res2-2)+nu, N1+1,N0+res1*(res2-2)+uu));
120 colors.insert(CImg<tc>::vector(r,g,b));
121 colors.insert(CImg<tc>::vector(r,g,b));
122 }
123 }
125 // Insert a fiber in a CImg 3D scene
126 //-----------------------------------
127 template<typename T,typename te,typename tp, typename tf, typename tc>
128 void insert_fiber(const CImg<T>& fiber, const CImg<te>& eigen, const CImg<tc>& palette,
129 const int xm, const int ym, const int zm,
130 const float vx, const float vy, const float vz,
131 CImgList<tp>& points, CImgList<tf>& primitives, CImgList<tc>& colors) {
132 const int N0 = points.size;
133 float x0 = fiber(0,0), y0 = fiber(0,1), z0 = fiber(0,2), fa0 = eigen.linear_atXYZ(x0,y0,z0,12);
134 points.insert(CImg<>::vector(vx*(x0-xm),vy*(y0-ym),vz*(z0-zm)));
135 for (int l=1; l<fiber.dimx(); l++) {
136 float x1 = fiber(l,0), y1 = fiber(l,1), z1 = fiber(l,2), fa1 = eigen.linear_atXYZ(x1,y1,z1,12);
137 points.insert(CImg<tp>::vector(vx*(x1-xm),vy*(y1-ym),vz*(z1-zm)));
138 primitives.insert(CImg<tf>::vector(N0+l-1,N0+l));
139 const unsigned char
140 icol = (unsigned char)(fa0*255),
141 r = palette(icol,0),
142 g = palette(icol,1),
143 b = palette(icol,2);
144 colors.insert(CImg<unsigned char>::vector(r,g,b));
145 x0=x1; y0=y1; z0=z1; fa0=fa1;
146 }
147 }
149 // Compute fiber tracking using 4th-order Runge Kutta integration
150 //-----------------------------------------------------------------
151 template<typename T>
152 CImg<> get_fibertrack(CImg<T>& eigen,
153 const int X0, const int Y0, const int Z0, const float lmax=100,
154 const float dl = 0.1f, const float FAmin=0.7f, const float cmin=0.5f) {
156 #define align_eigen(i,j,k) \
157 { T &u = eigen(i,j,k,3), &v = eigen(i,j,k,4), &w = eigen(i,j,k,5); \
158 if (u*cu+v*cv+w*cw<0) { u=-u; v=-v; w=-w; }}
160 CImgList<> resf;
162 // Forward tracking
163 float normU = 0, normpU = 0, l = 0, X = (float)X0, Y = (float)Y0, Z = (float)Z0;
164 T
165 pu = eigen(X0,Y0,Z0,3),
166 pv = eigen(X0,Y0,Z0,4),
167 pw = eigen(X0,Y0,Z0,5);
168 normpU = (float)std::sqrt(pu*pu+pv*pv+pw*pw);
169 bool stopflag = false;
171 while (!stopflag) {
172 if (X<0 || X>eigen.dimx()-1 || Y<0 || Y>eigen.dimy()-1 || Z<0 || Z>eigen.dimz()-1 ||
173 eigen((int)X,(int)Y,(int)Z,12)<FAmin || l>lmax) stopflag = true;
174 else {
175 resf.insert(CImg<>::vector(X,Y,Z));
177 const int
178 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>=eigen.dimx())?eigen.dimx()-1:cx+1,
179 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>=eigen.dimy())?eigen.dimy()-1:cy+1,
180 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>=eigen.dimz())?eigen.dimz()-1:cz+1;
181 const T cu = eigen(cx,cy,cz,3), cv = eigen(cx,cy,cz,4), cw = eigen(cx,cy,cz,5);
183 align_eigen(px,py,pz); align_eigen(cx,py,pz); align_eigen(nx,py,pz);
184 align_eigen(px,cy,pz); align_eigen(cx,cy,pz); align_eigen(nx,cy,pz);
185 align_eigen(px,ny,pz); align_eigen(cx,ny,pz); align_eigen(nx,ny,pz);
186 align_eigen(px,py,cz); align_eigen(cx,py,cz); align_eigen(nx,py,cz);
187 align_eigen(px,cy,cz); align_eigen(nx,cy,cz);
188 align_eigen(px,ny,cz); align_eigen(cx,ny,cz); align_eigen(nx,ny,cz);
189 align_eigen(px,py,nz); align_eigen(cx,py,nz); align_eigen(nx,py,nz);
190 align_eigen(px,cy,nz); align_eigen(cx,cy,nz); align_eigen(nx,cy,nz);
191 align_eigen(px,ny,nz); align_eigen(cx,ny,nz); align_eigen(nx,ny,nz);
193 const T
194 u0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,3),
195 v0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,4),
196 w0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,5),
197 u1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,3),
198 v1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,4),
199 w1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,5),
200 u2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,3),
201 v2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,4),
202 w2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,5),
203 u3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,3),
204 v3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,4),
205 w3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,5);
206 T
207 u = u0/3 + 2*u1/3 + 2*u2/3 + u3/3,
208 v = v0/3 + 2*v1/3 + 2*v2/3 + v3/3,
209 w = w0/3 + 2*w1/3 + 2*w2/3 + w3/3;
210 if (u*pu+v*pv+w*pw<0) { u=-u; v=-v; w=-w; }
211 normU = (float)std::sqrt(u*u+v*v+w*w);
212 const float scal = (u*pu+v*pv+w*pw)/(normU*normpU);
213 if (scal<cmin) stopflag=true;
215 X+=(pu=u); Y+=(pv=v); Z+=(pw=w);
216 normpU = normU;
217 l+=dl;
218 }
219 }
221 // Backward tracking
222 l = dl; X = (float)X0; Y = (float)Y0; Z = (float)Z0;
223 pu = eigen(X0,Y0,Z0,3);
224 pv = eigen(X0,Y0,Z0,4);
225 pw = eigen(X0,Y0,Z0,5);
226 normpU = (float)std::sqrt(pu*pu+pv*pv+pw*pw);
227 stopflag = false;
229 while (!stopflag) {
230 if (X<0 || X>eigen.dimx()-1 || Y<0 || Y>eigen.dimy()-1 || Z<0 || Z>eigen.dimz()-1 ||
231 eigen((int)X,(int)Y,(int)Z,12)<FAmin || l>lmax) stopflag = true;
232 else {
234 const int
235 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>=eigen.dimx())?eigen.dimx()-1:cx+1,
236 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>=eigen.dimy())?eigen.dimy()-1:cy+1,
237 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>=eigen.dimz())?eigen.dimz()-1:cz+1;
238 const T cu = eigen(cx,cy,cz,3), cv = eigen(cx,cy,cz,4), cw = eigen(cx,cy,cz,5);
240 align_eigen(px,py,pz); align_eigen(cx,py,pz); align_eigen(nx,py,pz);
241 align_eigen(px,cy,pz); align_eigen(cx,cy,pz); align_eigen(nx,cy,pz);
242 align_eigen(px,ny,pz); align_eigen(cx,ny,pz); align_eigen(nx,ny,pz);
243 align_eigen(px,py,cz); align_eigen(cx,py,cz); align_eigen(nx,py,cz);
244 align_eigen(px,cy,cz); align_eigen(nx,cy,cz);
245 align_eigen(px,ny,cz); align_eigen(cx,ny,cz); align_eigen(nx,ny,cz);
246 align_eigen(px,py,nz); align_eigen(cx,py,nz); align_eigen(nx,py,nz);
247 align_eigen(px,cy,nz); align_eigen(cx,cy,nz); align_eigen(nx,cy,nz);
248 align_eigen(px,ny,nz); align_eigen(cx,ny,nz); align_eigen(nx,ny,nz);
250 const T
251 u0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,3),
252 v0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,4),
253 w0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,5),
254 u1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,3),
255 v1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,4),
256 w1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,5),
257 u2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,3),
258 v2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,4),
259 w2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,5),
260 u3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,3),
261 v3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,4),
262 w3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,5);
263 T
264 u = u0/3 + 2*u1/3 + 2*u2/3 + u3/3,
265 v = v0/3 + 2*v1/3 + 2*v2/3 + v3/3,
266 w = w0/3 + 2*w1/3 + 2*w2/3 + w3/3;
267 if (u*pu+v*pv+w*pw<0) { u=-u; v=-v; w=-w; }
268 normU = (float)std::sqrt(u*u+v*v+w*w);
269 const float scal = (u*pu+v*pv+w*pw)/(normU*normpU);
270 if (scal<cmin) stopflag=true;
272 X-=(pu=u); Y-=(pv=v); Z-=(pw=w);
273 normpU=normU;
274 l+=dl;
276 resf.insert(CImg<>::vector(X,Y,Z),0);
277 }
278 }
280 return resf.get_append('x');
281 }
283 // Main procedure
284 //----------------
285 int main(int argc,char **argv) {
287 // Read and init data
288 //--------------------
289 cimg_usage("A viewer of Diffusion-Tensor MRI volumes.");
290 const char *file_i = cimg_option("-i",(char*)0,"Input : Filename of tensor field (volume wxhxdx6)");
291 const char* vsize = cimg_option("-vsize","1x1x1","Input : Voxel aspect");
292 const bool normalize = cimg_option("-normalize",true,"Input : Enable tensor normalization");
293 const char *file_f = cimg_option("-f",(char*)0,"Input : Input fibers\n");
294 const float dl = cimg_option("-dl",0.5f,"Fiber computation : Integration step");
295 const float famin = cimg_option("-famin",0.3f,"Fiber computation : Fractional Anisotropy threshold");
296 const float cmin = cimg_option("-cmin",0.2f,"Fiber computation : Curvature threshold");
297 const float lmin = cimg_option("-lmin",10.0f,"Fiber computation : Minimum length\n");
298 const float lmax = cimg_option("-lmax",1000.0f,"Fiber computation : Maximum length\n");
299 const float tfact = cimg_option("-tfact",1.2f,"Display : Tensor size factor");
300 const char *bgcolor = cimg_option("-bg","0,0,0","Display : Background color");
301 unsigned int bgr=0, bgg=0, bgb=0;
302 std::sscanf(bgcolor,"%u%*c%u%*c%u",&bgr,&bgg,&bgb);
304 CImg<> tensors;
305 if (file_i) {
306 std::fprintf(stderr,"\n- Loading tensors '%s'",cimg::basename(file_i));
307 tensors.load(file_i);
308 } else {
309 // Create a synthetic tensor field here
310 std::fprintf(stderr,"\n- No input files : Creating a synthetic tensor field");
311 tensors.assign(32,32,32,6);
312 const CImg<> Id = CImg<>::diagonal(0.3f,0.3f,0.3f);
313 cimg_forXYZ(tensors,x,y,z) {
314 const float
315 u = x-tensors.dimx()/2.0f,
316 v = y-tensors.dimy()/2.0f,
317 w = z-tensors.dimz()/2.0f,
318 norm = (float)std::sqrt(1e-5f+u*u+v*v+w*w),
319 nu = u/norm, nv = v/norm, nw = w/norm;
320 const CImg<>
321 dir1 = CImg<>::vector(nu,nv,nw),
322 dir2 = CImg<>::vector(-nv,nu,nw),
323 dir3 = CImg<>::vector(nw*(nv-nu),-nw*(nu+nv),nu*nu+nv*nv);
324 tensors.set_tensor_at(2.0*dir1*dir1.get_transpose() +
325 1.0*dir2*dir2.get_transpose() +
326 0.7*dir3*dir3.get_transpose(),
327 x,y,z);
328 }
329 }
330 float voxw=1,voxh=1,voxd=1;
331 std::sscanf(vsize,"%f%*c%f%*c%f",&voxw,&voxh,&voxd);
333 std::fprintf(stderr," : %ux%ux%u image, voxsize=%gx%gx%g.",
334 tensors.dimx(),tensors.dimy(),tensors.dimz(),
335 voxw,voxh,voxd);
338 CImgList<> fibers;
339 if (file_f) {
340 std::fprintf(stderr,"\n- Loading fibers '%s'.",cimg::basename(file_f));
341 fibers.load(file_f);
342 }
344 const CImg<unsigned char> fiber_palette =
345 CImg<>(2,1,1,3).fill(200,255,0,255,0,200).RGBtoHSV().resize(256,1,1,3,3).HSVtoRGB();
347 // Compute eigen elements
348 //------------------------
349 std::fprintf(stderr,"\n- Compute eigen elements.");
350 CImg<unsigned char> coloredFA(tensors.dimx(),tensors.dimy(),tensors.dimz(),3);
351 CImg<> eigen(tensors.dimx(),tensors.dimy(),tensors.dimz(),13);
352 CImg<> val,vec;
353 float eigmax = 0;
354 cimg_forXYZ(tensors,x,y,z) {
355 tensors.get_tensor_at(x,y,z).symmetric_eigen(val,vec);
356 eigen(x,y,z,0) = val[0]; eigen(x,y,z,1) = val[1]; eigen(x,y,z,2) = val[2];
357 if (val[0]<0) val[0]=0;
358 if (val[1]<0) val[1]=0;
359 if (val[2]<0) val[2]=0;
360 if (val[0]>eigmax) eigmax = val[0];
361 eigen(x,y,z,3) = vec(0,0); eigen(x,y,z,4) = vec(0,1); eigen(x,y,z,5) = vec(0,2);
362 eigen(x,y,z,6) = vec(1,0); eigen(x,y,z,7) = vec(1,1); eigen(x,y,z,8) = vec(1,2);
363 eigen(x,y,z,9) = vec(2,0); eigen(x,y,z,10) = vec(2,1); eigen(x,y,z,11) = vec(2,2);
364 const float fa = get_FA(val[0],val[1],val[2]);
365 eigen(x,y,z,12) = fa;
366 const int
367 r = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,0))),
368 g = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,1))),
369 b = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,2)));
370 coloredFA(x,y,z,0) = (unsigned char)r;
371 coloredFA(x,y,z,1) = (unsigned char)g;
372 coloredFA(x,y,z,2) = (unsigned char)b;
373 }
374 tensors.assign();
375 std::fprintf(stderr,"\n- Maximum diffusivity = %g, Maximum FA = %g",eigmax,eigen.get_shared_channel(12).max());
376 if (normalize) {
377 std::fprintf(stderr,"\n- Normalize tensors.");
378 eigen.get_shared_channels(0,2)/=eigmax;
379 }
381 // Init display and begin user interaction
382 //-----------------------------------------
383 std::fprintf(stderr,"\n- Open user window.");
384 CImgDisplay disp(256,256,"DTMRI Viewer",0);
385 CImgDisplay disp3d(800,600,"3D Local View",0,false,true);
386 unsigned int XYZ[3];
387 XYZ[0] = eigen.dimx()/2; XYZ[1] = eigen.dimy()/2; XYZ[2] = eigen.dimz()/2;
389 while (!disp.is_closed && disp.key!=cimg::keyQ && disp.key!=cimg::keyESC) {
390 const CImg<int> s = coloredFA.get_select(disp,2,XYZ);
391 if (!disp.is_closed) switch (disp.key) {
393 // Open 3D visualization window
394 //-----------------------------
395 case cimg::keyA:
396 case 0: {
397 unsigned char white[1] = { 255 };
398 disp3d.display(CImg<unsigned char>(disp3d.dimx(),disp3d.dimy(),1,1,0).draw_text(10,10,"Please wait...",white)).show();
399 int xm,ym,zm,xM,yM,zM;
400 if (!disp.key) { xm=s[0]; ym=s[1]; zm=s[2]; xM=s[3]; yM=s[4]; zM=s[5]; }
401 else { xm=ym=zm=0; xM=eigen.dimx()-1; yM=eigen.dimy()-1; zM=eigen.dimy()-1; }
402 const CImg<> img = eigen.get_crop(xm,ym,zm,xM,yM,zM);
403 CImgList<> points;
404 CImgList<unsigned int> primitives;
405 CImgList<unsigned char> colors;
407 // Add ellipsoids to the 3D scene
408 int X = img.dimx()/2, Y = img.dimy()/2, Z = img.dimz()/2;
409 { cimg_forXY(img,x,y) insert_ellipsoid(img.get_vector_at(x,y,Z),(float)x,(float)y,(float)Z,tfact,voxw,voxh,voxd,points,primitives,colors,10,6); }
410 { cimg_forXZ(img,x,z) insert_ellipsoid(img.get_vector_at(x,Y,z),(float)x,(float)Y,(float)z,tfact,voxw,voxh,voxd,points,primitives,colors,10,6); }
411 { cimg_forYZ(img,y,z) insert_ellipsoid(img.get_vector_at(X,y,z),(float)X,(float)y,(float)z,tfact,voxw,voxh,voxd,points,primitives,colors,10,6); }
413 // Add computed fibers to the 3D scene
414 const CImg<> veigen = eigen.get_crop(xm,ym,zm,xM,yM,zM);
415 cimglist_for(fibers,l) {
416 const CImg<>& fiber = fibers[l];
417 if (fiber.dimx()) insert_fiber(fiber,eigen,fiber_palette,
418 xm,ym,zm,voxw,voxh,voxd,
419 points,primitives,colors);
420 }
422 // Display 3D object
423 CImg<unsigned char> visu = CImg<unsigned char>(3,disp3d.dimx(),disp3d.dimy(),1,0).fill(bgr,bgg,bgb).permute_axes("yzvx");
424 bool stopflag = false;
425 while (!disp3d.is_closed && !stopflag) {
426 visu.display_object3d(disp3d,points,primitives,colors,true,4,-1,false,800,0.05f,1.0f);
427 switch (disp3d.key) {
428 case cimg::keyM: { // Create movie
429 std::fprintf(stderr,"\n- Movie mode.\n");
430 const unsigned int N = 256;
431 CImg<> pts = points.get_append('x');
432 CImgList<> cpoints(points);
433 CImg<> x = pts.get_shared_line(0), y = pts.get_shared_line(1), z = pts.get_shared_line(2);
434 float
435 xm, xM = x.maxmin(xm),
436 ym, yM = y.maxmin(ym),
437 zm, zM = z.maxmin(zm),
438 ratio = 2.0f*cimg::min(visu.dimx(),visu.dimy())/(3.0f*cimg::max(xM-xm,yM-ym,zM-zm)),
439 dx = 0.5f*(xM+xm), dy = 0.5f*(yM+ym), dz = 0.5f*(zM+zm);
440 cimglist_for(points,l) {
441 cpoints(l,0) = (float)((points(l,0)-dx)*ratio);
442 cpoints(l,1) = (float)((points(l,1)-dy)*ratio);
443 cpoints(l,2) = (float)((points(l,2)-dz)*ratio);
444 }
446 for (unsigned int i=0; i<N; i++) {
447 std::fprintf(stderr,"\r- Frame %u/%u.",i,N);
448 const float alpha = (float)(i*2*cimg::valuePI/N);
449 const CImg<> rot = CImg<>::rotation_matrix(0,1,0,alpha)*CImg<>::rotation_matrix(1,0,0,1.30f);
450 CImgList<> rotated(cpoints);
451 cimglist_for(rotated,l) rotated[l] = rot*cpoints[l];
452 visu.fill(0).draw_object3d(visu.dimx()/2.0f,visu.dimy()/2.0f,-500.0f,rotated,primitives,colors,
453 4,false,800.0f,visu.dimx()/2.0f,visu.dimy()/2.0f,-800.0f,0.05f,1.0f).display(disp3d);
454 visu.save("frame.png",i);
455 }
456 visu.fill(0);
457 } break;
458 default: stopflag = true;
459 }
460 }
461 if (disp3d.is_fullscreen) disp3d.toggle_fullscreen().resize(800,600).close();
462 } break;
464 // Compute region statistics
465 //---------------------------
466 case cimg::keyR: {
467 std::fprintf(stderr,"\n- Statistics computation. Select region."); std::fflush(stderr);
468 const CImg<int> s = coloredFA.get_select(disp,2,XYZ);
469 int xm,ym,zm,xM,yM,zM;
470 if (!disp.key) { xm=s[0]; ym=s[1]; zm=s[2]; xM=s[3]; yM=s[4]; zM=s[5]; }
471 else { xm=ym=zm=0; xM=eigen.dimx()-1; yM=eigen.dimy()-1; zM=eigen.dimy()-1; }
472 const CImg<> img = eigen.get_crop(xm,ym,zm,xM,yM,zM);
473 std::fprintf(stderr,"\n- Mean diffusivity = %g, Mean FA = %g\n",
474 eigen.get_shared_channel(0).mean(),
475 eigen.get_shared_channel(12).mean());
476 } break;
478 // Track fiber bundle (single region)
479 //----------------------------------
480 case cimg::keyF: {
481 std::fprintf(stderr,"\n- Tracking mode (single region). Select starting region.\n"); std::fflush(stderr);
482 const CImg<int> s = coloredFA.get_select(disp,2,XYZ);
483 const unsigned int N = fibers.size;
484 for (int z=s[2]; z<=s[5]; z++)
485 for (int y=s[1]; y<=s[4]; y++)
486 for (int x=s[0]; x<=s[3]; x++) {
487 const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin);
488 if (fiber.dimx()>lmin) {
489 std::fprintf(stderr,"\rFiber %u : Starting from (%d,%d,%d)\t\t",fibers.size,x,y,z);
490 fibers.insert(fiber);
491 }
492 }
493 std::fprintf(stderr,"\n- %u fiber(s) added (total %u).",fibers.size-N,fibers.size);
494 } break;
496 // Track fiber bundle (double regions)
497 //------------------------------------
498 case cimg::keyG: {
499 std::fprintf(stderr,"\n- Tracking mode (double region). Select starting region."); std::fflush(stderr);
500 const CImg<int> s = coloredFA.get_select(disp,2,XYZ);
501 std::fprintf(stderr," Select ending region."); std::fflush(stderr);
502 const CImg<int> ns = coloredFA.get_select(disp,2,XYZ);
503 const unsigned int N = fibers.size;
505 // Track from start to end
506 for (int z=s[2]; z<=s[5]; z++)
507 for (int y=s[1]; y<=s[4]; y++)
508 for (int x=s[0]; x<=s[3]; x++) {
509 const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin);
510 if (fiber.dimx()>lmin) {
511 bool valid_fiber = false;
512 cimg_forX(fiber,k) {
513 const int fx = (int)fiber(k,0), fy = (int)fiber(k,1), fz = (int)fiber(k,2);
514 if (fx>=ns[0] && fx<=ns[3] &&
515 fy>=ns[1] && fy<=ns[4] &&
516 fz>=ns[2] && fz<=ns[5]) valid_fiber = true;
517 }
518 if (valid_fiber) fibers.insert(fiber);
519 }
520 }
522 // Track from end to start
523 { for (int z=ns[2]; z<=ns[5]; z++)
524 for (int y=ns[1]; y<=ns[4]; y++)
525 for (int x=ns[0]; x<=ns[3]; x++) {
526 const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin);
527 if (fiber.dimx()>lmin) {
528 bool valid_fiber = false;
529 cimg_forX(fiber,k) {
530 const int fx = (int)fiber(k,0), fy = (int)fiber(k,1), fz = (int)fiber(k,2);
531 if (fx>=s[0] && fx<=s[3] &&
532 fy>=s[1] && fy<=s[4] &&
533 fz>=s[2] && fz<=s[5]) valid_fiber = true;
534 }
535 if (valid_fiber) {
536 std::fprintf(stderr,"\rFiber %u : Starting from (%d,%d,%d)\t\t",fibers.size,x,y,z);
537 fibers.insert(fiber);
538 }
539 }
540 }}
542 std::fprintf(stderr," %u fiber(s) added (total %u).",fibers.size-N,fibers.size);
543 } break;
545 // Clear fiber bundle
546 //-------------------
547 case cimg::keyC: {
548 std::fprintf(stderr,"\n- Fibers removed.");
549 fibers.assign();
550 } break;
552 // Save fibers
553 //-------------
554 case cimg::keyS: {
555 fibers.save("fibers.cimg");
556 std::fprintf(stderr,"\n- Fibers saved.");
557 } break;
559 }
560 }
562 std::fprintf(stderr,"\n- Exit.\n\n\n");
563 return 0;
564 }