矩阵类定义与操作

在图像处理中经常用矩阵，抽时间看了矩阵的内容，还是有很多数学应用和解决方案的，整理了一下。

1.行主映射，列主映射

在二维数组中，可以把二维映射到一维。

行主映射：从第一行开始，一次对每行引索从左至右编号，把二维数组映射为[0,n-1].

列主映射：从第一列开始，一次对每列引索从左至右编号

如：

2.矩阵中的转置，相加（两个矩阵行数和列数分别相等），相乘（a的列数等于b的行数），运算符重载

例：matrix.h

#ifndef matrix_
#define matrix_

#include "myExceptions.h"

using namespace std;
template<class T>
class matrix 
{
   friend ostream& operator<<(ostream&, const matrix<T>&);
   public:
      matrix(int theRows = 0, int theColumns = 0);
      matrix(const matrix<T>&);
      ~matrix() {delete [] element;}
      int rows() const {return theRows;}
      int columns() const {return theColumns;}
      T& operator()(int i, int j) const;
      matrix<T>& operator=(const matrix<T>&);
      matrix<T> operator+() const; // unary +
      matrix<T> operator+(const matrix<T>&) const;
      matrix<T> operator-() const; // unary minus
      matrix<T> operator-(const matrix<T>&) const;
      matrix<T> operator*(const matrix<T>&) const;
      matrix<T>& operator+=(const T&);
   private:
       int theRows,    // number of rows in matrix
           theColumns; // number of columns in matrix
       T *element;     // element array
};  

template<class T>
matrix<T>::matrix(int theRows, int theColumns)
{// matrix constructor.
   // validate theRows and theColumns
   if (theRows < 0 || theColumns < 0)
      throw illegalParameterValue("Rows and columns must be >= 0");
   if ((theRows == 0 || theColumns == 0)
                && (theRows != 0 || theColumns != 0))
      throw illegalParameterValue
      ("Either both or neither rows and columns should be zero");

   // create the matrix
   this->theRows = theRows;
   this->theColumns = theColumns;
   element = new T [theRows * theColumns];
}

template<class T>
matrix<T>::matrix(const matrix<T>& m)
{// Copy constructor for matrices.
   // create matrix
   theRows = m.theRows;
   theColumns = m.theColumns;
   element = new T [theRows * theColumns];

   // copy each element of m
   copy(m.element,
        m.element + theRows * theColumns,
        element);
}

template<class T>
matrix<T>& matrix<T>::operator=(const matrix<T>& m)
{// Assignment. (*this) = m.
   if (this != &m)
   {// not copying to self
      delete [] element;
      theRows = m.theRows;
      theColumns = m.theColumns;
      element = new T [theRows * theColumns];
      // copy each element
      copy(m.element,
           m.element + theRows * theColumns,
           element);
   }
   return *this;
}

template<class T>
T& matrix<T>::operator()(int i, int j) const
{// Return a reference to element (i,j).
   if (i < 1 || i > theRows
       || j < 1 || j > theColumns)
   throw matrixIndexOutOfBounds();
   return element[(i - 1) * theColumns + j - 1];
}

template<class T>
matrix<T> matrix<T>::operator+(const matrix<T>& m) const
{// Return w = (*this) + m.
   if (theRows != m.theRows
       || theColumns != m.theColumns)
      throw matrixSizeMismatch();

   // create result matrix w
   matrix<T> w(theRows, theColumns);
   for (int i = 0; i < theRows * theColumns; i++)
      w.element[i] = element[i] + m.element[i];

   return w;
}

template<class T>
matrix<T> matrix<T>::
          operator-(const matrix<T>& m) const
{// Return (*this) - m.
   if (theRows != m.theRows 
       || theColumns != m.theColumns)
      throw matrixSizeMismatch();

   // create result matrix w
   matrix<T> w(theRows, theColumns);
   for (int i = 0; i < theRows * theColumns; i++)
      w.element[i] = element[i] - m.element[i];

   return w;
}

template<class T>
matrix<T> matrix<T>::operator-() const
{// Return w = -(*this).

   // create result matrix w
   matrix<T> w(theRows, theColumns);
   for (int i = 0; i < theRows * theColumns; i++)
      w.element[i] = -element[i];
   return w;

}

template<class T>
matrix<T> matrix<T>::operator*(const matrix<T>& m) const
{// matrix multiply.  Return w = (*this) * m.
   if (theColumns != m.theRows) 
      throw matrixSizeMismatch();

   matrix<T> w(theRows, m.theColumns);  // result matrix

   // define cursors for *this, m, and w
   // and initialize to location of (1,1) element
   int ct = 0, cm = 0, cw = 0;

   // compute w(i,j) for all i and j
   for (int i = 1; i <= theRows; i++)
   {// compute row i of result
      for (int j = 1; j <= m.theColumns; j++)
      { // compute first term of w(i,j)
          T sum =  element[ct] * m.element[cm];

          // add in remaining terms
          for (int k = 2; k <= theColumns; k++)
          {
             ct++;  // next term in row i of *this
             cm += m.theColumns;  // next in column j of m
             sum += element[ct] * m.element[cm];
          }
          w.element[cw++] = sum;  // save w(i,j)
 
          // reset to start of row and next column
          ct -= theColumns - 1;
          cm = j;
      }

      // reset to start of next row and first column
      ct += theColumns;
      cm = 0;
   }

   return w;
}

