analysis.cpp

图像分割算法

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



#include "./Analysis.h"
#include "mkl.h"


namespace Analysis
{

	void StatusDisplay(long status)
	{
		long classError;                                        
		
		classError = DftiErrorClass(status, DFTI_ERROR_CLASS);  
		if(! classError)
		{
			cerr << "Error Status is not a member of Predefined Error Class\n"; 
		}
		else
		{
			char* errorMessage = DftiErrorMessage(status);
			cerr << "Error_message = " << errorMessage << endl;
		}
	}


	Vector<ComplexFloat> FFT(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}

	Vector<ComplexFloat>& FFTI(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		//Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}




	Vector<ComplexDouble> FFT(Vector<ComplexDouble>& x)
	{
		int n = x.Length();
		Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}

		return y;
	}

	Vector<ComplexDouble>& FFTI(Vector<ComplexDouble>& x)
	{
		int n = x.Length();
		//Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}

		return x;
	}


	Vector<ComplexFloat> FFT(Vector<float>& x)
	{
		int n = x.Length();
		Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_REAL, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}
		
		for(int i=1;i<(n+1)/2;i++)
		{
			y[n-i] = conj(y[i]);
		}

		return y;
	}



	Vector<ComplexDouble> FFT(Vector<double>& x)
	{
		int n = x.Length();
		Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_REAL, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}
		
		for(int i=1;i<(n+1)/2;i++)
		{
			y[n-i] = conj(y[i]);
		}

		return y;
	}








	// ////////////////////
	// Backward transform
	// ////////////////////


	Vector<ComplexFloat> IFFT(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Backward transform
		status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}

	Vector<ComplexFloat>& IFFTI(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		//Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Backward transform
		status = DftiComputeBackward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)