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/LocalStatsRegionModule.h

Go to the documentation of this file.

#ifndef LOCALSTATSREGIONMODULE_H
#define LOCALSTATSREGIONMODULE_H

//#define SECOND_TERM

#include <SegmentationModule/RegionModule.h>
#include <DericheFilter.h>

namespace segmentation {


        template <unsigned DIM, typename Pixel = float, typename PixelImage = float>
        class LocalStatsRegionModule :  public RegionModule<DIM, Pixel, PixelImage> 
        {
        public:

                typedef          RegionModule<DIM, Pixel, PixelImage>    base;
                typedef typename base::LevelSet                          LevelSet;
                typedef typename base::Index                             Index;
                typedef typename base::Histogram                         Histogram;
                typedef typename base::ImageType                         ImageType;

                typedef          BaseImage<DIM, float>                   ImageFloat;
                typedef          DericheFilter<DIM, PixelImage>          DericheFilter;



//TODO garder local variance ?

                LocalStatsRegionModule(ImageType * _data,                    
                                       float       _weight = 1.f,            
                                       float       _sigma = 20.0f,           
                                       bool        _localVariance = false    
                                       );

                virtual ~LocalStatsRegionModule();



                void set_sigma(float _sigma);
                bool & use_local_variance();



                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);
        


        protected:

                float          m_sigma;         
                DericheFilter  m_deriche;       
                ImageFloat     m_chi_pos;       
                ImageType      m_data_pos;      
                ImageType      m_data_neg;      

//TODO garder variance ???
                ImageType  m_var_pos;          
                ImageType  m_var_neg;          
                Pixel      m_image_variance;   
                bool       m_local_variance;   

#ifdef SECOND_TERM
// debug : variables en plus pr calculer terme en + ds derivee nrj
                ImageType  m_term1_pos;
                ImageType  m_term1_neg;
                ImageType  m_term2_pos;
                ImageType  m_term2_neg;
// debug : variables en plus pr calculer terme en + ds derivee nrj
#endif


        };


        // Constructors - Destructor


        template <unsigned DIM, typename Pixel, typename PixelImage>
                LocalStatsRegionModule<DIM, Pixel, PixelImage>::LocalStatsRegionModule(ImageType * _data, float _weight, float _sigma, bool _local_variance) 
                : RegionModule<DIM, Pixel, PixelImage>(_data, _weight, 0.0f, 0.0f, true), m_sigma(_sigma), m_deriche(_sigma), m_chi_pos(_data->shape()), m_data_pos(_data->shape()), m_data_neg(_data->shape()), m_var_pos(_data->shape()), m_var_neg(_data->shape()), m_image_variance(-1.0f), m_local_variance(_local_variance)
#ifdef SECOND_TERM
                ,m_term1_pos(_data->shape()),m_term2_pos(_data->shape()),m_term1_neg(_data->shape()),m_term2_neg(_data->shape())
