nddirectionalfilterbank.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
//
//////////////////////////////////////////////////////////////////////////
//
//	NdDirectionalFilterBank.cpp
//	
//	First created: 04-02-06
//	Last modified: 04-11-06
//
//////////////////////////////////////////////////////////////////////////

#include "NdDirectionalFilterBank.h"
#include "math.h"
#include <sstream>
#include <iostream>
#include <fstream>
#include <cassert>
#include "fftw3.h"

using namespace std;

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

//////////////////////////////////////////////////////////////////////////
//	Constructor
//////////////////////////////////////////////////////////////////////////

NdDirectionalFilterBank::NdDirectionalFilterBank()
{
	// do nothing
}


//////////////////////////////////////////////////////////////////////////
//	Destructor
//////////////////////////////////////////////////////////////////////////

NdDirectionalFilterBank::~NdDirectionalFilterBank()
{
	// do nothing
}


//////////////////////////////////////////////////////////////////////////
//  NDFB Decomposition
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  InArray
//      Input array
//
//  OutArrays
//      List of output subbands
//
//  OutputInFourierDomain
//      Specifying which domain the output subband should be in.
//
//  Levels
//      Decomposition level at each dimension
//
//  filter0, filter1
//      The two decomposition filters
//
//	center0, center1
//		The centers of the decomposition filters
//  
//  mSize
//      Diamond filter size
//
//  beta
//		For Kaiser window
//
//  lambda
//		See paper for more detail

void NdDirectionalFilterBank::GetDecomposition(SurfArray &InArray, SurfArray *OutArrays[], bool OutputInFourierDomain, 
    int Levels[], SurfMatrix &filter0, int center0[], SurfMatrix &filter1, int center1[], int mSize /* = 15 */, double beta /* = 2.5 */, double lambda /* = 4.0 */)
{
    int nDims = InArray.GetRank();
    
    assert(nDims >= 2);

    int *pDims = new int [nDims];
    InArray.GetDims(pDims);

    // Hourglass filter bank decomposition
    SurfArray *subbands = new SurfArray [nDims];
    hourglass_fb.GetDecomposition(InArray, subbands, true, mSize, beta, lambda);

	int i, j, full_dim, nSubs;
    int *pLevel = Levels;
    double *pSubbandData;

	// Note: we scale the filter coefficients (a small 2-D array)
	// by 0.5. By doing this, we avoid having to scale the
	// large N-D array by 0.5 in the downsampling stage.
	double *pf;
	int nCoeffs;
	pf = filter0.GetPointer();
	nCoeffs = filter0.nx * filter0.ny;
	for (i = 0; i < nCoeffs; i++)
	{
		*pf /= 2.0;
		pf++;
	}

	pf = filter1.GetPointer();
	nCoeffs = filter1.nx * filter1.ny;
	for (i = 0; i < nCoeffs; i++)
	{
		*pf /= 2.0;
		pf++;
	}

    for (i = 0; i < nDims; i++)
    {
        // number of output subbands
		nSubs = 1;
		for (j = 0; j < nDims; j++)
			nSubs *= ((j == i)? 1 : (int)pow(2.0, (double)pLevel[j]));

		OutArrays[i] = new SurfArray [nSubs];

		if (nSubs > 1)		
		{
			// pSubbandData will point to the newly-allocated memory space
			pSubbandData = subbands[i].NewHermitianSymmetryAxis(i);
	        
			// One of the dimensions is nearly reduced by half
			full_dim = pDims[i];
			pDims[i] = full_dim / 2 + 1;

			// release memory space
			subbands[i].Reset();

			// Iteratively resampled checkerboard filter bank decomposition
			IrcDecomposition(nDims, pDims, pSubbandData, i, full_dim, pLevel, filter0, center0, filter1, center1);

			// extract all the subbands
			BoxToCell(nDims, pDims, pSubbandData, OutArrays[i], pLevel, i, full_dim);

			delete [] pSubbandData;

			// restore the full dimension
			pDims[i] = full_dim;
		}
		else
		{
			assert(nSubs == 1);
			*OutArrays[i] = subbands[i];
		}

		pLevel += nDims;

        if (!OutputInFourierDomain)
        {
            for (j = 0; j < nSubs; j++)
                OutArrays[i][j].GetInPlaceBackwardFFT();
        }

    }
    
	// restore the filter coefficients
	pf = filter0.GetPointer();
	nCoeffs = filter0.nx * filter0.ny;
	for (i = 0; i < nCoeffs; i++)
	{
		*pf *= 2.0;
		pf++;
	}

	pf = filter1.GetPointer();
	nCoeffs = filter1.nx * filter1.ny;
	for (i = 0; i < nCoeffs; i++)
	{
		*pf *= 2.0;
		pf++;
	}


    FilterSlice0.Reset();
    FilterSlice1.Reset();
    delete [] pDims;
    delete [] subbands;

}


