Convolution.hxx

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_hxx
#define _Convolution_hxx

template <class ImageType>
inline void
IP3D::Convolution(const ImageType& imToConvolute, const Filter& imFilter, ImageType& imResult)
{
    ImageProcessorArgs<ImageType> args(ImArgs(imToConvolute),imResult);
    args.SameSize();
    ImageType &imDest=args.GetDest();

    typename ImageType::const_iteratorXYZ pIm;
    typename ImageType::iteratorFast pDest;
    Filter::const_iteratorXYZ pFilter;
    
    int rx = imFilter.GetCenter().x;
    int ry = imFilter.GetCenter().y;
    int rz = imFilter.GetCenter().z;

    for (pDest = imDest.begin(), pIm = imToConvolute.begin(); pDest!=imDest.end(); pDest++, pIm++)
    {
        (*pDest) = Zero<typename ImageType::value_type>();
        for (pFilter = imFilter.begin(); pFilter !=imFilter.end(); pFilter++)
        {
            (*pDest) += (typename ImageType::value_type) 
                ((*pFilter)*imToConvolute.SafeValue(pIm.x - (pFilter.x - rx) , 
                                                         pIm.y - (pFilter.y - ry) , 
                                                         pIm.z - (pFilter.z - rz) ) );
        }
//          if (!(rand()%1000)){
//              cout << "Value: " << *pDest << endl;
//              getchar();
//          }
    }
    //args.Finished();
}

// gcc3 complex doen't support float * comlex<double>
//#ifdef _FFT_hpp
namespace IP3D
{
template <>
inline void
Convolution<Image3Dcomplex >(const Image3Dcomplex& imToConvolute, const Filter& imFilter, Image3Dcomplex& imResult)
{
    ImageProcessorArgs<Image3Dcomplex> args(ImArgs(imToConvolute),imResult);
    args.SameSize();
    Image3Dcomplex &imDest=args.GetDest();

    Image3Dcomplex::const_iteratorXYZ pIm;
    Image3Dcomplex::iteratorFast pDest;
    Filter::const_iteratorXYZ pFilter;
    
    int rx = imFilter.GetCenter().x;
    int ry = imFilter.GetCenter().y;
    int rz = imFilter.GetCenter().z;

    for (pDest = imDest.begin(), pIm = imToConvolute.begin(); pDest!=imDest.end(); pDest++, pIm++)
    {
        (*pDest) = Zero<Image3Dcomplex::value_type>();
        for (pFilter = imFilter.begin(); pFilter !=imFilter.end(); pFilter++)
        {
            (*pDest) += (Image3Dcomplex::value_type) 
                (((double)(*pFilter))*imToConvolute.SafeValue(pIm.x - (pFilter.x - rx) , 
                                                         pIm.y - (pFilter.y - ry) , 
                                                         pIm.z - (pFilter.z - rz) ) );
        }
//          if (!(rand()%1000)){
//              cout << "Value: " << *pDest << endl;
//              getchar();
//          }
    }
    //args.Finished();
}
//#endif //_FFT_hpp
}

template <class ImageType>
void
IP3D::ConvolutionFFT(const ImageType& imToConvolute, const Filter& imFilter, ImageType& imResult)
{
    ImageProcessorArgs<ImageType> args(ImArgs(imToConvolute),imResult);
    args.SameSize();
//      ImageType &imDest=args.GetDest();

    ThrowError("IP3D::ConvolutionFFT Error: not yet implemented");
    
    //args.Finished();

}

template <class ImageType>
void
IP3D::SeparableConvolution(const ImageType& imToConvolute, const SeparableFilter& filter, ImageType& imResult)
{
    ImageProcessorArgs<ImageType> args(ImArgs(imToConvolute),imResult);
    args.SameSize();
    ImageType &imDest=args.GetDest();

    typename ImageType::value_type zero;
    SetZero<typename ImageType::value_type>::Set(zero);
    imDest.Fill(zero);

    int rx = filter.iCenter.x;
    int ry = filter.iCenter.y;
    int rz = filter.iCenter.z;

    typename ImageType::iteratorXYZ p;
    ImageProgress tracker("SeparableConvolution (3x)",p);

    // Convolute X direction
//  ImageType xRes(imDest.Size());
    ImageType &xRes=imDest;
//  cout << "Convolution in X direction:" << endl;
    for (p = xRes.begin(); p != xRes.end(); p++)
    {
        *p = Zero<typename ImageType::value_type>();

        for (int n=0; n< filter.size.width; n++)
        {
            (*p) += (typename ImageType::value_type) 
                (filter.xcomponent[n]*imToConvolute.SafeValue(p.x - (n - rx), 
                                                              p.y, 
                                                              p.z ));
        }
    }
//  cout << "Convolution in Y direction:" << endl;
    // Convolute Y direction
    ImageType yRes(imDest.Size());
    for (p = yRes.begin(); p != yRes.end(); p++)
    {
        *p = Zero<typename ImageType::value_type>();
        for (int n=0; n< filter.size.height; n++)
        {
            (*p) += (typename ImageType::value_type) 
                (filter.ycomponent[n]*xRes.SafeValue(p.x, 
                                                     p.y - (n - ry), 
                                                     p.z ));
        }
    }

//      cout << "Convolution in Z direction:" << endl;
    // Convolute Z direction
    for (p = imDest.begin(); p != imDest.end(); p++)
    {
        *p = Zero<typename ImageType::value_type>();
        for (int n=0; n< filter.size.depth; n++)
        {
            (*p) += (typename ImageType::value_type) 
                (filter.zcomponent[n]*yRes.SafeValue(p.x, 
                                                     p.y, 
                                                     p.z  - (n - rz)));
        }
    }
}


#endif // _Convolution_hxx

---