template<class T>
matrix<T>& matrix<T>::operator+=(const T& x)
{// Increment all elements of *this by x.
   for (int i = 0; i < theRows * theColumns; i++)
       element[i] += x;
   return *this;
}

template<class T>
ostream& operator<<(ostream& out, const matrix<T>& m)
{// Put matrix m into the stream out.
 // One row per line.
   int k = 0;  // index into element array
   for (int i = 0; i < m.theRows; i++)
   {// do row i
      for (int j = 0; j < m.theColumns; j++)
         out << m.element[k++] << "  ";

      // row i finished
      out << endl;
   }

   return out;
}

// for some reason compiler can't create this on its own
ostream& operator<<(ostream& out, const matrix<int>& m)
{// Put matrix m into the stream out.
 // One row per line.
   int k = 0;  // index into element array
   for (int i = 0; i < m.theRows; i++)
   {// do row i
      for (int j = 0; j < m.theColumns; j++)
         out << m.element[k++] << "  ";

      // row i finished
      out << endl;
   }

   return out;
}

#endif

matrix.cpp

// test matrix class

#include <iostream>
#include "matrix.h"

using namespace std;
int main(void)
{
   try
   {
      matrix<int> x(3,2), y, z;
      int i, j;
      for (i = 1; i <= 3; i++)
         for (j = 1; j <= 2; j++)
            x(i,j) = 2*i + j;
      cout << "Initialized x(i,j) = 2*i + j" << endl;
      cout << "x(3,1) = " << x(3,1) << endl;
      cout << "x is" << endl;;
      cout << x << endl;

      y = x;
      cout << "Assigned y = x" << endl;
      cout << "y is" << endl;
      cout << y << endl;

      x += 2;
      cout << "x incremented by 2 is" << endl;
      cout << x << endl;

      z = y + x;
      cout << "y + x is" << endl;
      cout << z << endl;

      cout << "-(y + x) is" << endl;
      cout << -z << endl;

      matrix<int> w(2,3);
      for (i = 1; i <= 2; i++)
         for (j = 1; j <= 3; j++)
            w(i,j) = i + j;
      cout << "Initialized w(i,j) = i + j" << endl;
      cout << "w is" << endl;
      cout << w << endl;

      z = y * w;
      cout << "y * w is" << endl;
      cout << z << endl;
   }
   catch (...) {
      cerr << "An exception has occurred" << endl;}

   return 0;
}

3.特殊矩阵中又有：对角矩阵，三对角矩阵，下三角矩阵，上三角矩阵，对称矩阵等

4稀疏矩阵：

在m*n 矩阵中，如果大多数元素都是0，就称为稀疏矩阵

而稀疏矩阵可以映射到一个线性表中，如果按行排列：

下面是稀疏矩阵重载操作符，矩阵转置，add：

sparseMatrix.h

// sparse matrix using an extended array list

#ifndef sparseMatrix_
#define sparseMatrix_

#include <iostream>
#include "matrixTerm.h"
#include "extendedArrayList.h"
#include "myExceptions.h"

template<class T>
class sparseMatrix
{
   friend ostream& operator<<
          (ostream&, sparseMatrix<T>&);
   friend istream& operator>>
          (istream&, sparseMatrix<T>&);
   public:
      void transpose(sparseMatrix<T> &b);
      void add(sparseMatrix<T> &b, sparseMatrix<T> &c);
   private:
      int rows,    // number of rows in matrix
          cols;    // number of columns in matrix
      arrayList<matrixTerm<T> > terms;
                   // list of nonzero terms
};

// overload <<
template <class T>
ostream& operator<<(ostream& out, sparseMatrix<T>& x)
{// Put x in output stream.

   // put matrix characteristics
   out << "rows = " << x.rows << " columns = "
       << x.cols  << endl;
   out << "nonzero terms = " << x.terms.size() << endl;

   // put terms, one per line
   for (arrayList<matrixTerm<T> >::iterator i = x.terms.begin();
        i != x.terms.end(); i++)
      out << "a(" << (*i).row << ',' << (*i).col
          << ") = " << (*i).value << endl;

   return out;
}

// overload >>
template<class T>
istream& operator>>(istream& in, sparseMatrix<T>& x)
{// Input a sparse matrix.

   // input matrix characteristics
   int numberOfTerms;
   cout << "Enter number of rows, columns, and #terms"
        << endl;
   in >> x.rows >> x.cols >> numberOfTerms;

   // set size of x.terms and ensure enough capacity
   x.terms.reSet(numberOfTerms);

   // input terms
   matrixTerm<T> mTerm;
   for (int i = 0; i < numberOfTerms; i++)
   {
      cout << "Enter row, column, and value of term " 
           << (i + 1) << endl;
      in >> mTerm.row >> mTerm.col >> mTerm.value;
      x.terms.set(i, mTerm);
   }

   return in;
}


