/* 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   Arithmetic.hxx
  \brief  Function definitions for classes declared in .hpp 


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


#ifndef _Arithmetic_hxx
#define _Arithmetic_hxx


//! Apply calculates imDest = im1 + im2
/*!
  This is a function that calculate im1 + im2 voxel by voxel
 */

template <class ImageType>
void
IP3D::Addition(const ImageType& im1, const ImageType& im2, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1,im2),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p1;
	typename ImageType::const_iteratorFast p2;
	typename ImageType::iteratorFast pDest;
	
	for (p1 = im1.begin(), p2 = im2.begin(), pDest = imDest.begin(); pDest != imDest.end(); p1++, p2++, pDest++)
	{
		(*pDest) = (*p1) + (*p2);
	}

}

//! Adds a scalar value to each element of Imagetype
/*!
  <long-description>

  \param const Imagetype& im1    
  \param const typename ImageType::value_type&scalarValue    
  \param Imagetype& imDest 
*/
template <class ImageType>
void
IP3D::AdditionWithConstant(const ImageType& im1, const typename ImageType::value_type& scalarValue, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p;
	typename ImageType::iteratorFast pDest;
	for(p = im1.begin(), pDest = imDest.begin(); pDest != imDest.end(); p++, pDest++)
	{
		(*pDest) = (*p) + scalarValue;
	}
}

//! Apply calculates imDest = im1 - im2
/*!
  
 */

template <class ImageType>
void
IP3D::Difference(const ImageType& im1, const ImageType& im2, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1,im2),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p1;
	typename ImageType::const_iteratorFast p2;
	typename ImageType::iteratorFast pDest;
	
	for (p1 = im1.begin(), p2 = im2.begin(), pDest = imDest.begin(); pDest != imDest.end(); p1++, p2++, pDest++)
	{
		(*pDest) = (*p1) - (*p2);
	}
}

//! Substracts a scalar value from each element of Imagetype
/*!
  <long-description>

  \param const Imagetype& im1    
  \param const typename ImageType::value_type&scalarValue    
  \param Imagetype& imDest 
*/
template <class ImageType>
void
IP3D::DifferenceWithConstant(const ImageType& im1, const typename ImageType::value_type& scalarValue, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p;
	typename ImageType::iteratorFast pDest;
	for(p = im1.begin(), pDest = imDest.begin(); pDest != imDest.end(); p++, pDest++)
	{
		(*pDest) = (*p) - scalarValue;
	}

}

//! Multiply imDest = im1 * im2, voxel by voxel
/*!
  
 */
template <class ImageType>
void
IP3D::Multiplication(const ImageType& im1, const ImageType& im2, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1,im2),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p1;
	typename ImageType::const_iteratorFast p2;
	typename ImageType::iteratorFast pDest;
	
	for (p1 = im1.begin(), p2 = im2.begin(), pDest = imDest.begin(); pDest != imDest.end(); p1++, p2++, pDest++)
	{
		(*pDest) = (*p1) * (*p2);
	}
}

//! Divide imDest = im1 * im2, voxel by voxel
/*!
  
 */
template <class ImageType>
void
IP3D::Division(const ImageType& im1, const ImageType& im2, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1,im2),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p1;
	typename ImageType::const_iteratorFast p2;
	typename ImageType::iteratorFast pDest;
	
	for (p1 = im1.begin(), p2 = im2.begin(), pDest = imDest.begin(); pDest != imDest.end(); p1++, p2++, pDest++)
	{
		(*pDest) = (*p1) / (*p2);
	}
}

//! Multiplies each voxel with the same scalar
/*!
  <long-description>

  \param const Imagetype& im1    
  \param const typename ImageType::value_type&scalarValue    
  \param Imagetype& imDest 
*/
template <class ImageType>
void
IP3D::MultiplicationWithConstant(const ImageType& im1, double v, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p;
	typename ImageType::iteratorFast pDest;
	for(p = im1.begin(), pDest = imDest.begin(); pDest != imDest.end(); p++, pDest++)
	{
		*pDest=(typename ImageType::value_type)(v*(*p));
	}
}

template <class ImageType>
void
IP3D::MaxImage(const ImageType& im1, const ImageType& im2, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1,im2),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p1;
	typename ImageType::const_iteratorFast p2;
	typename ImageType::iteratorFast pDest;
	
	for (p1 = im1.begin(), p2 = im2.begin(), pDest = imDest.begin(); pDest != imDest.end(); p1++, p2++, pDest++)
	{
		(*pDest) = max((*p1) , (*p2));
	}

}
template <class ImageType>
void
IP3D::Abs(const ImageType& im1, ImageType& _imDest)
{
	ImageProcessorArgs<ImageType> args(ImArgs(im1),_imDest);
	args.SameSize();
	args.InPlace();
	ImageType &imDest=args.GetDest();

	typename ImageType::const_iteratorFast p1;
	typename ImageType::iteratorFast pDest;
	
	for (p1 = im1.begin(), pDest = imDest.begin(); pDest != imDest.end(); p1++, pDest++)
	{
		(*pDest) = AbsoluteValue(*p1);
	}

}


// **************************** Image3Dlinear : Operators********************
template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator+=(const ImageType &other)
{
	IP3D::Addition(*image(), other, *image());
	return *image();
}

template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator+=(const Im3DValue& scalarValue)
{
	IP3D::AdditionWithConstant(*image(), scalarValue, *image());
	return *image();
}

template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator-=(const ImageType &other)
{
	IP3D::Difference(*image(), other, *image());
	return *image();
}

template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator-=(const Im3DValue& scalarValue)
{
	IP3D::DifferenceWithConstant(*image(), scalarValue, *image());
	return *image();
}

template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator*=(const ImageType &other)
{
	IP3D::Multiplication(*image(), other, *image());
	return *image();
}

template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator*=(double scalarValue)
{
	IP3D::MultiplicationWithConstant(*image(), scalarValue, *image());
	return *image();
}

template<class ImageType,class Im3DValue>
ImageType &ArithmeticOperators<ImageType,Im3DValue>::operator/=(const ImageType &other)
{
	IP3D::Division(*image(), other, *image());
	return *image();
}


#endif // _Arithmetic_hxx
