/* 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.
 */
/*! \file Convolution.hxx */
#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) = imDest.GetZero();
		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
template <>
inline void
IP3D::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) = imDest.GetZero();
		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

#ifdef HAS_FFTW

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();

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

}
#endif // HAS_FFTW

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();

	imDest.Fill(0);

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

	typename ImageType::iteratorXYZ 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=xRes.GetZero();
		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=yRes.GetZero();
  		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) = imDest.GetZero();
  		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