//////////////////////////////////////////////////////////////////////////
//  NDFB Reconstruction
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  InArrays
//      List of input arrays
//
//  OutArray
//      The output array
//
//  OutputInFourierDomain
//      Specifying which domain the output subband should be in.
//
//  Levels
//      Decomposition level at each dimension
//
//  filter0, filter1
//      The two reconstruction filters
//
//	center0, center1
//		The centers of the reconstruction filters
//  
//  mSize
//      Diamond filter size
//
//  beta
//		For Kaiser window
//
//  lambda
//		See paper for more details

void NdDirectionalFilterBank::GetReconstruction(SurfArray *InArrays[], SurfArray &OutArray, bool OutputInFourierDomain, 
	int Levels[], SurfMatrix &filter0, int center0[], SurfMatrix &filter1, int center1[], int mSize /* = 15 */, double beta /* = 2::5 */, double lambda /* = 4::0 */)
{
	int nDims = InArrays[0][0].GetRank();

	assert(nDims >= 2);

	SurfArray *subbands = new SurfArray [nDims];
    int *pDims = new int [nDims];
	InArrays[0][0].GetDims(pDims);
        
    int i, j, full_dim, nSubs;
	int *pLevel = Levels;
    double *pSubbandData;

	for (i = 0; i < nDims; i++)
	{
		pDims[i] *= ((i == 0)? 1 : ((int)pow(2.0, (double)pLevel[i])));
	}

	for (i = 0; i < nDims; i++)
	{
		nSubs = 1;
		for (j = 0; j < nDims; j++)
			nSubs *= ((j == i)? 1 : ((int)pow(2.0, (double)pLevel[j])));

        for (j = 0; j < nSubs; j++)
        {
            if (InArrays[i][j].IsRealValued())
                InArrays[i][j].GetInPlaceForwardFFT();
        }

		if (nSubs > 1)
		{
			pSubbandData = CellToBox(InArrays[i], i, pLevel);
	        
			full_dim = pDims[i];
			pDims[i] = full_dim / 2 + 1;

			IrcReconstruction(nDims, pDims, pSubbandData, i, full_dim, pLevel, filter0, center0, filter1, center1);

			// restore the dimension
			pDims[i] = full_dim;

			subbands[i].AllocateSpace(nDims, pDims);
			subbands[i].SetIsRealValued(false);
			subbands[i].RestoreHermianSymmetryAxis(pSubbandData, i);

			delete [] pSubbandData;
		}
		else
		{
			assert(nSubs == 1);
			subbands[i] = *(InArrays[i]);
		}

        pLevel += nDims;
	
	}

    FilterSlice0.Reset();
    FilterSlice1.Reset();
	
	// Hourglass filter bank reconstruction
	hourglass_fb.GetReconstruction(subbands, OutArray, OutputInFourierDomain, mSize, beta, lambda);	

	// release the temporary memory
    delete [] pDims;
	delete [] subbands;
    
}


//////////////////////////////////////////////////////////////////////////
//	Read the checkerboard filter bank coefficients from files
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	bo:
//		Order of the checkerboard filters
//
//	IsDecomposition:
//		Decomposition filters or reconstruction filters
//
//	filter0:
//		The imported first filter, stored in a SurfMatrix object
//
//	center0:
//		The center of the first filter
//
//	filter1:
//		The imported second filter, stored in a SurfMatrix object
//
//	center1:
//		The center of the second filter
//
//	dir_info:
//		A string containing the directory of the filter files
//
//	NOTE:
//
//	Each coefficient file are named as "cbd_coeffs_bo_XX.surf", where
//	XX is the bo number. Within each file, the coefficient info is stored as
//
//	dim[0] dim[1] center0[0] center0[1] coeffs | dim[0] dim[1] center1[0] center1[1] coeffs | dim[0] dim[1] center0[0] center0[1] coeffs | dim[0] dim[1] center1[0] center1[1] coeffs
//	**        Decompsition filter 0         **   **       Reconstruction filter 0        **   **       Decompsition filter 0          **   **       Reconstruction filter 0        **


