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/SegmentationModule/CircularityModule.h

Go to the documentation of this file.

#ifndef CIRCULARITYMODULE_H
#define CIRCULARITYMODULE_H

#include <SegmentationModule/SegmentationModule.h>
#include <image/RecFilters.H>
#include <image/IIR.H>
#include <DericheFilter.h> 


namespace segmentation {


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

                typedef           SegmentationModule<DIM, Pixel>  base;
                typedef typename  base::LevelSet                  LevelSet;
                typedef typename  base::Index                     Index;
                typedef typename  base::Schemes                   Schemes;
                
                typedef           BaseImage<DIM, float>           ImageFloat;
                typedef           Histogram<int>                  Histogram;



                CircularityModule(float       _weight = 1.f,      
                                  bool        _var=false          
                                 );

                virtual ~CircularityModule();



                void use_variance(bool _var=true);
                Histogram * hist_x();
                Histogram * hist_y();



                float value(const LevelSet & phi, const Index & ind);
                float energy(const LevelSet & phi, const Index & ind);
                void init(const LevelSet & phi);
                void update(const LevelSet & phi);
                void save_result(const LevelSet & phi, std::string evodir, int iter);

                void update_params(const LevelSet & phi);
                void init_histograms(const LevelSet & phi);
                void update_histograms(const LevelSet & phi);
                void save_histograms(std::string evodir, int iter);
        

        protected:


