1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/PTdecode/CImg-1.3.0/examples/mcf_levelsets.cpp Mon Aug 03 14:09:20 2009 +0100 1.3 @@ -0,0 +1,124 @@ 1.4 +/* 1.5 + # 1.6 + # File : mcf_levelsets.cpp 1.7 + # ( C++ source file ) 1.8 + # 1.9 + # Description : Implementation of the Mean Curvature Flow (classical 2d curve evolution), 1.10 + # using the framework of Level Sets. 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 +// get_level0() : Retrieve the curve corresponding to the zero level set of the distance function 1.61 +//------------- 1.62 +CImg<unsigned char> get_level0(const CImg<>& img) { 1.63 + CImg<unsigned char> dest(img); 1.64 + CImg_2x2(I,float); Inn = 0; 1.65 + cimg_for2x2(img,x,y,0,0,I) if (Icc*Inc<0 || Icc*Icn<0) dest(x,y) = 255; else dest(x,y) = Icc<0?100:0; 1.66 + return dest; 1.67 +} 1.68 + 1.69 +//----------------- 1.70 +// Main procedure 1.71 +//----------------- 1.72 +int main(int argc,char **argv) { 1.73 + cimg_usage("Perform a Mean Curvature Flow on closed curves, using Level Sets"); 1.74 + const float dt = cimg_option("-dt",0.8f,"PDE time step"); 1.75 + const unsigned int nb_iter = cimg_option("-iter",10000,"Number of iterations"); 1.76 + 1.77 + // Create a user-defined closed curve 1.78 + CImg<unsigned char> curve(256,256,1,2,0); 1.79 + unsigned char col1[2]={0,255}, col2[2]={200,255}, col3[2]={255,255}; 1.80 + curve.draw_grid(20,20,0,0,false,false,col1,0.4f,0xCCCCCCCC,0xCCCCCCCC). 1.81 + draw_text(5,5,"Please draw your curve\nin this window\n(Use your mouse)",col1); 1.82 + CImgDisplay disp(curve,"Mean curvature flow",0); 1.83 + int xo=-1,yo=-1,x0=-1,y0=-1,x1=-1,y1=-1; 1.84 + while (!disp.is_closed && (x0<0 || disp.button)) { 1.85 + if (disp.button && disp.mouse_x>=0 && disp.mouse_y>=0) { 1.86 + if (x0<0) { xo = x0 = disp.mouse_x; yo = y0 = disp.mouse_y; } else { 1.87 + x1 = disp.mouse_x; y1 = disp.mouse_y; 1.88 + curve.draw_line(x0,y0,x1,y1,col2).display(disp); 1.89 + x0 = x1; y0 = y1; 1.90 + } 1.91 + } 1.92 + disp.wait(); 1.93 + if (disp.is_resized) disp.resize(disp); 1.94 + } 1.95 + curve.draw_line(x1,y1,xo,yo,col2).channel(0).draw_fill(0,0,col3); 1.96 + CImg<> img = CImg<>(curve.get_shared_channel(0)).normalize(-1,1); 1.97 + 1.98 + // Perform the "Mean Curvature Flow" 1.99 + img.distance_hamilton(10); 1.100 + CImg_3x3(I,float); 1.101 + for (unsigned int iter=0; iter<nb_iter && !disp.is_closed && !disp.is_keyQ; iter++) { 1.102 + CImg<> veloc(img.dimx(),img.dimy(),img.dimz(),img.dimv()); 1.103 + cimg_for3x3(img,x,y,0,0,I) { 1.104 + const float 1.105 + ix = 0.5f*(Inc-Ipc), 1.106 + iy = 0.5f*(Icn-Icp), 1.107 + ixx = Inc+Ipc-2*Icc, 1.108 + iyy = Icn+Icp-2*Icc, 1.109 + ixy = 0.25f*(Ipp+Inn-Inp-Ipn), 1.110 + ngrad = ix*ix+iy*iy, 1.111 + iee = (ngrad>1e-5)?(( iy*iy*ixx - 2*ix*iy*ixy + ix*ix*iyy )/ngrad):0; 1.112 + veloc(x,y) = iee; 1.113 + } 1.114 + float m, M = veloc.maxmin(m); 1.115 + const double xdt = dt/cimg::max(cimg::abs(m),cimg::abs(M)); 1.116 + img+=xdt*veloc; 1.117 + if (!(iter%10)) { 1.118 + get_level0(img).resize(disp.dimx(),disp.dimy()).draw_grid(20,20,0,0,false,false,col3,0.4f,0xCCCCCCCC,0xCCCCCCCC). 1.119 + draw_text(5,5,"Iteration %d",col3,0,1,11,iter).display(disp); 1.120 + } 1.121 + if (!(iter%30)) img.distance_hamilton(1,3); 1.122 + if (disp.is_resized) disp.resize(); 1.123 + } 1.124 + 1.125 + // End of program 1.126 + return 0; 1.127 +}