1.1 diff -r 5edfbd3e7a46 -r 1204ebf9340d PTdecode/CImg-1.3.0/examples/dtmri_view.cpp 1.2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 +++ b/PTdecode/CImg-1.3.0/examples/dtmri_view.cpp Mon Aug 03 14:09:20 2009 +0100 1.4 @@ -0,0 +1,564 @@ 1.5 +/* 1.6 + # 1.7 + # File : dtmri_view.cpp 1.8 + # ( C++ source file ) 1.9 + # 1.10 + # Description : A viewer of Diffusion-Tensor MRI volumes (medical imaging). 1.11 + # This file is a part of the CImg Library project. 1.12 + # ( http://cimg.sourceforge.net ) 1.13 + # 1.14 + # Copyright : David Tschumperle 1.15 + # ( http://www.greyc.ensicaen.fr/~dtschump/ ) 1.16 + # 1.17 + # License : CeCILL v2.0 1.18 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.19 + # 1.20 + # This software is governed by the CeCILL license under French law and 1.21 + # abiding by the rules of distribution of free software. You can use, 1.22 + # modify and/ or redistribute the software under the terms of the CeCILL 1.23 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.24 + # "http://www.cecill.info". 1.25 + # 1.26 + # As a counterpart to the access to the source code and rights to copy, 1.27 + # modify and redistribute granted by the license, users are provided only 1.28 + # with a limited warranty and the software's author, the holder of the 1.29 + # economic rights, and the successive licensors have only limited 1.30 + # liability. 1.31 + # 1.32 + # In this respect, the user's attention is drawn to the risks associated 1.33 + # with loading, using, modifying and/or developing or reproducing the 1.34 + # software by the user in light of its specific status of free software, 1.35 + # that may mean that it is complicated to manipulate, and that also 1.36 + # therefore means that it is reserved for developers and experienced 1.37 + # professionals having in-depth computer knowledge. Users are therefore 1.38 + # encouraged to load and test the software's suitability as regards their 1.39 + # requirements in conditions enabling the security of their systems and/or 1.40 + # data to be ensured and, more generally, to use and operate it in the 1.41 + # same conditions as regards security. 1.42 + # 1.43 + # The fact that you are presently reading this means that you have had 1.44 + # knowledge of the CeCILL license and that you accept its terms. 1.45 + # 1.46 +*/ 1.47 + 1.48 +#include "CImg.h" 1.49 +using namespace cimg_library; 1.50 + 1.51 +// The lines below are necessary when using a non-standard compiler as visualcpp6. 1.52 +#ifdef cimg_use_visualcpp6 1.53 +#define std 1.54 +#endif 1.55 +#ifdef min 1.56 +#undef min 1.57 +#undef max 1.58 +#endif 1.59 + 1.60 +// Compute fractional anisotropy (FA) of a tensor 1.61 +//------------------------------------------- 1.62 +template<typename T> float get_FA(const T& val1, const T& val2, const T& val3) { 1.63 + const float 1.64 + l1 = val1>0?val1:0, l2 = val2>0?val2:0, l3 = val3>0?val3:0, 1.65 + lm = (l1+l2+l3)/3, 1.66 + tr2 = 2*( l1*l1 + l2*l2 + l3*l3 ), 1.67 + ll1 = l1-lm, 1.68 + ll2 = l2-lm, 1.69 + ll3 = l3-lm; 1.70 + if (tr2>0) return (float)std::sqrt( 3*(ll1*ll1 + ll2*ll2 + ll3*ll3)/tr2 ); 1.71 + return 0; 1.72 +} 1.73 + 1.74 +// Insert an ellipsoid in a CImg 3D scene 1.75 +//---------------------------------------- 1.76 +template<typename t,typename tp,typename tf,typename tc> 1.77 +void insert_ellipsoid(const CImg<t>& tensor,const float X,const float Y,const float Z,const float tfact, 1.78 + const float vx, const float vy, const float vz, 1.79 + CImgList<tp>& points, CImgList<tf>& faces, CImgList<tc>& colors, 1.80 + const unsigned int res1 = 20, const unsigned int res2 = 20) { 1.81 + 1.82 + // Compute eigen elements 1.83 + const float l1 = tensor[0], l2 = tensor[1], l3 = tensor[2], fa = get_FA(l1,l2,l3); 1.84 + 1.85 + CImg<> vec = CImg<>::matrix(tensor[3],tensor[6],tensor[9], 1.86 + tensor[4],tensor[7],tensor[10], 1.87 + tensor[5],tensor[8],tensor[11]); 1.88 + const int 1.89 + r = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[3]),255.0f), 1.90 + g = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[4]),255.0f), 1.91 + b = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[5]),255.0f); 1.92 + 1.93 + // Define mesh points 1.94 + const unsigned int N0 = points.size; 1.95 + for (unsigned int v=1; v<res2; v++) 1.96 + for (unsigned int u=0; u<res1; u++) { 1.97 + const float 1.98 + alpha = (float)(u*2*cimg::valuePI/res1), 1.99 + beta = (float)(-cimg::valuePI/2 + v*cimg::valuePI/res2), 1.100 + x = (float)(tfact*l1*std::cos(beta)*std::cos(alpha)), 1.101 + y = (float)(tfact*l2*std::cos(beta)*std::sin(alpha)), 1.102 + z = (float)(tfact*l3*std::sin(beta)); 1.103 + points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(x,y,z)).mul(CImg<tp>::vector(vx,vy,vz))); 1.104 + } 1.105 + const unsigned int N1 = points.size; 1.106 + points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(0,0,-l3*tfact))); 1.107 + points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(0,0,l3*tfact))); 1.108 + points[points.size-2](0)*=vx; points[points.size-2](1)*=vy; points[points.size-2](2)*=vz; 1.109 + points[points.size-1](0)*=vx; points[points.size-1](1)*=vy; points[points.size-1](2)*=vz; 1.110 + 1.111 + // Define mesh triangles 1.112 + for (unsigned int vv=0; vv<res2-2; vv++) 1.113 + for (unsigned int uu=0; uu<res1; uu++) { 1.114 + const int nv = (vv+1)%(res2-1), nu = (uu+1)%res1; 1.115 + faces.insert(CImg<tf>::vector(N0+res1*vv+nu,N0+res1*nv+uu,N0+res1*vv+uu)); 1.116 + faces.insert(CImg<tf>::vector(N0+res1*vv+nu,N0+res1*nv+nu,N0+res1*nv+uu)); 1.117 + colors.insert(CImg<tc>::vector(r,g,b)); 1.118 + colors.insert(CImg<tc>::vector(r,g,b)); 1.119 + } 1.120 + for (unsigned int uu=0; uu<res1; uu++) { 1.121 + const int nu = (uu+1)%res1; 1.122 + faces.insert(CImg<tf>::vector(N0+nu,N0+uu,N1)); 1.123 + faces.insert(CImg<tf>::vector(N0+res1*(res2-2)+nu, N1+1,N0+res1*(res2-2)+uu)); 1.124 + colors.insert(CImg<tc>::vector(r,g,b)); 1.125 + colors.insert(CImg<tc>::vector(r,g,b)); 1.126 + } 1.127 +} 1.128 + 1.129 +// Insert a fiber in a CImg 3D scene 1.130 +//----------------------------------- 1.131 +template<typename T,typename te,typename tp, typename tf, typename tc> 1.132 +void insert_fiber(const CImg<T>& fiber, const CImg<te>& eigen, const CImg<tc>& palette, 1.133 + const int xm, const int ym, const int zm, 1.134 + const float vx, const float vy, const float vz, 1.135 + CImgList<tp>& points, CImgList<tf>& primitives, CImgList<tc>& colors) { 1.136 + const int N0 = points.size; 1.137 + float x0 = fiber(0,0), y0 = fiber(0,1), z0 = fiber(0,2), fa0 = eigen.linear_atXYZ(x0,y0,z0,12); 1.138 + points.insert(CImg<>::vector(vx*(x0-xm),vy*(y0-ym),vz*(z0-zm))); 1.139 + for (int l=1; l<fiber.dimx(); l++) { 1.140 + float x1 = fiber(l,0), y1 = fiber(l,1), z1 = fiber(l,2), fa1 = eigen.linear_atXYZ(x1,y1,z1,12); 1.141 + points.insert(CImg<tp>::vector(vx*(x1-xm),vy*(y1-ym),vz*(z1-zm))); 1.142 + primitives.insert(CImg<tf>::vector(N0+l-1,N0+l)); 1.143 + const unsigned char 1.144 + icol = (unsigned char)(fa0*255), 1.145 + r = palette(icol,0), 1.146 + g = palette(icol,1), 1.147 + b = palette(icol,2); 1.148 + colors.insert(CImg<unsigned char>::vector(r,g,b)); 1.149 + x0=x1; y0=y1; z0=z1; fa0=fa1; 1.150 + } 1.151 +} 1.152 + 1.153 +// Compute fiber tracking using 4th-order Runge Kutta integration 1.154 +//----------------------------------------------------------------- 1.155 +template<typename T> 1.156 +CImg<> get_fibertrack(CImg<T>& eigen, 1.157 + const int X0, const int Y0, const int Z0, const float lmax=100, 1.158 + const float dl = 0.1f, const float FAmin=0.7f, const float cmin=0.5f) { 1.159 + 1.160 +#define align_eigen(i,j,k) \ 1.161 + { T &u = eigen(i,j,k,3), &v = eigen(i,j,k,4), &w = eigen(i,j,k,5); \ 1.162 + if (u*cu+v*cv+w*cw<0) { u=-u; v=-v; w=-w; }} 1.163 + 1.164 + CImgList<> resf; 1.165 + 1.166 + // Forward tracking 1.167 + float normU = 0, normpU = 0, l = 0, X = (float)X0, Y = (float)Y0, Z = (float)Z0; 1.168 + T 1.169 + pu = eigen(X0,Y0,Z0,3), 1.170 + pv = eigen(X0,Y0,Z0,4), 1.171 + pw = eigen(X0,Y0,Z0,5); 1.172 + normpU = (float)std::sqrt(pu*pu+pv*pv+pw*pw); 1.173 + bool stopflag = false; 1.174 + 1.175 + while (!stopflag) { 1.176 + if (X<0 || X>eigen.dimx()-1 || Y<0 || Y>eigen.dimy()-1 || Z<0 || Z>eigen.dimz()-1 || 1.177 + eigen((int)X,(int)Y,(int)Z,12)<FAmin || l>lmax) stopflag = true; 1.178 + else { 1.179 + resf.insert(CImg<>::vector(X,Y,Z)); 1.180 + 1.181 + const int 1.182 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>=eigen.dimx())?eigen.dimx()-1:cx+1, 1.183 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>=eigen.dimy())?eigen.dimy()-1:cy+1, 1.184 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>=eigen.dimz())?eigen.dimz()-1:cz+1; 1.185 + const T cu = eigen(cx,cy,cz,3), cv = eigen(cx,cy,cz,4), cw = eigen(cx,cy,cz,5); 1.186 + 1.187 + align_eigen(px,py,pz); align_eigen(cx,py,pz); align_eigen(nx,py,pz); 1.188 + align_eigen(px,cy,pz); align_eigen(cx,cy,pz); align_eigen(nx,cy,pz); 1.189 + align_eigen(px,ny,pz); align_eigen(cx,ny,pz); align_eigen(nx,ny,pz); 1.190 + align_eigen(px,py,cz); align_eigen(cx,py,cz); align_eigen(nx,py,cz); 1.191 + align_eigen(px,cy,cz); align_eigen(nx,cy,cz); 1.192 + align_eigen(px,ny,cz); align_eigen(cx,ny,cz); align_eigen(nx,ny,cz); 1.193 + align_eigen(px,py,nz); align_eigen(cx,py,nz); align_eigen(nx,py,nz); 1.194 + align_eigen(px,cy,nz); align_eigen(cx,cy,nz); align_eigen(nx,cy,nz); 1.195 + align_eigen(px,ny,nz); align_eigen(cx,ny,nz); align_eigen(nx,ny,nz); 1.196 + 1.197 + const T 1.198 + u0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,3), 1.199 + v0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,4), 1.200 + w0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,5), 1.201 + u1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,3), 1.202 + v1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,4), 1.203 + w1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,5), 1.204 + u2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,3), 1.205 + v2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,4), 1.206 + w2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,5), 1.207 + u3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,3), 1.208 + v3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,4), 1.209 + w3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,5); 1.210 + T 1.211 + u = u0/3 + 2*u1/3 + 2*u2/3 + u3/3, 1.212 + v = v0/3 + 2*v1/3 + 2*v2/3 + v3/3, 1.213 + w = w0/3 + 2*w1/3 + 2*w2/3 + w3/3; 1.214 + if (u*pu+v*pv+w*pw<0) { u=-u; v=-v; w=-w; } 1.215 + normU = (float)std::sqrt(u*u+v*v+w*w); 1.216 + const float scal = (u*pu+v*pv+w*pw)/(normU*normpU); 1.217 + if (scal<cmin) stopflag=true; 1.218 + 1.219 + X+=(pu=u); Y+=(pv=v); Z+=(pw=w); 1.220 + normpU = normU; 1.221 + l+=dl; 1.222 + } 1.223 + } 1.224 + 1.225 + // Backward tracking 1.226 + l = dl; X = (float)X0; Y = (float)Y0; Z = (float)Z0; 1.227 + pu = eigen(X0,Y0,Z0,3); 1.228 + pv = eigen(X0,Y0,Z0,4); 1.229 + pw = eigen(X0,Y0,Z0,5); 1.230 + normpU = (float)std::sqrt(pu*pu+pv*pv+pw*pw); 1.231 + stopflag = false; 1.232 + 1.233 + while (!stopflag) { 1.234 + if (X<0 || X>eigen.dimx()-1 || Y<0 || Y>eigen.dimy()-1 || Z<0 || Z>eigen.dimz()-1 || 1.235 + eigen((int)X,(int)Y,(int)Z,12)<FAmin || l>lmax) stopflag = true; 1.236 + else { 1.237 + 1.238 + const int 1.239 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>=eigen.dimx())?eigen.dimx()-1:cx+1, 1.240 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>=eigen.dimy())?eigen.dimy()-1:cy+1, 1.241 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>=eigen.dimz())?eigen.dimz()-1:cz+1; 1.242 + const T cu = eigen(cx,cy,cz,3), cv = eigen(cx,cy,cz,4), cw = eigen(cx,cy,cz,5); 1.243 + 1.244 + align_eigen(px,py,pz); align_eigen(cx,py,pz); align_eigen(nx,py,pz); 1.245 + align_eigen(px,cy,pz); align_eigen(cx,cy,pz); align_eigen(nx,cy,pz); 1.246 + align_eigen(px,ny,pz); align_eigen(cx,ny,pz); align_eigen(nx,ny,pz); 1.247 + align_eigen(px,py,cz); align_eigen(cx,py,cz); align_eigen(nx,py,cz); 1.248 + align_eigen(px,cy,cz); align_eigen(nx,cy,cz); 1.249 + align_eigen(px,ny,cz); align_eigen(cx,ny,cz); align_eigen(nx,ny,cz); 1.250 + align_eigen(px,py,nz); align_eigen(cx,py,nz); align_eigen(nx,py,nz); 1.251 + align_eigen(px,cy,nz); align_eigen(cx,cy,nz); align_eigen(nx,cy,nz); 1.252 + align_eigen(px,ny,nz); align_eigen(cx,ny,nz); align_eigen(nx,ny,nz); 1.253 + 1.254 + const T 1.255 + u0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,3), 1.256 + v0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,4), 1.257 + w0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,5), 1.258 + u1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,3), 1.259 + v1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,4), 1.260 + w1 = 0.5f*dl*eigen.linear_atXYZ(X+u0,Y+v0,Z+w0,5), 1.261 + u2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,3), 1.262 + v2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,4), 1.263 + w2 = 0.5f*dl*eigen.linear_atXYZ(X+u1,Y+v1,Z+w1,5), 1.264 + u3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,3), 1.265 + v3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,4), 1.266 + w3 = 0.5f*dl*eigen.linear_atXYZ(X+u2,Y+v2,Z+w2,5); 1.267 + T 1.268 + u = u0/3 + 2*u1/3 + 2*u2/3 + u3/3, 1.269 + v = v0/3 + 2*v1/3 + 2*v2/3 + v3/3, 1.270 + w = w0/3 + 2*w1/3 + 2*w2/3 + w3/3; 1.271 + if (u*pu+v*pv+w*pw<0) { u=-u; v=-v; w=-w; } 1.272 + normU = (float)std::sqrt(u*u+v*v+w*w); 1.273 + const float scal = (u*pu+v*pv+w*pw)/(normU*normpU); 1.274 + if (scal<cmin) stopflag=true; 1.275 + 1.276 + X-=(pu=u); Y-=(pv=v); Z-=(pw=w); 1.277 + normpU=normU; 1.278 + l+=dl; 1.279 + 1.280 + resf.insert(CImg<>::vector(X,Y,Z),0); 1.281 + } 1.282 + } 1.283 + 1.284 + return resf.get_append('x'); 1.285 +} 1.286 + 1.287 +// Main procedure 1.288 +//---------------- 1.289 +int main(int argc,char **argv) { 1.290 + 1.291 + // Read and init data 1.292 + //-------------------- 1.293 + cimg_usage("A viewer of Diffusion-Tensor MRI volumes."); 1.294 + const char *file_i = cimg_option("-i",(char*)0,"Input : Filename of tensor field (volume wxhxdx6)"); 1.295 + const char* vsize = cimg_option("-vsize","1x1x1","Input : Voxel aspect"); 1.296 + const bool normalize = cimg_option("-normalize",true,"Input : Enable tensor normalization"); 1.297 + const char *file_f = cimg_option("-f",(char*)0,"Input : Input fibers\n"); 1.298 + const float dl = cimg_option("-dl",0.5f,"Fiber computation : Integration step"); 1.299 + const float famin = cimg_option("-famin",0.3f,"Fiber computation : Fractional Anisotropy threshold"); 1.300 + const float cmin = cimg_option("-cmin",0.2f,"Fiber computation : Curvature threshold"); 1.301 + const float lmin = cimg_option("-lmin",10.0f,"Fiber computation : Minimum length\n"); 1.302 + const float lmax = cimg_option("-lmax",1000.0f,"Fiber computation : Maximum length\n"); 1.303 + const float tfact = cimg_option("-tfact",1.2f,"Display : Tensor size factor"); 1.304 + const char *bgcolor = cimg_option("-bg","0,0,0","Display : Background color"); 1.305 + unsigned int bgr=0, bgg=0, bgb=0; 1.306 + std::sscanf(bgcolor,"%u%*c%u%*c%u",&bgr,&bgg,&bgb); 1.307 + 1.308 + CImg<> tensors; 1.309 + if (file_i) { 1.310 + std::fprintf(stderr,"\n- Loading tensors '%s'",cimg::basename(file_i)); 1.311 + tensors.load(file_i); 1.312 + } else { 1.313 + // Create a synthetic tensor field here 1.314 + std::fprintf(stderr,"\n- No input files : Creating a synthetic tensor field"); 1.315 + tensors.assign(32,32,32,6); 1.316 + const CImg<> Id = CImg<>::diagonal(0.3f,0.3f,0.3f); 1.317 + cimg_forXYZ(tensors,x,y,z) { 1.318 + const float 1.319 + u = x-tensors.dimx()/2.0f, 1.320 + v = y-tensors.dimy()/2.0f, 1.321 + w = z-tensors.dimz()/2.0f, 1.322 + norm = (float)std::sqrt(1e-5f+u*u+v*v+w*w), 1.323 + nu = u/norm, nv = v/norm, nw = w/norm; 1.324 + const CImg<> 1.325 + dir1 = CImg<>::vector(nu,nv,nw), 1.326 + dir2 = CImg<>::vector(-nv,nu,nw), 1.327 + dir3 = CImg<>::vector(nw*(nv-nu),-nw*(nu+nv),nu*nu+nv*nv); 1.328 + tensors.set_tensor_at(2.0*dir1*dir1.get_transpose() + 1.329 + 1.0*dir2*dir2.get_transpose() + 1.330 + 0.7*dir3*dir3.get_transpose(), 1.331 + x,y,z); 1.332 + } 1.333 + } 1.334 + float voxw=1,voxh=1,voxd=1; 1.335 + std::sscanf(vsize,"%f%*c%f%*c%f",&voxw,&voxh,&voxd); 1.336 + 1.337 + std::fprintf(stderr," : %ux%ux%u image, voxsize=%gx%gx%g.", 1.338 + tensors.dimx(),tensors.dimy(),tensors.dimz(), 1.339 + voxw,voxh,voxd); 1.340 + 1.341 + 1.342 + CImgList<> fibers; 1.343 + if (file_f) { 1.344 + std::fprintf(stderr,"\n- Loading fibers '%s'.",cimg::basename(file_f)); 1.345 + fibers.load(file_f); 1.346 + } 1.347 + 1.348 + const CImg<unsigned char> fiber_palette = 1.349 + CImg<>(2,1,1,3).fill(200,255,0,255,0,200).RGBtoHSV().resize(256,1,1,3,3).HSVtoRGB(); 1.350 + 1.351 + // Compute eigen elements 1.352 + //------------------------ 1.353 + std::fprintf(stderr,"\n- Compute eigen elements."); 1.354 + CImg<unsigned char> coloredFA(tensors.dimx(),tensors.dimy(),tensors.dimz(),3); 1.355 + CImg<> eigen(tensors.dimx(),tensors.dimy(),tensors.dimz(),13); 1.356 + CImg<> val,vec; 1.357 + float eigmax = 0; 1.358 + cimg_forXYZ(tensors,x,y,z) { 1.359 + tensors.get_tensor_at(x,y,z).symmetric_eigen(val,vec); 1.360 + eigen(x,y,z,0) = val[0]; eigen(x,y,z,1) = val[1]; eigen(x,y,z,2) = val[2]; 1.361 + if (val[0]<0) val[0]=0; 1.362 + if (val[1]<0) val[1]=0; 1.363 + if (val[2]<0) val[2]=0; 1.364 + if (val[0]>eigmax) eigmax = val[0]; 1.365 + eigen(x,y,z,3) = vec(0,0); eigen(x,y,z,4) = vec(0,1); eigen(x,y,z,5) = vec(0,2); 1.366 + eigen(x,y,z,6) = vec(1,0); eigen(x,y,z,7) = vec(1,1); eigen(x,y,z,8) = vec(1,2); 1.367 + eigen(x,y,z,9) = vec(2,0); eigen(x,y,z,10) = vec(2,1); eigen(x,y,z,11) = vec(2,2); 1.368 + const float fa = get_FA(val[0],val[1],val[2]); 1.369 + eigen(x,y,z,12) = fa; 1.370 + const int 1.371 + r = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,0))), 1.372 + g = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,1))), 1.373 + b = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,2))); 1.374 + coloredFA(x,y,z,0) = (unsigned char)r; 1.375 + coloredFA(x,y,z,1) = (unsigned char)g; 1.376 + coloredFA(x,y,z,2) = (unsigned char)b; 1.377 + } 1.378 + tensors.assign(); 1.379 + std::fprintf(stderr,"\n- Maximum diffusivity = %g, Maximum FA = %g",eigmax,eigen.get_shared_channel(12).max()); 1.380 + if (normalize) { 1.381 + std::fprintf(stderr,"\n- Normalize tensors."); 1.382 + eigen.get_shared_channels(0,2)/=eigmax; 1.383 + } 1.384 + 1.385 + // Init display and begin user interaction 1.386 + //----------------------------------------- 1.387 + std::fprintf(stderr,"\n- Open user window."); 1.388 + CImgDisplay disp(256,256,"DTMRI Viewer",0); 1.389 + CImgDisplay disp3d(800,600,"3D Local View",0,false,true); 1.390 + unsigned int XYZ[3]; 1.391 + XYZ[0] = eigen.dimx()/2; XYZ[1] = eigen.dimy()/2; XYZ[2] = eigen.dimz()/2; 1.392 + 1.393 + while (!disp.is_closed && disp.key!=cimg::keyQ && disp.key!=cimg::keyESC) { 1.394 + const CImg<int> s = coloredFA.get_select(disp,2,XYZ); 1.395 + if (!disp.is_closed) switch (disp.key) { 1.396 + 1.397 + // Open 3D visualization window 1.398 + //----------------------------- 1.399 + case cimg::keyA: 1.400 + case 0: { 1.401 + unsigned char white[1] = { 255 }; 1.402 + disp3d.display(CImg<unsigned char>(disp3d.dimx(),disp3d.dimy(),1,1,0).draw_text(10,10,"Please wait...",white)).show(); 1.403 + int xm,ym,zm,xM,yM,zM; 1.404 + if (!disp.key) { xm=s[0]; ym=s[1]; zm=s[2]; xM=s[3]; yM=s[4]; zM=s[5]; } 1.405 + else { xm=ym=zm=0; xM=eigen.dimx()-1; yM=eigen.dimy()-1; zM=eigen.dimy()-1; } 1.406 + const CImg<> img = eigen.get_crop(xm,ym,zm,xM,yM,zM); 1.407 + CImgList<> points; 1.408 + CImgList<unsigned int> primitives; 1.409 + CImgList<unsigned char> colors; 1.410 + 1.411 + // Add ellipsoids to the 3D scene 1.412 + int X = img.dimx()/2, Y = img.dimy()/2, Z = img.dimz()/2; 1.413 + { 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); } 1.414 + { 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); } 1.415 + { 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); } 1.416 + 1.417 + // Add computed fibers to the 3D scene 1.418 + const CImg<> veigen = eigen.get_crop(xm,ym,zm,xM,yM,zM); 1.419 + cimglist_for(fibers,l) { 1.420 + const CImg<>& fiber = fibers[l]; 1.421 + if (fiber.dimx()) insert_fiber(fiber,eigen,fiber_palette, 1.422 + xm,ym,zm,voxw,voxh,voxd, 1.423 + points,primitives,colors); 1.424 + } 1.425 + 1.426 + // Display 3D object 1.427 + CImg<unsigned char> visu = CImg<unsigned char>(3,disp3d.dimx(),disp3d.dimy(),1,0).fill(bgr,bgg,bgb).permute_axes("yzvx"); 1.428 + bool stopflag = false; 1.429 + while (!disp3d.is_closed && !stopflag) { 1.430 + visu.display_object3d(disp3d,points,primitives,colors,true,4,-1,false,800,0.05f,1.0f); 1.431 + switch (disp3d.key) { 1.432 + case cimg::keyM: { // Create movie 1.433 + std::fprintf(stderr,"\n- Movie mode.\n"); 1.434 + const unsigned int N = 256; 1.435 + CImg<> pts = points.get_append('x'); 1.436 + CImgList<> cpoints(points); 1.437 + CImg<> x = pts.get_shared_line(0), y = pts.get_shared_line(1), z = pts.get_shared_line(2); 1.438 + float 1.439 + xm, xM = x.maxmin(xm), 1.440 + ym, yM = y.maxmin(ym), 1.441 + zm, zM = z.maxmin(zm), 1.442 + ratio = 2.0f*cimg::min(visu.dimx(),visu.dimy())/(3.0f*cimg::max(xM-xm,yM-ym,zM-zm)), 1.443 + dx = 0.5f*(xM+xm), dy = 0.5f*(yM+ym), dz = 0.5f*(zM+zm); 1.444 + cimglist_for(points,l) { 1.445 + cpoints(l,0) = (float)((points(l,0)-dx)*ratio); 1.446 + cpoints(l,1) = (float)((points(l,1)-dy)*ratio); 1.447 + cpoints(l,2) = (float)((points(l,2)-dz)*ratio); 1.448 + } 1.449 + 1.450 + for (unsigned int i=0; i<N; i++) { 1.451 + std::fprintf(stderr,"\r- Frame %u/%u.",i,N); 1.452 + const float alpha = (float)(i*2*cimg::valuePI/N); 1.453 + const CImg<> rot = CImg<>::rotation_matrix(0,1,0,alpha)*CImg<>::rotation_matrix(1,0,0,1.30f); 1.454 + CImgList<> rotated(cpoints); 1.455 + cimglist_for(rotated,l) rotated[l] = rot*cpoints[l]; 1.456 + visu.fill(0).draw_object3d(visu.dimx()/2.0f,visu.dimy()/2.0f,-500.0f,rotated,primitives,colors, 1.457 + 4,false,800.0f,visu.dimx()/2.0f,visu.dimy()/2.0f,-800.0f,0.05f,1.0f).display(disp3d); 1.458 + visu.save("frame.png",i); 1.459 + } 1.460 + visu.fill(0); 1.461 + } break; 1.462 + default: stopflag = true; 1.463 + } 1.464 + } 1.465 + if (disp3d.is_fullscreen) disp3d.toggle_fullscreen().resize(800,600).close(); 1.466 + } break; 1.467 + 1.468 + // Compute region statistics 1.469 + //--------------------------- 1.470 + case cimg::keyR: { 1.471 + std::fprintf(stderr,"\n- Statistics computation. Select region."); std::fflush(stderr); 1.472 + const CImg<int> s = coloredFA.get_select(disp,2,XYZ); 1.473 + int xm,ym,zm,xM,yM,zM; 1.474 + if (!disp.key) { xm=s[0]; ym=s[1]; zm=s[2]; xM=s[3]; yM=s[4]; zM=s[5]; } 1.475 + else { xm=ym=zm=0; xM=eigen.dimx()-1; yM=eigen.dimy()-1; zM=eigen.dimy()-1; } 1.476 + const CImg<> img = eigen.get_crop(xm,ym,zm,xM,yM,zM); 1.477 + std::fprintf(stderr,"\n- Mean diffusivity = %g, Mean FA = %g\n", 1.478 + eigen.get_shared_channel(0).mean(), 1.479 + eigen.get_shared_channel(12).mean()); 1.480 + } break; 1.481 + 1.482 + // Track fiber bundle (single region) 1.483 + //---------------------------------- 1.484 + case cimg::keyF: { 1.485 + std::fprintf(stderr,"\n- Tracking mode (single region). Select starting region.\n"); std::fflush(stderr); 1.486 + const CImg<int> s = coloredFA.get_select(disp,2,XYZ); 1.487 + const unsigned int N = fibers.size; 1.488 + for (int z=s[2]; z<=s[5]; z++) 1.489 + for (int y=s[1]; y<=s[4]; y++) 1.490 + for (int x=s[0]; x<=s[3]; x++) { 1.491 + const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin); 1.492 + if (fiber.dimx()>lmin) { 1.493 + std::fprintf(stderr,"\rFiber %u : Starting from (%d,%d,%d)\t\t",fibers.size,x,y,z); 1.494 + fibers.insert(fiber); 1.495 + } 1.496 + } 1.497 + std::fprintf(stderr,"\n- %u fiber(s) added (total %u).",fibers.size-N,fibers.size); 1.498 + } break; 1.499 + 1.500 + // Track fiber bundle (double regions) 1.501 + //------------------------------------ 1.502 + case cimg::keyG: { 1.503 + std::fprintf(stderr,"\n- Tracking mode (double region). Select starting region."); std::fflush(stderr); 1.504 + const CImg<int> s = coloredFA.get_select(disp,2,XYZ); 1.505 + std::fprintf(stderr," Select ending region."); std::fflush(stderr); 1.506 + const CImg<int> ns = coloredFA.get_select(disp,2,XYZ); 1.507 + const unsigned int N = fibers.size; 1.508 + 1.509 + // Track from start to end 1.510 + for (int z=s[2]; z<=s[5]; z++) 1.511 + for (int y=s[1]; y<=s[4]; y++) 1.512 + for (int x=s[0]; x<=s[3]; x++) { 1.513 + const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin); 1.514 + if (fiber.dimx()>lmin) { 1.515 + bool valid_fiber = false; 1.516 + cimg_forX(fiber,k) { 1.517 + const int fx = (int)fiber(k,0), fy = (int)fiber(k,1), fz = (int)fiber(k,2); 1.518 + if (fx>=ns[0] && fx<=ns[3] && 1.519 + fy>=ns[1] && fy<=ns[4] && 1.520 + fz>=ns[2] && fz<=ns[5]) valid_fiber = true; 1.521 + } 1.522 + if (valid_fiber) fibers.insert(fiber); 1.523 + } 1.524 + } 1.525 + 1.526 + // Track from end to start 1.527 + { for (int z=ns[2]; z<=ns[5]; z++) 1.528 + for (int y=ns[1]; y<=ns[4]; y++) 1.529 + for (int x=ns[0]; x<=ns[3]; x++) { 1.530 + const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin); 1.531 + if (fiber.dimx()>lmin) { 1.532 + bool valid_fiber = false; 1.533 + cimg_forX(fiber,k) { 1.534 + const int fx = (int)fiber(k,0), fy = (int)fiber(k,1), fz = (int)fiber(k,2); 1.535 + if (fx>=s[0] && fx<=s[3] && 1.536 + fy>=s[1] && fy<=s[4] && 1.537 + fz>=s[2] && fz<=s[5]) valid_fiber = true; 1.538 + } 1.539 + if (valid_fiber) { 1.540 + std::fprintf(stderr,"\rFiber %u : Starting from (%d,%d,%d)\t\t",fibers.size,x,y,z); 1.541 + fibers.insert(fiber); 1.542 + } 1.543 + } 1.544 + }} 1.545 + 1.546 + std::fprintf(stderr," %u fiber(s) added (total %u).",fibers.size-N,fibers.size); 1.547 + } break; 1.548 + 1.549 + // Clear fiber bundle 1.550 + //------------------- 1.551 + case cimg::keyC: { 1.552 + std::fprintf(stderr,"\n- Fibers removed."); 1.553 + fibers.assign(); 1.554 + } break; 1.555 + 1.556 + // Save fibers 1.557 + //------------- 1.558 + case cimg::keyS: { 1.559 + fibers.save("fibers.cimg"); 1.560 + std::fprintf(stderr,"\n- Fibers saved."); 1.561 + } break; 1.562 + 1.563 + } 1.564 + } 1.565 + 1.566 + std::fprintf(stderr,"\n- Exit.\n\n\n"); 1.567 + return 0; 1.568 +}