Jerome Piovano

Research in Medical Imaging and Computer Vision

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

/home/jpiovano/src/Odyssee++/trunk/Libs/LevelSet/src/Segmentation.h

Go to the documentation of this file.

#ifndef SEGMENTATION_H
#define SEGMENTATION_H

#include <sstream>
#include <vector>

#include <Globals.h>
#include <Schemes.h>
#include <LevelSet.h>
#include <Initializer.h>
#include <DericheFilter.h>

#include <SegmentationModule/SegmentationModule.h>

using namespace Images;
using namespace levelset;


namespace segmentation {

        template <class _Modules,class _LevelSet,class _Index>
        float sum(_Modules & modules, const _LevelSet & phi, const _Index & ind, int& indice)
        {
                //TODO automatic estimation via CFL conditions
                //TODO normaliser a 1 la somme des weights => dt <= 1/v_max
                float result = 0;
                for (int i=0; i<modules.size(); i++){
                        const float value  = modules[i]->value(phi, ind);
                        const float weight = modules[i]->weight(ind);
                        result += value*weight;
                }
                return result;
        }

        template<unsigned DIM, typename Pixel = float>
        class Segmentation 
        {
        public: 

                typedef           BaseImage<DIM, Pixel>              ImageType;
                typedef           Initializer<DIM, float>            Initializer;
                typedef           LevelSet<DIM, float>               LevelSet;
                typedef typename  LevelSet::Index                    Index;
                typedef           SegmentationModule<DIM,float>      SegmentationModule;
                typedef           std::vector<SegmentationModule *>  Modules;

                typedef           float (* ModuleAssembler)(Modules &, const LevelSet &, const Index &, int &);

                typedef           Schemes<DIM,float>                 Schemes;
                typedef           Pixels::RGB<unsigned char>         RGBPixel;
                typedef           BaseImage<DIM, RGBPixel>           ColorImageType;


                Segmentation(LevelSet *              _phi,              
                             Initializer *           _init,             
                             ImageType *             _data  = NULL,     
                             BaseImage<DIM, bool> *  _mask  = NULL,     
                             float                   _dt = 0.5f,        
                             int                     _reinit_freq = 1   
                             ); 

                virtual ~Segmentation();



                void add_module(SegmentationModule & _sm);
                void remove_module(SegmentationModule & _sm);
                void set_module_assembler(ModuleAssembler _assembler);
                void set_level_set(LevelSet & _phi);
                void set_initializer(Initializer & _init);
                void set_data_term(ImageType & _data);
                void set_mask(BaseImage<DIM, bool>& _mask);
                float & dt();
                int & reinit_freq();
        

                
                ColorImageType get_result(RGBPixel neg = RGBPixel(0,0,0), RGBPixel pos = RGBPixel(255,255,255));
                ImageType get_binary_result();
                ImageType get_speed(int _module = 0);
                float get_energy();



                void segment(int nbiter = 30, std::string evodir="", RGBPixel neg = RGBPixel(0,0,0), RGBPixel pos = RGBPixel(255,255,255));
                bool stop_condition(int iter, int nbiter);
                void init();
                void iterate(int iter, int nbiter);

        protected:



                LevelSet *                            m_phi;                
                BaseImage<DIM, float>                 m_dphi;               
                Initializer *                         m_init;               
                ImageType *                           m_data;               
                std::vector<SegmentationModule *>     m_sm;                 
                ModuleAssembler                       m_assembler;          
                BaseImage<DIM, bool> *                m_mask;               
                int                                   m_reinit_freq;        
                float                                 m_dt;                 
                bool                                  m_stop_condition;     

        };



        // Constructors - Destructor - Initializers


        template<unsigned DIM, typename Pixel>
                Segmentation<DIM, Pixel>::Segmentation(LevelSet * _phi,  Initializer * _init, ImageType * _data, BaseImage<DIM, bool> * _mask, float _dt, int _reinit_freq)
                :m_phi(_phi), m_dphi(_phi->shape()), m_init(_init), m_data(_data), m_assembler(&sum), m_mask(_mask), m_reinit_freq(_reinit_freq), m_dt(_dt), m_stop_condition(false)
        { 
                if (_data != NULL && _data->shape() != _phi->shape()) {
                        std::cerr << "Error : Shape of the data term and shape of the LevelSet don't match"  << std::endl;
                        exit(1);
                }
                if (_mask != NULL && _mask->shape() != _phi->shape()) {
                        std::cerr << "Error : Shape of the mask and shape of the LevelSet don't match"  << std::endl;
                        exit(1);
                }
        }

