/*-----------------------------------------------------------------------------*
 * File Name: pro_mat_utils	 												   *
 * Creation: Derek 07/14/2010												   *
 * Purpose: implementation of general utility functions for matrix calculation *
 * Copyright (c) OriginLab Corp.2010                                     	   *
 * All Rights Reserved														   *
 * 																			   *
 * Modification Log:														   *
 * Derek 07/14/2010 2D_FFT_FILTERING                    					   *
 * Derek 07/28/2010 ADD_2D_THRESHOLD_FILTER                                    *
 *-----------------------------------------------------------------------------*/



#include <origin.h>
#include <..\OriginLab\pro_mat_utils.h>
#include <fft_utils.h>


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

// Derek 07/14/2010 2D_FFT_FILTERING
int fft_band_filter_2D(const matrix& mi, const matrix& mf, matrix& mo)
{

	// size mismatch
    if( (mi.GetNumRows() != mf.GetNumRows()) || (mi.GetNumCols() != mf.GetNumCols()) || 
		(mi.GetNumRows() != mo.GetNumRows()) || (mi.GetNumCols() != mo.GetNumCols()) )
	{
        return MISMATCHED_SIZE;
	}

    matrix<complex> mFFT, mIFFT;

    // Compute 2D FFT
    if( FFT2(mi, mFFT) != 0 )
        return FFT_ERROR;
	
    // Shift FFT result to place DC component in the center
    if( mFFT.FFTShift(mFFT) != 0)
        return FFT_SHIFT_ERROR;
	
    // Filtering with the filter
    if( mFFT.DotMultiply(mf) == false )
        return MATRIX_DOT_MULTIPLY_ERROR;
	
    // Inverse FFT Shift
    if( mFFT.IFFTShift(mFFT) != 0 )
        return IFFT_SHIFT_ERROR;
	
    // Compute 2D inverse FFT
    if( IFFT2(mFFT, mIFFT) != 0 )
        return IFFT_ERROR;
	
    // Get the amplitude of inverse FFT
    if( mIFFT.GetAmplitude(mo) == false )
        return GET_AMPLITUDE_ERROR;
	
    return NO_ERROR;
}
// End of 2D_FFT_FILTERING



// Derek 07/28/2010 ADD_2D_THRESHOLD_FILTER
int fft_threshold_filter_2D(const matrix& mi, const double fLowerThresh, const double fUpperThresh, matrix& mo)
{

	// size mismatch
    if( (mi.GetNumRows() != mo.GetNumRows()) || (mi.GetNumCols() != mo.GetNumCols()))
	{
        return MISMATCHED_SIZE;
	}
	
	if( fLowerThresh < 0 || fLowerThresh > 1 || fUpperThresh < 0 || fUpperThresh > 1 || fUpperThresh < fLowerThresh )
	{
		return INVALID_THRESHOLD_VALUE;
	}
	
	matrix<complex> mFFT, mIFFT;
	matrix mFFTAmp;
	
    if( FFT2(mi, mFFT) != 0 )
        return FFT_ERROR;
    
    if( mFFT.GetAmplitude(mFFTAmp) == false )
        return GET_AMPLITUDE_ERROR;
    
    double fMax = mFFTAmp.GetMax();
    
    mFFTAmp.ReplaceLessThan(fMax*fLowerThresh,0);
    mFFTAmp.ReplaceGreaterThan(fMax*fUpperThresh,0);
    mFFTAmp.Replace(fMax*fLowerThresh,1, MATREPL_TEST_GREATER | MATREPL_TEST_EQUAL);
    
    if( mFFT.DotMultiply(mFFTAmp) == false )
        return MATRIX_DOT_MULTIPLY_ERROR; 

    if( IFFT2(mFFT, mIFFT) != 0 )
        return IFFT_ERROR;
	
    if( mIFFT.GetAmplitude(mo) == false )
        return GET_AMPLITUDE_ERROR;    
}
// End ADD_2D_THRESHOLD_FILTER

////////////////////////////////////////////////////////////////////////////////
//////////////////////////////End of file///////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////