matrix.h

图像分割算法

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



#pragma once
#ifndef MATRIX_H
#define MATRIX_H



#include <iostream>

using std::cout;
using std::cerr;
using std::endl;
using std::ostream;
using std::right;
using std::fixed;

#include <iomanip>

using std::setw;

#include <cstdlib>
#include <ctime>
#include <cmath>

#include <string>

using std::string;

#include <sstream>

using std::ostringstream;

#include <typeinfo>

#include <limits>
using std::numeric_limits;


using namespace std;


#include "cimpl.h"


namespace CIMPL
{


// forward declaration
//template< class T > class Array;
template< class T > class Matrix;
template< class T > class SubMatrix;
template< class T > class Vector;

/// \brief 2-D Matrix class.
template< class T >
class Matrix 
{

	//friend class Array<T>;
	friend class Vector<T>;
	friend class SubMatrix<T>;


protected:
	T *data;
	int ndims;  // always 2
	int length; // xDim*yDim
	int xDim;
	int yDim;
	Cleaner<T> *clean;   // Reference counting Garbage collector.

	Vector<T> *columns;
	
public:
	Matrix();
	Matrix(const int rows, const int columns);
	Matrix(const int rows, const int columns, T init);
	
	Matrix(string str);

	Matrix(Matrix<T> &m); // copy constructor

	~Matrix();
	
	void Clean();

	/// Make sure enough memory is already allocated.
	/// Memory is from now on handled by the matrix... Don't free it manually.
	void Set(const int rows, const int columns, T *_data); 

	const T* DataPtr();
	T* Data();

	Matrix<T> Clone() const;
	
	void SwitchColumns(int i, int j);

	/// \brief reorganizes data array such that matrix elements are 
	/// ordered in column major ordering. Calling this function will 
	/// allocate new memory with correct ordering of elements and old memory is freed.
	///
	/// If you created external (manual) pointers to matrix data, you might experience 
	/// memory/pointer errors after the data is rearranged.
	void SyncData2Columns();
	
//	SubMatrix<T> Slice(int row1, int row2, int col1, int col2);
//	SubMatrix<T> Slice(string str1, string str2);

	Matrix<T> Slice(int row1, int row2, int col1, int col2);
	Matrix<T> Slice(string str1, string str2);

	Vector<T> Slice(string str);


	const int Columns() const;
	const int Rows() const;
	const int XDim() const;
	const int YDim() const;
	const int* Dims() const;
	const int Length() const;
	const int Numel() const;
	const int NDims() const;
	void Init(const T init);

	Matrix<T>&  Rand(const double max);
	static Matrix<T> Rand(const int rows, const int cols, const double max);

	void ReadFromMatrix(const Matrix<T>& m);
	void ReadFromMatrix(const Matrix<T>& m, const int rowStart, const int colStart);

