surfarray.cpp

A set of C++ and Matlab routines implementing the surfacelet transform surfacelet的一个非常好用的工具箱

//////////////////////////////////////////////////////////////////////////
//	SurfBox-C++ (c) 
//////////////////////////////////////////////////////////////////////////
//
//	Yue Lu and Minh N. Do
//
//	Department of Electrical and Computer Engineering
//	Coordinated Science Laboratory
//	University of Illinois at Urbana-Champaign
//
//////////////////////////////////////////////////////////////////////////
//
//	SurfArray.cpp
//	
//	First created: 03-13-06
//	Last modified: 04-10-06
//
//////////////////////////////////////////////////////////////////////////

#include <cassert>
#include "SurfArray.h"
#include "fftw3.h"
#include <math.h>

extern void MessageLog(const char*, const char*, const char*);


//////////////////////////////////////////////////////////////////////////
//  Class constructor
//////////////////////////////////////////////////////////////////////////

SurfArray::SurfArray(void)
{
	// Initialize the data structures
    pData = NULL;
	pDims_logical = pDims_real = NULL;

    nDims = 0;

	// Default: the array constains real-valued data;
	IsReal = true;
}


//////////////////////////////////////////////////////////////////////////
//  Class destructor
//////////////////////////////////////////////////////////////////////////

SurfArray::~SurfArray(void)
{
	// pData is allocated by fftw, so it should be de-allocated by fftw too.
    if (pData) fftw_free(pData);

	if (pDims_real) delete [] pDims_real;

	if (pDims_logical) delete [] pDims_logical;
}


//////////////////////////////////////////////////////////////////////////
//  Reset the current array object, free resources
//////////////////////////////////////////////////////////////////////////

void SurfArray::Reset()
{
    if (pData) fftw_free(pData);

    if (pDims_real) delete [] pDims_real;

    if (pDims_logical) delete [] pDims_logical;

    pData = NULL;
    pDims_logical = pDims_real = NULL;

    // Default: the array constains real-valued data;
    IsReal = true;

    nDims = 0;
}


//////////////////////////////////////////////////////////////////////////
//  Get the pointer to the data array
//////////////////////////////////////////////////////////////////////////
//
// PARAMETERS:
//
// rPadding: 
//      return the internal padding value

double* SurfArray::GetPointer(int& rPadding)
{
	rPadding = padding;

	return pData;
}


//////////////////////////////////////////////////////////////////////////
//  Tell whether the current object is real-valued
//////////////////////////////////////////////////////////////////////////

bool SurfArray::IsRealValued()
{
	return IsReal;
}


//////////////////////////////////////////////////////////////////////////
//  Set the value domain
//////////////////////////////////////////////////////////////////////////

void SurfArray::SetIsRealValued(bool RealValue)
{
	IsReal = RealValue;
}


//////////////////////////////////////////////////////////////////////////
//  Get the number of dimension of the current array
//////////////////////////////////////////////////////////////////////////

int SurfArray::GetRank()
{
	return nDims;
}


//////////////////////////////////////////////////////////////////////////
//  Get the logical dimensions
//////////////////////////////////////////////////////////////////////////

void SurfArray::GetDims(int pDims[])
{
	for (register int i = 0; i < nDims; i++)
	{
		pDims[i] =  pDims_logical[i];
	}
}


//////////////////////////////////////////////////////////////////////////
//  Allocate necessary memory spaces
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  N:
//      dimension of the array
//
//  dims:
//      an array of length N indicating the size of the array along each dimension

void SurfArray::AllocateSpace(int N, int* dims)
{
	nDims = N;

	// tricky: we use integer division here
	padding = 2 * (dims[nDims - 1] / 2 + 1) - dims[nDims - 1];

	int nPoints, i;

    // total number of data points
	nPoints = 1;
	for (i = 0; i <= nDims - 2; i++)
	{
		nPoints *= dims[i];
	}
	nPoints *= (dims[nDims - 1] + padding);
	
	// Allocate memory
	try
	{
		pDims_real = new int [nDims];
		pDims_logical = new int [nDims];
		// This thing is what we are really worrying about, considering its size ...
		pData = (double *)fftw_malloc(sizeof(double) * nPoints);

		if (!pData)
			throw 1;
	}
	catch (...)
	{
		MessageLog("SurfArray", "AllocateSpace", "Insufficient memory!");
		
		if (pDims_real)
		{
			delete [] pDims_real;
			pDims_real = NULL;
		}
		if (pDims_logical)
		{
			delete [] pDims_logical;
			pDims_logical = NULL;
		}
		if (pData) 
		{
			fftw_free(pData);
			pData = NULL;
		}	

		throw;
	}	

	// Note: we adopt the "row-major" convension in addressing 
	// multidimensional arrays.
	for (i = 0; i < nDims; i++)
	{
		pDims_real[i] = pDims_logical[i] = dims[i];
	}
	pDims_real[nDims - 1] += padding;

}


//////////////////////////////////////////////////////////////////////////
//  Apply forward FFT on the current array (in-place version)
//////////////////////////////////////////////////////////////////////////
//
//  NOTE:
//
//  We assume the input array is real-valued, so the output of FFT has Hermitian symmetry.
//  We only keep nearly half of the FFT samples, along the last (continuous) dimension.
//
//  PARAMETERS:
//
//  mode:
//      0: do not normalize
//      1: normalize by 1 / sqrt(n_total_points)

