matrix.c

图像处理的压缩算法

/*------------------------------------------------------------------------------*
 * File Name:matrix.c				 											*
 * Creation:July 22th, 2002														*
 * Purpose: OriginC Source C file												*
 * Copyright (c) OriginLab Corp.2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 * TCZ 07/29/02 QA70-2497	ADD_GLOBAL_STATS_SUMMARY_FUNCTION					*
 * ER 01/29/03 QA70_3802	ADD_VECTOR_FFT_IFFT
 * CPY 7/7/03 changed to include only origin.h to take advantage of PCH			
 *------------------------------------------------------------------------------*/
 
#include <origin.h> // main Origin C header that is precompiled and already include most headers 

// this file include most of the other header files except the NAG header, which takes longer to compile
// NAG routines
#include <NAG\OCN_c06.h> // this contains the FFT NAG headers, 
#include <NAG\OCN_f.h>//this contain the linear algebra function
#include <NAG\OCN_g01.h>//this file contain the simple statistics function
////////////////////////////////////////////////////////////////////////////////////
#include <matrix.h>
// start your functions here


//local functions for the internal use

 
// This function will do the 2d fft with the data form required by nag funtion.
static int fft2_for_complex_matrix(int nNumRows, int nNumCols,  matrix& matReal, matrix& matImage, matrix<complex>& matFFT2Result);


//  	This function will take a matrix and set the result matrix accoring to the parameter values.
//	in other words, it will do the neccessory padding with zero or truncting.
static int set_matrix_with_padding_truncting(matrix & matSource, matrix & matResult, int nRowProcess , int nColProcess);


//   This function will take care of the complex data matrix and sperate it into two double matrix.
static int seperate_complex_mat_to_real_and_imag_mat(matrix<complex>& matSource, matrix& matReal, matrix&  matImage, int nRowProcess, int nColProcess);

/**

	int FFT2(matrix& matSource, matrix<complex>& matFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
	int FFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
			These two Functions will do the 2 dimension Fast Fourier Transform.
	Remarks:
			Right now, we only support the data matrix in double and complex form. If user have the data
			matrix in integer form, user should first transform it to double. It is a member function of matrix.
			Also, if user has a vector data, it should firstly be transformed into a matrix and then do the 
			FFT.
	Example:
	Parameters:
			matSource: The data source matrix
			matFFT2Result: The matrix containing the 2 dimension FFT results
			nRowProcess:
			nColProcess: The number of rows or columns user provided. If they are -1, there is no padding
	 					or truncating. If they are more than the sizes of the matrix, the function will pad zero to meet the size
						If less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number. Must be -1 or positive integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The column or the row size of matrix is less then 1
 			 73: System failed to allocate memory

*/
int FFT2(matrix& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess)
{
	matrix matReal;
	matrix matTemp(matSource);
	int nRet;
	if(nRet = set_matrix_with_padding_truncting(matTemp, matReal, nRowProcess, nColProcess))
		return nRet;
	matrix matImage;
		
	int nNumCols = matReal.GetNumCols();
	int nNumRows = matReal.GetNumRows();
	matImage.SetSize(nNumRows, nNumCols);
	matImage = 0;	
	if(nRet = fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImage, matFFT2Result))
		return nRet;
	return nRet;
}

static int fft2_for_complex_matrix(int nNumRows, int nNumCols,  matrix& matReal, matrix& matImage, matrix<complex>& matFFT2Result)
{
	vector vecTrigRows;
	vecTrigRows.SetSize(2 * nNumRows);
	//cal nage function to prepare the trig value
	int nRet;
	if( nRet = nag_fft_init_trig(nNumRows, vecTrigRows))
		return nRet;
	
	vector vecTrigCols;
	vecTrigCols.SetSize(2 * nNumCols);
	if( nRet = nag_fft_init_trig(nNumCols, vecTrigCols))
		return nRet;
	
	//call the nag function to do the fft2
	if(nRet = nag_fft_2d_complex(nNumRows, nNumCols, matReal, matImage, vecTrigRows, vecTrigCols))
		return nRet;
	//set the data back
	if( nRet = matFFT2Result.MakeComplex(matReal, matImage))
		return nRet;
	
}
		
int FFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess )
{
	int nRet;
	matrix matReal, matImage;
	if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImage, nRowProcess, nColProcess))
		return nRet;
	
	int nNumCols = matImage.GetNumCols();
	int nNumRows = matImage.GetNumRows();
	if(nRet = fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImage, matFFT2Result))
		return nRet;
	
}
		
static int set_matrix_with_padding_truncting(matrix & matSource, matrix & matResult, int nRowProcess , int nColProcess)
{
		
	
	if(nRowProcess < -1 || nRowProcess == 0 ||nColProcess < -1 ||nColProcess == 0)
		return FFT_ERROR_COULUMN_ROW_NUMBER_NEGATIVE;
	int nNumCols = matSource.GetNumCols();
	int nNumRows = matSource.GetNumRows();
	if(nNumCols == 0 || nNumRows == 0)
		return FFT_ERROR_SOURCE_MATRIX_EMPTY;
		
	int nIntendNumCols = nNumCols;
	int nIntendNumRows = nNumRows;
	if(nRowProcess != -1) 
		nIntendNumRows = nRowProcess;
	if(nColProcess != -1)
		nIntendNumCols = nColProcess;
	
	int nRet;
	matResult = matSource;
	
	if(nIntendNumRows > nNumRows)//so we should increase the dimension and do the padding
	{
		matrix matTemp;
		matTemp.SetSize(nIntendNumRows - nNumRows, nNumCols );
		matTemp = 0;
		if(nRet = matSource.Concatenate(1,matTemp, matResult))
			return nRet;
	}
	BOOL bOK = TRUE;
	if (nIntendNumRows < nNumRows)
    	bOK = matSource.GetSubMatrix(matResult, 0, -1, 0, nIntendNumRows - 1);
	if(!bOK)
		return FFT_ERROR_GETSUBMATRIX_FUNCTION_ERROR;
    		
	    
    if(nIntendNumCols > nNumCols)
	{
		matrix matTemp;
		matTemp.SetSize(nIntendNumRows, nIntendNumCols - nNumCols );
		matTemp = 0;
		if(nRet = matResult.Concatenate(2, matTemp, matResult))
			return nRet;
		
	}
	
	if(nIntendNumCols < nNumCols)
	{
		matrix matResultTemp(matResult); 
		bOK = matResultTemp.GetSubMatrix(matResult, 0, nIntendNumCols -1, 0, -1);
		if(!bOK)
			return FFT_ERROR_GETSUBMATRIX_FUNCTION_ERROR;
	}
	return 0;
}

