1.1 diff -r 5edfbd3e7a46 -r 1204ebf9340d PTdecode/CImg-1.3.0/examples/mcf_levelsets3d.cpp 1.2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 +++ b/PTdecode/CImg-1.3.0/examples/mcf_levelsets3d.cpp Mon Aug 03 14:09:20 2009 +0100 1.4 @@ -0,0 +1,183 @@ 1.5 +/* 1.6 + # 1.7 + # File : mcf_levelsets3d.cpp 1.8 + # ( C++ source file ) 1.9 + # 1.10 + # Description : Implementation of the Mean Curvature Flow on Surfaces 1.11 + # using the framework of Level Sets 3D. 1.12 + # This file is a part of the CImg Library project. 1.13 + # ( http://cimg.sourceforge.net ) 1.14 + # 1.15 + # Copyright : David Tschumperle 1.16 + # ( http://www.greyc.ensicaen.fr/~dtschump/ ) 1.17 + # 1.18 + # License : CeCILL v2.0 1.19 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.20 + # 1.21 + # This software is governed by the CeCILL license under French law and 1.22 + # abiding by the rules of distribution of free software. You can use, 1.23 + # modify and/ or redistribute the software under the terms of the CeCILL 1.24 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.25 + # "http://www.cecill.info". 1.26 + # 1.27 + # As a counterpart to the access to the source code and rights to copy, 1.28 + # modify and redistribute granted by the license, users are provided only 1.29 + # with a limited warranty and the software's author, the holder of the 1.30 + # economic rights, and the successive licensors have only limited 1.31 + # liability. 1.32 + # 1.33 + # In this respect, the user's attention is drawn to the risks associated 1.34 + # with loading, using, modifying and/or developing or reproducing the 1.35 + # software by the user in light of its specific status of free software, 1.36 + # that may mean that it is complicated to manipulate, and that also 1.37 + # therefore means that it is reserved for developers and experienced 1.38 + # professionals having in-depth computer knowledge. Users are therefore 1.39 + # encouraged to load and test the software's suitability as regards their 1.40 + # requirements in conditions enabling the security of their systems and/or 1.41 + # data to be ensured and, more generally, to use and operate it in the 1.42 + # same conditions as regards security. 1.43 + # 1.44 + # The fact that you are presently reading this means that you have had 1.45 + # knowledge of the CeCILL license and that you accept its terms. 1.46 + # 1.47 +*/ 1.48 + 1.49 +#include "CImg.h" 1.50 +using namespace cimg_library; 1.51 + 1.52 +// The lines below are necessary when using a non-standard compiler as visualcpp6. 1.53 +#ifdef cimg_use_visualcpp6 1.54 +#define std 1.55 +#endif 1.56 +#ifdef min 1.57 +#undef min 1.58 +#undef max 1.59 +#endif 1.60 + 1.61 +// Apply the Mean curvature flow PDE 1.62 +//----------------------------------- 1.63 +template<typename T> CImg<T>& mcf_PDE(CImg<T>& img, const unsigned int nb_iter, 1.64 + const float dt=0.25f, const float narrow=4.0f) { 1.65 + CImg<T> veloc(img.dimx(),img.dimy(),img.dimz(),img.dimv()); 1.66 + CImg_3x3x3(I,float); 1.67 + for (unsigned int iter=0; iter<nb_iter; iter++) { 1.68 + cimg_for3x3x3(img,x,y,z,0,I) if (cimg::abs(Iccc)<narrow) { 1.69 + const float 1.70 + ix = 0.5f*(Incc-Ipcc), 1.71 + iy = 0.5f*(Icnc-Icpc), 1.72 + iz = 0.5f*(Iccn-Iccp), 1.73 + norm = (float)std::sqrt(1e-5f+ix*ix+iy*iy+iz*iz), 1.74 + ixx = Incc+Ipcc-2*Iccc, 1.75 + ixy = 0.25f*(Ippc+Innc-Inpc-Ipnc), 1.76 + ixz = 0.25f*(Ipcp+Incn-Incp-Ipcn), 1.77 + iyy = Icnc+Icpc-2*Iccc, 1.78 + iyz = 0.25f*(Icpp+Icnn-Icnp-Icpn), 1.79 + izz = Iccn+Iccp-2*Iccc, 1.80 + a = ix/norm, 1.81 + b = iy/norm, 1.82 + c = iz/norm, 1.83 + inn = a*a*ixx + b*b*iyy + c*c*izz + 2*a*b*ixy + 2*a*c*ixz + 2*b*c*iyz; 1.84 + veloc(x,y,z) = ixx+iyy+izz-inn; 1.85 + } else veloc(x,y,z) = 0; 1.86 + float m, M = veloc.maxmin(m); 1.87 + const double xdt = dt/cimg::max(cimg::abs(m),cimg::abs(M)); 1.88 + img+=xdt*veloc; 1.89 + } 1.90 + return img; 1.91 +} 1.92 + 1.93 +// Main procedure 1.94 +//---------------- 1.95 +int main(int argc,char **argv) { 1.96 + cimg_usage("Mean curvature flow of a surface, using 3D level sets"); 1.97 + const char *file_i = cimg_option("-i",(char*)0,"Input image"); 1.98 + const float dt = cimg_option("-dt",0.05f,"PDE Time step"); 1.99 + const float narrow = cimg_option("-band",5.0f,"Size of the narrow band"); 1.100 + const bool both = cimg_option("-both",false,"Show both evolving and initial surface"); 1.101 + 1.102 + // Define the signed distance map of the initial surface 1.103 + CImg<> img; 1.104 + if (file_i) { 1.105 + const float sigma = cimg_option("-sigma",1.2f,"Segmentation regularity"); 1.106 + const float alpha = cimg_option("-alpha",5.0f,"Region growing tolerance"); 1.107 + img.load(file_i).channel(0); 1.108 + CImg<int> s; 1.109 + CImgDisplay disp(img,"Please select a starting point"); 1.110 + while (!s || s[0]<0) s = img.get_select(0,disp); 1.111 + CImg<> region; 1.112 + float tmp[1] = { 0 }; 1.113 + img.draw_fill(s[0],s[1],s[2],tmp,1,region,alpha); 1.114 + ((img = region.normalize(-1,1))*=-1).blur(sigma); 1.115 + 1.116 + } 1.117 + else { // Create synthetic implicit function 1.118 + img.assign(60,60,60); 1.119 + const float exte[1]={1}, inte[1]={-1}; 1.120 + img.fill(*exte).draw_rectangle(15,15,15,45,45,45,inte).draw_rectangle(25,25,0,35,35,img.dimz()-1,exte). 1.121 + draw_rectangle(0,25,25,img.dimx()-1,35,35,exte).draw_rectangle(25,0,25,35,img.dimy()-1,35,exte); 1.122 + } 1.123 + img.distance_hamilton(10,0,0.1f); 1.124 + 1.125 + // Compute corresponding surface triangularization by the marching cube algorithm (isovalue 0) 1.126 + CImg<> points0; 1.127 + CImgList<unsigned int> faces0; 1.128 + if (both) points0 = img.get_isovalue3d(faces0,0); 1.129 + const CImgList<unsigned char> colors0(faces0.size,CImg<unsigned char>::vector(100,200,255)); 1.130 + const CImgList<> opacities0(faces0.size,1,1,1,1,0.2f); 1.131 + 1.132 + // Perform MCF evolution 1.133 + CImgDisplay disp(256,256,"",1), disp3d(512,512,"",0); 1.134 + float alpha = 0, beta = 0; 1.135 + for (unsigned int iter=0; !disp.is_closed && !disp3d.is_closed && !disp.is_keyESC && !disp3d.is_keyESC && 1.136 + !disp.is_keyQ && !disp3d.is_keyQ; iter++) { 1.137 + disp.set_title("3D implicit Function (iter. %u)",iter); 1.138 + disp3d.set_title("Mean curvature flow 3D - Isosurface (iter. %u)",iter); 1.139 + 1.140 + // Apply PDE on the distance function 1.141 + mcf_PDE(img,1,dt,narrow); // Do one iteration of mean curvature flow 1.142 + if (!(iter%10)) img.distance_hamilton(1,narrow,0.5f); // Every 10 steps, do one iteration of distance function re-initialization 1.143 + 1.144 + // Compute surface triangularization by the marching cube algorithm (isovalue 0) 1.145 + CImgList<unsigned int> faces; 1.146 + CImg<> points = img.get_isovalue3d(faces,0); 1.147 + CImgList<unsigned char> colors(faces.size,CImg<unsigned char>::vector(200,128,100)); 1.148 + CImgList<> opacities(faces.size,CImg<>::vector(1.0f)); 1.149 + const float fact = 3*cimg::max(disp3d.dimx(),disp3d.dimy())/(4.0f*cimg::max(img.dimx(),img.dimy())); 1.150 + 1.151 + // Append initial object if necessary. 1.152 + if (both) { 1.153 + points.append_object3d(faces,points0,faces0); 1.154 + colors.insert(colors0); 1.155 + opacities.insert(opacities0); 1.156 + } 1.157 + 1.158 + // center and rescale the objects 1.159 + cimg_forX(points,l) { 1.160 + points(l,0)=(points(l,0)-img.dimx()/2)*fact; 1.161 + points(l,1)=(points(l,1)-img.dimy()/2)*fact; 1.162 + points(l,2)=(points(l,2)-img.dimz()/2)*fact; 1.163 + } 1.164 + 1.165 + // Display 3D object on the display window. 1.166 + CImg<unsigned char> visu(disp3d.dimx(),disp3d.dimy(),1,3,0); 1.167 + const CImg<> rot = CImg<>::rotation_matrix(1,0,0,(beta+=0.01f))*CImg<>::rotation_matrix(0,1,1,(alpha+=0.05f)); 1.168 + if (points.size()) { 1.169 + visu.draw_object3d(visu.dimx()/2.0f,visu.dimy()/2.0f,0.0f, 1.170 + rot*points,faces,colors,opacities,3, 1.171 + false,500.0,0.0f,0.0f,-8000.0f).display(disp3d); 1.172 + } else visu.fill(0).display(disp3d); 1.173 + img.display(disp.wait(20)); 1.174 + 1.175 + if ((disp3d.button || disp3d.key) && points.size()) { 1.176 + unsigned char white[3]={ 255,255,255 }; 1.177 + visu.fill(0).draw_text(10,10,"Time stopped, press any key to start again",white). 1.178 + display_object3d(disp3d,points,faces,colors,opacities,true,4,3,false,500,0.4f,0.3f); 1.179 + disp3d.key = 0; 1.180 + } 1.181 + if (disp.is_resized) disp.resize(false); 1.182 + if (disp3d.is_resized) disp3d.resize(false); 1.183 + } 1.184 + 1.185 + // Exit 1.186 + return 0; 1.187 +}