Convolution.cpp

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.
 */

#include<ImLib3D/FFT.hpp>
#include<ImLib3D/Convolution.hpp>
#include<ImLib3D/Transformation.hpp>
#include<ImLib3D/Arithmetic.hpp>

void
IP3D::FFTLowPassFilterApodizedIdeal(const Image3Df& src,float cutoffFreq,Image3Df &res,float rolloffFactor)
{
    Image3Dcomplexf fft;
    IP3D::FFT(src,fft,false);
    Vect3Di center=fft.SizeV()/2;
    center.x=0;// half fft trf
    // cuttoff value in each direction
    Vect3Df cutoff=cutoffFreq*.5*(Vect3Df)fft.SizeV();
    cutoff.x=cutoffFreq*fft.SizeV().x;// half fft trf
    cutoff.Show("cutoff:","\n");
    Vect3Df rolloff=cutoffFreq*rolloffFactor*(Vect3Df)cutoff;
    rolloff.Show("rolloff:","\n");
    for(Image3Dcomplexf::iteratorXYZ p=fft.begin();p!=fft.end();++p)
    {
        Vect3Di f=p.Pos();
        if((f.y+center.y)>=fft.Height()){f.y-=fft.Height();}
        if((f.z+center.z)>=fft.Depth ()){f.z-=fft.Depth();}

        if(!rolloffFactor)
        {
            if(abs(f.x)>cutoff.x || abs(f.y)>cutoff.y || abs(f.z)>cutoff.z){*p=complexf(0);}
        }
        else
        {
            for(uint i=0;i<3;i++)
            {
                float u=(abs(f(i))-cutoff(i))/rolloff(i);   
                if(u> 1){*p =complexf(0);break;}
                if(u>-1){double m=cos((1+u)*M_PI/4);*p*=complexf(m);}
            }
        }
    }
    IP3D::FFTInverse(fft,res,src.Width(),true);
    res*=1.0/res.GetNVoxels();
}

void
IP3D::FFTLowPassFilterButterworth(const Image3Df& src,float cutoffFreq,float order,Image3Df &res)
{
    Image3Dcomplexf fft;
    IP3D::FFT(src,fft,false);
    Vect3Di center=fft.SizeV()/2;
    center.x=0;// half fft trf
    for(Image3Dcomplexf::iteratorXYZ p=fft.begin();p!=fft.end();++p)
    {
        Vect3Df f=p.Pos();
        if((f.y+center.y)>=fft.Height()){f.y-=fft.Height();}
        if((f.z+center.z)>=fft.Depth ()){f.z-=fft.Depth();}
        f.x/=   fft.Width();
        f.y/=.5*fft.Height();
        f.z/=.5*fft.Depth();
        float r=f.Norm();
        (*p)*=1/(1+pow(r/cutoffFreq,2*order));
    }
    IP3D::FFTInverse(fft,res,src.Width(),true);
    res*=1.0/res.GetNVoxels();
}

void
IP3D::FFTLowPassFilterChebyscheff(const Image3Df& src,float cutoffFreq,float order,Image3Df &res)
{
    Image3Dcomplexf fft;
    IP3D::FFT(src,fft,false);
    Vect3Di center=fft.SizeV()/2;
    center.x=0;// half fft trf
    for(Image3Dcomplexf::iteratorXYZ p=fft.begin();p!=fft.end();++p)
    {
        Vect3Df f=p.Pos();
        if((f.y+center.y)>=fft.Height()){f.y-=fft.Height();}
        if((f.z+center.z)>=fft.Depth ()){f.z-=fft.Depth();}
        f.x/=   fft.Width();
        f.y/=.5*fft.Height();
        f.z/=.5*fft.Depth();
        float w=f.Norm()/cutoffFreq;
        float t;
        if(w<1){t=cos (order*acos (w));}
        else{   t=cosh(order*acosh(w));}
        (*p)*=1/(1+t*t);
    }
    IP3D::FFTInverse(fft,res,src.Width(),true);
    res*=1.0/res.GetNVoxels();
}

---