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

Go to the documentation of this file.

#ifndef DERICHEFILTER_H
#define DERICHEFILTER_H

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

#include <Image.H>

using namespace Images;

namespace Images {

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

                        typedef          BaseImage<DIM, Pixel>             Image;
                        typedef typename Image::Index                      Index;
                        typedef          BaseImage<DIM, float>             ImageVariance;
                        typedef          blitz::RectDomain<DIM>            RectDomain;



                        DericheFilter(float _sigma,               
                                      int   _order = 0,           
                                      bool  _border_cond = true   
                                      );

                        DericheFilter(float _sigma[DIM],          
                                      int   _order = 0,           
                                      bool  _border_cond = true   
                                      );

                        DericheFilter(ImageVariance& _sigma,               
                                      int            _order = 0,           
                                      bool           _border_cond = true   
                                      );

                        ~DericheFilter();



                        void filter(Image& image);
                        void filter_axis(Image& image, const int axis);



                        void set_filter_type(int _filter_type);
                        void set_order(int _order);
                        void set_border_cond(bool _border_cond);

                        void set_isotropic_sigma(float _sigma);
                        void set_anisotropic_sigma(float _sigma[DIM]);
                        void set_image_sigma(const ImageVariance& _sigma);

                protected:



                        void _init_coefs();     
                        void _coefs_sigma(float sigma, float coefs[6],  float border_coefs[4]);

                        void _deriche_map(Image line, const int axis, Pixel* tmp_buf, Index lbound, Index ubound);




                        int              m_filter_type;                   
                        int              m_order;                         
                        bool             m_border_cond;                   

                        float            m_sigma[DIM];                    
                        ImageVariance *  m_sigma_image;                   

                        float            m_coefs[6*DIM];                  
                        ImageVariance *  m_variant_coefs[6];              