#endif
        { }

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


        // Modify/Access member data

        
        template <unsigned DIM, typename Pixel, typename PixelImage>
                void LocalStatsRegionModule<DIM, Pixel, PixelImage>::set_sigma(float _sigma)
        { 
                m_sigma = _sigma;
                m_deriche.set_isotropic_sigma(_sigma);
        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                bool & LocalStatsRegionModule<DIM, Pixel, PixelImage>::use_local_variance()
        {
                return m_local_variance;
        }


        // Member functions

        
        template <unsigned DIM, typename Pixel, typename PixelImage>
                float LocalStatsRegionModule<DIM, Pixel, PixelImage>::value(const LevelSet & phi, const Index & ind)
        { 
                PixelImage val = (* base::m_data)(ind);

                //TODO garder local variance ?
                float var_pos = m_local_variance?m_var_pos(ind)/m_chi_pos(ind):m_image_variance;
                float var_neg = m_local_variance?m_var_neg(ind)/(1-m_chi_pos(ind)):m_image_variance;

                float ll_pos = -0.5f*log(2*M_PI*var_pos) - 0.5f*Maths::sqr(val-m_data_pos(ind)/   m_chi_pos(ind) )/var_pos;
                float ll_neg = -0.5f*log(2*M_PI*var_neg) - 0.5f*Maths::sqr(val-m_data_neg(ind)/(1-m_chi_pos(ind)))/var_pos;

#ifdef SECOND_TERM
                ll_pos += 2*val*m_term1_pos(ind)/var1 - 2*m_term2_pos(ind)/var_pos;
                ll_neg += 2*val*m_term1_neg(ind)/var2 - 2*m_term2_neg(ind)/var_neg;
#endif

                return -(ll_neg - ll_pos)/3.0;
        }       

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

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

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

         template <unsigned DIM, typename Pixel, typename PixelImage>
                 void LocalStatsRegionModule<DIM, Pixel, PixelImage>::update_params(const LevelSet & phi) 
        {
                // computation of the variance of the image
                if(!m_local_variance && (m_image_variance<0)) {
                        PixelImage min;
                        PixelImage max;
                        PixelImage mean;
                        
                        Misc<DIM,PixelImage>::stats(*base::m_data, min, max, mean, &m_image_variance);
                }


                // computation of m_chi_pos, m_data_pos, m_data_neg
                for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
                        const Pixel      _phi = phi(ind);
                        const PixelImage _val = (*base::m_data)(ind);

                        if (_phi>=0) {
                                m_chi_pos(ind)   = 1.0f;
                                m_data_pos(ind)  = _val;
                                m_data_neg(ind) = 0.0f;
                                if(m_local_variance) {
                                        m_var_pos(ind) = _val;
                                        m_var_neg(ind) = 0.0f;
                                }
                        }
                        else {
                                m_chi_pos(ind)   = 0.0f;
                                m_data_pos(ind)  = 0.0f;
                                m_data_neg(ind) = _val;
                                if(m_local_variance) {
                                        m_var_pos(ind) = 0.0f;
                                        m_var_neg(ind) = _val;
                                }
                        }
                }

#pragma omp parallel sections
{
        #pragma omp section
                m_deriche.filter(m_chi_pos);
        #pragma omp section
                m_deriche.filter(m_data_pos);
        #pragma omp section
                m_deriche.filter(m_data_neg);
}


#ifdef SECOND_TERM
                // computation of m_term1_pos, m_term2_pos, m_term1_neg, m_term2_neg
                for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
                        const Pixel      _phi = phi(ind);
                        const PixelImage _val = (*base::m_data)(ind);

                        if (_phi>=0) {
                                m_term1_pos(ind) = ((*m_data)(ind) - (m_data_pos(ind)/m_chi_pos(ind))) / m_chi_pos(ind);
                                m_term2_pos(ind) = ((*m_data)(ind) - (m_data_pos(ind)/m_chi_pos(ind))) * (m_data_pos(ind)/m_chi_pos(ind))  / m_chi_pos(ind);
                                m_term1_neg(ind) = 0.0f;
                                m_term2_neg(ind) = 0.0f;
                        }
                        else {
                                m_term1_pos(ind) = 0.0f;
                                m_term2_pos(ind) = 0.0f;
                                m_term1_neg(ind) = ((*m_data)(ind) - (m_data_neg(ind)/(1 - m_chi_pos(ind))))  / (1 - m_chi_pos(ind));
                                m_term2_neg(ind) = ((*m_data)(ind) - (m_data_neg(ind)/(1 - m_chi_pos(ind)))) * (m_data_neg(ind)/(1 - m_chi_pos(ind)))  / (1 - m_chi_pos(ind));
                        }
                }
                
#pragma omp parallel sections
{
        #pragma omp section
                m_deriche.filter(m_term1_pos);
        #pragma omp section
                m_deriche.filter(m_term2_pos);
        #pragma omp section
                m_deriche.filter(m_term1_neg);
        #pragma omp section
                m_deriche.filter(m_term2_neg);
}
#endif

                //TODO garder variance ?
                if(m_local_variance) {
                                
                        for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
                                const Pixel      _phi = phi(ind);
                                const PixelImage _val = (*base::m_data)(ind);   
                                if (_phi>=0) {
                                        m_var_pos(ind) = Maths::sqr(m_var_pos(ind) - m_data_pos(ind)/m_chi_pos(ind));
                                        m_var_neg(ind) = 0.0f;
                                }
                                else {
                                        m_var_pos(ind) = 0.0f;
                                        m_var_neg(ind) = Maths::sqr(m_var_neg(ind) - m_data_neg(ind)/(1-m_chi_pos(ind)));
                                }
                        }
#pragma omp parallel sections
{
        #pragma omp section
                        m_deriche.filter(m_var_pos);
        #pragma omp section
                        m_deriche.filter(m_var_neg);    
}
                }
        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                void LocalStatsRegionModule<DIM, Pixel, PixelImage>::save_result(const LevelSet & phi, std::string evodir, int iter)
        {
                base::save_result(phi, evodir, iter);

                ImageType tmp(phi.shape());
                for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
                        if(phi(i) > 0) tmp(i) = m_data_pos(i)/m_chi_pos(i);
                        else           tmp(i) = m_data_neg(i)/(1-m_chi_pos(i));
                }

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

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

#endif // LOCALSTATSREGIONMODULE_H

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