void SurfArray::GetInPlaceForwardFFT(int mode /* = 0 */)
{
	assert((pData != NULL) && (pDims_real != NULL) && (pDims_logical != NULL));
	
	assert(IsReal);

	fftw_plan planFFT;
	planFFT = fftw_plan_dft_r2c(nDims, pDims_logical, pData, (fftw_complex *)pData, FFTW_ESTIMATE);
	
	fftw_execute(planFFT);
	
	fftw_destroy_plan(planFFT);

	IsReal = false;

    double *ptr;
    int nPoints;
    double normalizer;

    switch(mode)
    {
    case 0:
        // do nothing
    	break;
    
    case 1:
        register int i;
        normalizer = 1.0;
        nPoints = 1;
        for (i = 0; i < nDims; i++)
        {
            normalizer /= pDims_logical[i];
            nPoints *= pDims_real[i];
        }
        normalizer = sqrt(normalizer);
        ptr = pData;
        for (i = 0; i < nPoints; i++)
            *(ptr++) *= normalizer;    
        
        break;

    default:
        // this should not happen
        assert(false);
        break;
    }

	return;
}


//////////////////////////////////////////////////////////////////////////
//  Apply backward FFT on the current array
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  mode:
//      0: normalize by 1 / n_total_points
//      1: normalize by 1 / sqrt(n_total_points)
//      2: do not normalize

void SurfArray::GetInPlaceBackwardFFT(int mode /* = 0 */)
{
	assert((pData != NULL) && (pDims_real != NULL) && (pDims_logical != NULL));
	
	assert(!IsReal);
		
	fftw_plan planFFT;
	planFFT = fftw_plan_dft_c2r(nDims, pDims_logical, (fftw_complex *)pData, pData, FFTW_ESTIMATE);

	fftw_execute(planFFT);

	fftw_destroy_plan(planFFT);

	IsReal = true;

    if ((mode == 0) || (mode == 1))
    {     
        // Normalize the FFT result;
        register int i;
        double normalizer = 1.0;
        int nPoints = 1;
        for (i = 0; i < nDims; i++)
        {
            normalizer /= pDims_logical[i];
            nPoints *= pDims_real[i];
        }

        if (mode == 1)
            normalizer = sqrt(normalizer);

        double *ptr = pData;
        for (i = 0; i < nPoints; i++)
            *(ptr++) *= normalizer;    
    }
    else if ( mode != 2)
    {
        assert(false);
    }

	return;
}


//////////////////////////////////////////////////////////////////////////
//  Copy an array object
//////////////////////////////////////////////////////////////////////////
//
//  NOTE:
//
//  The array being copied to will allocate its own memory space.

SurfArray& SurfArray::operator = (SurfArray& Src)
{
    // free the current resource if there is any ...
    Reset();

    int *pDims;
    int N;
    try
    {
        N = Src.nDims;
        pDims = new int [N];
        Src.GetDims(pDims);
        AllocateSpace(N, pDims);
    }
    catch (std::bad_alloc)
    {
        MessageLog("SurfArray", "Operator =", "Out of memory");
        throw;
    }

    int i, n_total_points = 1;
    for (i = 0; i < N; i++)
        n_total_points *= pDims_real[i];

    // copy data values
    int padding;
    memcpy(pData, Src.GetPointer(padding), n_total_points * sizeof(double));

    // copy array status
    IsReal = Src.IsRealValued();

    delete [] pDims;
    return *this;
}


//////////////////////////////////////////////////////////////////////////
// Fill the internal data buffer with data pointed to by pIn.
// Appropriate padding will be taken care of.
//////////////////////////////////////////////////////////////////////////

void SurfArray::ImportRealValues(double* pIn)
{
	assert(IsReal);
	
	register int i, j;
	int nRows, lengthRow;

	nRows = 1;
	for (i = 0; i <= nDims - 2; i++)
	{
		nRows *= pDims_logical[i];
	}

	double *pVal = pData;
	lengthRow = pDims_logical[nDims - 1];

	for (i = 0; i < nRows; i++)
	{
		for (j = 0; j < lengthRow; j++)
			*(pVal++) = *(pIn++);

		pVal += padding;
	}
}


//////////////////////////////////////////////////////////////////////////
//	Export the internal data values.
//////////////////////////////////////////////////////////////////////////
//
//	Appropriate padding will be taken care of.

void SurfArray::ExportRealValues(double* pOut)
{
	assert(IsReal);
	
	register int i, j;
	int nRows, lengthRow;
		
	nRows = 1;
	for (i = 0; i < nDims - 1; i++)
	{
		nRows *= pDims_logical[i];
	}

	double *pVal = pData;
	lengthRow = pDims_logical[nDims - 1];
	
	for (i = 0; i < nRows; i++)
	{
		for (j = 0; j < lengthRow; j++)
			*(pOut++) = *(pVal++);

		pVal += padding;
	}
}


//////////////////////////////////////////////////////////////////////////
//	Get the magnitude frequency response
//////////////////////////////////////////////////////////////////////////

void SurfArray::GetMagnitudeResponse()
{
	// The array must be in the Fourier domain
	assert(!IsReal);

	register int i, j;
	int nRows, lengthRow;

	// Total number of rows
	nRows = 1;
	for (i = 0; i < nDims - 1; i++)
	{