#ifndef IM3DMATHS_HEADER
#define IM3DMATHS_HEADER 1


/////////////////////////////////////////////////////////////////////
//! userlevel = Normal
//! docentry="Image3D"
//! rcsid="$Id: Im3dMaths.hh,v 1.2 1999/11/05 15:33:54 ees1cg Exp $"
//! file="amma/Image3D/Im3dMaths.hh"
//! lib=Image3D
//! author="Manuel Segovia"
//! date= "11/11/98"

#include "amma/Image3d.hh"
#include "amma/SArr3Iter.hh"
#include "amma/SArr3Iter2.hh"
#include "amma/SArr3Iter3.hh"

template<class DataT>
Image3dC<DataT> operator+(const Image3dC<DataT> & a, const Image3dC<DataT> & b );
//: overload of the adding function

template<class DataT>
Image3dC<DataT> operator-(const Image3dC<DataT> & a, const Image3dC<DataT> & b);
//: overload of the substracting function

template<class DataT>
Image3dC<DataT>  operator*(const Image3dC<DataT> & a, const Image3dC<DataT> & b);
//: overload of the multiplication function

template<class DataT>
Image3dC<DataT>  operator/(const Image3dC<DataT> & a, const Image3dC<DataT> & b);
//: overload of the division function

template<class DataT>
Image3dC<DataT> sqr(const Image3dC<DataT> & in);
//: overload of sqr

template<class DataT>
Image3dC<DataT> sqrt(const Image3dC<DataT> & in);
//: overload of sqrt

template<class DataT>
Image3dC<DataT> operator+(const Image3dC<DataT> & a, const Image3dC<DataT> & b )
{
  Image3dC<DataT> result(a.Dim1(),a.Dim2(),a.Dim3(),a.SubSamplingZ(),a.SubSamplingY());
  SArray3dIter3C< DataT, DataT, DataT> it(const_cast<Image3dC<DataT> &>(a), const_cast<Image3dC<DataT> &>(b), result);
  for( it.First() ;it.IsElm();it.Next() )
    it.Data3() =  it.Data1() + it.Data2();
  return result;
}


template<class DataT>
Image3dC<DataT> operator-(const Image3dC<DataT> & a, const Image3dC<DataT> & b)
{
  Image3dC<DataT> result(a.Dim1(),a.Dim2(),a.Dim3(),a.SubSamplingZ(),a.SubSamplingY());
  SArray3dIter3C< DataT, DataT, DataT> it(const_cast<Image3dC<DataT> &>(a), const_cast<Image3dC<DataT> &>(b), result);
  for( it.First() ;it.IsElm();it.Next() )
    it.Data3() =  it.Data1() - it.Data2();
  return result;
}


template<class DataT>
Image3dC<DataT>  operator*(const Image3dC<DataT> & a, const Image3dC<DataT> & b)
{
  Image3dC<DataT> result(a.Dim1(),a.Dim2(),a.Dim3(),a.SubSamplingZ(),a.SubSamplingY());
  SArray3dIter3C< DataT, DataT, DataT> it(const_cast<Image3dC<DataT> &>(a), const_cast<Image3dC<DataT> &>(b), result);
  for( it.First() ;it.IsElm();it.Next() )
    it.Data3() =  it.Data1() * it.Data2();
  return result;
}

template<class DataT>
Image3dC<DataT>  operator/(const Image3dC<DataT> & a, const Image3dC<DataT> & b)
{
  Image3dC<DataT> result(a.Dim1(),a.Dim2(),a.Dim3(),a.SubSamplingZ(),a.SubSamplingY());
  SArray3dIter3C< DataT, DataT, DataT> it(const_cast<Image3dC<DataT> &>(a), const_cast<Image3dC<DataT> &>(b), result);
  for( it.First() ;it.IsElm();it.Next() )
    it.Data3() =  it.Data1() / it.Data2();
  return result;
}


template<class DataT>
Image3dC<DataT> sqr(const Image3dC<DataT> & in)
{
  Image3dC<DataT> result(in.Dim1(), in.Dim2(), in.Dim3(), in.SubSamplingZ(), in.SubSamplingY());
  for(SArray3dIter2C<DataT,DataT> it(const_cast<Image3dC<DataT> &>(in), result); it.IsElm(); it.Next())
    it.Data2() = it.Data1()*it.Data1();
  return result;
}

template<class DataT>
Image3dC<DataT> sqrt(const Image3dC<DataT> & in)
{
  Image3dC<DataT> result(in.Dim1(), in.Dim2(), in.Dim3(), in.SubSamplingZ(), in.SubSamplingY());
  for(SArray3dIter2C<DataT,DataT> it(const_cast<Image3dC<DataT> &>(in), result); it.IsElm(); it.Next())
    it.Data2() = sqrt(it.Data1());
  return result;
}



#endif