analysis.cpp

图像分割算法

		{
			temp1 = input;
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}


		if(b == Symmetric)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexFloat> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
				temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexFloat> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
				temp1.ReadFromMatrix(c2, 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexFloat> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexFloat> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r2, input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}
		else if(b == Replicate)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexFloat> c1 = input.Slice(0, input.Rows()-1, 0, 0);
				temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexFloat> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexFloat> r1 = input.Slice(0, 0, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexFloat> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q11( hei1, wid1, input.ElemNC(0,0) );
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}

		
		// Allocate area for the filter
		Matrix<ComplexFloat> temp2;
		if(b == ZeroPad || b == Symmetric || b == Replicate)
		{
			temp2 = Matrix<ComplexFloat>(fftrows, fftcols,0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp2 = Matrix<ComplexFloat>(input.Rows(), input.Columns(),0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}


		// FFT and inverse FFT
		Matrix<ComplexFloat> o = IFFT2( FFT2(temp1)*FFT2(temp2) );

		// Based on the filter area, crop the border.
		if(fa == Same && b != Circular)
		{
			o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
		}
		else if(fa == Valid)
		{
			o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
		}
		else if(fa == Full && PowerOfTwo)
		{
			o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2);
		}

		//else if(fa != Full)
		//{
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Invalid FilterArea type!");
		//}
		
		return o;
		
	}


	
	Matrix<ComplexFloat> Conv2(Vector<ComplexFloat>& filter1, Vector<ComplexFloat>& filter2, Matrix<ComplexFloat>& input, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter sizes is smaller than the input size
		if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
		}
		
		Vector<ComplexFloat> *cols = new (std::nothrow) Vector<ComplexFloat>[input.Columns()];
		Utility::CheckPointer(cols);
		
		for(int i=0; i<input.Columns(); i++)
		{
			cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo);
		}
		
		int rlength = cols[0].Length();
		
		Vector<ComplexFloat> *rows = new (std::nothrow) Vector<ComplexFloat>[rlength];
		Utility::CheckPointer(rows);
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Vector<ComplexFloat>(input.Columns());
			for(int j=0; j<input.Columns(); j++)
			{
				rows[i][j] = cols[j][i];
			}
		}
		
		delete [] cols;
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo);
		}
		
		int clength = rows[0].Length();
		Matrix<ComplexFloat> temp(rlength, clength);
		for(int i=0; i<rlength; i++)
		{
			for(int j=0; j<clength; j++)
			{
				temp.ElemNC(i,j) = rows[i][j];
			}
		}
		
		delete [] rows;
		
		return temp;
	}




	Matrix<ComplexDouble> Conv2(Matrix<ComplexDouble>& input, Matrix<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter size is smaller than the input size
		if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
		}
		
		
		// Allocate area for input using the border
		
		int fftrows = input.Rows()+filter.Rows()-1;
		int fftcols = input.Columns()+filter.Columns()-1;
		if(PowerOfTwo)
		{
			fftrows = NextPow2(fftrows);
			fftcols = NextPow2(fftcols);
		}
		
		Matrix<ComplexDouble> temp1;
		if(b == ZeroPad || b == Symmetric || b == Replicate) 
		{
			temp1 = Matrix<ComplexDouble>(fftrows, fftcols,0);
			for(int j=0; j<input.Columns(); j++)
			{
				for(int i=0; i<input.Rows(); i++)
				{
					temp1.ElemNC(i,j) = input.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp1 = input;
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}
		
		if(b == Symmetric)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexDouble> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
				temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexDouble> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
				temp1.ReadFromMatrix(c2, 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexDouble> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexDouble> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r2, input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}
		else if(b == Replicate)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexDouble> c1 = input.Slice(0, input.Rows()-1, 0, 0);
				temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexDouble> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexDouble> r1 = input.Slice(0, 0, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexDouble> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q11( hei1, wid1, input.ElemNC(0,0) );
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}

		
		// Allocate area for the filter
		Matrix<ComplexDouble> temp2;
		if(b == ZeroPad || b == Symmetric || b == Replicate)
		{
			temp2 = Matrix<ComplexDouble>(fftrows, fftcols,0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp2 = Matrix<ComplexDouble>(input.Rows(), input.Columns(),0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}


		// FFT and inverse FFT
		Matrix<ComplexDouble> o = IFFT2( FFT2(temp1) * FFT2(temp2) );
		
		// Based on the filter area, crop the border.
		if(fa == Same && b != Circular)
		{
			o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
		}
		else if(fa == Valid)
		{
			o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
		}
		else if(fa == Full && PowerOfTwo)
		{
			o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2);
		}
		//else if(fa != Full)
		//{
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Invalid FilterArea type!");
		//}
		
		return o;		
	}
	
	Matrix<ComplexDouble> Conv2(Vector<ComplexDouble>& filter1, Vector<ComplexDouble>& filter2, Matrix<ComplexDouble>& input, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter sizes is smaller than the input size
		if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
		}
		
		Vector<ComplexDouble> *cols = new (std::nothrow) Vector<ComplexDouble>[input.Columns()];
		Utility::CheckPointer(cols);
		
		for(int i=0; i<input.Columns(); i++)
		{
			cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo);
		}
		
		int rlength = cols[0].Length();
		
		Vector<ComplexDouble> *rows = new (std::nothrow) Vector<ComplexDouble>[rlength];
		Utility::CheckPointer(rows);
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Vector<ComplexDouble>(input.Columns());
			for(int j=0; j<input.Columns(); j++)
			{
				rows[i][j] = cols[j][i];
			}
		}
		
		delete [] cols;
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo);
		}
		
		int clength = rows[0].Length();
		Matrix<ComplexDouble> temp(rlength, clength);
		for(int i=0; i<rlength; i++)
		{
			for(int j=0; j<clength; j++)
			{
				temp.ElemNC(i,j) = rows[i][j];
			}
		}
		
		delete [] rows;
		
		return temp;
	}





};