/*
	int IFFT2(matrix<complex>& matSource, matrix<complex>& matIFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
	Remarks:
			Right now, we only support the data matrix in complex form. If user have the data
			matrix in double form, user should add the imaginary part as zero to make a complex matrix.
			Matrix class has a member function to make the complex matrix. Also, if user has a vector data,
			it should firstly be transformed into a matrix and then one can do the inverse FFT.
	Example:
	Parameters:
			matSource: The data source matrix
			matIFFT2Result: The matrix containing the 2 dimension inverse FFT results
			nRowProcess:
			nColProcess:The number of rows or columns user provided. if they are -1, there are no padding
	 					or truncating. if they are larger than the size of the matrix, the function will pad zero to meet
						the size. If less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number. Must be -1 or positive integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The column or the row size of matrix is less then 1
 			 73: System failed to allocate memory

  
*/
int IFFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess)
{

	int nRet;
	matrix<complex> matSourceCopy(matSource);
	if(nRet = matSourceCopy.Conjugate())
		return nRet;
	
	if(nRet = FFT2(matSourceCopy, matFFT2Result,nRowProcess, nColProcess ))
		return nRet;
	
	if(nRet = matFFT2Result.Conjugate())
		return nRet;
}		
/*
	int FFT(matrix& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
	int FFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);

				These two Functions will do the 1 dimension Fast Fourier Transform.
  	Remarks:
			Right now, we only support the data matrix in double and complex form. If user has the data
			matrix in integer form, user should firstly transform it to double. It is a member function of matrix.
			Also, if user has a vector data, he should firstly transform data into a matrix and then one can do the 
			FFT. Unlike the matlab, this function will do the FFT on each row (Matlab do FFT on each column).
	Example:
	Parameters:
			matSource: The data source matrix
			matFFTResult: The matrix containing the 1 dimension FFT results
			nColProcess: The number of columns user provided. If it is -1, function will not do any padding
	 					or truncating. But if it is larger than the size of the  matrix, it will pad zero to meet the size
						if less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number. Must be -1 or positive integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The column or the row size of matrix is less then 1
 			 73: System failed to allocate memory

*/

int FFT(matrix& matSource, matrix<complex>& matFFTResult, int nColProcess)
{
	matrix matTemp;
	int nRet;
	if(nRet = set_matrix_with_padding_truncting(matSource, matTemp, -1, nColProcess))
		return nRet;
	int nNumCols = matTemp.GetNumCols();
	int nNumRows = matTemp.GetNumRows();
	vector vecTrig;
	vecTrig.SetSize(2 * nNumCols); 
	if(nRet = nag_fft_init_trig(nNumCols, vecTrig))
		return nRet;
	if(nRet = nag_fft_multiple_real(nNumRows, nNumCols, matTemp, vecTrig))
		return nRet;
	matrix matReal, matImage;
	matReal.SetSize(nNumRows, nNumCols);
	matImage.SetSize(nNumRows, nNumCols);
	
	if(nRet = nag_multiple_hermitian_to_complex(nNumRows, nNumCols, matTemp, matReal, matImage))
		return nRet;
	if(nRet = matFFTResult.MakeComplex(matReal, matImage))
		return nRet;
	return nRet;
}
		
int FFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess)
{
	int nRet;
	matrix matReal, matImage;
	if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImage, -1, nColProcess))
		return nRet;
	int nNumCols = matImage.GetNumCols();
	int nNumRows = matImage.GetNumRows();
	vector vecTrig;
	vecTrig.SetSize(2 * nNumCols); 
	if(nRet = nag_fft_init_trig(nNumCols, vecTrig))
		return nRet;
	if(nRet = nag_fft_multiple_complex(nNumRows, nNumCols, matReal, matImage, vecTrig ))
		return nRet;
	if(nRet = matFFTResult.MakeComplex(matReal, matImage))
		return nRet;
	return nRet;
}
/*	  
	int IFFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
	
			This Function will do the 1 dimension inverse Fast Fourier Transform.
	Remarks:
			Right now, we only support the data matrix in complex form. If user has the data
			matrix in double form, user should add the imaginary part as zero to make a complex matrix.
			Matrix class has a member function to make the complex matrix. also, if user has a vector data,
			it should be firstly transformed into a matrix and then one can do the inverse FFT.

	Example:
	Parameters:
			matSource:The data source matrix
			matIFFT2Result: The matrix containing the 1 dimension inverse FFT results
			nColProcess:The number of columns user provided. if it is -1, there is no padding
	 					or truncting. if it is more than the size of the  matrix, it will pad zero to meet the size
						if less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number must be -1 or positve integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The coloumn or the row size of matrix is less then 1;
 			 73: System failed to allocate memory

  
*/


		
int IFFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess)
{
	int nRet;
	matrix<complex> matSourceCopy(matSource);
	if(nRet = matSourceCopy.Conjugate())
		return nRet;