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

Go to the documentation of this file.

#ifndef LEVELSET_H
#define LEVELSET_H

#include <cmath>
#include <Image.H>
#include <Initializer.h>

#include <Fluxes.h>

#define _USE_MATH_DEFINES 
 
using namespace Images;
using namespace levelset;


namespace levelset {

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

                typedef          BaseImage<DIM, Pixel>             base;
                typedef typename base::Shape                       Shape;
                typedef typename base::Index                       Index;
                typedef          Initializer<DIM, Pixel>           Initializer;

//TODO type de l'approximateur fixe en dur dans le code . autre idee ???

                typedef          Weno3<DIM, Pixel>                 Approximator;
                typedef          Fluxes<DIM, Approximator, Pixel>  Fluxes;
                typedef          Schemes<DIM, Pixel>               Schemes;


                LevelSet();

                LevelSet(const Shape &  _shape,              
                         float          _thickness = 3.0f    
                         );

                template <typename Pixel2>
                LevelSet(const LevelSet<DIM, Pixel2>& ls):base(ls), m_thickness(ls.thickness()), m_mask(ls.shape()), m_nbpts(0), m_pos_to_neg(0), m_neg_to_pos(0)
                {
                        m_mask = -1;
                        m_band = new Index[ls.size()];
                        for (int i=0 ; i<ls.nb_points() ; i++)
                                add_point(ls.get_point(i));
                }

                template <typename Pixel2>
                LevelSet operator=(const LevelSet<DIM, Pixel2>& ls)
                {
                        if (base::shape != ls.shape())
                                resize(ls.shape());
                        
                        base::operator=(ls);
                        m_thickness = ls.thickness();
                        m_nbpts = 0;
                        for (int i=0 ; i<ls.nb_points() ; i++)
                                add_point(ls.get_point(i));
                        return *this;
                }

                virtual ~LevelSet();

                using base::operator=;

                void resize(const Shape& s);



                void update(const BaseImage<DIM, Pixel> &dphi, const float dt);



                bool is_in_band(const Index& ind) const;
                bool add_point(const Index& ind);
                bool remove_point(const Index& ind);
                const Index& get_point(unsigned i) const;
                int nb_points() const;
                int pos_to_neg_index() const;
                int neg_to_pos_index() const;
                float & thickness();
                void clear_band();



                void init_sphere(const Index& center, float radius, Initializer& init, bool add=true);
                void init_small_spheres(float radius, float window, Initializer& init, bool add=true);
                void init_border(const unsigned border, Initializer& init, bool add=true);
                void init_rectangle(const Index& tl, const Index& br, Initializer& init, bool add=true);
                void init_image(const BaseImage<DIM, Pixel>& image, Initializer& init, const Pixel& pix = 0);



//TODO accelerer reinit pde . globalement reinit pde assez lent ... gprof dis que ca viens de blitz tiny vector .... peut etre essayer d implementer avec des tableaux plutot.

                void reinit_PDE(int iter = 2, float dt = 0.4);
                void extend_data(BaseImage<DIM, Pixel>& u, short type = 0, int iter = 20, float dt = 0.4) const;
                Pixel band_average(const BaseImage<DIM, Pixel>& u, const float epsilon = 2.0) const;
                Pixel band_integral(const BaseImage<DIM, Pixel>& u, const float epsilon = 2.0) const;
                Pixel area(const float epsilon = 2.0) const;

        protected:



                float                            m_thickness;   
                Index *                          m_band;        
                BaseImage<DIM, int>              m_mask;        
                unsigned                         m_nbpts;       
                unsigned                         m_pos_to_neg;  
                unsigned                         m_neg_to_pos;  

        };



        // Constructors - Destructor - Copy

        template<unsigned DIM, typename Pixel>
                        LevelSet<DIM, Pixel>::LevelSet()
        : base(), m_thickness(0), m_band(NULL), m_mask(), m_nbpts(0), m_pos_to_neg(0), m_neg_to_pos(0)
        { }

        template<unsigned DIM, typename Pixel>
                        LevelSet<DIM, Pixel>::LevelSet(const Shape &  _shape, float _thickness)
        : base(_shape), m_thickness(_thickness), m_mask(_shape), m_nbpts(0), m_pos_to_neg(0), m_neg_to_pos(0)
        { 
                m_mask = -1;
                base::operator=(-m_thickness);
                m_band = new Index[base::size()];
        }

        template<unsigned DIM, typename Pixel>
                LevelSet<DIM, Pixel>::~LevelSet() 
        { 
                if (m_band != NULL) delete[] m_band;
        }

        template<unsigned DIM, typename Pixel>
                void LevelSet<DIM, Pixel>::resize(const Shape& s)
        {
                base::resize(s);
                m_mask.resize(s);
                base::operator=(-m_thickness);
                m_mask = -1;
                if (m_band != NULL) delete[] m_band;
                m_band = new Index[base::size()];
                m_nbpts = 0;
        }


        // Update of the narrow band


        // note : after the update, the array 'm_band' is organized as follow :
        //
        //          --------------------------------------------------------------------------------------
        //  m_band  |      pos->neg      |                  other                   |     neg->pos       |
        //          --------------------------------------------------------------------------------------
        //          ^                    ^                                         ^                    ^
        //          0               m_pos_to_neg                             m_neg_to_pos            m_nbpts
        //
        // where,
        //    - pos->neg [0,            m_pos_to_neg-1] : points where the sign of phi changed from positive to negative.
        //    - neg->pos [m_pos_to_neg, m_neg_to_pos-1] : points where the sign of phi didn't change
        //    - other    [m_neg_to_pos,      m_nbpts-1] : points where the sign of phi changed from negative to positive
        // 
        // cf : "The Problem of the Dutch National Flag", Dijkstra

        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::update(const BaseImage<DIM, Pixel> &dphi, const float dt)
        {
                m_pos_to_neg = 0;
                m_neg_to_pos = m_nbpts;

                unsigned int i = 0;
                while (i < m_neg_to_pos)
                {
                        const Index& ind = get_point(i);

                        Pixel tmp = base::operator()(ind);

                        base::operator()(ind) += dt * dphi(ind);

                        // keep the zero level set in the band
                        if (m_thickness > 0) {
                                base::operator()(ind) = (base::operator()(ind) >  m_thickness) ?  m_thickness : base::operator()(ind);
                                base::operator()(ind) = (base::operator()(ind) < -m_thickness) ? -m_thickness : base::operator()(ind);
                        }

                        // sign changes
                        if(tmp * base::operator()(ind) < 0)
                        {  
                                //pos->neg
                                if (tmp>0){
                                        std::swap(m_band[m_pos_to_neg], m_band[i]);
                                        std::swap(m_mask(m_band[m_pos_to_neg]), m_mask(m_band[i]));
                                        m_pos_to_neg++;
                                        i++;
                                }
                                //neg->pos
                                else if(tmp<0) {
                                        m_neg_to_pos--;
                                        std::swap(m_band[m_neg_to_pos], m_band[i]);
                                        std::swap(m_mask(m_band[m_neg_to_pos]), m_mask(m_band[i]));
                                }
                        }
                        //others
                        else i++;
                }
        }
        

        //  Access/Modif of the narrow band 
        

        template<unsigned DIM, typename Pixel>
        bool LevelSet<DIM, Pixel>::is_in_band(const Index& ind) const
        {
                return (m_mask(ind) != -1);
        }
                
        template<unsigned DIM, typename Pixel>
        bool LevelSet<DIM, Pixel>::add_point(const Index& ind)
        {
                if (!is_in_band(ind)){

                        // we put the indice of the point in the mask of indice
                        m_mask(ind) = m_nbpts;

                        // we add the point at the end of m_band
                        m_band[m_nbpts++] = ind;
                        return true;
                }
                return false;
        }

        template<unsigned DIM, typename Pixel>
        bool LevelSet<DIM, Pixel>::remove_point(const Index& ind)
        {
                if (is_in_band(ind)){
                        const int position = m_mask(ind);

                        // we remove the point from the band
                        m_mask(ind) = -1;
                        m_nbpts--;

                        // we put the last point at its place
                        if (position != m_nbpts){
                                m_band[position] = m_band[m_nbpts];
                                m_mask(m_band[position]) = position;
                        }
                        return true;
                }
                return false;
        }
                
        template<unsigned DIM, typename Pixel>
        const typename LevelSet<DIM, Pixel>::Index& LevelSet<DIM, Pixel>::get_point(unsigned i) const
        {
                return m_band[i];
        }
                
        template<unsigned DIM, typename Pixel>
        int LevelSet<DIM, Pixel>::nb_points() const
        {
                return m_nbpts;
        }       

        template<unsigned DIM, typename Pixel>
        int LevelSet<DIM, Pixel>::pos_to_neg_index() const
        {
                return m_pos_to_neg;
        }       
                
        template<unsigned DIM, typename Pixel>
        int LevelSet<DIM, Pixel>::neg_to_pos_index() const
        {
                return m_neg_to_pos;
        }       
                
        template<unsigned DIM, typename Pixel>
        float & LevelSet<DIM, Pixel>::thickness() 
        {
                return m_thickness;
        }
        
        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::clear_band()
        {
                m_nbpts = 0;
                m_mask = -1;
        }

        
        // Initializers
        

        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::init_sphere(const Index& center, float radius, Initializer& init, bool add)
        {
                for (typename base::template iterator<domain> i=this->begin() ; i!=this->end() ; ++i ){

                        //compute the squared distance of the point 'i' to the point 'center'
                        float dist = 0;
                        for (unsigned dim=1 ; dim<=base::dimension() ; dim++)
                                dist += Maths::sqr(i(dim)-center(dim));

                        dist = radius - sqrt(dist);

                        if (add) base::operator()(i) = ((dist>=0) || (base::operator()(i)>0))? 1 : -1;
                        else     base::operator()(i) = ((dist< 0) && (base::operator()(i)>0))? 1 : -1;
                }

                init.init_band(*this);
        }

        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::init_small_spheres(float radius, float window, Initializer& init, bool add)
        {
                for (typename base::template iterator<domain> i=this->begin() ; i!=this->end() ; ++i ){

                        //compute the squared distance of the point 'i % windows' to the point of coordinate 'windows / 2'
                        float dist = 0;
                        for (unsigned dim=1 ; dim<=base::dimension() ; dim++)
                                dist += Maths::sqr(i(dim)%int(window)-window/2);

                        dist = radius - sqrt(dist);

                        if (add) base::operator()(i) = ((dist>=0) || (base::operator()(i)>0))? 1 : -1;
                        else     base::operator()(i) = ((dist< 0) && (base::operator()(i)>0))? 1 : -1;
                }

                init.init_band(*this);
        }
        
        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::init_border(const unsigned border, Initializer& init, bool add)
        {
                for (typename base::template iterator<domain> i=this->begin() ; i!=this->end() ; ++i ){

                        //inside=true if the point 'i' inside the region
                        bool inside = true;
                        for (unsigned dim=1 ; dim<=base::dimension() ; dim++)
                                inside &= (i(dim) >= border) && (i(dim) <= base::extent(dim-1) - border);

                        if (add) base::operator()(i) = ( inside || (base::operator()(i)>0))? 1 : -1;
                        else     base::operator()(i) = (!inside && (base::operator()(i)>0))? 1 : -1;
                }

                init.init_band(*this);
        }
        
        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::init_rectangle(const Index& tl, const Index& br, Initializer& init, bool add)
        {
                for (typename base::template iterator<domain> i=this->begin() ; i!=this->end() ; ++i ){

                        //inside=true if the point 'i' inside the region
                        bool inside = true;
                        for (unsigned dim=1 ; dim<=base::dimension() ; dim++)
                                inside &= (i(dim) >= tl(dim)) && (i(dim) <= br(dim));

                        if (add) base::operator()(i) = ( inside || (base::operator()(i)>0))? 1 : -1;
                        else     base::operator()(i) = (!inside && (base::operator()(i)>0))? 1 : -1;
                }

                init.init_band(*this);
        }

        template<unsigned DIM, typename Pixel>          
        void LevelSet<DIM, Pixel>::init_image(const BaseImage<DIM, Pixel>& image, Initializer& init, const Pixel& pix)
        {
                resize(image.shape());
                
                operator=(image);
                
                for (typename base::template iterator<domain> i=this->begin() ; i!=this->end() ; ++i ){
                        base::operator()(i) = (base::operator()(i) > pix)? 1 : -1;
                }

                init.init_band(*this);
        }
        

        // Auxiliary functions


        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::reinit_PDE(int iter, float dt)
        {
                // note : pb si sur stack et que temp trop gros, foire .. donc alloc/desalloc explicite.
                Pixel* temp = new Pixel[nb_points()];
                        
                for (int t=0;t<iter;t++) {      
                        #pragma omp parallel for
                        for (int i = 0 ; i < nb_points() ; i++){
                                const Index& ind = get_point(i);
                                temp[i] = Fluxes::reinitialization(*this,ind);
                        }

                        #pragma omp parallel for
                        for (int i = 0 ; i < nb_points() ; i++){
                                const Index& ind = get_point(i);
                                base::operator()(ind) += dt*temp[i];
                        }
                }

                delete[] temp;
        }

        template<unsigned DIM, typename Pixel>
        void LevelSet<DIM, Pixel>::extend_data(BaseImage<DIM, Pixel>& u, short type, int iter, float dt) const
        {
                Pixel temp[nb_points()];

                for (int t=0;t<iter;t++) {
                        #pragma omp parallel for
                        for (int i = 0 ; i < nb_points() ; i++){
                                const Index& ind = get_point(i);
                                switch(type){
                                        case  1 : temp[i] = (base::operator()(ind) >  0) ? Fluxes::extension(*this,u,ind) : 0; break;
                                        case -1 : temp[i] = (base::operator()(ind) <= 0) ? Fluxes::extension(*this,u,ind) : 0; break;
                                        default : temp[i] = Fluxes::extension(*this,u,ind); break;
                                }
                        }
                                
                        #pragma omp parallel for
                        for (int i = 0 ; i < nb_points() ; i++){
                                const Index& ind = get_point(i);
                                u(ind)  += dt * temp[i];
                        }
                }
        }

        template<unsigned DIM, typename Pixel>
        Pixel LevelSet<DIM, Pixel>::band_average(const BaseImage<DIM, Pixel>& u, const float epsilon) const
        {
                Pixel num = 0;
                Pixel den = 0;

                for (int i = 0 ; i < nb_points() ; i++){
                        const Index& ind = get_point(i);
                        const Pixel  val = base::operator()(ind);

                        if (fabs(val) < epsilon) {
                                const Pixel tmp = Globals::delta(val, epsilon) * Schemes::norm_gradient(*this, ind);
                                den += tmp;
                                num += tmp * u(ind);
                        }
                }

                if (den < Globals::EPS) return 0;

                return num / den;
        }

        template<unsigned DIM, typename Pixel>
        Pixel LevelSet<DIM, Pixel>::band_integral(const BaseImage<DIM, Pixel>& u, const float epsilon) const
        {
                Pixel sum = 0;

                for (int i = 0 ; i < nb_points() ; i++){
                        const Index& ind = get_point(i);
                        const Pixel  val = base::operator()(ind);

                        if (fabs(val) < epsilon)
                                sum += u(ind) * Globals::delta(val, epsilon) * Schemes::norm_gradient(*this, ind);
                }

                return sum;
        }

        template<unsigned DIM, typename Pixel>
        Pixel LevelSet<DIM, Pixel>::area(const float epsilon) const
        {
                Pixel sum = 0;
                
                for (int i = 0 ; i < nb_points() ; i++){
                        const Index& ind = get_point(i);
                        const Pixel  val = base::operator()(ind);

                        if (fabs(val) < epsilon) 
                                sum += Globals::delta(val, epsilon) * Schemes::norm_gradient(*this, ind);
                }

                return sum;
        }
}


#endif  // LEVELSET_H

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