        template<unsigned DIM, typename Pixel>
                Segmentation<DIM, Pixel>::~Segmentation()
        { }


        // Modify segmentation information


        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::add_module(SegmentationModule & _sm) 
        {
                m_sm.push_back(&_sm);
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::remove_module(SegmentationModule & _sm)  
        {
                for (int i=0 ; i<m_sm.size() ; i++)
                        if (m_sm[i] == &_sm){
                                std::swap(m_sm[i], *m_sm.end());
                                m_sm.pop_back();
                                return;
                        }
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::set_module_assembler(ModuleAssembler _assembler)
        {
                m_assembler = _assembler;
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::set_mask( BaseImage<DIM, bool>& _mask) 
        {
                if (_mask.shape() != m_phi->shape()) {
                        std::cerr << "Error : Shape of the mask and shape of the LevelSet don't match"  << std::endl;
                        exit(1);
                }
                m_mask = &_mask;
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::set_initializer(Initializer & _init)
        {
                m_init = &_init;
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::set_level_set(LevelSet & _phi)
        {
                m_phi = &_phi;
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::set_data_term(ImageType & _data)
        {
                if (_data.shape() != m_phi->shape()) {
                        std::cerr << "Error : Shape of the data term and shape of the LevelSet don't match"  << std::endl;
                        exit(1);
                }
                m_data = &_data;
        }

        template<unsigned DIM, typename Pixel>
                float & Segmentation<DIM, Pixel>::dt()
        {
                return m_dt;
        }

        template<unsigned DIM, typename Pixel>
                int & Segmentation<DIM, Pixel>::reinit_freq()
        {
                return m_reinit_freq;
        }


        // Results of the segmentation


        template<unsigned DIM, typename Pixel>
                typename Segmentation<DIM, Pixel>::ColorImageType Segmentation<DIM, Pixel>::get_result(RGBPixel neg, RGBPixel pos)
        {
                ColorImageType result;

                result.resize(m_phi->shape());

                if (m_data == NULL) 
                        result = 255;
                else 
                        result = blitz::cast<unsigned char>(*m_data);

                #pragma omp parallel for
                for (int i = 0 ; i < m_phi->nb_points() ; i++){
                        const Index& ind = m_phi->get_point(i);

                        const float _phi = (*m_phi)(ind);

                        int inf_coords[DIM];
                        int sup_coords[DIM];
                        Schemes::find_neighbours(*m_phi, ind, inf_coords, sup_coords);

                        for (int i=0 ; i<DIM ; i++){
                                Index inf = ind;
                                Index sup = ind;
                                inf(i+1) = inf_coords[i];
                                sup(i+1) = sup_coords[i];
                                if ((_phi * (*m_phi)(inf) <= 0) || (_phi * (*m_phi)(sup) <= 0))
                                        result(ind) = (_phi>0)? pos: neg;
                        }
                }

                return result;
        }

        template<unsigned DIM, typename Pixel>
                typename Segmentation<DIM, Pixel>::ImageType Segmentation<DIM, Pixel>::get_binary_result() 
        {
                ImageType result(m_phi->shape());
                for (typename ImageType::template iterator<domain> ind=result.begin() ; ind!=result.end() ; ++ind )
                        result(ind) = ((*m_phi)(ind) < 0)? 0 : 255;

                return result;
        }

        template<unsigned DIM, typename Pixel>
                typename Segmentation<DIM, Pixel>::ImageType Segmentation<DIM, Pixel>::get_speed(int _module)
        {
                ImageType result(m_phi->shape());
                result = 0;

                if (_module < 0 ||  _module > m_sm.size()){
                        std::cerr << " Module does not exist !" << std::endl;
                        return result;
                }

                if (_module == 0){
                        int indice;
                        #pragma omp parallel for
                        for (int i = 0 ; i < m_phi->nb_points() ; i++){
                                const Index& ind = m_phi->get_point(i);
                                if (m_mask==NULL || (*m_mask)(ind)==true)
                                        result(ind) = m_assembler(m_sm, *m_phi, ind, indice);
                        }
                }
                else {
                        #pragma omp parallel for
                        for (int i = 0 ; i < m_phi->nb_points() ; i++){
                                const Index& ind = m_phi->get_point(i);
                                if (m_mask==NULL || (*m_mask)(ind)==true){
                                        const float value  = m_sm[_module-1]->value(*m_phi, ind);
                                        const float weight = m_sm[_module-1]->weight(ind);
                                        result(ind) = value*weight;
                                }
                        }
                }
                
                return result;
        }

        template<unsigned DIM, typename Pixel>
                float Segmentation<DIM, Pixel>::get_energy() 
        {
                float total_energy = 0;
        
                for (unsigned int i = 0 ; i < m_phi->nb_points() ; i++){

                        const Index& ind = m_phi->get_point(i);

                        if (m_mask==NULL || (*m_mask)(ind)!=0) 
                                for (int j=0; j<m_sm.size(); j++){
                                        const float energy = m_sm[j]->energy(*m_phi, ind);
                                        const float weight = m_sm[j]->weight(ind);
                                        total_energy += energy*weight;
                                }
                }

                return total_energy;
        }


        // Launch the segmentation

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::segment(int nbiter, std::string evodir, RGBPixel neg, RGBPixel pos) 
        {
                // initialization
                init();

                // main loop
                for (int i=0 ; !m_stop_condition ; i++){

                        std::cout << i << std::endl;

                        #ifdef SAVE_EVOLUTION
                        if((i%int(SAVE_EVOLUTION))==0) {

                                std::ostringstream name;

                                //TODO oblige de faire 3 objets ofstream sinon foire 
                                std::ofstream file1;
                                std::ofstream file2;
                                std::ofstream file3;

                                #if 1
                                //TODO energy (?)
                                std::cout <<  " - Energy : " << get_energy() << std::endl;

                                name.str(""); 
                                name << evodir << "phi_" << std::setw(6) << std::setfill('0') << i << ".inr";
                                file1.open(name.str().c_str());
                                file1 << format("inrimage") << *m_phi;
                                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                                #endif

                                #if 1
                                if (DIM == 2){
                                        // colored ppm only in 2D
                                        name.str("");
                                        name << evodir << "result_" << std::setw(6) << std::setfill('0') << i << ".ppm";
                                        file2.open(name.str().c_str());
                                        file2 << format("rawppm") << get_result(neg, pos);
                                        std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                                }
                                else {
                                        // colored inrimage-5 in 2D/3D
                                        name.str("");
                                        name << evodir << "result_" << std::setw(6) << std::setfill('0') << i << ".inr";
                                        file2.open(name.str().c_str());
                                        ImageType result(m_phi->shape());
                                        result = blitz::cast<float>(get_result(neg, pos));
                                        file2 << format("Inrimage-5") << result; //blitz::cast<float>(get_result(neg, pos)); 
                                        std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                                }
                                #endif

                                #if 0  // speed
                                for (int j=1; j<=m_sm.size(); j++){
                                        std::ofstream file3;
                                        name.str("");
                                        name << evodir << "speed_module" << j << "_" << std::setw(6) << std::setfill('0') << i << ".inr";
                                        file3.open(name.str().c_str());
                                        file3 << format("inrimage") << get_speed(j) ;
                                        std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                                }
                                name.str("");
                                name << evodir << "speed_" << std::setw(6) << std::setfill('0') << i << ".inr";
                                file3.open(name.str().c_str());
                                file3 << format("inrimage") << get_speed(0) ;
                                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                                #endif

                                #if 1
                                for (int j=0; j<m_sm.size(); j++)
                                        m_sm[j]->save_result(*m_phi, evodir, i);
                                #endif
                        }
                        #endif

                        iterate(i, nbiter);
                }
        }

        template<unsigned DIM, typename Pixel>
                bool Segmentation<DIM, Pixel>::stop_condition(int iter, int nbiter)
        {
                // stop at convergence
                if (nbiter < 0){

                        // static object that keep traces of variation of the contour.
                        // stop if no variations (but maybe oscillation of max period 50) of the contour during 10 iter, or iter > 2000 (moche mais comment ?)
                        static BaseImage<DIM, short> sign_change(m_phi->shape());
                        static int count = 0;

                        if (iter == 0) {
                                sign_change = 0;
                                count = 0;
                        }

                        int nb_sign_change = 0;
                        for (int i = 0 ; i < m_phi->pos_to_neg_index() ; i++){
                                const Index& ind = m_phi->get_point(i);
                                if (iter - sign_change(ind) > 50) // no oscillation
                                        nb_sign_change++;
                                sign_change(ind) = iter;
                        }

                        for (int i = m_phi->neg_to_pos_index() ; i < m_phi->nb_points() ; i++){
                                const Index& ind = m_phi->get_point(i);
                                if (iter - sign_change(ind) > 50) // no oscillation
                                        nb_sign_change++;
                                sign_change(ind) = iter;
                        }

                        if (nb_sign_change > 0) count = 0;
                        return ((nb_sign_change == 0) && (++count > 20)) || (iter >= 2000);
                }

                // stop when (iter > iter max)
                else return (iter >= nbiter);
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::init()
        {
                // resize of d_phi
                if (m_dphi.shape() != m_phi->shape())
                        m_dphi.resize(m_phi->shape());
                        
                m_dphi=0.0f;

                // initialization of the modules
                for (int i=0; i<m_sm.size(); i++)
                        m_sm[i]->init(*m_phi);

                m_stop_condition = false;
        }

        template<unsigned DIM, typename Pixel>
                void Segmentation<DIM, Pixel>::iterate(int iter, int nbiter) 
        {
                #ifdef PRINT_TIME
                timeval _time;
                Globals::initChrono(_time);
                #endif

                m_dphi=0.0f;

//ImageType result(m_phi->shape());
// result = 0;
int indice;

                // speed computation
                #pragma omp parallel for
                for (int i = 0 ; i < m_phi->nb_points() ; i++){

                        const Index& ind = m_phi->get_point(i);
                        m_dphi(ind)=0.0f;

                        if (m_mask==NULL || (*m_mask)(ind)==true)
                                m_dphi(ind) = m_assembler(m_sm, *m_phi, ind, indice);
                                //result(ind) = indice;
                }

//std::ostringstream name;
//name.str("");
//std::ofstream file3;
//name << std::setw(6) << std::setfill('0') << iter << ".inr";
//file3.open(name.str().c_str());
//file3 << format("inrimage") << result;
//std::cout <<  "       - Saved " << name.str() << " ..."  << std::endl;

                #ifdef PRINT_TIME
                Globals::printChrono("  => compute speed : ", _time);
                Globals::initChrono(_time);
                #endif

                // phi += dt * d_phi
                m_phi->update(m_dphi, m_dt);

                #ifdef PRINT_TIME
                Globals::printChrono("  => update Level set : ", _time);
                Globals::initChrono(_time);
                #endif

                // stop condition
                m_stop_condition = stop_condition(iter, nbiter);

                #ifdef PRINT_TIME
                Globals::printChrono("  => stop condition : ", _time);
                Globals::initChrono(_time);
                #endif

                // update segmentation modules
                for (int i=0; i<m_sm.size(); i++)
                        if (m_sm[i]->adaptive()) m_sm[i]->update(*m_phi);

                #ifdef PRINT_TIME
                Globals::printChrono("  => total update modules : ", _time);
                Globals::initChrono(_time);
                #endif

                // level-set reinitialization
                if(iter%m_reinit_freq == 0) 
                        m_init->reinit_band(*m_phi);

                #ifdef PRINT_TIME
                Globals::printChrono("  => total reinit band : ", _time);
                Globals::initChrono(_time);
                #endif
        }
}

#endif // SEGMENTATION_H

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