dtmri_view.cpp

this a image processing program

  const char *file_i   = cimg_option("-i",(char*)0,"Input : Filename of tensor field (volume wxhxdx6)");
  const char* vsize    = cimg_option("-vsize","1x1x1","Input : Voxel aspect");
  const bool normalize = cimg_option("-normalize",true,"Input : Enable tensor normalization");
  const char *file_f   = cimg_option("-f",(char*)0,"Input : Input fibers\n");
  const float dl       = cimg_option("-dl",0.5f,"Fiber computation : Integration step");
  const float famin    = cimg_option("-famin",0.3f,"Fiber computation : Fractional Anisotropy threshold");
  const float cmin     = cimg_option("-cmin",0.2f,"Fiber computation : Curvature threshold");
  const float lmin     = cimg_option("-lmin",10.0f,"Fiber computation : Minimum length\n");
  const float lmax     = cimg_option("-lmax",1000.0f,"Fiber computation : Maximum length\n");
  const float tfact    = cimg_option("-tfact",1.2f,"Display : Tensor size factor");

  CImg<> tensors;
  if (file_i) {
    std::fprintf(stderr,"\n- Loading tensors '%s'",cimg::basename(file_i));
    tensors.load(file_i);
  } else {
    // Create a synthetic tensor field here
    std::fprintf(stderr,"\n- No input files : Creating a synthetic tensor field");
    tensors.assign(32,32,32,6);
    const CImg<> Id = CImg<>::diagonal(0.3f,0.3f,0.3f);
    cimg_forXYZ(tensors,x,y,z) {
      const float
	u = x-tensors.dimx()/2.0f,
	v = y-tensors.dimy()/2.0f,
	w = z-tensors.dimz()/2.0f,
	norm = (float)std::sqrt(1e-5f+u*u+v*v+w*w),
	nu = u/norm, nv = v/norm, nw = w/norm;
      const CImg<>
	dir1 = CImg<>::vector(nu,nv,nw),
	dir2 = CImg<>::vector(-nv,nu,nw),
	dir3 = CImg<>::vector(nw*(nv-nu),-nw*(nu+nv),nu*nu+nv*nv);
      tensors.set_tensor_at(2.0f*dir1*dir1.get_transpose() +
			    1.0f*dir2*dir2.get_transpose() +
			    0.7f*dir3*dir3.get_transpose(),
			    x,y,z);
    }
  }
  float voxw=1,voxh=1,voxd=1;
  std::sscanf(vsize,"%f%*c%f%*c%f",&voxw,&voxh,&voxd);

  std::fprintf(stderr," : %ux%ux%u image, voxsize=%gx%gx%g.",
	       tensors.width,tensors.height,tensors.depth,
	       voxw,voxh,voxd);


  CImgList<> fibers;
  if (file_f) {
    std::fprintf(stderr,"\n- Loading fibers '%s'.",cimg::basename(file_f));
    fibers.load(file_f);
  }

  const CImg<unsigned char> fiber_palette =
    CImg<>(2,1,1,3).fill(200,255,0,255,0,200).RGBtoHSV().resize(256,1,1,3,3).HSVtoRGB();

  // Compute eigen elements
  //------------------------
  std::fprintf(stderr,"\n- Compute eigen elements.");
  CImg<unsigned char> coloredFA(tensors.dimx(),tensors.dimy(),tensors.dimz(),3);
  CImg<> eigen(tensors.dimx(),tensors.dimy(),tensors.dimz(),13);
  CImg<> val,vec;
  float eigmax = 0;
  cimg_forXYZ(tensors,x,y,z) {
    tensors.get_tensor_at(x,y,z).symmetric_eigen(val,vec);
    eigen(x,y,z,0) = val[0]; eigen(x,y,z,1) = val[1]; eigen(x,y,z,2) = val[2];
    if (val[0]<0) val[0]=0;
    if (val[1]<0) val[1]=0;
    if (val[2]<0) val[2]=0;
    if (val[0]>eigmax) eigmax = val[0];
    eigen(x,y,z,3) = vec(0,0); eigen(x,y,z,4)  = vec(0,1); eigen(x,y,z,5)  = vec(0,2);
    eigen(x,y,z,6) = vec(1,0); eigen(x,y,z,7)  = vec(1,1); eigen(x,y,z,8)  = vec(1,2);
    eigen(x,y,z,9) = vec(2,0); eigen(x,y,z,10) = vec(2,1); eigen(x,y,z,11) = vec(2,2);
    const float fa = get_FA(val[0],val[1],val[2]);
    eigen(x,y,z,12) = fa;
    const int
      r = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,0))),
      g = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,1))),
      b = (int)cimg::min(255.0f,1.5f*cimg::abs(255*fa*vec(0,2)));
    coloredFA(x,y,z,0) = (unsigned char)r;
    coloredFA(x,y,z,1) = (unsigned char)g;
    coloredFA(x,y,z,2) = (unsigned char)b;
  }
  tensors.assign();
  CImgStats stats(eigen.get_shared_channel(12),false);
  std::fprintf(stderr,"\n- Maximum diffusivity = %g, Maximum FA = %g",eigmax,stats.max);
  if (normalize) {
    std::fprintf(stderr,"\n- Normalize tensors.");
    eigen.get_shared_channels(0,2)/=eigmax;
  }

  // Init display and begin user interaction
  //-----------------------------------------
  std::fprintf(stderr,"\n- Open user window.");
  CImgDisplay disp(256,256,"DTMRI Viewer",0);
  CImgDisplay disp3d(800,600,"3D Local View",0,3,false,true);
  int s[6];
  unsigned int XYZ[3];
  XYZ[0] = eigen.dimx()/2; XYZ[1] = eigen.dimy()/2; XYZ[2] = eigen.dimz()/2;

  while (!disp.is_closed) {
    coloredFA.feature_selection(s,2,disp,XYZ);
    if (!disp.is_closed) switch (disp.key) {

      // Open 3D visualization window
      //-----------------------------
    case cimg::keyA:
    case 0: {
      unsigned char white[1] = { 255 };
      disp3d.display(CImg<unsigned char>(disp3d.dimx(),disp3d.dimy(),1,1,0).draw_text("Please wait...",10,10,white)).show();
      int xm,ym,zm,xM,yM,zM;
      if (!disp.key) { xm=s[0]; ym=s[1]; zm=s[2]; xM=s[3]; yM=s[4]; zM=s[5]; }
      else { xm=ym=zm=0; xM=eigen.dimx()-1; yM=eigen.dimy()-1; zM=eigen.dimy()-1; }
      const CImg<> img = eigen.get_crop(xm,ym,zm,xM,yM,zM);
      CImgList<> points;
      CImgList<unsigned int> primitives;
      CImgList<unsigned char> colors;

      // Add ellipsoids to the 3D scene
      int X = img.dimx()/2, Y = img.dimy()/2, Z = img.dimz()/2;
      { 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); }
      { 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); }
      { 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); }

      // Add computed fibers to the 3D scene
      const CImg<> veigen = eigen.get_crop(xm,ym,zm,xM,yM,zM);
      cimglist_for(fibers,l) {
	const CImg<>& fiber = fibers[l];
	if (fiber.dimx()) insert_fiber(fiber,eigen,fiber_palette,
				       xm,ym,zm,voxw,voxh,voxd,
				       points,primitives,colors);
      }

      // Display 3D object
      CImg<unsigned char> visu(disp3d.dimx(),disp3d.dimy(),1,3,0);
      cimg_forXY(visu,x,y) visu(x,y,2) = y*32/visu.dimy();
      bool stopflag = false;
      while (!disp3d.is_closed && !stopflag) {
	visu.display_object3d(points,primitives,colors,disp3d,true,4,-1,false,800,0.2f);
	switch (disp3d.key) {
	case cimg::keyM: { // Create movie
	  std::fprintf(stderr,"\n- Movie mode.\n");
	  const unsigned int N = 256;
	  CImg<> pts = points.get_append('x');
	  CImgList<> cpoints(points);
	  CImg<> x = pts.get_shared_line(0), y = pts.get_shared_line(1), z = pts.get_shared_line(2);
	  const CImgStats sx(x,false), sy(y,false), sz(z,false);
	  const float
	    xm = (float)sx.min, xM = (float)sx.max,
	    ym = (float)sy.min, yM = (float)sy.max,
	    zm = (float)sz.min, zM = (float)sz.max,
	    ratio = 2.0f*cimg::min(visu.dimx(),visu.dimy())/(3.0f*cimg::max(xM-xm,yM-ym,zM-zm)),
	    dx = 0.5f*(xM+xm), dy = 0.5f*(yM+ym), dz = 0.5f*(zM+zm);
	  cimglist_for(points,l) {
	    cpoints(l,0) = (float)((points(l,0)-dx)*ratio);
	    cpoints(l,1) = (float)((points(l,1)-dy)*ratio);
	    cpoints(l,2) = (float)((points(l,2)-dz)*ratio);
	  }

	  for (unsigned int i=0; i<N; i++) {
	    std::fprintf(stderr,"\r- Frame %u/%u.",i,N);
	    const float alpha = (float)(i*2*cimg::PI/N);
	    const CImg<> rot = CImg<>::rotation_matrix(0,1,0,alpha)*CImg<>::rotation_matrix(1,0,0,1.30f);
	    CImgList<> rotated(cpoints);
	    cimglist_for(rotated,l) rotated[l] = rot*cpoints[l];
	    visu.fill(0).draw_object3d(visu.dimx()/2.0f,visu.dimy()/2.0f,-500.0f,rotated,primitives,colors,
				       4,false,800.0f,visu.dimx()/2.0f,visu.dimy()/2.0f,-800.0f,0.2f).display(disp3d);
	    visu.save("frame.png",i);
	  }
	  visu.fill(0);
	} break;
	default: stopflag = true;
	}
      }
      if (disp3d.is_fullscreen) disp3d.toggle_fullscreen().resize(800,600).close();
    } break;

    // Compute region statistics
    //---------------------------
    case cimg::keyR: {
      std::fprintf(stderr,"\n- Statistics computation. Select region."); std::fflush(stderr);
      coloredFA.feature_selection(s,2,disp,XYZ);
      int xm,ym,zm,xM,yM,zM;
      if (!disp.key) { xm=s[0]; ym=s[1]; zm=s[2]; xM=s[3]; yM=s[4]; zM=s[5]; }
      else { xm=ym=zm=0; xM=eigen.dimx()-1; yM=eigen.dimy()-1; zM=eigen.dimy()-1; }
      const CImg<> img = eigen.get_crop(xm,ym,zm,xM,yM,zM);
      CImgStats stats(eigen.get_shared_channel(0),false), stats2(eigen.get_shared_channel(12),false);
      std::fprintf(stderr,"\n- Mean diffusivity = %g, Mean FA = %g\n",stats.mean,stats2.mean);
    } break;

    // Track fiber bundle (single region)
    //----------------------------------
    case cimg::keyF: {
      std::fprintf(stderr,"\n- Tracking mode (single region). Select starting region.\n"); std::fflush(stderr);
      coloredFA.feature_selection(s,2,disp,XYZ);
      const unsigned int N = fibers.size;
      for (int z=s[2]; z<=s[5]; z++)
	for (int y=s[1]; y<=s[4]; y++)
	  for (int x=s[0]; x<=s[3]; x++) {
	    const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin);
	    if (fiber.width>lmin) {
	      std::fprintf(stderr,"\rFiber %u : Starting from (%d,%d,%d)\t\t",fibers.size,x,y,z);
	      fibers.insert(fiber);
	    }
	  }
      std::fprintf(stderr,"\n- %u fiber(s) added (total %u).",fibers.size-N,fibers.size);
    } break;

    // Track fiber bundle (double regions)
    //------------------------------------
    case cimg::keyG: {
      std::fprintf(stderr,"\n- Tracking mode (double region). Select starting region."); std::fflush(stderr);
      int ns[6];
      coloredFA.feature_selection(s,2,disp,XYZ);
      std::fprintf(stderr," Select ending region."); std::fflush(stderr);
      coloredFA.feature_selection(ns,2,disp,XYZ);
      const unsigned int N = fibers.size;

      // Track from start to end
      for (int z=s[2]; z<=s[5]; z++)
	for (int y=s[1]; y<=s[4]; y++)
	  for (int x=s[0]; x<=s[3]; x++) {
	    const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin);
	    if (fiber.width>lmin) {
	      bool valid_fiber = false;
	      cimg_forX(fiber,k) {
		const int fx = (int)fiber(k,0), fy = (int)fiber(k,1), fz = (int)fiber(k,2);
		if (fx>=ns[0] && fx<=ns[3] &&
		    fy>=ns[1] && fy<=ns[4] &&
		    fz>=ns[2] && fz<=ns[5]) valid_fiber = true;
	      }
	      if (valid_fiber) fibers.insert(fiber);
	    }
	  }

      // Track from end to start
      { for (int z=ns[2]; z<=ns[5]; z++)
	for (int y=ns[1]; y<=ns[4]; y++)
	  for (int x=ns[0]; x<=ns[3]; x++) {
	    const CImg<> fiber = get_fibertrack(eigen,x,y,z,lmax,dl,famin,cmin);
	    if (fiber.width>lmin) {
	      bool valid_fiber = false;
	      cimg_forX(fiber,k) {
		const int fx = (int)fiber(k,0), fy = (int)fiber(k,1), fz = (int)fiber(k,2);
		if (fx>=s[0] && fx<=s[3] &&
		    fy>=s[1] && fy<=s[4] &&
		    fz>=s[2] && fz<=s[5]) valid_fiber = true;
	      }
	      if (valid_fiber) {
		std::fprintf(stderr,"\rFiber %u : Starting from (%d,%d,%d)\t\t",fibers.size,x,y,z);
		fibers.insert(fiber);
	      }
	    }
	  }}

      std::fprintf(stderr," %u fiber(s) added (total %u).",fibers.size-N,fibers.size);
    } break;

    // Clear fiber bundle
    //-------------------
    case cimg::keyC: {
      std::fprintf(stderr,"\n- Fibers removed.");
      fibers.assign();
    } break;

    // Save fibers
    //-------------
    case cimg::keyS: {
      fibers.save("fibers.cimg");
      std::fprintf(stderr,"\n- Fibers saved.");
    } break;

    }
  }

  std::fprintf(stderr,"\n- Exit.\n\n\n");
  return 0;
}