core.inl

图像分割算法

	template <class T>
	Matrix<T> Rot90(Matrix<T>& m)
	{
		return Rot90(m,1);
	}

	
	template <class T>
	Matrix<T> Rot90(Matrix<T>& m, int k)
	{
		int rot = k % 4;
		if(rot<0)
		{
			rot += 4;
		}
		
		if(rot==0)
		{
			return m.Clone();
		}
		else if(rot==1)
		{
			return FlipUD(~m);
		}
		else if(rot==2)
		{
			return FlipLR(FlipUD(m));
		}
		else if(rot==3)
		{
			return FlipLR(~m);
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Something went wrong when trying to rotate the matrix!");
		}
	}
	
	template <class T>
	Vector<T> Maximum(Matrix<T>& m)
	{
		Vector<T> maxes(m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			maxes[j] = Maximum(m[j]);
		}
		return maxes;
	}
	
	template <class T>
	T Maximum(Vector<T>& m)
	{
		T max = m.ElemNC(0);
		for(int i=1; i<m.Length(); i++)
		{
			max = max >= m.ElemNC(i) ? max : m.ElemNC(i);
		}
		return max;
	}
	
	template <class T>
	Vector<T> Minimum(Matrix<T>& m)
	{
		Vector<T> mins(m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			mins[j] = Minimum(m[j]);
		}
		return mins;
	}
	
	template <class T>
	T Minimum(Vector<T>& m)
	{
		T min = m.ElemNC(0);
		for(int i=1; i<m.Length(); i++)
		{
			min = min <= m.ElemNC(i) ? min : m.ElemNC(i);
		}
		return min;
	}

	
	template <class T>
	Vector<double> Mean(Matrix<T>& m)
	{
		Vector<double> means(m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			means[j] = Mean(m[j]);
		}
		return means;
	}
	
	template <class T>
	double Mean(Vector<T>& m)
	{
		double mean = 0;
		for(int i=0; i<m.Length(); i++)
		{
			mean += m.ElemNC(i);
		}
		mean /= m.Length();
		return mean;
	}
	
	template <class T>
	Vector<T> Median(Matrix<T>& m)
	{
		Vector<T> medians(m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			medians[j] = Median(m[j]);
		}
		return medians;
	}
	
	template <class T>
	T Median(Vector<T>& m)
	{
		Vector<T> temp = m.Clone();
		T *arr = temp.Data();
		int n = temp.Length();
		int low, high;
		int median;
		int middle, ll, hh;

		low = 0;
		high = m.Length()-1;
		median = (low + high)/2;
		for(;;)
		{
			if(high <= low)
			{
				return arr[median];
			}
			if(high == low + 1)
			{
				if(arr[low] > arr[high])
				{
					Swap(arr[low], arr[high]);
				}
				return arr[median];
			}

			middle = (low + high)/2;
			if(arr[middle] > arr[high])
			{
				Swap(arr[middle], arr[high]);
			}
			if(arr[low] > arr[high])
			{
				Swap(arr[low], arr[high]);
			}
			if(arr[middle] > arr[low])
			{
				Swap(arr[middle], arr[low]);
			}

			Swap(arr[middle], arr[low+1]) ;

			ll = low + 1;
			hh = high;
			for(;;)
			{
				do ll++; while(arr[low] > arr[ll]);
				do hh--; while(arr[hh]  > arr[low]);

				if (hh < ll)
				{
					break;
				}

				Swap(arr[ll], arr[hh]);
			}

			Swap(arr[low], arr[hh]);

			if (hh <= median)
			{
				low = ll;
			}
			if (hh >= median)
			{
				high = hh - 1;
			}
		}
	}

	
	template <class T>
	Vector<double> Var(Matrix<T>& m)
	{
		return Var(m,true);
	}
	
	template <class T>
	double Var(Vector<T>& m)
	{
		return Var(m,true);
	}

	
	template <class T>
	Vector<double> Var(Matrix<T>& m, bool unbiased)
	{
		Vector<double> vars(m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			vars[j] = Var(m[j], unbiased);
		}
		return vars;
	}
	
	template <class T>
	double Var(Vector<T>& m, bool unbiased)
	{
		double mean = Mean(m);
		double var = 0;
		for(int i=0; i<m.Length(); i++)
		{
			var += (m.ElemNC(i)-mean)*(m.ElemNC(i)-mean);
		}
		
		if(unbiased)
		{
			return var/(m.Length()-1);
		}
		else
		{
			return var/m.Length();
		}
	}


	
	template <class T>
	Vector<double> Std(Matrix<T>& m)
	{
		return Std(m,true);
	}
	
	template <class T>
	double Std(Vector<T>& m)
	{
		return Std(m,true);
	}

	
	template <class T>
	Vector<double> Std(Matrix<T>& m, bool unbiased)
	{
		Vector<double> stds(m.Columns());
		for(int j=0; j<m.Columns(); j++)
		{
			stds[j] = Std(m[j], unbiased);
		}
		return stds;
	}
	
	template <class T>
	double Std(Vector<T>& m, bool unbiased)
	{
		return sqrt(Var(m,unbiased));
	}



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

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



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

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

	
	template <class T>
	Vector<T> Diff(Vector<T>& m)
	{
		if(m.Length() == 1)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Diff can only be called on vectors of length 2 or more!");
		}
		Vector<T> temp(m.Length()-1);
		for(int i=0;i<temp.Length(); i++)
		{
			temp.ElemNC(i) = m.ElemNC(i+1) - m.ElemNC(i);
		}
		return temp;
	}
	
	template <class T>
	Matrix<T> Diff(Matrix<T>& m)
	{
		if(m.Rows() == 1)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Diff can only be called on matrices with 2 or more rows!");
		}
		Matrix<T> temp(m.Rows()-1, m.Columns());
		for(int j=0; j<temp.Columns(); j++)
		{
			for(int i=0; i<temp.Rows(); i++)
			{
				
				temp.ElemNC(i,j) = m.ElemNC(i+1,j) - m.ElemNC(i,j);
			}
		}
		return temp;
	}



	
	template <class T>
	Matrix<T> RepMat(Matrix<T>& m, int M, int N)
	{
		Matrix<T> temp(M*m.Rows(),N*m.Columns());
		for(int k=0; k<M; k++)
		{
			for(int l=0; l<N; l++)
			{
				int s = k*m.Rows();
				int t = l*m.Columns();
				for(int i=0; i<m.Rows(); i++)
				{
					for(int j=0; j<m.Columns(); j++)
					{
						temp.ElemNC(s+i,t+j) = m.ElemNC(i,j);
					}
				}		
			}
		}
		return temp;
	}


	
	template <class T>
	Matrix<T> Reshape(Matrix<T>& m, int M, int N)
	{
		if(m.Length() != M*N)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Number of elements should be the same when using Reshape()!");
		}
		Matrix<T> temp(M,N);
		memcpy(temp.Data(), m.Data(), sizeof(T)*m.Length());
		return temp;
	}

	template <class T>
	Matrix<T> Reshape(Vector<T>& m, int M, int N)
	{
		if(m.Length() != M*N)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Number of elements should be the same when using Reshape()!");
		}
		Matrix<T> temp(M,N);
		memcpy(temp.Data(), m.Data(), sizeof(T)*m.Length());
		return temp;
	}
	
	template <class T>
	Vector<T>& QuickSort(Vector<T>& m)
	{
		int M = 7;
		int NSTACK = 50;
		
		int i, j, k; 
		int l=1;
		int ir = m.Length();
		Vector<int> stack(m.Length());
		int *istack = stack.Data()-1;
		int jstack = 0;
		
		T a;
		T* arr = m.Data()-1;
		
		for (;;) 
		{
			if (ir-l < M) 
			{
				for (j=l+1;j<=ir;j++) 
				{
					a=arr[j];
					for (i=j-1;i>=l;i--) 
					{
						if (arr[i] <= a) 
						{
							break;
						}
						arr[i+1]=arr[i];
					}
					arr[i+1]=a;
				}
				if (jstack == 0) 
				{
					break;
				}
				ir=istack[jstack--];
				l=istack[jstack--];
			} 
			else 
			{
				k=(l+ir) >> 1;
				Swap(arr[k],arr[l+1]);
				if (arr[l] > arr[ir]) 
				{
					Swap(arr[l],arr[ir]);
				}
				if (arr[l+1] > arr[ir]) 
				{
					Swap(arr[l+1],arr[ir]);
				}
				if (arr[l] > arr[l+1]) 
				{
					Swap(arr[l],arr[l+1]);
				}
				i=l+1; 
				j=ir;
				a=arr[l+1];
				for (;;) 
				{
					do i++; while (arr[i] < a); 
					do j--; while (arr[j] > a); 
					if (j < i) 
					{
						break; 
					}
					Swap(arr[i],arr[j]);
				}
				arr[l+1]=arr[j];
				arr[j]=a;
				jstack += 2;
				if (jstack > NSTACK) 
				{
					cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
					Utility::RunTimeError("Stack size is too small in quick sort!");
				}
				if (ir-i+1 >= j-l)
				{
					istack[jstack]=ir;
					istack[jstack-1]=i;
					ir=j-1;
				} 
				else
				{
					istack[jstack]=j-1;
					istack[jstack-1]=l;
					l=i;
				}
			}
		}
		return m;
	}
	
		
	template <class T>
	Vector<T> Sort(Vector<T>& m)
	{
		Vector<T> sorted = m.Clone();
		QuickSort(sorted);
		return sorted;
	}
	
	template <class T>
	Matrix<T> Sort(Matrix<T>& m, int dimension)
	{
		Matrix<T> temp;
		if(dimension == 1)
		{
			temp = m.Clone();
			for(int j=0; j<temp.Columns(); j++)
			{
				QuickSort(temp[j]);
			}
		}
		else if(dimension == 2)
		{
			temp = Matrix<T>(m.Rows(), m.Columns());
			for(int i=0; i<temp.Rows(); i++)
			{
				Vector<T> dummy(temp.Columns());
				for(int j=0; j<temp.Columns(); j++)
				{
					dummy.ElemNC(j) = m.ElemNC(i,j);
				}
				QuickSort(dummy);
				for(int j=0; j<temp.Columns(); j++)
				{
					temp.ElemNC(i,j) = dummy.ElemNC(j);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
		}
		return temp;
	}
	
	template <class T>
	Matrix<T> Sort(Matrix<T>& m)
	{
		return Sort(m,1);
	}




	
	template <class T>
	Vector<T>& SortI(Vector<T>& m)
	{
		return QuickSort(m);
	}
	
	template <class T>
	Matrix<T>& SortI(Matrix<T>& m, int dimension)
	{
		if(dimension == 1)
		{
			for(int j=0; j<m.Columns(); j++)
			{
				QuickSort(m[j]);
			}
		}
		else if(dimension == 2)
		{