void NdDirectionalFilterBank::GetCheckerboardFilters(int bo, bool IsDecomposition, 
	SurfMatrix &filter0, int center0[], SurfMatrix& filter1, int center1[], string& dir_info)
{
	// bo is an even number and between 4 and 18
	assert((bo >= 4) && (bo <= 18) && (bo % 2 == 0));

	// Get the file name for the filter coefficients
	stringstream convert;
	convert << bo;
	string filter_filename = dir_info + "cbd_coeffs_bo_" + convert.str() + ".surf";

	int dims[2], center[2];

	try
	{
		ifstream coeff_file;
		coeff_file.open(filter_filename.c_str(), ios::in | ios::binary);
		
		if (!coeff_file.is_open())
			throw 1;
			
		
		if (!IsDecomposition)
		{
			coeff_file.read((char *)dims, 2 * sizeof(int));
			coeff_file.read((char *)center, 2 * sizeof(int));
			coeff_file.seekg(dims[0] * dims[1] * sizeof(double), ios_base::cur);
			
			coeff_file.read((char *)dims, 2 * sizeof(int));
			coeff_file.read((char *)center, 2 * sizeof(int));
			coeff_file.seekg(dims[0] * dims[1] * sizeof(double), ios_base::cur);
		}

		coeff_file.read((char *)dims, 2 * sizeof(int));
		coeff_file.read((char *)center0, 2 * sizeof(int));
		filter0.AllocateSpace(dims[0], dims[1]);
		coeff_file.read((char *)filter0.GetPointer(), dims[0] * dims[1] * sizeof(double));

		coeff_file.read((char *)dims, 2 * sizeof(int));
		coeff_file.read((char *)center1, 2 * sizeof(int));
		filter1.AllocateSpace(dims[0], dims[1]);
		coeff_file.read((char *)filter1.GetPointer(), dims[0] * dims[1] * sizeof(double));

		coeff_file.close();

        //// modify the center values to be zero-based
        //center0[0] -= 1;
        //center0[1] -= 1;
        //center1[0] -= 1;
        //center1[1] -= 1;
	}
	catch (std::bad_alloc)
	{
		MessageLog("NdDirectionalFilterBank", "GetCheckerboardFilters", "Out of memory!");
		throw;
	}
	catch (...)
	{
		MessageLog("NdDirectionalFilterBank", "GetCheckerboardFilters", "Filter coefficients file is corrupted!");
		throw 1;
	}
}


//////////////////////////////////////////////////////////////////////////
//	Iteratively resampled checkerboard filter bank decomposition
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nDims:
//		Number of dimensions
//
//	pDims:
//		An array specifying the dimension of the input (complex-valued) array
//
//	pData:
//		Pointer to the input array
//
//	dim_axis:
//		The dimension along which we utilize the Hermitian symmetry
//
//	levels:
//		Levels of 2-D decomposition along each pair of dimensions
//
//	filter0, filter1:
//		The two decomposition filters
//
//	center0, center1:
//		The centers of the decomposition filters

void NdDirectionalFilterBank::IrcDecomposition(int nDims, int pDims[], double *pData, 
	int dim_axis, int full_dim, int levels[], SurfMatrix &filter0, int center0[], SurfMatrix &filter1, int center1[])
{	
	assert(levels[dim_axis] == -1);

    int m, iLevel, nSubbands = 1;

	try
	{
		for (m = nDims - 1; m >= 0; m--)
		{
			if (m == dim_axis) continue;

			for (iLevel = 1; iLevel <= levels[m]; iLevel++)
			{
				// Get the 2-D filter slices
				GetIrcFilterSlices(full_dim, pDims[m], iLevel, filter0, center0, FilterSlice0);
				GetIrcFilterSlices(full_dim, pDims[m], iLevel, filter1, center1, FilterSlice1);
				
				FilterDownsampleF(nDims, pDims, pData, dim_axis, m, iLevel);
			}
		}
	}
	catch(...)
	{
		// Filter coefficient file not found.
		throw 1;
	}
}


//////////////////////////////////////////////////////////////////////////
//	Iteratively resampled checkerboard filter bank reconstruction
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nDims:
//		Number of dimensions
//
//	pDims:
//		An array specifying the dimension of the input (complex-valued) array
//