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