/* 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
/*!
  \file   Convolution.hpp
  \brief  Convolution of two images 

  <long description>

  See AUTHORS and COPYING files for more information
  \date   2001-01-24
*/

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



namespace IP3D
{
/*! #&Convolution&#     Convolution of an image by a user defined filter.
        The filter class is just like a float valued image, with a center point 
        added. Ignores masks.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Filter        type:Filter
      <LI>Result image  type:ImageType
     </UL>
~*/
	template <class ImageType>
	void Convolution(const ImageType& imToConvolve, const Filter& imFilter, ImageType& imResult);
};

namespace IP3D
{
/*! #&SeparableConvolution&#     Convolution of an image by a user defined filter.
        The filter class is just like a float valued image, with a center point 
        added. Ignores masks.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Filter        type:Filter
      <LI>Result image  type:ImageType
     </UL>
~*/
	template <class ImageType>
	void SeparableConvolution(const ImageType& imToConvolve, const SeparableFilter& filter, ImageType& imResult);
};


#ifdef HAS_FFTW
namespace IP3D
{
/*! #&ConvolutionFFT&#     Convolution of an image by a user defined filter using FFT.
        It uses  multiplication in fourier domain with forward and inverse fft. This 
        is faster than normal convolution for lmarge filter sizes. Ignores masks. 
        UNIMPLEMENTED!
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Filter        type:Filter
      <LI>Result image  type:ImageType
     </UL>
~*/
	template <class ImageType>
	void ConvolutionFFT(const ImageType& imToConvolve, const Filter& imFilter, ImageType& imResult);
};
#endif // HAS_FFTW


namespace IP3D
{
/*! #&PartialSum&#     Computes the sum of an image in a rectangular zone. It uses an image 
    computed using CumulativeSum.
        This is a fast way of doing filtering with big square convolution filters. 
        Result is the sum of all points in the rectzone: the rectzone *includes 
        both* of its defining points.
     <UL>Args:
      <LI>CumulativeSum image  type:ImageType
      <LI>zone                 type:RectZone3Di
      <LI>result               type:Im3DValue
     </UL>
~*/
	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
{
/*! #&CumulativeSum&#     Computes an image that I(i,j,k)=Sum(0<=r<=i, 0<=s<=j , 0<=t<=k : S(r,s,t) )
        This can be used by PartialSum for fast filtering.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Result image  type:ImageType
     </UL>
~*/
	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
{
/*! #&BoxFilter&#     Very fast convolution by an integer sized box signal
        Ignores masks.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>size          type:int
      <LI>Result image  type:ImageType
     </UL>
~*/
	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=src.GetZero();
				for(int z=0;z<src.Depth()+center-1;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=src.GetZero();
					for(int y=0;y<src.Height()+center-1;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=src.GetZero();
					for(int x=0;x<src.Width()+center-1;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
{
/*! #&GaussianApproxFilter&#     Fast convolution by an approximate "integer" guassian filter/
        This just applies several boxfilters (boxsize=size 
        and order times), and is only efficient for low orders. Ignores masks.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>order         type:int
      <LI>size          type:int
      <LI>Result image  type:ImageType
     </UL>
~*/
	template <class ImageType>
	void GaussianApproxFilter(const ImageType& src,int order,int size,ImageType& res)
	{
		res=src;
		for(int i=0;i<order;i++)
		{
			IP3D::BoxFilter(res,size,res);
		}		
	}

};

#ifdef HAS_FFTW
namespace IP3D
{
/*! #&CFFTFilterLowPass&# ~*/
		void FFTFilterLowPass(const Image3Df& src,float cutoffFreq,Image3Df &res,float rolloffFactor=0);
};
#endif // HAS_FFTW

#include<ImLib3D/Convolution.hxx>

#endif //_Convolution_hpp