/****************************************************************/
// explict code tooverload with T = int for test as compiler
// unable to generate

// overload <<
ostream& operator<<(ostream& out, sparseMatrix<int>& x)
{// Put x in output stream.

   // put matrix characteristics
   out << "rows = " << x.rows << " columns = "
       << x.cols  << endl;
   out << "nonzero terms = " << x.terms.size() << endl;

   // put terms, one per line
   for (arrayList<matrixTerm<int> >::iterator i = x.terms.begin();
        i != x.terms.end(); i++)
      out << "a(" << (*i).row << ',' << (*i).col
          << ") = " << (*i).value << endl;

   return out;
}

// overload >>
istream& operator>>(istream& in, sparseMatrix<int>& x)
{// Input a sparse matrix.

   // input matrix characteristics
   int numberOfTerms;
   cout << "Enter number of rows, columns, and #terms"
        << endl;
   in >> x.rows >> x.cols >> numberOfTerms;

   // set size of x.terms and ensure enough capacity
   x.terms.reSet(numberOfTerms);

   // input terms
   matrixTerm<int> mTerm;
   for (int i = 0; i < numberOfTerms; i++)
   {
      cout << "Enter row, column, and value of term " 
           << (i + 1) << endl;
      in >> mTerm.row >> mTerm.col >> mTerm.value;
      x.terms.set(i, mTerm);
   }

   return in;
}
/****************************************************************/
template<class T>
void sparseMatrix<T>::transpose(sparseMatrix<T> &b)
{// Return transpose of *this in b.

   // set transpose characteristics
   b.cols = rows;
   b.rows = cols;
   b.terms.reSet(terms.size());

   // initialize to compute transpose
   int* colSize = new int[cols + 1];
   int* rowNext = new int[cols + 1];
   
   // find number of entries in each column of *this
   for (int i = 1; i <= cols; i++) // initialize
      colSize[i] = 0;  
   for (arrayList<matrixTerm<T> >::iterator i = terms.begin();
        i != terms.end(); i++)
      colSize[(*i).col]++;  
   
   // find the starting point of each row of b
   rowNext[1] = 0;
   for (int i = 2; i <= cols; i++)
      rowNext[i] = rowNext[i - 1] + colSize[i - 1];  
   
   // perform the transpose copying from *this to b
   matrixTerm<T> mTerm;
   for (arrayList<matrixTerm<T> >::iterator i = terms.begin();
        i != terms.end(); i++)
   {
      int j = rowNext[(*i).col]++; // position in b
      mTerm.row = (*i).col;
      mTerm.col = (*i).row;
      mTerm.value = (*i).value;
      b.terms.set(j, mTerm);
   }
}

template<class T>
void sparseMatrix<T>::add(sparseMatrix<T> &b, sparseMatrix<T> &c)
{// Compute c = (*this) + b.

   // verify compatibility
   if (rows != b.rows || cols != b.cols)
     throw matrixSizeMismatch(); // incompatible matrices

   // set characteristics of result c
   c.rows = rows;
   c.cols = cols;
   c.terms.clear();
   int cSize = 0;

   // define iterators for *this and b
   arrayList<matrixTerm<T> >::iterator it = terms.begin();
   arrayList<matrixTerm<T> >::iterator ib = b.terms.begin();
   arrayList<matrixTerm<T> >::iterator itEnd = terms.end();
   arrayList<matrixTerm<T> >::iterator ibEnd = b.terms.end();

   // move through *this and b adding like terms
   while (it != itEnd && ib != ibEnd)
   {
     // row-major index plus cols of each term
     int tIndex = (*it).row * cols + (*it).col;
     int bIndex = (*ib).row * cols + (*ib).col;

     if (tIndex < bIndex)
     {// b term comes later
   	 c.terms.insert(cSize++, *it);
         it++;
     }
     else {if (tIndex == bIndex)
           {// both in same position

              // append to c only if sum not zero
              if ((*it).value + (*ib).value != 0)
              {
                 matrixTerm<T> mTerm;
                 mTerm.row = (*it).row;
                 mTerm.col = (*it).col;
                 mTerm.value = (*it).value + (*ib).value;
                 c.terms.insert(cSize++, mTerm);
              }

              it++; 
              ib++;
           }
           else
           {// a term comes later
              c.terms.insert(cSize++, *ib);
              ib++;
           }
          }
     }

   // copy over remaining terms
   for (; it != itEnd; it++)
      c.terms.insert(cSize++, *it);
   for (; ib != ibEnd; ib++)
      c.terms.insert(cSize++, *ib);
}

#endif

sparseMatrix.cpp:

// test array based sparse matrix class

#include <iostream>
#include "sparseMatrix.h"

int main(void)
{
   sparseMatrix<int> a, b, c;

   // test input and output
   cin >> a;
   cout << "Matrix a is" << endl << a;
   cin >> b;
   cout << "Matrix b is" << endl << b;

   // test transpose
   a.transpose(c);
   cout << "The transpose of a is" << endl << c;

   // test add
   a.add(b,c);
   cout << "The sum of a and b is" << endl << c;

   return 0;
}