1 /*
2 #
3 # File : mcf_levelsets2d.cpp
4 # ( C++ source file )
5 #
6 # Description : Implementation of the Mean Curvature Flow on a 2D curve,
7 # using the framework of Level Sets.
8 # This file is a part of the CImg Library project.
9 # ( http://cimg.eu )
10 #
11 # Copyright : David Tschumperlé
12 # ( http://tschumperle.users.greyc.fr/ )
13 #
14 # License : CeCILL v2.0
15 # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
16 #
17 # This software is governed by the CeCILL license under French law and
18 # abiding by the rules of distribution of free software. You can use,
19 # modify and/ or redistribute the software under the terms of the CeCILL
20 # license as circulated by CEA, CNRS and INRIA at the following URL
21 # "http://www.cecill.info".
22 #
23 # As a counterpart to the access to the source code and rights to copy,
24 # modify and redistribute granted by the license, users are provided only
25 # with a limited warranty and the software's author, the holder of the
26 # economic rights, and the successive licensors have only limited
27 # liability.
28 #
29 # In this respect, the user's attention is drawn to the risks associated
30 # with loading, using, modifying and/or developing or reproducing the
31 # software by the user in light of its specific status of free software,
32 # that may mean that it is complicated to manipulate, and that also
33 # therefore means that it is reserved for developers and experienced
34 # professionals having in-depth computer knowledge. Users are therefore
35 # encouraged to load and test the software's suitability as regards their
36 # requirements in conditions enabling the security of their systems and/or
37 # data to be ensured and, more generally, to use and operate it in the
38 # same conditions as regards security.
39 #
40 # The fact that you are presently reading this means that you have had
41 # knowledge of the CeCILL license and that you accept its terms.
42 #
43 */
44
45 #include "CImg.h"
46 using namespace cimg_library;
47 #undef min
48 #undef max
49
50 // Retrieve the curve corresponding to the zero level set of the distance function.
51 template<typename T>
get_level0(const CImg<T> & img)52 CImg<unsigned char> get_level0(const CImg<T>& img) {
53 CImg<unsigned char> dest(img);
54 CImg_2x2(I,T); Inn = 0;
55 cimg_for2x2(img,x,y,0,0,I,T) if (Icc*Inc<0 || Icc*Icn<0) dest(x,y) = 255; else dest(x,y) = Icc<0?100:0;
56 return dest;
57 }
58
59 /*--------------------
60
61 Main procedure
62
63 ----------------------*/
main(int argc,char ** argv)64 int main(int argc,char **argv) {
65 cimg_usage("Perform a Mean Curvature Flow on closed curves, using Level Sets");
66 const float dt = cimg_option("-dt",0.8f,"PDE time step");
67 const unsigned int nb_iterations = cimg_option("-iter",10000,"Number of iterations");
68
69 // Create a user-defined closed curve.
70 CImg<unsigned char> curve(256,256,1,2,0);
71 unsigned char col1[] = {0,255}, col2[] = {200,255}, col3[] = {255,255};
72 curve.draw_grid(20,20,0,0,false,false,col1,0.4f,0xCCCCCCCC,0xCCCCCCCC).
73 draw_text(5,5,"Please draw your curve\nin this window\n(Use your mouse)",col1);
74 CImgDisplay disp(curve,"Mean curvature flow",0);
75 int xo = -1, yo = -1, x0 = -1, y0 = -1, x1 = -1, y1 = -1;
76 while (!disp.is_closed() && (x0<0 || disp.button())) {
77 if (disp.button() && disp.mouse_x()>=0 && disp.mouse_y()>=0) {
78 if (x0<0) { xo = x0 = disp.mouse_x(); yo = y0 = disp.mouse_y(); } else {
79 x1 = disp.mouse_x(); y1 = disp.mouse_y();
80 curve.draw_line(x0,y0,x1,y1,col2).display(disp);
81 x0 = x1; y0 = y1;
82 }
83 }
84 disp.wait();
85 if (disp.is_resized()) disp.resize(disp);
86 }
87 curve.draw_line(x1,y1,xo,yo,col2).channel(0).draw_fill(0,0,col3);
88 CImg<> img = CImg<>(curve.get_shared_channel(0)).normalize(-1,1);
89
90 // Perform the "Mean Curvature Flow".
91 img.distance_eikonal(10);
92 CImg_3x3(I,float);
93 for (unsigned int iteration = 0; iteration<nb_iterations && !disp.is_closed() &&
94 !disp.is_keyQ() && !disp.is_keyESC(); ++iteration) {
95 CImg<float> velocity(img.width(),img.height(),img.depth(),img.spectrum());
96 float *ptrd = velocity.data(), veloc_max = 0;
97 cimg_for3x3(img,x,y,0,0,I,float) {
98 const float
99 ix = (Inc - Ipc)/2,
100 iy = (Icn - Icp)/2,
101 ixx = Inc + Ipc - 2*Icc,
102 iyy = Icn + Icp - 2*Icc,
103 ixy = (Ipp + Inn - Inp - Ipn)/4,
104 ngrad = ix*ix + iy*iy,
105 iee = (ngrad>1e-5)?((iy*iy*ixx - 2*ix*iy*ixy + ix*ix*iyy)/ngrad):0;
106 *(ptrd++) = iee;
107 if (iee>veloc_max) veloc_max = iee; else if (-iee>veloc_max) veloc_max = -iee;
108 }
109 if (veloc_max>0) img+=(velocity*=dt/veloc_max);
110 if (!(iteration%10)) {
111 get_level0(img).resize(disp.width(),disp.height()).
112 draw_grid(20,20,0,0,false,false,col3,0.4f,0xCCCCCCCC,0xCCCCCCCC).
113 draw_text(5,5,"Iteration %d",col3,0,1,13,iteration).display(disp);
114 }
115 if (!(iteration%60)) img.distance_eikonal(1,3);
116 if (disp.is_resized()) disp.resize();
117 }
118
119 return 0;
120 }
121