core.inl

图像分割算法

//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.


namespace MathCore
{


	template< class T >
	void Plot(Vector<T>& v, int height)
	{
		int bufferWidth = 80;
		Vector<T> ind(bufferWidth-2);
		for(int i=0; i<ind.Numel(); i++)
		{
			int newInd = int(v.Numel()*(i+0.5)/ind.Numel()+0.5);
			if(newInd > v.Numel()-1)
			{
				newInd = v.Numel()-1;
			}
			ind(i) = v(newInd);
		}
		T max = Maximum(ind);
		T min = Minimum(ind);

		ind -= min;
		ind *= (T)((height-1.1)/(max-min));
		Matrix<int> m(height, ind.Numel(), 0);
		for(int i=0; i<ind.Numel(); i++)
		{
			m(int(ind(i)+0.5), i) = 1;
		}


		cout << endl << endl;
		for(int i=0; i<bufferWidth; i++)
		{
			cout << "-";
		}
		cout << endl << max << endl;
		for(int i=0; i<m.Rows(); i++)
		{
			cout << " ";
			for(int j=0; j<m.Columns(); j++)
			{
				if(m(i,j) == 1)
				{
					cout << "*";
				}
				else
				{
					cout << " ";
				}
			}
			cout << " ";
			cout << endl;
		}
		cout << min;
		cout << endl << endl;
		for(int i=0; i<bufferWidth; i++)
		{
			cout << "-";
		}
		cout << endl << endl;

	}

	template <class T>
	Vector<int> Find(Vector<T>& m)
	{
		std::vector<int> ind;
		for(int i=0; i<m.Length(); i++)
		{
			if(m.ElemNC(i) != 0)
			{
				ind.push_back(i);
			}
		}
		Vector<int> temp((int)ind.size());
		memcpy(temp.Data(), &ind[0], sizeof(int)*ind.size() );
		//std::vector<int>::const_iterator constIterator;
		//int j = 0;
		//for(constIterator = ind.begin(); constIterator != ind.end(); constIterator++)
		//{
		//	temp[j] = *constIterator;
		//	j++;
		//}
		return temp;
	}
	
	template <class T>
	Vector<int> Find(Matrix<T>& m)
	{
		std::vector<int> ind;
		for(int i=0; i<m.Length(); i++)
		{
			if(m.ElemNC(i) != 0)
			{
				ind.push_back(i);
			}
		}
		Vector<int> temp((int)ind.size());
		memcpy(temp.Data(), &ind[0], sizeof(int)*ind.size() );
		return temp;
	}
	
