Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

Convolution.hpp

Go to the documentation of this file.
00001 /* ImLib3D
00002  * Copyright (c) 2001, ULP-IPB Strasbourg.
00003  *
00004  * This program is free software; you can redistribute it and/or modify
00005  * it under the terms of the GNU General Public License as published by
00006  * the Free Software Foundation; either version 2 of the License, or (at
00007  * your option) any later version.
00008  * 
00009  * This program is distributed in the hope that it will be useful, but
00010  * WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * General Public License for more details.
00013  * 
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software
00016  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00017  */
00018 #ifndef _Convolution_hpp
00019 #define _Convolution_hpp
00020 
00030 #include<ImLib3D/ImageProcessor.hpp>
00031 #include<ImLib3D/Display.hpp>
00032 #include<ImLib3D/TestPatterns.hpp>
00033 #include<ImLib3D/Filter.hpp>
00034 #include<ImLib3D/Arithmetic.hpp>
00035 #include<ImLib3D/ImageProgress.hpp>
00036 #include<ImLib3D/Transformation.hpp>
00037 
00038 namespace IP3D
00039 {
00049     template <class ImageType>
00050     void Convolution(const ImageType& imToConvolve, const Filter& imFilter, ImageType& imResult);
00051 };
00052 
00053 namespace IP3D
00054 {
00064     template <class ImageType>
00065     void SeparableConvolution(const ImageType& imToConvolve, const SeparableFilter& filter, ImageType& imResult);
00066 };
00067 
00068 
00069 namespace IP3D
00070 {
00081     template <class ImageType>
00082     void ConvolutionFFT(const ImageType& imToConvolve, const Filter& imFilter, ImageType& imResult);
00083 };
00084 
00085 
00086 namespace IP3D
00087 {
00099     template <class ImageType>
00100     void PartialSum(const ImageType& cumSum,const RectZone3Di &r0,typename ImageType::value_type& res)
00101     {       
00102         RectZone3Di r=r0;
00103         r.x0--;
00104         r.y0--;
00105         r.z0--;
00106         res=cumSum(r.x1,r.y1,r.z1);
00107 
00108         if(                      r.z0>=0){res+=-cumSum(r.x1,r.y1,r.z0);}
00109         if(           r.y0>=0           ){res+=-cumSum(r.x1,r.y0,r.z1);}
00110         if(           r.y0>=0 && r.z0>=0){res+= cumSum(r.x1,r.y0,r.z0);}
00111         if(r.x0>=0                      ){res+=-cumSum(r.x0,r.y1,r.z1);}
00112         if(r.x0>=0 &&            r.z0>=0){res+= cumSum(r.x0,r.y1,r.z0);}
00113         if(r.x0>=0 && r.y0>=0           ){res+= cumSum(r.x0,r.y0,r.z1);}
00114         if(r.x0>=0 && r.y0>=0 && r.z0>=0){res+=-cumSum(r.x0,r.y0,r.z0);}
00115     }
00116 };
00117 
00118 namespace IP3D
00119 {
00127     template <class ImageType>
00128     void CumulativeSum(const ImageType& src, ImageType& res0)
00129     {
00130         ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
00131         args.SameSize();
00132         ImageType &res=args.GetDest();
00133         for(typename ImageType::iteratorXYZ p=res.begin();p!=res.end();++p)
00134         {
00135             const Vect3Di &v=p.Pos();
00136             *p=src(v);
00137             typename ImageType::value_type psr;
00138             
00139             if(v.z>0){IP3D::PartialSum(res,RectZone3Di(  0,  0,  0,  v.x  ,v.y  ,v.z-1),psr);*p+=psr;}
00140             if(v.y>0){IP3D::PartialSum(res,RectZone3Di(  0,  0,v.z,  v.x  ,v.y-1,v.z  ),psr);*p+=psr;}
00141             if(v.x>0){IP3D::PartialSum(res,RectZone3Di(  0,v.y,v.z,  v.x-1,v.y  ,v.z  ),psr);*p+=psr;}
00142         }
00143     }
00144 };
00145 
00146 namespace IP3D
00147 {
00156     template <class ImageType>
00157     void BoxFilter(const ImageType& src, int size, ImageType& res0)
00158     {
00159         ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
00160         args.SameSize();
00161         ImageType &res=args.GetDest();
00162 
00163         int center=size/2;
00164         typedef typename ImageType::value_type Im3DValue;
00165         for(int x=0;x<src.Width() ;x++)
00166             for(int y=0;y<src.Height();y++)
00167             {       
00168                 Im3DValue sum= Zero<typename ImageType::value_type>();
00169 
00170                 for(int z=0;z<src.Depth()+center;z++)
00171                 {
00172                     if(z<src.Depth()){sum+=src(x,y,z);}
00173                     if(z>=size)      {sum-=src(x,y,z-size);}
00174                     if(z>=center)res(x,y,z-center)=sum;
00175                 }
00176             }
00177         {
00178             ImageType tmp=res;
00179             for(int x=0;x<src.Width() ;x++)
00180                 for(int z=0;z<src.Depth();z++)
00181                 {       
00182                     Im3DValue sum= Zero<typename ImageType::value_type>();
00183                     for(int y=0;y<src.Height()+center;y++)
00184                     {
00185                         if(y<src.Height()){sum+=tmp(x,y,z);}
00186                         if(y>=size)       {sum-=tmp(x,y-size,z);}
00187                         if(y>=center)res(x,y-center,z)=sum;
00188                     }
00189                 }
00190         }
00191         {
00192             ImageType tmp=res;
00193             for(int y=0;y<src.Height() ;y++)
00194                 for(int z=0;z<src.Depth();z++)
00195                 {       
00196                     Im3DValue sum= Zero<typename ImageType::value_type>();
00197                     for(int x=0;x<src.Width()+center;x++)
00198                     {
00199                         if(x<src.Width()){sum+=tmp(x,y,z);}
00200                         if(x>=size)      {sum-=tmp(x-size,y,z);}
00201                         if(x>=center)res(x-center,y,z)=sum;
00202                     }
00203                 }
00204         }
00205         res*=1.0/pow((double)size,3);
00206     }
00207 };
00208 
00209 namespace IP3D
00210 {
00224     template <class ImageType>
00225     void GaussianApproxFilter(const ImageType& src,int order,int size,ImageType& res)
00226     {
00227         res=src;
00228         //      if(!(size%2)){ThrowError("GaussianApproxFilter: size (%d) must be unpair (3,5,7...)",size);}
00229         for(int i=0;i<order;i++)
00230         {
00231             IP3D::BoxFilter(res,size,res);
00232             // recenter image by 1 voxel for every second application of the bosfilter when the size is pair
00233             if (i>0 && (i%2) && !(size%2))
00234             {
00235                 IP3D::WrapTranslate(res,Vect3Di(1,1,1),res);
00236             }
00237         }       
00238     }
00239 
00252     template <class ImageType>
00253     void BorderCorrectedGaussianApproxFilter(const ImageType &src,int size,int nbFilterIter,ImageType &res)
00254     {
00255         ImageType lowFreq=src;
00256         if(lowFreq.HasMask())
00257         {
00258             lowFreq.Mask().Normalize();
00259             lowFreq.Mask().ApplyOn(lowFreq);
00260         }
00261 
00262         // apply filter in normal fashion
00263         IP3D::GaussianApproxFilter(lowFreq,nbFilterIter,size,lowFreq);
00264     
00265         // compensate for unmasked and border regions in filtering 
00266         {
00267             Image3Df smask(src.Size());
00268             if(src.HasMask()){smask=src.Mask();}
00269             else
00270             {smask.Fill(1);}
00271             IP3D::GaussianApproxFilter(smask,nbFilterIter,size,smask);
00272             for(typename ImageType::iteratorXYZ p=lowFreq.begin();p!=lowFreq.end();p++)
00273             {
00274                 if(smask(p.Pos())){*p/=smask(p.Pos());}
00275             }
00276         }
00277         if(src.HasMask()){src.Mask().ApplyOn(lowFreq);}
00278 
00279         res=lowFreq;
00280     }
00281 
00282 };
00283 
00284 namespace IP3D
00285 {
00298     template <class ImageType>
00299     void SymmetricBinomialFilter(const ImageType& src,int order,ImageType& res)
00300     {
00301         BinomialFilter filter(order);
00302 //      cout << "Filtering with filter:" << endl << filter;
00303         res=src;
00304         IP3D::SeparableConvolution(src,filter,res);
00305     }
00306 
00307 };
00308 
00309 namespace IP3D
00310 {
00322     void FFTLowPassFilterApodizedIdeal(const Image3Df& src,float cutoffFreq,Image3Df &res,float rolloffFactor=0);
00323 };
00324 
00325 namespace IP3D
00326 {
00338     void FFTLowPassFilterButterworth(const Image3Df& src,float cutoffFreq,float order,Image3Df &res);
00339 };
00340 
00341 namespace IP3D
00342 {
00354     void FFTLowPassFilterChebyscheff(const Image3Df& src,float cutoffFreq,float order,Image3Df &res);
00355 };
00356 
00357 
00358 #include<ImLib3D/Convolution.hxx>
00359 
00360 #endif //_Convolution_hpp

Generated on Fri Jun 17 13:36:00 2005 for ImLib3D by  doxygen 1.4.2