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

Go to the documentation of this file.

#ifndef INITIALIZER_H
#define INITIALIZER_H

#ifdef _USE_INTEL_COMPILER
#include <omp.h>
#endif

#include <LevelSet.h>
#include <FastMarching.h>
#include <Schemes.h>
#include <Globals.h>

namespace levelset {

        template <unsigned DIM, typename Pixel> class LevelSet;

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

                typedef  LevelSet<DIM, Pixel>      LevelSetType;


                Initializer(){}
                virtual ~Initializer(){}



                virtual void init_band(LevelSetType & phi) = 0 ;
                virtual void reinit_band(LevelSetType & phi) = 0 ;
        };

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

                typedef           Initializer<DIM, Pixel>       base;
                typedef typename  base::LevelSetType            LevelSetType;
                typedef typename  LevelSetType::Index           Index;
                typedef           Schemes<DIM,Pixel>            Schemes;


                FastMarchingInitializer(int   _iter = 2,   // # of iterations of the reinitialization PDE
                                        float _dt = 0.4    // time step of the reinitialization PDE
                                       );
                virtual ~FastMarchingInitializer();



                void init_band(LevelSetType & phi);
                void reinit_band(LevelSetType & phi);
                int& iter();
                float& dt();

        protected:

                void _add_alive_point(LevelSetType & phi, const Index & ind);
                void _reinit_fast_marching(LevelSetType & phi);
        private:



                int                            m_iter;     
                float                          m_dt;       
                Eikonal<DIM, Pixel, +1>      * m_reinit1;  
                Eikonal<DIM, Pixel, -1>      * m_reinit2;  

        };



        // FastMarching Initializer 


        template<unsigned DIM, typename Pixel>
        FastMarchingInitializer<DIM, Pixel>::FastMarchingInitializer(int _iter, float _dt)
                : base(), m_iter(_iter), m_dt(_dt)
        {
                m_reinit1 = new Eikonal<DIM, Pixel, +1>();
                m_reinit2 = new Eikonal<DIM, Pixel, -1>();
        }

        template<unsigned DIM, typename Pixel>
        FastMarchingInitializer<DIM, Pixel>::~FastMarchingInitializer()
        {
                delete m_reinit1;
                delete m_reinit2;
        }

        template<unsigned DIM, typename Pixel>
        void FastMarchingInitializer<DIM, Pixel>::init_band(LevelSetType & phi)
        {
                // Initialization of the fast marching datas
                m_reinit1->init(phi);
                m_reinit2->init(phi);

                // clear the band
                phi.clear_band();

                #ifdef PRINT_TIME
                timeval _time;
                Globals::initChrono(_time);
                #endif

                // we add the points neighboring the zero level set to the alive points
                int inf_coords[DIM];
                int sup_coords[DIM];
                for (typename LevelSetType::template iterator<domain> ind=phi.begin() ; ind!=phi.end() ; ++ind ){

                        const Pixel _phi = phi(ind);
                        bool  add        = false;

                        // find the coordinates neighboring the point 'ind'
                        Schemes::find_neighbours(phi, ind.position(), inf_coords, sup_coords);

                        //for each point 'inf' and 'sup' neighboring 'ind' ...
                        for (int i=0 ; i<DIM ; i++){
                                Index inf = ind.position();
                                Index sup = ind.position();
                                inf(i+1) = inf_coords[i];
                                sup(i+1) = sup_coords[i];
                                if (_phi * phi(inf) <= 0) { add = true; _add_alive_point(phi,inf); }
                                if (_phi * phi(sup) <= 0) { add = true; _add_alive_point(phi,sup); }
                        }

                        if (add) { _add_alive_point(phi,ind.position());}
                }

                #ifdef PRINT_TIME
                Globals::printChrono("          - points neighboring the zero level set : ", _time);
                Globals::initChrono(_time);
                #endif

                // We add to the band the points alive (for the reinit EDP which is done first).
                // we set the value of the others points to '1' ( border condition ) (2 if alive+trial too)
                for (typename LevelSetType::template iterator<domain> ind=phi.begin() ; ind!=phi.end() ; ++ind ){
                        const Pixel _phi = phi(ind);
                        if (_phi>0) {
                                if ((*m_reinit1)(ind) == 2) phi(ind) =  1;
                                else phi.add_point(ind.position());
                        }
                        else {
                                if ((*m_reinit2)(ind) == 2) phi(ind) = -1;
                                else phi.add_point(ind.position());
                        }
                }

                #ifdef PRINT_TIME
                Globals::printChrono("          - add to the band the points alive: ", _time);
                Globals::initChrono(_time);
                #endif

                // reinitialization PDE (2*iter because first initialisation)
                phi.reinit_PDE(m_iter,m_dt);

                #ifdef PRINT_TIME
                Globals::printChrono("          - reinitialization PDE: ", _time);
                Globals::initChrono(_time);
                #endif

                // init of the trial points from the alive points
                m_reinit1->init_trial_from_alive();
                m_reinit2->init_trial_from_alive();

                #ifdef PRINT_TIME
                Globals::printChrono("          - init of the trial points from the alive points : ", _time);
                Globals::initChrono(_time);
                #endif

                // Fast marching
                _reinit_fast_marching(phi);

                #ifdef PRINT_TIME
                Globals::printChrono("          - FastMarching : ", _time);
                Globals::initChrono(_time);
                #endif

                // we put the other points at the max value
                for (typename LevelSetType::template iterator<domain> ind=phi.begin() ; ind!=phi.end() ; ++ind ){
                        const Pixel _phi = phi(ind);
                        if      (_phi>0 && (*m_reinit1)(ind)==2) phi(ind) =  phi.thickness();
                        else if (_phi<0 && (*m_reinit2)(ind)==2) phi(ind) = -phi.thickness();
                }

                #ifdef PRINT_TIME
                Globals::printChrono("          - put the other points at the max value : ", _time);
                #endif

        }

        template<unsigned DIM, typename Pixel>
        void FastMarchingInitializer<DIM, Pixel>::reinit_band(LevelSetType & phi)
        {

                #ifdef PRINT_TIME
                timeval _time;
                Globals::initChrono(_time);
                #endif

                // reinitialization EDP 
                phi.reinit_PDE(m_iter,m_dt);

                // Initialization of the fast marching datas
                m_reinit1->init(phi);
                m_reinit2->init(phi);

                #ifdef PRINT_TIME
                Globals::printChrono("          - reinitialization PDE : ", _time);
                Globals::initChrono(_time);
                #endif

                // we add the points neighboring the zero level set to the alive points
                int inf_coords[DIM];
                int sup_coords[DIM];
                for (unsigned int i = 0 ; i < phi.nb_points() ; i++){

                        const Index& ind = phi.get_point(i);
                        const Pixel _phi = phi(ind);
                        bool  add        = false;

                        // find the coordinates neighboring the point 'ind'
                        Schemes::find_neighbours(phi, ind, inf_coords, sup_coords);

                        //for each point 'inf' and 'sup' neighboring 'ind' ...
                        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 * phi(inf) <= 0) { add = true; _add_alive_point(phi,inf); }
                                if (_phi * phi(sup) <= 0) { add = true; _add_alive_point(phi,sup); }
                        }

                        if (add) { _add_alive_point(phi,ind);}
                }

                #ifdef PRINT_TIME
                Globals::printChrono("          - points neighboring the zero level set : ", _time);
                Globals::initChrono(_time);
                #endif

                // trial points are neighbor of alive points
                m_reinit1->init_trial_from_alive();
                m_reinit2->init_trial_from_alive();

                #ifdef PRINT_TIME
                Globals::printChrono("          - init of the trial points from the alive points : ", _time);
                Globals::initChrono(_time);
                #endif

                // we remove from the band the points that are not trial nor alive
                // we put their value at max
                int i = 0;
                while (i<phi.nb_points()){
                        const Index& ind = phi.get_point(i);
                        const Pixel _phi = phi(ind);

                        if (_phi>0 && (*m_reinit1)(ind) == 2){
                                phi(ind) = phi.thickness();
                                phi.remove_point(ind);
                        }
                        else if  (_phi<=0 && (*m_reinit2)(ind) == 2) {
                                phi(ind) = -phi.thickness();
                                phi.remove_point(ind);
                        }
                        else i++;
                }

                #ifdef PRINT_TIME
                Globals::printChrono("          - remove from the band the points not alive or trial : ", _time);
                #endif

                // Fast marching
                _reinit_fast_marching(phi);

                #ifdef PRINT_TIME
                Globals::printChrono("          - FastMarching : ", _time);
                #endif
        }

        template<unsigned DIM, typename Pixel>
        int& FastMarchingInitializer<DIM, Pixel>::iter()
        {
                return m_iter;
        }

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

        template<unsigned DIM, typename Pixel>
        void FastMarchingInitializer<DIM, Pixel>::_add_alive_point(LevelSetType & phi, const Index & ind)
        {
                if (phi(ind)>0) {
                        m_reinit1->add_alive_point(ind);
                        m_reinit2->add_forbidden_point(ind);
                }
                else {
                        m_reinit1->add_forbidden_point(ind);
                        m_reinit2->add_alive_point(ind);
                }       
        }

        template<unsigned DIM, typename Pixel>
        void FastMarchingInitializer<DIM, Pixel>::_reinit_fast_marching(LevelSetType & phi) 
        {
                // we launch the two fast marching
                if ( phi.thickness()>0 ) {
#pragma omp parallel sections
{
        #pragma omp section
                        m_reinit1->run( phi.thickness());
        #pragma omp section
                        m_reinit2->run(-phi.thickness());
}
                }
                else {
#pragma omp parallel sections
{
        #pragma omp section
                        m_reinit1->run();
        #pragma omp section
                        m_reinit2->run();
}
                }
                
                // we add the alive points to the band
                for (unsigned int i = 0 ; i < m_reinit1->nb_alive_points() ; i++)
                        phi.add_point(m_reinit1->get_alive_point(i));

                for (unsigned int i = 0 ; i < m_reinit2->nb_alive_points() ; i++)
                        phi.add_point(m_reinit2->get_alive_point(i));
        }
}

#endif // INITIALIZER_H

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