                        float            m_border_coefs[4*DIM];           


                };



                // Constructors - Destructor


                template<unsigned DIM, typename Pixel>
                        DericheFilter<DIM, Pixel>::DericheFilter(float _sigma, int _order, bool _border_cond)
                        :m_filter_type(0), m_order(_order), m_sigma_image(NULL), m_border_cond(_border_cond)
                {
                        for (int i=0 ; i<DIM; i++) 
                                m_sigma[i] = _sigma;

                        for (int i=0 ; i<6; i++) 
                                m_variant_coefs[i] = NULL;

                        _init_coefs();
                }

                template<unsigned DIM, typename Pixel>
                        DericheFilter<DIM, Pixel>::DericheFilter(float _sigma[DIM], int _order, bool _border_cond)
                        :m_filter_type(1), m_order(_order), m_sigma_image(NULL), m_border_cond(_border_cond)
                { 
                        for (int i=0 ; i<DIM; i++) 
                                m_sigma[i] = _sigma[i];

                        for (int i=0 ; i<6; i++) 
                                m_variant_coefs[i] = NULL;

                        _init_coefs();
                }

                template<unsigned DIM, typename Pixel>
                        DericheFilter<DIM, Pixel>::DericheFilter(ImageVariance& _sigma, int _order, bool _border_cond)
                        :m_filter_type(2), m_order(_order), m_sigma_image(&_sigma), m_border_cond(_border_cond)
                {
                        for (int i=0 ; i<DIM; i++)
                                m_sigma[i] = 0.0;

                        for (int i=0 ; i<6; i++) 
                                m_variant_coefs[i] = NULL;

                        _init_coefs();
                }

                template<unsigned DIM, typename Pixel>
                        DericheFilter<DIM, Pixel>::~DericheFilter()
                {
                        for (int i=0 ; i<6; i++)
                                if (m_variant_coefs[i] != NULL)
                                         delete m_variant_coefs[i];
                }


                // Main Method


                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::filter(Image& image)
                {
                        for (int i=1 ; i<=DIM ; i++)
                                filter_axis(image, i);
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::filter_axis(Image& image, const int axis)
                {
                        // get the first slice in the given axis
                        Index tl = image.lbound();
                        Index br = image.ubound(); 
                        tl(axis) = 0;
                        br(axis) = 0;
                                
                        blitz::RectDomain<DIM> rect_domain(tl, br);
                        Image slice = image(rect_domain);
                        int ubound = image.ubound(axis-1);

                #if 0   

                        Pixel* tmp_buf = new Pixel[image.extent(axis-1)];
                
                        // on each pixel of the slice, get its "perpendicular line" in the axis direction, and make a deriche filter on it
                        for (typename Image::template iterator<domain> i=slice.begin() ; i!=slice.end() ; ++i){
                                rect_domain.lbound() = i.position();
                                rect_domain.ubound() = i.position();
                                rect_domain.ubound(axis-1) = ubound;

                                _deriche_map(image(rect_domain), axis, tmp_buf, rect_domain.lbound(), rect_domain.ubound());
                        }

                        delete[] tmp_buf;

                #else   

                        Index tmp[slice.size()];
                        int i = 0;
                        for (typename Image::template iterator<domain> it=slice.begin() ; it!=slice.end() ; ++it)
                                tmp[i++] = it.position();

                        #pragma omp parallel 
                        {
                        // ATTENTION un Ybuf par thread.
                        Pixel* tmp_buf = new Pixel[image.extent(axis-1)];
                        #pragma omp for
                        for (int i = 0 ; i < slice.size() ; i++){
                                blitz::RectDomain<DIM> rect_domain(tmp[i],tmp[i]);
                                rect_domain.ubound(axis-1) = ubound;
                                _deriche_map(image(rect_domain), axis, tmp_buf, rect_domain.lbound(), rect_domain.ubound());
                        }
                        delete[] tmp_buf;
                        }
                #endif

                }


                // Modify/Access member data 


                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::set_filter_type(int _filter_type)
                {
                        m_filter_type = _filter_type;
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::set_order(int _order)
                {
                        m_order = _order;
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::set_border_cond(bool _border_cond)
                {
                        m_border_cond = _border_cond;
                }       

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::set_isotropic_sigma(float _sigma)
                {
                        m_filter_type = 0;
                        for (int i=0 ; i<DIM; i++)
                                m_sigma[i] = _sigma;
                        _init_coefs();
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::set_anisotropic_sigma(float _sigma[DIM])
                {
                        m_filter_type = 1;
                        for (int i=0 ; i<DIM; i++)
                                m_sigma[i] = _sigma[i];
                        _init_coefs();
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::set_image_sigma(const ImageVariance& _sigma)
                {
                        m_filter_type = 2;
                        m_sigma_image = &_sigma;
                        _init_coefs();
                }


                // Internal method 


                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::_init_coefs()
                {

                        float coefs[6];
                        float border_coefs[4];

                        switch(m_filter_type){

                                case 0 :                                     // isotropic filter
                                        _coefs_sigma(m_sigma[0], coefs, border_coefs);
                                        for (int i = 0 ; i < 6*DIM ; i++)
                                                m_coefs[i] = coefs[i%6];
                                        for (int i = 0 ; i < 4*DIM ; i++)
                                                m_border_coefs[i] = border_coefs[i%4];
                                        break;

                                case 1 :                                     // anisotropic filter
                                        for (int i = 0 ; i < DIM ; i++){
                                                _coefs_sigma(m_sigma[i], coefs, border_coefs);
                                                for (int j = 0 ; j < 6 ; j++)
                                                        m_coefs[i*6 + j] = coefs[j];
                                                for (int j = 0 ; j < 4 ; j++)
                                                        m_border_coefs[i*4 + j] = border_coefs[j];
                                        }
                                        break;

                                case 2 :                                     // space variant filter
                                        for (int i=0 ; i<6 ; i++){
                                                if (m_variant_coefs[i] != NULL) delete m_variant_coefs[i];
                                                m_variant_coefs[i] = new ImageVariance(m_sigma_image->shape());
                                        }
                                        for (typename ImageVariance::template iterator<domain> i=m_sigma_image->begin() ; i!=m_sigma_image->end() ; ++i ){
                                                _coefs_sigma((*m_sigma_image)(i), coefs, border_coefs);
                                                for (int j = 0 ; j < 6 ; j++)
                                                        (*m_variant_coefs[j])(i) = coefs[j];    
                                        }
                                        break;
                        }
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::_coefs_sigma(float sigma, float coefs[6],  float border_coefs[4])
                {
                
                        float a1,a2,a3,a4,b1,b2;         // coefs
                        float parity, g0, sumg1, sumg0;  // border coefs

                        // by default, dirichlet border condition
                        g0 = sumg0 = sumg1 = 0;

                        // precomputation of some values
                        const float alpha=sigma>0?1.695f/sigma:0;
                        const float ea=(float)exp(alpha);
                        const float ema=(float)exp(-alpha);
                        const float em2a=ema*ema;

                        // b1, b2 do not depend on the order of derivation
                        b1=2*ema;
                        b2=-em2a;

                        float ek,ekn;

                        // a1, a2, a3, a4 depend on the order of derivation
                        switch(m_order)
                        {
                                case 1:                 // first derivative
                                        ek = -(1-ema)*(1-ema)*(1-ema)/(2*(ema+1)*ema);
                                        a1 = 0;
                                        a2 = ek*ema;
                                        a3 = -ek*ema;
                                        a4 = 0;
                                        parity = 1;
                                        if (m_border_cond) { 
                                                sumg1 = (ek*ea)/((ea-1)*(ea-1));
                                                g0    = 0;
                                                sumg0 = g0+sumg1; 
                                        }
                                        break;
                                case 2:               // second derivative
                                        ekn = ( -2*(-1+3*ea-3*ea*ea+ea*ea*ea)/(3*ea+1+3*ea*ea+ea*ea*ea) );
                                        ek = -(em2a-1)/(2*alpha*ema);
                                        a1 = ekn;
                                        a2 = -ekn*(1+ek*alpha)*ema;
                                        a3 = ekn*(1-ek*alpha)*ema;
                                        a4 = -ekn*em2a;
                                        parity = -1;
                                        if (m_border_cond) {
                                                sumg1 = ekn/2;
                                                g0    = ekn;
                                                sumg0 = g0+sumg1;
                                        }
                                        break;
                                default:              // smoothing
                                        ek = (1-ema)*(1-ema) / (1+2*alpha*ema - em2a);
                                        a1 = ek;
                                        a2 = ek*ema*(alpha-1);
                                        a3 = ek*ema*(alpha+1);
                                        a4 = -ek*em2a;
                                        parity = 1;
                                        if (m_border_cond) {
                                                sumg1 = ek*(alpha*ea+ea-1)/((ea-1)*(ea-1));
                                                g0    = ek;
                                                sumg0 = g0+sumg1;
                                        }
                                        break;
                        }

                        coefs[0] = a1;
                        coefs[1] = a2;
                        coefs[2] = a3;
                        coefs[3] = a4;
                        coefs[4] = b1;
                        coefs[5] = b2;

                        border_coefs[0] = parity;
                        border_coefs[1] = g0;
                        border_coefs[2] = sumg1;
                        border_coefs[3] = sumg0;
                }

                template<unsigned DIM, typename Pixel>
                        void DericheFilter<DIM, Pixel>::_deriche_map(Image line, const int axis, Pixel* temp_buf, Index lbound, Index ubound)
                {
                        float b1, b2, a1, a2, a3, a4;    // coefs
                        float parity, g0, sumg1, sumg0;  // border coefs

                        if (m_filter_type != 2){
                                a1 = m_coefs[(axis-1)*6];
                                a2 = m_coefs[(axis-1)*6 + 1];
                                a3 = m_coefs[(axis-1)*6 + 2];
                                a4 = m_coefs[(axis-1)*6 + 3];
                                b1 = m_coefs[(axis-1)*6 + 4];
                                b2 = m_coefs[(axis-1)*6 + 5];

                                parity = m_border_coefs[(axis-1)*4];
                                g0     = m_border_coefs[(axis-1)*4 + 1];
                                sumg1  = m_border_coefs[(axis-1)*4 + 2];
                                sumg0  = m_border_coefs[(axis-1)*4 + 3];
                        }
                        
                        typename Image::template iterator<domain> i = line.begin();
                        Pixel Y0, I0, Y1, I1, Y2, I2;

                        // forward passing border condition

                        if (m_filter_type != 2){ 
                                I1 = line(i); ++i;
                                I2 = line(i); ++i;
                                Y2 = *(temp_buf++) = sumg0*I2;
                                Y1 = *(temp_buf++) = g0*I1 + sumg1*I2;
                        }
                        else{
                                I1 = line(i);
                                I2 = I1;
                                Y1 = I1/2.0;
                                Y2 = I1/2.0;
                        }
                        
                        // forward passing

                        for (; i != line.end() ; ++i){  
                                if (m_filter_type == 2){
                                        a1 = (*m_variant_coefs[0])(i.position() + lbound);
                                        a2 = (*m_variant_coefs[1])(i.position() + lbound);
                                        b1 = (*m_variant_coefs[4])(i.position() + lbound);
                                        b2 = (*m_variant_coefs[5])(i.position() + lbound);
                                }
                                I1 = line(i);
                                Y0 = *(temp_buf++) = a1*I1 + a2*I2 + b1*Y1 + b2*Y2; 
                                I2=I1;
                                Y2=Y1;
                                Y1=Y0;
                        }

                        // backward passing border condition
                        
                        line.reverseSelf(axis-1);
                        i = line.begin();
                        
                        if (m_filter_type != 2){ 
                                I2 = line(i);
                                Y2 = Y1 = parity*sumg1*I2;
                                line(i) = *(--temp_buf)+Y2; ++i;
                                I0 = line(i);
                                line(i) = *(--temp_buf)+Y1; ++i;
                        }
                        else{
                                I1 = line(i);
                                I2 = I1;
                                Y1 = I1/2.0;
                                Y2 = I1/2.0;
                        }
                        
                        // backward passing
                        
                        for (; i != line.end() ; ++i){
                                if (m_filter_type == 2){
                                        a3 = (*m_variant_coefs[2])(ubound - i.position());
                                        a4 = (*m_variant_coefs[3])(ubound - i.position());
                                        b1 = (*m_variant_coefs[4])(ubound - i.position());
                                        b2 = (*m_variant_coefs[5])(ubound - i.position());
                                }
                                Y0= a3*I1 + a4*I2 + b1*Y1 + b2*Y2;
                                I2=I1;
                                Y2=Y1;
                                Y1=Y0;
                                I1=line(i);
                                line(i)= *(--temp_buf) + Y0;  //stocking the new results!
                        }
                }
}

#endif // DERICHEFILTER_H

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