	template <class T>
	void Find(Matrix<T>& m, Vector<int>& I, Vector<int>& J)
	{
		std::vector<int> indI;
		std::vector<int> indJ;
		for(int j=0; j<m.Columns(); j++)
		{
			for(int i=0; i<m.Rows(); i++)
			{
				if(m.ElemNC(i,j) != 0)
				{
					indI.push_back(i);
					indJ.push_back(j);
				}
			}
		}
		I = Vector<int>((int)indI.size());
		memcpy(I.Data(), &indI[0], sizeof(int)*indI.size() );
		J = Vector<int>((int)indJ.size());
		memcpy(J.Data(), &indJ[0], sizeof(int)*indJ.size() );
	}


	
	// Vectors need to be of length 3
	template <class T>
	Vector<T> Cross(Vector<T>& x, Vector<T>& y)
	{
		Vector<T> temp(3);
		temp[0] = x[1]*y[2]-x[2]*y[1];
		temp[1] = x[2]*y[0]-x[0]*y[2];
		temp[2] = x[0]*y[1]-x[1]*y[0];
		
		return temp;
	}

	
	template <class T>
	Matrix<T> Cross(Matrix<T>& x, Matrix<T>& y)
	{
		if(x.Rows() == 3)
		{
			return Cross(x,y,1);
		}
		else if(x.Columns() == 3)
		{
			return Cross(x,y,2);
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("One of the dimension of the matrices x and y should be of length 3!");
		}
	}

	
	template <class T>
	Matrix<T> Cross(Matrix<T>& x, Matrix<T>& y, int dimension)
	{
		if(x.Rows() != y.Rows() || x.Columns() != y.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrices are not of the same size!");
		}
		
		Matrix<T> temp(x.Rows(), x.Columns());
		
		if(dimension == 1)
		{
			if(temp.Rows() != 3)
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Both matrices x and y should be of length 3 along the dimension on which the cross product is taken!");
			}
			for(int z=0; z<temp.Columns(); z++)
			{
				temp.ElemNC(0,z) = x.ElemNC(1,z)*y.ElemNC(2,z) - x.ElemNC(2,z)*y.ElemNC(1,z);
				temp.ElemNC(1,z) = x.ElemNC(2,z)*y.ElemNC(0,z) - x.ElemNC(0,z)*y.ElemNC(2,z);
				temp.ElemNC(2,z) = x.ElemNC(0,z)*y.ElemNC(1,z) - x.ElemNC(1,z)*y.ElemNC(0,z);
			}
		}
		else if(dimension == 2)
		{
			if(temp.Columns() != 3)
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Both matrices x and y should be of length 3 along the dimension on which the cross product is taken!");
			}
			for(int z=0; z<temp.Rows(); z++)
			{
				temp.ElemNC(z,0) = x.ElemNC(z,1)*y.ElemNC(z,2) - x.ElemNC(z,2)*y.ElemNC(z,1);
				temp.ElemNC(z,1) = x.ElemNC(z,2)*y.ElemNC(z,0) - x.ElemNC(z,0)*y.ElemNC(z,2);
				temp.ElemNC(z,2) = x.ElemNC(z,0)*y.ElemNC(z,1) - x.ElemNC(z,1)*y.ElemNC(z,0);
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
		}
		return temp;
	}


	
	template <class T>
	Vector<T> CumSum(Vector<T>& m)
	{
		Vector<T> temp(m.Length());
		temp.ElemNC(0) = m.ElemNC(0);
		for(int i=1; i<m.Length(); i++)
		{
			temp.ElemNC(i) = temp.ElemNC(i-1) + m.ElemNC(i);
		}
		return temp;
	}
	
	template <class T>
	Matrix<T> CumSum(Matrix<T>& m)
	{
		return CumSum(m,1);
	}
	
	template <class T>
	Matrix<T> CumSum(Matrix<T>& m, int dimension)
	{
		Matrix<T> temp(m.Rows(), m.Columns());
		if(dimension == 1)
		{
			for(int z=0; z<m.Columns(); z++)
			{
				temp.ElemNC(0,z) = m.ElemNC(0,z);
				for(int i=1; i<m.Rows(); i++)
				{
					temp.ElemNC(i,z) = temp.ElemNC(i-1,z) + m.ElemNC(i,z);
				}
			}
		}
		else if(dimension == 2)
		{
			for(int z=0; z<m.Rows(); z++)
			{
				temp.ElemNC(z,0) = m.ElemNC(z,0);
				for(int i=1; i<m.Columns(); i++)
				{
					temp.ElemNC(z,i) = temp.ElemNC(z,i-1) + m.ElemNC(z,i);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
		}	
		return temp;
	}


	
	
	template <class T>
	Vector<T>& CumSumI(Vector<T>& m)
	{
		for(int i=1; i<m.Length(); i++)
		{
			m.ElemNC(i) = m.ElemNC(i-1) + m.ElemNC(i);
		}
		return m;
	}
	
	template <class T>
	Matrix<T>& CumSumI(Matrix<T>& m)
	{
		return CumSum(m,1);
	}
	
	template <class T>
	Matrix<T>& CumSumI(Matrix<T>& m, int dimension)
	{
		if(dimension == 1)
		{
			for(int z=0; z<m.Columns(); z++)
			{
				for(int i=1; i<m.Rows(); i++)
				{
					m.ElemNC(i,z) = m.ElemNC(i-1,z) + m.ElemNC(i,z);
				}
			}
		}
		else if(dimension == 2)
		{
			for(int z=0; z<m.Rows(); z++)
			{
				for(int i=1; i<m.Columns(); i++)
				{
					m.ElemNC(z,i) = m.ElemNC(z,i-1) + m.ElemNC(z,i);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
		}	
		return m;
	}


	
	template <class T>
	Vector<T> CumProd(Vector<T>& m)
	{
		Vector<T> temp(m.Length());
		temp.ElemNC(0) = m.ElemNC(0);
		for(int i=1; i<m.Length(); i++)
		{
			temp.ElemNC(i) = temp.ElemNC(i-1) * m.ElemNC(i);
		}
		return temp;
	}
	
	template <class T>
	Matrix<T> CumProd(Matrix<T>& m)
	{
		return CumProd(m,1);
	}
	
	template <class T>
	Matrix<T> CumProd(Matrix<T>& m, int dimension)
	{
		Matrix<T> temp(m.Rows(), m.Columns());
		if(dimension == 1)
		{
			for(int z=0; z<m.Columns(); z++)
			{
				temp.ElemNC(0,z) = m.ElemNC(0,z);
				for(int i=1; i<m.Rows(); i++)
				{
					temp.ElemNC(i,z) = temp.ElemNC(i-1,z) * m.ElemNC(i,z);
				}
			}
		}
		else if(dimension == 2)
		{
			for(int z=0; z<m.Rows(); z++)
			{
				temp.ElemNC(z,0) = m.ElemNC(z,0);
				for(int i=1; i<m.Columns(); i++)
				{
					temp.ElemNC(z,i) = temp.ElemNC(z,i-1) * m.ElemNC(z,i);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
		}	
		return temp;
	}

	
	template <class T>
	Vector<T>& CumProdI(Vector<T>& m)
	{
		for(int i=1; i<m.Length(); i++)
		{
			m.ElemNC(i) = m.ElemNC(i-1) * m.ElemNC(i);
		}
		return m;
	}
	
	template <class T>
	Matrix<T>& CumProdI(Matrix<T>& m)
	{
		return CumProd(m,1);
	}
	
	template <class T>
	Matrix<T>& CumProdI(Matrix<T>& m, int dimension)
	{
		if(dimension == 1)
		{
			for(int z=0; z<m.Columns(); z++)
			{
				for(int i=1; i<m.Rows(); i++)
				{
					m.ElemNC(i,z) = m.ElemNC(i-1,z) * m.ElemNC(i,z);
				}
			}
		}
		else if(dimension == 2)
		{
			for(int z=0; z<m.Rows(); z++)
			{
				for(int i=1; i<m.Columns(); i++)
				{
					m.ElemNC(z,i) = m.ElemNC(z,i-1) * m.ElemNC(z,i);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
		}	
		return m;
	}

	
	template <class T>
	Matrix<T> Diag(Vector<T>& m)
	{
		return Diag(m,0);
	}
	
	template <class T>
	Matrix<T> Diag(Vector<T>& m, int k)
	{
		int absk = abs(k);
		int dim = m.Length() + absk;
		Matrix<T> temp(dim,dim,0);
		if(k>=0)
		{
			for(int i=0; i<m.Length(); i++)
			{
				temp.ElemNC(i,i+absk) = m[i];
			}
		}
		else
		{
			for(int i=0; i<m.Length(); i++)
			{
				temp.ElemNC(i+absk,i) = m[i];
			}
		}
		return temp;
	}

	
	template <class T>
	Vector<T> Diag(Matrix<T>& m)
	{
		return Diag(m,0);
	}

	
	template <class T>
	Vector<T> Diag(Matrix<T>& m, int k)
	{
		Vector<T> temp;
		int absk = abs(k);
		int dim;
		if(k>=0)
		{
			if(absk>=m.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("diagonal index k is outside the matrix dimensions!");
			}
			dim = m.Columns()-absk > m.Rows() ? m.Rows() : m.Columns()-absk;
			temp = Vector<T>(dim);
			for(int i=0; i<temp.Length(); i++)
			{
				temp[i] = m.ElemNC(i,i+absk);
			}
		}
		else
		{
			if(absk>=m.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("diagonal index k is outside the matrix dimensions!");
			}
			dim = m.Rows()-absk > m.Columns() ? m.Columns() : m.Rows()-absk;
			temp = Vector<T>(dim);
			for(int i=0; i<temp.Length(); i++)
			{
				temp[i] = m.ElemNC(i+absk,i);
			}
		}
		
		return temp;
	}


	
	template <class T>
	Matrix<T> FlipLR(Matrix<T>& m)
	{
		Matrix<T> temp(m.Rows(), m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			for(int i=0; i<m.Rows(); i++)
			{
				temp.ElemNC(i,m.Columns()-j-1) = m.ElemNC(i,j);
			}
		}
		return temp;
	}

	
	template <class T>
	Matrix<T>& FlipLRI(Matrix<T>& m)
	{
		for(int j=0; j<m.Columns()/2; j++)
		{
			for(int i=0; i<m.Rows(); i++)
			{
				 T dummy = m.ElemNC(i,j);
				 m.ElemNC(i,j) = m.ElemNC(i,m.Columns()-j-1);
				 m.ElemNC(i,m.Columns()-j-1) = dummy;
			}
		}
		return m;
	}


	
	template <class T>
	Matrix<T> FlipUD(Matrix<T>& m)
	{
		Matrix<T> temp(m.Rows(), m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			for(int i=0; i<m.Rows(); i++)
			{
				temp.ElemNC(m.Rows()-i-1,j) = m.ElemNC(i,j);
			}
		}
		return temp;
	}

	
	template <class T>
	Matrix<T>& FlipUDI(Matrix<T>& m)
	{
		for(int j=0; j<m.Columns(); j++)
		{
			for(int i=0; i<m.Rows()/2; i++)
			{
				T dummy = m.ElemNC(i,j);
				m.ElemNC(i,j) = m.ElemNC(m.Rows()-i-1,j);
				m.ElemNC(m.Rows()-i-1,j) = dummy;
			}
		}
		return m;
	}


	
	template <class T>
	Vector<T> Reverse(Vector<T>& m)
	{
		Vector<T> temp(m.Length());
		for(int i=0; i<m.Length(); i++)
		{
			temp.ElemNC(m.Length()-i-1) = m.ElemNC(i);
		}
		return temp;
	}

	
	template <class T>
	Vector<T>& ReverseI(Vector<T>& m)
	{
		for(int i=0; i<(m.Length()+1)/2; i++)
		{
			T dummy = m.ElemNC(i);
			m.ElemNC(i) = m.ElemNC(m.Length()-i-1);
			m.ElemNC(m.Length()-i-1) = dummy;
		}
		return m;
	}