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

Go to the documentation of this file.

#ifndef PARZENREGIONMODULE_H
#define PARZENREGIONMODULE_H

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

namespace segmentation {


        template <unsigned DIM, typename Pixel, typename PixelImage = float> 
        class ParzenRegionModule : 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          DericheFilter<1, float>                 DericheFilter;


                ParzenRegionModule(ImageType * _data,                
                                   float       _weight = 1.f,        
                                   float       _sigma = 4,           
                                   bool        _adaptive = true      
                                   );

                virtual ~ParzenRegionModule();



                void set_sigma(float _sigma);



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



        protected:

                float                 m_sigma;       
                DericheFilter         m_deriche;     
                std::vector<Pixel>    m_ll_pos;      
                std::vector<Pixel>    m_ll_neg;      
                Image1D<float>        m_parzen_pos;  
                Image1D<float>        m_parzen_neg;  
                
        };



        // Constructors - Destructor


        template <unsigned DIM, typename Pixel, typename PixelImage>
                ParzenRegionModule<DIM, Pixel, PixelImage>::ParzenRegionModule(ImageType * _data, float _weight, float _sigma, bool _adaptive)
                :RegionModule<DIM, Pixel, PixelImage>(_data, _weight, _adaptive), m_sigma(_sigma), m_deriche(_sigma)
        { }

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


        // Modify/Access member data 


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


        // Member functions


        template <unsigned DIM, typename Pixel, typename PixelImage>
                float ParzenRegionModule<DIM, Pixel, PixelImage>::value(const LevelSet & phi, const Index & ind)
        {
                PixelImage val = (* base::m_data)(ind);
                const Pixel ll_pos = m_ll_pos[base::m_hist_pos->get_index(val)];
                const Pixel ll_neg = m_ll_neg[base::m_hist_neg->get_index(val)];

                return  Globals::delta(phi(ind)) * (ll_pos - ll_neg);
        }       

        template <unsigned DIM, typename Pixel, typename PixelImage>
                float ParzenRegionModule<DIM, Pixel, PixelImage>::energy(const LevelSet & phi, const Index & ind)
        {
                PixelImage val = (* base::m_data)(ind);
                float ll_pos = m_ll_pos[base::m_hist_pos->get_index(val)];
                float ll_neg = m_ll_neg[base::m_hist_neg->get_index(val)];
                float h = Globals::heaviside(phi(ind));

                return -(1-h)*ll_pos - h*ll_neg;
        }

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

        template <unsigned DIM, typename Pixel, typename PixelImage>
                void ParzenRegionModule<DIM, Pixel, PixelImage>::update(const LevelSet & phi)
        { 
                base::update_histograms(phi);
                update_parzen_params();
        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                void ParzenRegionModule<DIM, Pixel, PixelImage>::update_parzen_params()
        { 
                // if the size of the histogram has changed, resize parzen densities and log-likelihood
                if (m_parzen_neg.size() != base::hist_neg()->size()){
                        m_parzen_neg.resize(base::hist_neg()->size());
                        m_ll_neg.resize(base::hist_neg()->size());
                }
                if (m_parzen_pos.size() != base::hist_pos()->size()){
                        m_parzen_pos.resize(base::hist_pos()->size());
                        m_ll_pos.resize(base::hist_pos()->size());
                }

                // initialize the parzen densities from the histograms
                for(int i=0; i < m_parzen_neg.size() ; i++) m_parzen_neg(i) = base::hist_neg()->data()[i];
                for(int i=0; i < m_parzen_pos.size() ; i++) m_parzen_pos(i) = base::hist_pos()->data()[i];

                // smooth the parzen densities
                m_deriche.filter(m_parzen_neg);
                m_deriche.filter(m_parzen_pos);

                // normalize the parzen densities
                float sum_neg=0.0f;
                float sum_pos=0.0f;
                for(int i=0; i < m_parzen_neg.size() ; i++) sum_neg += m_parzen_neg(i);
                for(int i=0; i < m_parzen_pos.size() ; i++) sum_pos += m_parzen_pos(i);
                for(int i=0; i < m_parzen_neg.size() ; i++) m_parzen_neg(i) /= sum_neg;
                for(int i=0; i < m_parzen_pos.size() ; i++) m_parzen_pos(i) /= sum_pos;

                // log likelihood computation
                for(int i=0; i < m_ll_neg.size(); i++)  m_ll_neg[i] = log(m_parzen_neg(i)+Globals::EPS);
                for(int i=0; i < m_ll_pos.size(); i++)  m_ll_pos[i] = log(m_parzen_pos(i)+Globals::EPS);
        }

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

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

                name.str("");
                name << evodir << "Parzen_pos_" << std::setw(6) << std::setfill('0') << iter << ".txt";
                file.open(name.str().c_str());
                for(int i=0; i < m_parzen_pos.size() ; i++) 
                        file << i + base::m_hist_pos->min() << " " << m_parzen_pos(i) << std::endl;
                file.close();
                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;

                name.str("");
                name << evodir << "Parzen_neg_" << std::setw(6) << std::setfill('0') << iter << ".txt";
                file.open(name.str().c_str());
                for(int i=0; i < m_parzen_neg.size() ; i++) 
                        file << i + base::m_hist_neg->min()<< " " << m_parzen_neg(i) << std::endl;
                file.close();
                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
        }
}

#endif // PARZENREGIONMODULE_H

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