/////////////////////////////////////////////////////
//This file defines operations to do math on matrices
/////////////////////////////////////////////////////
#ifndef MATRIX_ARITH_LIBRARY
#define MATRIX_ARITH_LIBRARY

//#include "cmatrix"
#define UNDEFINED 0

//////////////////////////
//Interface
//////////////////////////
///Checks to see if two matrices have the same dimensions
template <class T, class U>
bool samedimen(CMatrix<T> &LHS, CMatrix<U> &RHS);

///Convolves matrix 'I' with matrix 'Mask' and puts the results in 'Out.'  Note that this is a discrete approximation of an actual convolution.
template <class T, class U>
void convolve(CMatrix<T> &I, CMatrix<U> &Mask, CMatrix<T> & Out);

///Convovles matrix 'I' with matrix 'Mask' and puts the result in 'I'
template <class T, class U>
void convolve(CMatrix<T> &I, CMatrix<U> &Mask);

///Takes the discrete derivative in the x direction using a simple 2x2 window.
template <class T, class U>
void derivex(CMatrix<T> &In, CMatrix<U> &Out);

///Takes the discrete derivative in the y direction using a simple 2x2 window
template <class T, class U>
void derivey(CMatrix<T> &In, CMatrix<U> &Out);

///Takes the discrete derivative in t (subtracting the two images)
template <class T, class U>
void derivet(CMatrix<T> &In1, CMatrix<U> &In2, CMatrix<T> &Out);

///does Out = OP1 * OP2
template <class T, class U> //Performs matrix multiplication
T & dotprod(CMatrix<T> &OP1,CMatrix<U> &OP2, CMatrix<T> & Out);

///Transposes a matrix and puts the result in 'Out'
template <class T>
void transpose(CMatrix<T> &In, CMatrix<T> &Out);

///////////////////////////////////////////////////////////
//Implementation
///////////////////////////////////////////////////////////

template <class T, class U>
bool samedimen(CMatrix<T> &LHS, CMatrix<U> &RHS)
{
	return	(LHS.xmax()==RHS.xmax() & 
			LHS.xmin()==RHS.xmin() & 
			LHS.ymax()==RHS.ymax() &
			LHS.ymin()==RHS.ymin())
}

template <class T, class U>
void convolve(CMatrix<T> &I, CMatrix<U> &Mask, CMatrix<T> &Out)
{
double masksum=0,tempsum, realstorage;
for(i=0; i< mask.size(); i++)
	masksum+=mask[i];
if(I.size()!=Out.size())
	Out.SetRange(I);
for(i=I.xmin();i<I.xmax();i++)
   for(j=I.ymin();j<I.ymax();j++)
   {   Out(i,j)=0; tempsum=0;  realstorage=0;
       for(m1=Mask.xmin();m1<Mask.xmax();m1++)  //For even convolutions, the result has to run differently.
         for(m2=Mask.ymin();m2<Mask.ymax();m2++)//This is why there is a ternary operator.
		 { //Assume every part outside the boundary of the image is filled with zeros.         
	     if (
			 ( ((i+m1)<I.xmax()) &
			   ((j+m2)<I.ymax()) ) & 
			 ( ((j+m2)>=I.ymin())   & 
			   ((i+m1)>=I.xmin())  )  )
			   realstorage+=
			           (U)input(i+m1,j+m2)*mask(m1,m2);
		  else
			  tempsum+=mask(m1,m2);				
		 }
		 if(tempsum>0) out(i,j)=(T)realstorage*(1+tempsum);
   }					//Used instead of using 0's around the edges.
}

template <class T, class U>
void convolve(CMatrix<T> &I, CMatrix<U> &Mask)
{
	CMatrix<T> TEMP;
	TEMP.SetSize(I,0);
	convolve(I,Mask,TEMP);
	I=TEMP;
}

//There are two convolutions: one to subtract in the dimension in which the derivative
//is being taken, and one to average the result of this operation.
template <class T, class U>
void derivex(CMatrix<T> &In, CMatrix<T> &Out)
{
	CMatrix<T> Mask(1,0);
	Mask[0]=1;
	Mask[1]=-1;
	convolve(In,Mask,Out)
	Mask.SetSize(0,1);
	Mask[0]=.5;
	Mask[1]=.5;
	convolve(Out,Mask);
}

template <class T, class U>
void derivey(CMatrix<T> &In, CMatrix<T> &Out)
{
	CMatrix<T> Mask(0,1);
	Mask[0]=1;
	Mask[1]=-1;
	convolve(In,Mask,Out)
	Mask.SetSize(0,1);
	Mask[0]=.5;
	Mask[1]=.5;
	convolve(Out,Mask);
}

template <class T, class U>
void derivet(CMatrix<T> &In1, CMatrix<U> &In2, CMatrix<T> &Out)
{
	CMatrix<T> Mask(0,1);
	int	Xmin=In1.xmin()>In2.xmin()?In1.xmin():In2.xmin(),
		Xmax=In1.xmax()<In2.xmax()?In1.xmax():In2.xmax(),
		Ymin=In1.ymin()>In2.ymin()?In1.ymin():in2.ymin(),
		Ymax=In1.ymax()>In2.ymax()?In1.ymax():In2.ymax();
	for(i=Xmin;i<Xmax; i++)
		for(j=Ymin; j<Ymax; j++)
			Out=In1(i,j)-In2(i,j);
	//Usually, the differences are averaged slightly.  We'll do that.
	Mask[0]=.5;
	Mask[1]=.5;
	convolve(Out,Mask);
	Mask.SetSize(1,0);
	Mask[0]=.5;
	Mask[1]=.5;
	convolve(Out,Mask);
}

template <class T, class U> //Performs matrix multiplication
T & dotprod(CMatrix<T> &OP1,CMatrix<U> &OP2, CMatrix<T> & Out)
{
	int i=0,j=0,k;
	if(OP1.xmin()!=OP2.ymin()||OP1.xmax()!=OP2.ymax())
	{
		cerr << "Multiplication is not possible for these matrices.  Perhaps you must transpose one of them." << endl;
		exit(0);
	}
	if((OP1.y()!=Out.y()) || (OP2.x()!=Out.x()) )
		out.initialize(OP2.x(),OP1.y());
	for(k=0; k<rhs.x(); k++)
		for(j=0; j<lhs.y(); j++)
		{
			out(k,j)=0;
			for(i=0; i<lhs.x(); i++)
				Out(k,j)+=OP1(i,j)*OP2(k,i);
		};
		return out;
}

template <class T>
void transpose(CMatrix<T> &In,CMatrix<T> & Out)
{
	Out.SetSize(In.ymin(),In.xmin(),In.ymax(),In.xmax());
	int i=0, j=0;
	for(i=0; i<xmax; i++)
		for(j=0; j<ymax; j++)
			Out(j,i)=OP1(i,j);
		
}


#endif //MATRIX_ARITH_LIBRARY