Convolution.hpp

Go to the documentation of this file.

/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#ifndef _Convolution_hpp
#define _Convolution_hpp

#include<ImLib3D/ImageProcessor.hpp>
#include<ImLib3D/Display.hpp>
#include<ImLib3D/TestPatterns.hpp>
#include<ImLib3D/Filter.hpp>
#include<ImLib3D/Arithmetic.hpp>
#include<ImLib3D/ImageProgress.hpp>
#include<ImLib3D/Transformation.hpp>

namespace IP3D
{
    template <class ImageType>
    void Convolution(const ImageType& imToConvolve, const Filter& imFilter, ImageType& imResult);
};

namespace IP3D
{
    template <class ImageType>
    void SeparableConvolution(const ImageType& imToConvolve, const SeparableFilter& filter, ImageType& imResult);
};


namespace IP3D
{
    template <class ImageType>
    void ConvolutionFFT(const ImageType& imToConvolve, const Filter& imFilter, ImageType& imResult);
};


namespace IP3D
{
    template <class ImageType>
    void PartialSum(const ImageType& cumSum,const RectZone3Di &r0,typename ImageType::value_type& res)
    {       
        RectZone3Di r=r0;
        r.x0--;
        r.y0--;
        r.z0--;
        res=cumSum(r.x1,r.y1,r.z1);

        if(                      r.z0>=0){res+=-cumSum(r.x1,r.y1,r.z0);}
        if(           r.y0>=0           ){res+=-cumSum(r.x1,r.y0,r.z1);}
        if(           r.y0>=0 && r.z0>=0){res+= cumSum(r.x1,r.y0,r.z0);}
        if(r.x0>=0                      ){res+=-cumSum(r.x0,r.y1,r.z1);}
        if(r.x0>=0 &&            r.z0>=0){res+= cumSum(r.x0,r.y1,r.z0);}
        if(r.x0>=0 && r.y0>=0           ){res+= cumSum(r.x0,r.y0,r.z1);}
        if(r.x0>=0 && r.y0>=0 && r.z0>=0){res+=-cumSum(r.x0,r.y0,r.z0);}
    }
};

namespace IP3D
{
    template <class ImageType>
    void CumulativeSum(const ImageType& src, ImageType& res0)
    {
        ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
        args.SameSize();
        ImageType &res=args.GetDest();
        for(typename ImageType::iteratorXYZ p=res.begin();p!=res.end();++p)
        {
            const Vect3Di &v=p.Pos();
            *p=src(v);
            typename ImageType::value_type psr;
            
            if(v.z>0){IP3D::PartialSum(res,RectZone3Di(  0,  0,  0,  v.x  ,v.y  ,v.z-1),psr);*p+=psr;}
            if(v.y>0){IP3D::PartialSum(res,RectZone3Di(  0,  0,v.z,  v.x  ,v.y-1,v.z  ),psr);*p+=psr;}
            if(v.x>0){IP3D::PartialSum(res,RectZone3Di(  0,v.y,v.z,  v.x-1,v.y  ,v.z  ),psr);*p+=psr;}
        }
    }
};

namespace IP3D
{
    template <class ImageType>
    void BoxFilter(const ImageType& src, int size, ImageType& res0)
    {
        ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
        args.SameSize();
        ImageType &res=args.GetDest();

        int center=size/2;
        typedef typename ImageType::value_type Im3DValue;
        for(int x=0;x<src.Width() ;x++)
            for(int y=0;y<src.Height();y++)
            {       
                Im3DValue sum= Zero<typename ImageType::value_type>();

                for(int z=0;z<src.Depth()+center;z++)
                {
                    if(z<src.Depth()){sum+=src(x,y,z);}
                    if(z>=size)      {sum-=src(x,y,z-size);}
                    if(z>=center)res(x,y,z-center)=sum;
                }
            }
        {
            ImageType tmp=res;
            for(int x=0;x<src.Width() ;x++)
                for(int z=0;z<src.Depth();z++)
                {       
                    Im3DValue sum= Zero<typename ImageType::value_type>();
                    for(int y=0;y<src.Height()+center;y++)
                    {
                        if(y<src.Height()){sum+=tmp(x,y,z);}
                        if(y>=size)       {sum-=tmp(x,y-size,z);}
                        if(y>=center)res(x,y-center,z)=sum;
                    }
                }
        }
        {
            ImageType tmp=res;
            for(int y=0;y<src.Height() ;y++)
                for(int z=0;z<src.Depth();z++)
                {       
                    Im3DValue sum= Zero<typename ImageType::value_type>();
                    for(int x=0;x<src.Width()+center;x++)
                    {
                        if(x<src.Width()){sum+=tmp(x,y,z);}
                        if(x>=size)      {sum-=tmp(x-size,y,z);}
                        if(x>=center)res(x-center,y,z)=sum;
                    }
                }
        }
        res*=1.0/pow((double)size,3);
    }
};

namespace IP3D
{
    template <class ImageType>
    void GaussianApproxFilter(const ImageType& src,int order,int size,ImageType& res)
    {
        res=src;
        //      if(!(size%2)){ThrowError("GaussianApproxFilter: size (%d) must be unpair (3,5,7...)",size);}
        for(int i=0;i<order;i++)
        {
            IP3D::BoxFilter(res,size,res);
            // recenter image by 1 voxel for every second application of the bosfilter when the size is pair
            if (i>0 && (i%2) && !(size%2))
            {
                IP3D::WrapTranslate(res,Vect3Di(1,1,1),res);
            }
        }       
    }

    template <class ImageType>
    void BorderCorrectedGaussianApproxFilter(const ImageType &src,int size,int nbFilterIter,ImageType &res)
    {
        ImageType lowFreq=src;
        if(lowFreq.HasMask())
        {
            lowFreq.Mask().Normalize();
            lowFreq.Mask().ApplyOn(lowFreq);
        }

        // apply filter in normal fashion
        IP3D::GaussianApproxFilter(lowFreq,nbFilterIter,size,lowFreq);
    
        // compensate for unmasked and border regions in filtering 
        {
            Image3Df smask(src.Size());
            if(src.HasMask()){smask=src.Mask();}
            else
            {smask.Fill(1);}
            IP3D::GaussianApproxFilter(smask,nbFilterIter,size,smask);
            for(typename ImageType::iteratorXYZ p=lowFreq.begin();p!=lowFreq.end();p++)
            {
                if(smask(p.Pos())){*p/=smask(p.Pos());}
            }
        }
        if(src.HasMask()){src.Mask().ApplyOn(lowFreq);}

        res=lowFreq;
    }

};

namespace IP3D
{
    template <class ImageType>
    void SymmetricBinomialFilter(const ImageType& src,int order,ImageType& res)
    {
        BinomialFilter filter(order);
//      cout << "Filtering with filter:" << endl << filter;
        res=src;
        IP3D::SeparableConvolution(src,filter,res);
    }

};

namespace IP3D
{
    void FFTLowPassFilterApodizedIdeal(const Image3Df& src,float cutoffFreq,Image3Df &res,float rolloffFactor=0);
};

namespace IP3D
{
    void FFTLowPassFilterButterworth(const Image3Df& src,float cutoffFreq,float order,Image3Df &res);
};

namespace IP3D
{
    void FFTLowPassFilterChebyscheff(const Image3Df& src,float cutoffFreq,float order,Image3Df &res);
};


#include<ImLib3D/Convolution.hxx>

#endif //_Convolution_hpp

---