                std::vector<Pixel>          m_ll_x;       
                std::vector<Pixel>          m_ll_y;       
                Histogram *                 m_hist_x;     
                Histogram *                 m_hist_y;     
                Images::Index<DIM, float>   m_centroid;   
                Images::Index<DIM, float>   m_variance;   
                bool                        m_var;        


        };


        // Constructors - Destructor


        template <unsigned DIM, typename Pixel>
                CircularityModule<DIM, Pixel>::CircularityModule(float _weight, bool _var) 
                : SegmentationModule<DIM, Pixel>(_weight, true), m_hist_x(NULL), m_hist_y(NULL)
        { }

        template <unsigned DIM, typename Pixel>
                CircularityModule<DIM, Pixel>::~CircularityModule() 
        { 
                if (m_hist_x != NULL){
                        delete m_hist_x;
                        delete m_hist_y;
                }
        }


        //  Member functions


        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::use_variance(bool _var)
        {
                m_var = _var;
        }

        template <unsigned DIM, typename Pixel>
                typename CircularityModule<DIM, Pixel>::Histogram* CircularityModule<DIM, Pixel>::hist_x()
        {
                return m_hist_x;
        }

        template <unsigned DIM, typename Pixel>
                typename CircularityModule<DIM, Pixel>::Histogram* CircularityModule<DIM, Pixel>::hist_y()
        {
                return m_hist_y;
        }

        // Member functions


        template <unsigned DIM, typename Pixel>
                float CircularityModule<DIM, Pixel>::value(const LevelSet & phi, const Index & ind)
        { 
                const int _x = ind(1);
                const int _y = ind(2);

                return  - Globals::delta(phi(ind)) * (Maths::sqr(_x - m_centroid(1))  + Maths::sqr(_y - m_centroid(2)));

                
/*              
                const int _x = ind(1);
                const int _y = ind(2);

                const Pixel ll_x = m_ll_x[m_hist_x->get_index(_x)];
                const Pixel ll_y = m_ll_y[m_hist_y->get_index(_y)];

                //return Globals::delta(phi(ind)) * (ll_pos - ll_neg);
                return  Globals::delta(phi(ind)) * (ll_x + ll_y);
*/
        }       

        template <unsigned DIM, typename Pixel>
                float CircularityModule<DIM, Pixel>::energy(const LevelSet & phi, const Index & ind)
        {
                //TODO 
                return 0;
        }

        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::init(const LevelSet & phi)
        {
                base::init(phi);
                init_histograms(phi);
                update_params(phi);
        }

        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::update(const LevelSet & phi)
        {
                update_histograms(phi);
                update_params(phi);
        }

        template <unsigned DIM, typename Pixel>
                 void CircularityModule<DIM, Pixel>::update_params(const LevelSet & phi) 
        {
                m_centroid = 0.0f;
                m_variance = 0.0f;

                // compute the mean and variance of the shape in the x and y dimension
                m_hist_x->compute_statistics(m_centroid(1), m_variance(1), m_var);
                m_hist_y->compute_statistics(m_centroid(2), m_variance(2), m_var);

                // if the size of the histogram has changed, resize de log-likelihood vectors
                if (m_ll_x.size() != m_hist_x->size()) m_ll_x.resize(m_hist_x->size());
                if (m_ll_y.size() != m_hist_y->size()) m_ll_y.resize(m_hist_y->size());

                if (m_var) {
                        for(int i=0; i<m_ll_x.size(); i++)
                                m_ll_x[i] = -0.5f*log(2*M_PI*m_variance(1)) - 0.5f*Maths::sqr( i+m_hist_x->min() - m_centroid(1))/m_variance(1);
                        for(int i=0; i<m_ll_y.size(); i++)
                                m_ll_y[i] = -0.5f*log(2*M_PI*m_variance(2)) - 0.5f*Maths::sqr( i+m_hist_y->min() - m_centroid(2))/m_variance(2);
                }
                else{
                        for(int i=0; i<m_ll_x.size(); i++)
                                m_ll_x[i] =  - Maths::sqr(i+m_hist_x->min() - m_centroid(1));
                        for(int i=0; i<m_ll_y.size(); i++) 
                                m_ll_y[i] =  - Maths::sqr(i+m_hist_y->min() - m_centroid(2));
                }
        }

        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::init_histograms(const LevelSet & phi)
        {
                if (m_hist_x == NULL){
                        m_hist_x = new Histogram(0, phi.extent(0), true);
                        m_hist_y = new Histogram(0, phi.extent(1), true);
                }

                //computation of the first histograms

                m_hist_x->clear();
                m_hist_y->clear();

                for (typename LevelSet::template iterator<domain> ind=phi.begin() ; ind!=phi.end() ; ++ind ){
                        const Pixel _phi = phi(ind);
                        const int     _x = (ind.position())(1);
                        const int     _y = (ind.position())(2);

                        if (_phi>=0) {
                                m_hist_x->add_pixel(_x);
                                m_hist_y->add_pixel(_y);
                        }
                }
        }

        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::update_histograms(const LevelSet & phi)
        {
                 if (m_hist_x == NULL){
                        std::cerr << "Error : Histograms not initialized"  << std::endl;
                        exit(1);
                }

                //computation of the histograms using previous ones

                if (!m_hist_x->adaptative() && !m_hist_y->adaptative()) return;

                for (unsigned int i = 0 ; i < phi.pos_to_neg_index() ; i++){
                        const Index& ind = phi.get_point(i);
                        const int     _x = ind(1);
                        const int     _y = ind(2);

                        m_hist_x->remove_pixel(_x);
                        m_hist_y->remove_pixel(_y);
                }

                for (unsigned int i = phi.neg_to_pos_index() ; i < phi.nb_points() ; i++){
                        const Index& ind = phi.get_point(i);
                        const int     _x = ind(1);
                        const int     _y = ind(2);

                        m_hist_x->add_pixel(_x);
                        m_hist_y->add_pixel(_y);
                }
        }

        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::save_result(const LevelSet & phi, std::string evodir, int iter)
        {
                ImageFloat tmp(phi.shape());
                std::ostringstream name;

                for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
                        const Pixel _phi = phi(i);
                        if (_phi<=0) tmp(i) = 1;
                        else {
                                double coeff = exp(-0.5f*( Maths::sqr(i(1)-m_centroid(1))/Maths::sqr(m_variance(1)) + Maths::sqr(i(2)-m_centroid(2))/Maths::sqr(m_variance(2)) ))/2;
                                tmp(i) = coeff;
                        }
                }

                name.str("");
                name << evodir << "centroid_" << std::setw(6) << std::setfill('0') << iter << ".inr";
                std::ofstream file(name.str().c_str());
                file << format("inrimage") << tmp;
                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;

                save_histograms(evodir, iter);

                if (m_var) {
                        
                        std::ofstream file2;
                        // plot des gaussienes associes a mean et var

                        name.str("");
                        name << evodir << "Gaussian_x_" << std::setw(6) << std::setfill('0') << iter << ".txt";
                        file2.open(name.str().c_str());
                        for(int i=m_hist_x->min(); i < m_hist_x->max() ; i++) 
                                file2 << i << " " << Globals::gauss_density(m_centroid(1), m_variance(1), i) << std::endl;
                        file2.close();
                        std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
        
                        name.str("");
                        name << evodir << "Gaussian_y_" << std::setw(6) << std::setfill('0') << iter << ".txt";
                        file2.open(name.str().c_str());
                        for(int i=m_hist_y->min(); i < m_hist_y->max() ; i++) 
                                file2 << i << " " << Globals::gauss_density(m_centroid(2), m_variance(2), i) << std::endl;
                        file2.close();
                        std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                }
        }

        
        template <unsigned DIM, typename Pixel>
                void CircularityModule<DIM, Pixel>::save_histograms(std::string evodir, int iter) 
        {
                if (m_hist_x == NULL){
                        std::cerr << "Error : Histograms not initialized"  << std::endl;
                        exit(1);
                }

                std::ostringstream name;
                std::ofstream file;

                name.str("");
                name << evodir << "Hist_x_" << std::setw(6) << std::setfill('0') << iter << ".txt";
                file.open(name.str().c_str());
                for(int i=m_hist_x->min(); i < m_hist_x->max() ; i++) 
                        file << i << " " << (*m_hist_x)[i]  << std::endl;
                file.close();
                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;

                name.str("");
                name << evodir << "Hist_y_" << std::setw(6) << std::setfill('0') << iter << ".txt";
                file.open(name.str().c_str());
                for(int i=m_hist_y->min(); i < m_hist_y->max() ; i++) 
                        file << i << " " << (*m_hist_y)[i] << std::endl;
                file.close();
                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
        }

}

#endif // CIRCULARITYMODULE_H

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