	static Matrix<T> Cat(int dimension, Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<T> Cat(int dimension, Matrix<T>& m1, Vector<T>& m2);
	static Matrix<T> Cat(int dimension, Vector<T>& m1, Matrix<T>& m2);
	static Matrix<T> Cat(int dimension, Vector<T>& m1, Vector<T>& m2);

	static Matrix<T> Ones(int side);
	static Matrix<T> Ones(int rows, int cols);
	static Matrix<T> Zeros(int side);
	static Matrix<T> Zeros(int rows, int cols);


	static bool IsSquare(Matrix<T>& m);
	static bool IsM2MCompatible(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<T> MMultiply(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<T> MMultiply(Matrix<T>& m1, Vector<T>& m2);
	static Matrix<T> MMultiply(Vector<T>& m1, Matrix<T>& m2);

	static Matrix<T> Transpose(Matrix<T>& m);
	Matrix<T>& Transpose();

	//static Matrix<T> Hermitian(Matrix<T>& m);
	//Matrix<T>& Hermitian();

	static bool IsCompatible(Matrix<T>& m1, Matrix<T>& m2);
	
	// Boolean Operations...
	static Matrix<int> And(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Or(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Lt(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Gt(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Le(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Ge(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Eq(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<int> Ne(Matrix<T>& m1, Matrix<T>& m2);

	// Boolean Operations with value type...
	static Matrix<int> And(Matrix<T>& m, T v);
	static Matrix<int> Or(Matrix<T>& m, T v);
	static Matrix<int> Lt(Matrix<T>& m, T v);
	static Matrix<int> Gt(Matrix<T>& m, T v);
	static Matrix<int> Le(Matrix<T>& m, T v);
	static Matrix<int> Ge(Matrix<T>& m, T v);
	static Matrix<int> Eq(Matrix<T>& m, T v);
	static Matrix<int> Ne(Matrix<T>& m, T v);


// Add matrix to another matrix
	static Matrix<T> Add(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<T> Subtract(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<T> Multiply(Matrix<T>& m1, Matrix<T>& m2);
	static Matrix<T> Divide(Matrix<T>& m1, Matrix<T>& m2);

// Add value to another matrix
	static Matrix<T> Add(Matrix<T>& m1, T v2);
	static Matrix<T> Subtract(Matrix<T>& m1, T v2);
	static Matrix<T> Subtract(T v2, Matrix<T>& m1);
	static Matrix<T> Multiply(Matrix<T>& m1, T v2);
	static Matrix<T> Divide(Matrix<T>& m1, T v2);
	static Matrix<T> Divide(T v2, Matrix<T>& m1);


// Add matrix to "this" matrix
	Matrix<T>& Add(Matrix<T>& m);
	Matrix<T>& Subtract(Matrix<T>& m);
	Matrix<T>& Multiply(Matrix<T>& m);
	Matrix<T>& Divide(Matrix<T>& m);

// Add value to "this" matrix
	Matrix<T>& Add(T v);
	Matrix<T>& Subtract(T v);
	Matrix<T>& Multiply(T v);
	Matrix<T>& Divide(T v);


//OPERATORS
	//Matrix<T>& operator= (T init);
	Matrix<T>& operator= (Matrix<T>& m);
	//Matrix<T>& operator= (Array<T>& m);
	Matrix<T>& operator= (Vector<T>& m);
	Matrix<T>& operator= (SubMatrix<T>& m);

	Matrix<T>& operator= (string str);

	Matrix<T> operator+ ();
	Matrix<T> operator- ();
	Matrix<int> operator! ();
	
	friend Matrix<T> operator, (Matrix<T>& m1, Matrix<T>& m2)
	{
		return Matrix<T>::Cat(2, m1, m2);
	}
	
	friend Matrix<T> operator, (Matrix<T>& m1, Vector<T>& m2)
	{
		return Matrix<T>::Cat(2, m1, m2);
	}

	friend Matrix<T> operator, (Vector<T>& m1, Matrix<T>& m2)
	{
		return Matrix<T>::Cat(2, m1, m2);
	}



	// Transpose
	Matrix<T> operator~ ();

	friend ostream& operator<< (ostream& output, Matrix<T>& m)
	{
		int bufferWidth = 80;
		int lm = m.Rows();
		int rowNoWidth = (int)log10((double)(lm-1))+2;
		int infoWidth = 3+rowNoWidth+2;

		int maxLength  = 1;
		for(int i=0; i<m.Length(); i++)
		{
			ostringstream oS;
			oS << m.ElemNC(i);
			int l = (int)oS.str().length();
			if(l>maxLength)
			{
				maxLength = l;
			}
		}

		int columnWidth = maxLength+3;
		int columnMax = (bufferWidth-infoWidth)/columnWidth;
		
		output << typeid(m).name() << " of size " << m.Rows() << "x" << m.Columns() << endl;
		output << "=======================" << endl;
		if(m.Columns() <= columnMax)
		{
			for(int i=0;i<m.Rows();i++)