Jerome Piovano

Research in Medical Imaging and Computer Vision

• Home
• Research
• Publications
• Teaching
• Biography
• Shared code and data

2D_SoftShapePrior.C

Use of a 2D soft shape prior

#include <Segmentation.h>
#include <SegmentationModule/RegularizationModule.h>
#include <SegmentationModule/SoftShapeModule.h>
#include <SegmentationModule/BalloonModule.h>
#include <SegmentationModule/GaussianRegionModule.h>
#include <SegmentationModule/CircularityModule.h>

int main(int argc, char* argv[])
{       
        using namespace Images;
        using namespace levelset;
        using namespace segmentation;

#if 0 // creation sulci

        const Dimension nx = 300;
        const Dimension ny = 300;
        const Dimension dims[] = {nx, ny};
        Shape<2> shape(dims);
 
        // initialization of the segmentations objects
        FastMarchingInitializer<2,float> init;                         // Fast Marching Initializer
        LevelSet<2,float>                phi(shape, (double) 2.0);     // Level Set of band thickness = 2
        Segmentation<2,float>            segm(&phi, &init);            // Segmentation  

        LevelSet<2,float>::Index c(150, 150);
        phi.init_sphere(c, 70);
 
        // initialization of the modules 
        RegularizationModule<2,float>     regul(1.f);             // parameters : weight
        BalloonModule<2,float>            balloon(1.f);           // parameters : weight

        BaseImage<2,float> weight(shape);
        weight = 0;

        // bbox
        int height = 110;
        int width  = 30;
        int up    = 130;
        int down  = 240;
        int left  = 135; 
        int right = 165;

        for (BaseImage<2, float>::iterator<domain> ind=weight.begin() ; ind!=weight.end() ; ++ind){
                if (ind(1)>left && ind(1)<right && ind(2)>up && ind(2)<down){
                        int x = ind(1) - left;
                        int y = ind(2) - up;
                        if( fabs(x - width/2) <  y*(width/2)/height)
                                weight(ind) = -1.f;
                }
        }
                        
        balloon.set_image_weight(&weight);

        segm.add_module(regul);
        segm.add_module(balloon);

        timeval _time;
        Globals::initChrono(_time); 

        // launch the segmentation 
        segm.segment(1000, "result_softshapeprior/2D_brain_sulci/soft_shape_2D_");

        Globals::printChrono("\n        => Soft Shape Prior : ", _time);
        Globals::initChrono(_time);

        std::cout << "OK" << std::endl;

        system("./convert_cimg.sh png result_softshapeprior/2D_brain_sulci/soft_shape_2D_*ppm");
 

        return 0;
#endif

#if 1 // circle creation

        BaseImage<2, float> image;
        //std::ifstream f("image/test_circularite.inr");
        std::ifstream f("image/sylvain/S_test_G_1841_ROIs/S_test_G_1841_ROI1.inr");
        f >> format("inrimage") >> image;
        Shape<2> shape = image.shape();

        // initialization of the segmentations objects
        FastMarchingInitializer<2,float> init;                         // Fast Marching Initializer
        LevelSet<2,float>                phi(shape, (double) 3.0);     // Level Set of band thickness = 2
        Segmentation<2,float>            segm(&phi, &init, &image);    // Segmentation  

        // initalization of the level set
        //phi.init_border(40, init);
        phi.init_small_spheres(20, 50, init);
        //LevelSet<2,float>::Index tl(50, 50);
        //LevelSet<2,float>::Index br(250, 250);
        //phi.init_rectangle(tl, br, init, true);
        
        // initialization of the modules 
        RegularizationModule<2,float>        regul(1.f);                   // parameters : weight
        CircularityModule<2,float>           circular(0.00005f, true);        // parameters : weight
        //BalloonModule<2,float>               balloon(0.2f);                // parameters : weight
        GaussianRegionModule<2,float,float>  gauss(&image, 0.113f, true);  // parameters : image, weight

        segm.add_module(regul);
        segm.add_module(circular);
        //segm.add_module(balloon);
        segm.add_module(gauss);

        timeval _time;
        Globals::initChrono(_time);

        Segmentation<2,float>::RGBPixel neg(215,152,0); 
        Segmentation<2,float>::RGBPixel pos(255,252,0); 

        std::ostringstream name;
        name.str("");
        name << "mkdir result_softshapeprior/circle/";

        name.str("");
        name << "result_softshapeprior/circle/";
                
        segm.segment(3000,  name.str().c_str(), neg, pos);

        name.str("");
        name << "convert_cimg.sh png result_softshapeprior/circle/*ppm";
        system(name.str().c_str());

        Globals::printChrono("\n        => Circular : ", _time);
        Globals::initChrono(_time);
        
        std::cout << "OK" << std::endl;
 
        return 0;
#endif

#if 0 // 2D version

        const Dimension nx = 300;
        const Dimension ny = 300;
        const Dimension dims[] = {nx, ny};
        Shape<2> shape(dims);
 
        // initialization of the segmentations objects
        FastMarchingInitializer<2,float> init;                         // Fast Marching Initializer
        LevelSet<2,float>                phi(shape, (double) 2.0);     // Level Set of band thickness = 2
        Segmentation<2,float>            segm(&phi, &init);            // Segmentation  

        // initalization of the level set
        //phi.init_small_spheres(15, 40);    // parameters : radius, windows 
        //phi.init_small_spheres(10, 60);
        //LevelSet<2,float>::Index tl(100,100);
        //LevelSet<2,float>::Index br(140,140);
        //phi.init_rectangle(tl, br);
        LevelSet<2,float>::Index c(150, 150);
        phi.init_sphere(c, 70);
        //BaseImage<2, float> image;
        //std::ifstream f("soft_shape_2D_phi_001000.inr");
        //f >> format("inrimage") >> image;
        //phi.init_image(image, init);
 
        // initialization of the modules 
        RegularizationModule<2,float>     regul(1.f);             // parameters : weight
        SoftShapeModule<2,float>          soft_shape(1.f);        // parameters : weight, radius of tubular shape
        BalloonModule<2,float>            balloon(1.f);           // parameters : weight
        //BalloonModule<2,float>            balloon(0.75f);           // parameters : weight

#if 0
        float radius1 = 80;
        float ampl1   = 5;
        float freq1   = 8;
        float radius2 = 45;
        float ampl2   = 5;
        float freq2   = 5;

        BaseImage<2,bool> mask(shape);
        for (BaseImage<2, float>::iterator<domain> ind=mask.begin() ; ind!=mask.end() ; ++ind){
                double pol = atan2(ind(1) - 150, ind(2) - 150);
                if (   (sqrt( (double) Maths::sqr(ind(1) - 150) +  Maths::sqr( ind(2) - 150)) < ampl1*sin(freq1*pol) + radius1 + 3*sin(10*pol + 2))
                    && (sqrt( (double) Maths::sqr(ind(1) - 150) +  Maths::sqr( ind(2) - 150)) > ampl2*sin(freq2*pol) + radius2 + 3*sin(10*pol + 2)) )
                         mask(ind) = true;
                else    
                        mask(ind) = false;
        }
        segm.set_mask(mask);
#endif
        float radius = 65;
        float width  = 20;
        BaseImage<2,float> weight(shape);
        for (BaseImage<2, float>::iterator<domain> ind=weight.begin() ; ind!=weight.end() ; ++ind){
                float dist = sqrt( (double) Maths::sqr(ind(1) - 150) +  Maths::sqr( ind(2) - 150)) - radius;
                if ( fabs(dist) < width)
                         weight(ind) = (width - fabs(dist))/width;
                else    
                        weight(ind) = 0;
        }
        soft_shape.set_image_weight(&weight);
        balloon.set_image_weight(&weight);

        soft_shape.set_sigma(7.0f);
        soft_shape.p() = 2;

        segm.add_module(regul);
        segm.add_module(soft_shape); 
        segm.add_module(balloon);

        timeval _time;
        Globals::initChrono(_time); 

        // launch the segmentation 
        segm.segment(1000, "result_softshapeprior/2D_brain_weight/soft_shape_2D_");

        Globals::printChrono("\n        => Soft Shape Prior : ", _time);
        Globals::initChrono(_time);

        std::cout << "OK" << std::endl;

        system("./convert_cimg.sh png result_softshapeprior/2D_brain_weight/soft_shape_2D_*ppm");
 

        return 0;
#endif

#if 0 // 3D version 


        const Dimension nx = 300;
        const Dimension ny = 300; 
        const Dimension nz = 300;
        const Dimension dims[] = {nx, ny, nz};
        Shape<3> shape(dims);
 
        // initialization of the segmentations objects
        FastMarchingInitializer<3,float> init;                         // Fast Marching Initializer
        LevelSet<3,float>                phi(shape, (double) 2.0);     // Level Set of band thickness = 2
        Segmentation<3,float>            segm(&phi, &init);            // Segmentation

        // initalization of the level set
        //LevelSet<3,float>::Index c(150, 150, 150);
        //phi.init_sphere(c, 60);
        BaseImage<3, float> image;
        std::ifstream f("a.inr");
        f >> format("inrimage") >> image;
        phi.init_image(image, init);

        // initialization of the modules
        RegularizationModule<3,float>     regul(1.f);             // parameters : weight
        SoftShapeModule<3,float>          soft_shape(1.f);  // parameters : weight, radius of tubular shape
        BalloonModule<3,float>            balloon(1.f);         // parameters : weight

#if 0
        float radius1 = 80;
        float ampl1   = 5;
        float freq1   = 8;
        float radius2 = 45;
        float ampl2   = 5;
        float freq2   = 5;
        float radius3 = 45;
        //float ampl3   = 5;
        float freq3   = 5;
        BaseImage<3,bool> mask(shape);
        for (BaseImage<3, float>::iterator<domain> ind=mask.begin() ; ind!=mask.end() ; ++ind){
                double pol1 = atan2(ind(1) - 150, ind(2) - 150);
                double pol2 = atan2(ind(1) - 150, ind(3) - 150);
                if (   (sqrt( (double) Maths::sqr(ind(1) - 150) +  Maths::sqr( ind(2) - 150) + Maths::sqr( ind(3) - 150)) < radius1 + ampl1*sin(freq1*pol1)*sin(freq3*pol2))
                    && (sqrt( (double) Maths::sqr(ind(1) - 150) +  Maths::sqr( ind(2) - 150) + Maths::sqr( ind(3) - 150)) > radius2 + ampl2*sin(freq2*pol1)*sin(freq3*pol2)) )
                         mask(ind) = true;
                else    
                        mask(ind) = false;
        }
        segm.set_mask(mask);
#endif
        float radius = 65;
        //float width  = 20;
        BaseImage<3,float> weight(shape);
        for (BaseImage<3, float>::iterator<domain> ind=weight.begin() ; ind!=weight.end() ; ++ind){
                float dist = sqrt((double) Maths::sqr(ind(1) - 150) +  Maths::sqr( ind(2) - 150) + Maths::sqr( ind(3) - 150)) - radius;
                weight(ind) = exp(-0.5f*Maths::sqr(dist)/64);
                //if ( fabs(dist) < width)
                //       weight(ind) = (width - fabs(dist))/width;
                //else  
                //      weight(ind) = 0;
        }
        soft_shape.set_image_weight(&weight);
        balloon.set_image_weight(&weight);

        std::ofstream fw("weight.inr");
        fw << format("inrimage") << weight;
        fw.close();

        soft_shape.set_sigma(7.f);
        soft_shape.p() = 2;

        segm.add_module(regul);
        segm.add_module(soft_shape);
        segm.add_module(balloon);

        // launch the segmentation
        segm.segment(1000, "result_softshapeprior/3D_brain_weight/soft_shape_3D_");

        std::cout << "OK" << std::endl;


        return 0;

#endif

#if 0 // test changement topologie
 
        // initialization of the segmentations objects
        FastMarchingInitializer<2,float> init;                    // Fast Marching Initializer
        LevelSet<2,float>                phi;                     // Level Set of band thickness = 2
        Segmentation<2,float>            segm(&phi, &init);       // Segmentation

        phi.thickness() = 2;
        // binary image
        BaseImage<2, float> image;
        std::ifstream f("phi_000925.inr");
        f >> format("inrimage") >> image;
        phi.init_image(image, init);

        // initialization of the modules 
        RegularizationModule<2,float>     regul(1.f);             // parameters : weight
        SoftShapeModule<2,float>          soft_shape(1.f, 5.0f);  // parameters : weight, radius of tubular shape
        BalloonModule<2,float>            balloon(0.75f);         // parameters : weight

        soft_shape.set_sigma(5.0f); 
        soft_shape.p() = 2;

        segm.add_module(regul);
        //segm.add_module(soft_shape); 
        //segm.add_module(balloon);

        timeval _time;
        Globals::initChrono(_time); 

        // d aborsd une fois pour virer le bruit
        // launch the segmentation 
        segm.segment(10, "soft_shape_2D_");

        segm.add_module(soft_shape); 
        segm.add_module(balloon);
        
        // vraie segmentation
        segm.segment(300, "soft_shape_2D_");

        Globals::printChrono("\n        => Soft Shape Prior : ", _time);
        Globals::initChrono(_time); 

        std::cout << "OK" << std::endl;
 

        return 0;

#endif


}

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