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

Go to the documentation of this file.

#ifndef GEODESICACTIVECONTOURMODULE_H
#define GEODESICACTIVECONTOURMODULE_H

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

namespace segmentation {


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

                typedef          SegmentationModule<DIM, Pixel>  base;
                typedef typename base::LevelSet                  LevelSet;
                typedef typename base::Index                     Index;
                typedef          DericheFilter<DIM, PixelImage>  DericheFilter;

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


                GeodesicActiveContourModule(ImageType *  _data,             
                                            float        _weight = 1.f,     
                                            float        _v = 0.f,          
                                            float        _smooth = 3.0f     
                                            );

                virtual ~GeodesicActiveContourModule();



                void compute_metric(ImageType * _data, float _smooth = 3.0f);
                void set_metric(const ImageFloat & g, ImageFloat * g_deriv = NULL);
                float& balloon();
                const ImageFloat & metric();



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


        protected:      

                inline float _g(float v, float p)
                {
                        return 1/(1+pow(v,p));
                }


                ImageFloat   m_g;             
                ImageFloat   m_g_deriv[DIM];  
                float        m_v;             


        };



        // Constructors - Destructor


        template <unsigned DIM, typename Pixel, typename PixelImage>
                GeodesicActiveContourModule<DIM, Pixel, PixelImage>::GeodesicActiveContourModule(ImageType * _data, float _weight, float _v, float _smooth)
                :SegmentationModule<DIM, Pixel>(_weight, false), m_g(_data->shape()), m_v(_v)
        { 
                for (int i=0 ; i<DIM ; i++)
                        m_g_deriv[i].resize(_data->shape());

                compute_metric(_data, _smooth);
        }

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


        // Modify/Access member data


        template <unsigned DIM, typename Pixel, typename PixelImage>
                void GeodesicActiveContourModule<DIM, Pixel, PixelImage>::compute_metric(ImageType * _data, float _smooth)
        {
                // smoothed image 
                ImageType image_smooth(_data->shape());
                image_smooth = _data->copy();
                DericheFilter deriche(_smooth);
                deriche.filter(image_smooth);

                // metric computation : g = 1 / (1 + sqrt(|grad I|) )
                if (m_g.shape() != _data->shape()) m_g.resize(_data->shape());
                for (typename ImageFloat::template iterator<domain> ind=m_g.begin() ; ind!=m_g.end() ; ++ind )
                        m_g(ind) = _g(Schemes<DIM,PixelImage>::norm_gradient(image_smooth, ind.position()), 0.5f);

                // -0.2 to assure a range of value where the metric is 0
                Globals::normalize<DIM, PixelImage>(m_g, -0.2f, 1);
                Globals::cut<DIM, PixelImage>(m_g, 0, 1);

                // directional derivative of the metric
                for (int i=0 ; i<DIM ; i++)
                        if (m_g_deriv[i].shape() != _data->shape()) m_g_deriv[i].resize(_data->shape());

                PixelImage grad[DIM];
                for (typename ImageFloat::template iterator<domain> ind=m_g.begin() ; ind!=m_g.end() ; ++ind ){
                        Schemes<DIM, float>::gradient(m_g, ind.position(), grad);
                        for (int i=0 ; i<DIM ; i++)
                                m_g_deriv[i](ind.position()) = grad[i];
                }

        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                void GeodesicActiveContourModule<DIM, Pixel, PixelImage>::set_metric(const ImageFloat & g, ImageFloat * g_deriv)
        {
                m_g = g;

                if (g_deriv != NULL)
                        for (int i=0 ; i<DIM ; i++)
                                m_g_deriv[i] = g_deriv[i];
                        
                else {
                        // directional derivative of the metric
                        for (int i=0 ; i<DIM ; i++)
                                if (m_g_deriv[i].shape() != m_g.shape()) m_g_deriv[i].resize(m_g.shape());

                        PixelImage grad[DIM];
                        for (typename ImageFloat::template iterator<domain> ind=m_g.begin() ; ind!=m_g.end() ; ++ind ){
                                Schemes<DIM, float>::gradient(m_g, ind.position(), grad);
                                for (int i=0 ; i<DIM ; i++)
                                        m_g_deriv[i](ind.position()) = grad[i];
                        }
                }
        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                float& GeodesicActiveContourModule<DIM, Pixel, PixelImage>::balloon()
        {
                return m_v;
        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                const typename GeodesicActiveContourModule<DIM, Pixel, PixelImage>::ImageFloat & GeodesicActiveContourModule<DIM, Pixel, PixelImage>::metric()
        {
                return m_g;
        }


        // Member functions


        template <unsigned DIM, typename Pixel, typename PixelImage>
                float GeodesicActiveContourModule<DIM, Pixel, PixelImage>::value(const LevelSet & phi, const Index & ind)
        {
                Pixel grad[DIM];
                Pixel gnorm = 0;

                Schemes<DIM,Pixel>::gradient(phi, ind, grad);
                for (int i=0 ; i<DIM ; i++)
                        gnorm += Maths::sqr(grad[i]);
                gnorm = sqrt(gnorm);

                float result = m_g(ind)*gnorm*(Schemes<DIM,Pixel>::mean_curvature(phi,ind) - m_v);
                for (int i=0 ; i<DIM ; i++)
                        result += grad[i]*m_g_deriv[i](ind);

                return result;
        }        

        template <unsigned DIM, typename Pixel, typename PixelImage>
                float GeodesicActiveContourModule<DIM, Pixel, PixelImage>::energy(const LevelSet & phi, const Index & ind)
        {
                return m_g(ind) * Globals::delta(phi(ind));
        }

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

                if (m_g.shape() != phi.shape()) {
                        std::cerr << "Error : Shape of the Level Set and shape of the metric don't match"  << std::endl;
                        exit(1);
                }
        }

        template <unsigned DIM, typename Pixel, typename PixelImage>
                void GeodesicActiveContourModule<DIM, Pixel, PixelImage>::save_result(const LevelSet & phi, std::string evodir, int iter)
        {
                if (iter == 0){
                        std::ostringstream name;
                        name.str("");
                        name << evodir << "Metric.inr";
                        std::ofstream file(name.str().c_str());
                        file << format("inrimage") << m_g ;
                        std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                        for (int i=0 ; i<DIM ; i++){
                                name.str("");
                                name << evodir << "Metric_" << i << ".inr";
                                std::ofstream file(name.str().c_str());
                                file << format("inrimage") << m_g_deriv[i] ;
                                std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
                        }
                }
        }

}

#endif // GEODESICACTIVECONTOURMODULE_H

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