Jerome Piovano

Research in Medical Imaging and Computer Vision

2D_LocalGaussian.C

Local Gaussian Region Segmentation of a 2D Image

00001 #include <Segmentation.h>
00002 #include <SegmentationModule/RegularizationModule.h>
00003 #include <SegmentationModule/LocalGaussianRegionModule.h>
00004 #include <SegmentationModule/GaussianRegionModule.h>
00005 
00006 #include <SegmentationModule/SoftShapeModule.h>
00007 #include <SegmentationModule/BalloonModule.h>
00008 
00009 #include <cmath>
00010 #include <dirent.h>
00011 #include <DericheFilter.h>
00012 
00013         // module assembler qui prends en compte le module 0 + le premier module une valeur non nulle
00014         // ==> penser a ajouter les modules suivant une priorite 
00015         template <class _Modules, class _LevelSet, class _Index>
00016         float meancurv_depending_on_localstats(_Modules & modules, const _LevelSet & phi, const _Index & ind, int& indice)
00017         {
00018                 float value_regul  = modules[0]->value(phi, ind);
00019                 float weight_regul = modules[0]->weight(ind);
00020 
00021                 float value  = modules[1]->value(phi, ind);
00022                 float weight = modules[1]->weight(ind);
00023 
00024                 float result;
00025                 if (fabs(value) < 0.05) result = value_regul*weight_regul*3 + value*weight;
00026                 else                    result = value_regul*weight_regul   + value*weight;
00027                 
00028                 indice = 1;
00029                 return result;
00030         }
00031 
00034 int main(int argc, char* argv[])
00035 {
00036         using namespace Images;
00037         using namespace levelset;
00038         using namespace segmentation;   
00039 
00040 #if 1 // 2D
00041         // image to segment
00042         BaseImage<2, float> image;
00043         //std::ifstream f("image/tiger.inr");
00044         std::ifstream f("image/last.inr");
00045         //std::ifstream f("image/noise_plus.inr");
00046         //std::ifstream f("image/synth3/bruit/000475.inr");
00047         //std::ifstream f("image/presk_flou.inr");
00048         //std::ifstream f("image/line_presk_flou.inr");
00049         //std::ifstream f("image/head.inr");
00050         //std::ifstream f("image/head2.inr");
00051         //std::ifstream f("image/head3.inr");
00052         //std::ifstream f("image/head4.inr");
00053         //std::ifstream f("image/triangle.inr");
00054         //std::ifstream f("image/biais.inr");
00055         //std::ifstream f("image/nois_var.inr");
00056         //std::ifstream f("image/biais_gaussian.inr");
00057         //std::ifstream f("image/blur_gaussian.inr");
00058         //std::ifstream f("image/cercle_biais.inr");
00059         f >> format("inrimage") >> image;
00060 
00061         Globals::normalize(image, 0.0f, 250.0f);
00062 
00063         // initialization of the segmentations objects
00064         FastMarchingInitializer<2,float> init;                               // Fast Marching Initializer
00065         LevelSet<2,float>                phi(image.shape(), (double) 2.0);   // Level Set
00066         Segmentation<2,float>            segm(&phi, &init, &image);          // Segmentation object
00067         Segmentation<2,float>            segm_global(&phi, &init, &image);   // Segmentation object
00068 
00069         //segm.set_module_assembler(&meancurv_depending_on_localstats);
00070 
00071         // Initalization of the level set
00072         phi.init_small_spheres(10, 30,init);    // parameters : radius, windows
00073 /*
00074         LevelSet<2,float>::Index tl(100, 20);
00075         LevelSet<2,float>::Index br(200, 130);
00076         phi.init_rectangle(tl, br, init);
00077         LevelSet<2,float>::Index tl2(100, 170);
00078         LevelSet<2,float>::Index br2(200, 280);
00079         phi.init_rectangle(tl2, br2, init);
00080 */
00081         //LevelSet<2,float>::Index tl(20, 100);
00082         //LevelSet<2,float>::Index br(280, 200);
00083 
00084 
00085         //initialization of the modules
00086         RegularizationModule<2,float>             regul(0.8f);                        // parameters : weight
00087         LocalGaussianRegionModule<2,float,float>  local_gaussian( &image, 1.f, 2);    // parameters : image, weight, smooth
00088         GaussianRegionModule<2,float,float>       gaussian( &image, 1.f);             // parameters : image, weight
00089 
00090         //local_gaussian.use_variance(false);
00091         gaussian.use_variance(true);
00092 
00093         segm.add_module(regul);
00094         //segm.add_module(local_gaussian);
00095 
00096         segm_global.add_module(regul);
00097         segm_global.add_module(gaussian);
00098 
00099         // launch the segmentation
00100 /*
00101         Segmentation<2,float>::RGBPixel neg(40,90,220);
00102         Segmentation<2,float>::RGBPixel pos(249,249,0);
00103 */
00104         Segmentation<2,float>::RGBPixel neg(215,152,0); 
00105         Segmentation<2,float>::RGBPixel pos(255,252,0); 
00106 
00107         timeval _time;
00108         Globals::initChrono(_time);
00109         std::ostringstream name;
00110  
00111         name.str("");  
00112         name << "mkdir CVPR/";
00113         system (name.str().c_str()); 
00114 
00115 #if 1
00116 
00117         name.str("");  
00118         name << "mkdir CVPR/solution_coherence";
00119         system (name.str().c_str());
00120 
00121         for(int sigma = 2 ; sigma <= 64 ; sigma*=2 ){
00122                 sigma = 16;
00123 
00124                 name.str("");  
00125                 name << "mkdir CVPR/solution_coherence/"  << sigma << "/";
00126                 system (name.str().c_str());
00127                 name.str("");
00128                 name << "CVPR/solution_coherence/"  << sigma << "/";    
00129                 phi = 0;
00130                 //phi.init_rectangle(tl, br, init);
00131                 phi.init_small_spheres(10, 30,init);
00132                 LocalGaussianRegionModule<2,float,float>  local_gaussian( &image, 0.4f , sigma, false);    // parameters : image, weight, smooth
00133                 segm.add_module(local_gaussian); 
00134                 segm.segment(300, name.str().c_str() , neg, pos);
00135                 segm.remove_module(local_gaussian);
00136                 name.str("");
00137                 name << "convert_cimg.sh png CVPR/solution_coherence/"  << sigma <<  "/var/*ppm";
00138                 system (name.str().c_str());
00139 
00140                 return 0;
00141         }
00142 
00143         name.str("");  
00144         name << "mkdir CVPR/drawback_16/global/";
00145         system (name.str().c_str());
00146         name.str("");
00147         name << "CVPR/drawback_16/global/";     
00148         phi = 0;
00149         //phi.init_rectangle(tl, br, init);
00150         GaussianRegionModule<2,float,float>  gaussian1( &image, 2.f, true);    // parameters : image, weight, smooth
00151         segm.add_module(gaussian1);
00152         segm.segment(8000, name.str().c_str() , neg, pos);
00153         segm.remove_module(gaussian1);
00154         name.str("");
00155         name << "convert_cimg.sh png CVPR/drawback_16/global/var/*ppm";
00156         system (name.str().c_str());
00157 #endif
00158 
00159 #if 0 // 2eme image
00160 
00161         BaseImage<2, float> image2;
00162         std::ifstream f2("image/blur_gaussian.inr");
00163         f2 >> format("inrimage") >> image2;
00164         Globals::normalize(image2, 0.0f, 250.0f);
00165         segm.set_data_term(image2);
00166 
00167         name.str("");  
00168         name << "mkdir CVPR/blur_gaussian";
00169         system (name.str().c_str());
00170 
00171         for(int sigma = 2 ; sigma <= 64 ; sigma*=2 ){
00172 
00173                 name.str("");  
00174                 name << "mkdir CVPR/blur_gaussian/"  << sigma << "/";
00175                 system (name.str().c_str());
00176                 name.str("");
00177                 name << "CVPR/blur_gaussian/"  << sigma << "/"; 
00178                 phi = 0;
00179                 phi.init_rectangle(tl, br, init);phi.init_rectangle(tl2, br2, init);
00180                 LocalGaussianRegionModule<2,float,float>  local_gaussian( &image2, 1.f, sigma, false);    // parameters : image, weight, smooth
00181                 segm.add_module(local_gaussian);
00182                 segm.segment(8000, name.str().c_str() , neg, pos);
00183                 segm.remove_module(local_gaussian);
00184                 name.str("");
00185                 name << "convert_cimg.sh png CVPR/blur_gaussian/"  << sigma <<  "/var/*ppm";
00186                 system (name.str().c_str()); 
00187 
00188         }
00189         {
00190                 name.str("");  
00191                 name << "mkdir CVPR/blur_gaussian/global/";
00192                 system (name.str().c_str());
00193                 name.str("");
00194                 name << "CVPR/blur_gaussian/global/";   
00195                 phi = 0;
00196                 phi.init_rectangle(tl, br, init);phi.init_rectangle(tl2, br2, init);
00197                 GaussianRegionModule<2,float,float>  gaussian1( &image2, 2.f, true);    // parameters : image, weight, smooth
00198                 segm.add_module(gaussian1);
00199                 segm.segment(8000, name.str().c_str() , neg, pos);
00200                 segm.remove_module(gaussian1);
00201                 name.str("");
00202                 name << "convert_cimg.sh png CVPR/blur_gaussian/global/var/*ppm";
00203                 system (name.str().c_str());
00204         }
00205 #endif
00206 
00207 #if 0 // 3eme image
00208 
00209         BaseImage<2, float> image3;
00210         std::ifstream f3("image/mix_gaussian.inr");
00211         f3 >> format("inrimage") >> image3;
00212         Globals::normalize(image3, 0.0f, 250.0f);
00213         segm.set_data_term(image3);
00214 
00215         name.str("");  
00216         name << "mkdir CVPR/mix_gaussian";
00217         system (name.str().c_str());
00218 
00219         for(int sigma = 2 ; sigma <= 64 ; sigma*=2 ){
00220 
00221                 name.str("");  
00222                 name << "mkdir CVPR/mix_gaussian/"  << sigma << "/";
00223                 system (name.str().c_str());
00224                 name.str("");
00225                 name << "CVPR/mix_gaussian/"  << sigma << "/";  
00226                 phi = 0;
00227                 phi.init_rectangle(tl, br, init);phi.init_rectangle(tl2, br2, init);
00228                 LocalGaussianRegionModule<2,float,float>  local_gaussian( &image3, 1.f, sigma, false);    // parameters : image, weight, smooth
00229                 segm.add_module(local_gaussian);
00230                 segm.segment(8000, name.str().c_str() , neg, pos);
00231                 segm.remove_module(local_gaussian);
00232                 name.str("");
00233                 name << "convert_cimg.sh png CVPR/mix_gaussian/"  << sigma <<  "/var/*ppm";
00234                 system (name.str().c_str()); 
00235         }
00236         {
00237                 name.str("");  
00238                 name << "mkdir CVPR/mix_gaussian/global/";
00239                 system (name.str().c_str());
00240                 name.str("");
00241                 name << "CVPR/mix_gaussian/global/";    
00242                 phi = 0;
00243                 phi.init_rectangle(tl, br, init);phi.init_rectangle(tl2, br2, init);
00244                 GaussianRegionModule<2,float,float>  gaussian1( &image3, 1.f, true);    // parameters : image, weight, smooth
00245                 segm.add_module(gaussian1);
00246                 segm.segment(8000, name.str().c_str() , neg, pos);
00247                 segm.remove_module(gaussian1);
00248                 name.str("");
00249                 name << "convert_cimg.sh png CVPR/mix_gaussian/global/var/*ppm";
00250                 system (name.str().c_str());
00251         }
00252 #endif
00253 
00254 
00255         Globals::printChrono("\n        => LocalStats : ", _time);
00256         Globals::initChrono(_time); 
00257 
00258     #if 0 // + tests sigmas successifs init small sphere
00259 
00260         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00261         local_gaussian.set_sigma(5); 
00262         segm.segment(500, "result_synth_variant_blur/noise_var/var/5/localstats_", neg, pos);
00263         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/var/5/localstats_*ppm");
00264         
00265         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00266         local_gaussian.set_sigma(10);
00267         segm.segment(500, "result_synth_variant_blur/noise_var/var/10/localstats_", neg, pos);
00268         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/var/10/localstats_*ppm");
00269         
00270         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00271         local_gaussian.set_sigma(20);
00272         segm.segment(500, "result_synth_variant_blur/noise_var/var/20/localstats_", neg, pos);
00273         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/var/20/localstats_*ppm");
00274         
00275         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00276         local_gaussian.set_sigma(40);
00277         segm.segment(500, "result_synth_variant_blur/noise_var/var/40/localstats_", neg, pos);
00278         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/var/40/localstats_*ppm");
00279         
00280         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00281         local_gaussian.set_sigma(80);
00282         segm.segment(500, "result_synth_variant_blur/noise_var/var/80/localstats_", neg, pos);
00283         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/var/80/localstats_*ppm");
00284         
00285         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00286         segm_global.segment(500, "result_synth_variant_blur/noise_var/var/global/localstats_", neg, pos);
00287         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/var/global/localstats_*ppm");
00288 
00289     #endif
00290 
00291     #if 0 // + tests sigmas successifs init sphere
00292 
00293         LevelSet<2,float>::Index c(100, 100);
00294 
00295         phi.init_sphere(c, 90);
00296         local_gaussian.set_sigma(5);
00297         segm.segment(500, "test_variant_contour/var05_", neg, pos);
00298         system("./convert_cimg.sh png test_variant_contour/var05_*ppm");
00299         
00300         phi.init_sphere(c, 90);
00301         local_gaussian.set_sigma(10);
00302         segm.segment(500, "test_variant_contour/var10_", neg, pos);
00303         system("./convert_cimg.sh png test_variant_contour/var10_*ppm");
00304         
00305         phi.init_sphere(c, 90);
00306         local_gaussian.set_sigma(20);
00307         segm.segment(500, "test_variant_contour/var20_", neg, pos);
00308         system("./convert_cimg.sh png test_variant_contour/var20_*ppm");
00309         
00310         phi.init_sphere(c, 90);
00311         local_gaussian.set_sigma(40);
00312         segm.segment(500, "test_variant_contour/var40_", neg, pos);
00313         system("./convert_cimg.sh png test_variant_contour/var40_*ppm");
00314         
00315         phi.init_sphere(c, 90);
00316         local_gaussian.set_sigma(80);
00317         segm.segment(500, "test_variant_contour/var80_", neg, pos);
00318         system("./convert_cimg.sh png test_variant_contour/var80_*ppm");
00319 
00320         phi.init_sphere(c, 90);
00321         segm_global.segment(500, "test_variant_contour/global_", neg, pos);
00322         system("./convert_cimg.sh png test_variant_contour/global_*ppm");
00323 
00324     #endif
00325 
00326     #if 0 // second test autre image 2D
00327 
00328         BaseImage<2, float> image2;
00329         std::ifstream f2("image/head4.inr");
00330         f2 >> format("inrimage") >> image2;
00331 
00332         phi.resize(image2.shape());
00333         segm.set_data_term(image2);
00334         segm_global.set_data_term(image2);
00335         local_gaussian.set_data_term(image2);
00336         gaussian.set_data_term(image2);
00337 
00338         local_gaussian.use_variance(false);
00339         gaussian.use_variance(false);
00340 
00341         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00342         local_gaussian.set_sigma(2);
00343         segm_global.segment(2000, "result_synth_variant_blur/noise_var/no_var/2/localstats_", neg, pos);
00344         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/2/localstats_*ppm");
00345 
00346         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00347         local_gaussian.set_sigma(5);
00348         segm.segment(2000, "result_synth_variant_blur/noise_var/no_var/5/localstats_", neg, pos);
00349         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/5/localstats_*ppm");
00350         
00351         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00352         local_gaussian.set_sigma(10);
00353         segm.segment(2000, "result_synth_variant_blur/noise_var/no_var/10/localstats_", neg, pos);
00354         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/10/localstats_*ppm");
00355          
00356         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00357         local_gaussian.set_sigma(20);
00358         segm.segment(2000, "result_synth_variant_blur/noise_var/no_var/20/localstats_", neg, pos);
00359         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/20/localstats_*ppm");
00360         
00361         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00362         local_gaussian.set_sigma(40);
00363         segm.segment(2000, "result_synth_variant_blur/noise_var/no_var/40/localstats_", neg, pos);
00364         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/40/localstats_*ppm");
00365         
00366         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00367         local_gaussian.set_sigma(80);
00368         segm.segment(2000, "result_synth_variant_blur/noise_var/no_var/80/localstats_", neg, pos);
00369         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/80/localstats_*ppm");
00370         
00371         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00372         segm_global.segment(2000, "result_synth_variant_blur/noise_var/no_var/global/localstats_", neg, pos);
00373         system("./convert_cimg.sh png result_synth_variant_blur/noise_var/no_var/global/localstats_*ppm");
00374         
00375     #endif 
00376         std::cout << "OK" << std::endl;
00377 #endif 
00378 
00379 #if 0  // 3D
00380 
00381         //TODO traiter les images en mettant unsigned fixed a la place de fixed dans l entete inr
00382         
00383         // image to segment
00384         BaseImage<3, unsigned short> image;
00385         //BaseImage<3, unsigned char> image;
00386         //std::ifstream f("image/T1_512_msp.inr");
00387         //std::ifstream f("image/brain.inr");
00388         std::ifstream f("image/mt_diff_7_marche.inr");
00389 
00390         f >> format("inrimage") >> image;
00391 
00392         std::cout << image.shape() << std::endl; 
00393 
00394         // initialization of the segmentations objects
00395         FastMarchingInitializer<3, float> init;                               // Fast Marching Initializer
00396         LevelSet<3,float>                 phi(image.shape(), (double) 2.0);   // Level Set 
00397         Segmentation<3, unsigned short>    segm(&phi, &init, &image);         // Segmentation object
00398 
00399         // Initalization of the level set
00400         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00401 
00402         RegularizationModule<3,float>                     regul(0.75f);                         // parameters : weight
00403         LocalGaussianRegionModule<3,float, float>         local_gaussian( &image, 0.75f, 2);    // parameters : image, weight, smooth
00404         SoftShapeModule<3,float>                          soft_shape(0.75f);                    // parameters : weight, radius of tubular shape
00405         BalloonModule<3,float>                            balloon(0.5f);                        // parameters : weight
00406 
00407         local_gaussian.use_variance(false);
00408         //soft_shape.set_sigma(3.0f); 
00409 
00410         segm.add_module(regul);
00411         segm.add_module(local_gaussian);
00412         segm.add_module(balloon);
00413  
00414         // launch the segmentation
00415         Segmentation<3,float>::RGBPixel neg(215,152,0); 
00416         Segmentation<3,float>::RGBPixel pos(255,252,0); 
00417 
00418         timeval _time;
00419         Globals::initChrono(_time);
00420 
00421          segm.segment(1000, "result_3D/full_head_2/localstats_", neg, pos);
00422         //segm.segment(0, "localstats_", neg, pos); 
00423 
00424         Globals::printChrono("\n        => LocalStats : ", _time); 
00425         Globals::initChrono(_time); 
00426  
00427         //system("./convert_cimg.sh vti result_good_init/mt_diff7_10/localstats_phi*");
00428 
00429    #if 0 // test autre image variant sccessifs
00430         
00431         BaseImage<3, unsigned char> image2;
00432         std::ifstream f2("image/T1_512_msp.inr");
00433         f2 >> format("inrimage") >> image2;
00434         phi.resize(image2.shape());
00435         segm.set_data_term(image2);
00436         local_gaussian.set_data_term(image2);
00437         phi.init_small_spheres(10, 30);    // parameters : radius, windows
00438         segm.segment(1000, "result_good_init/T1_512_msp_10/localstats_", neg, pos);
00439         system("./convert_cimg.sh vti result_good_init/lT1_512_msp_10/localstats_phi*");
00440         
00441         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00442         local_gaussian.set_sigma(20);
00443         segm.segment(1000, "result_good_init/localstats_20/localstats_", neg, pos);
00444         system("./convert_cimg.sh png result_good_init/localstats_20/localstats_*ppm");
00445         
00446         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00447         local_gaussian.set_sigma(40);
00448         segm.segment(1000, "result_good_init/localstats_40/localstats_", neg, pos);
00449         system("./convert_cimg.sh png result_good_init/localstats_40/localstats_*ppm");
00450         
00451         phi.init_small_spheres(10, 25);    // parameters : radius, windows
00452         local_gaussian.set_sigma(80);
00453         segm.segment(1000, "result_good_init/localstats_80/localstats_", neg, pos);
00454         system("./convert_cimg.sh png result_good_init/localstats_80/localstats_*ppm");
00455    #endif
00456         std::cout << "OK" << std::endl;
00457 #endif
00458 }

For further information, please contact Jerome Piovano - Last update 2008-02-08