---

matrixarith

Go to the documentation of this file.

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

//#include "cmatrix"
#define UNDEFINED 0

//Interface
template <class T, class U>
bool samedimen(CMatrix<T> &LHS, CMatrix<U> &RHS);

template <class T, class U>
void convolve(CMatrix<T> &I, CMatrix<U> &Mask, CMatrix<T> & Out);

template <class T, class U>
void convolve(CMatrix<T> &I, CMatrix<U> &Mask);

template <class T, class U>
void derivex(CMatrix<T> &In, CMatrix<U> &Out);

template <class T, class U>
void derivey(CMatrix<T> &In, CMatrix<U> &Out);

template <class T, class U>
void derivet(CMatrix<T> &In1, CMatrix<U> &In2, CMatrix<T> &Out);

template <class T, class U> //Performs matrix multiplication
T & dotprod(CMatrix<T> &OP1,CMatrix<U> &OP2, CMatrix<T> & 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

---