nddirectionalfilterbank.cpp

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

//	pData:
//		Pointer to the input array
//
//	dim_axis:
//		The dimension along which we utilize the Hermitian symmetry
//
//	levels:
//		Levels of 2-D reconstruction along each pair of dimensions
//
//	filter0, filter1:
//		The two reconstruction filters
//
//	center0, center1:
//		The centers of the reconstruction filters

void NdDirectionalFilterBank::IrcReconstruction(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;

	for (m = 0; m < nDims; m++)
	{
		if (m == dim_axis) continue;

		for (iLevel = levels[m]; iLevel >= 1; 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);

			FilterUpsampleF(nDims, pDims, pData, dim_axis, m, iLevel);

		}
	}
}


//////////////////////////////////////////////////////////////////////////
//	Obtain the 2-D slice containing the frequency values of the irc filters
//////////////////////////////////////////////////////////////////////////
//	
//	PARAMETERS:
//	
//	dim_K
//		The major dimension
//	
//	dim_M
//		The minor dimension
//
//	nLevel
//		The pyramid decomposition level
//
//	filter
//		A 2-D matrix containing the filter coefficients
//
//	center
//		The center location of the 2-D filter
//
//	FilterSlice
//		The resulting 2-D slice
//
//	Note:
//	
//	FilterSlice will be a matrix of dimensions dim_M * (2 * dim_K);

void NdDirectionalFilterBank::GetIrcFilterSlices(int dim_K, int dim_M, int nLevel, SurfMatrix &filter, 
	int center[], SurfMatrix &FilterSlice)
{
	
	int nChannels = (int)pow(2.0,  (double)(nLevel - 1));
	
	assert((nLevel >= 1) && (dim_M % nChannels == 0));

	int nBlock = dim_M / nChannels;
	
	try
	{
		// complex-valued array of dimensions dim_M * dim_K
		if (FilterSlice.nx * FilterSlice.ny == dim_M * 2 * dim_K)
        {
            FilterSlice.nx = dim_M;
            FilterSlice.ny = 2 * dim_K;
        }
        else
            FilterSlice.AllocateSpace(dim_M, 2 * dim_K);

		FilterSlice.ZeroFill();

		int new_center[2];
		int shifting_factor;
		int ChannelIndex;
		

		SurfMatrix SmallFilterSlice;
		SmallFilterSlice.AllocateSpace(dim_M / nChannels, 2 * dim_K);
		fftw_plan planFFT;

		for (ChannelIndex = 0; ChannelIndex < nChannels; ChannelIndex++)
		{
			shifting_factor = nChannels - 1 - 2 * ChannelIndex;

			new_center[0] = center[0];
			if (shifting_factor >= 0)
				new_center[1] = center[1] + shifting_factor * center[0];
			else
				new_center[1] = dim_K - filter.ny + center[1] + shifting_factor * center[0];

			SmallFilterSlice.FillSubMatrixF(filter, shifting_factor >= 0);

            // make sure the resampled (sheared) filter can still fit in
            assert(filter.ny + (filter.nx - 1) * shifting_factor <= dim_K);

			// Note: the matrix is complex valued, so we need to double the shifting factor
			SmallFilterSlice.ResampleRow(2 * shifting_factor);

			SmallFilterSlice.CircularShift(- new_center[0], -2 * new_center[1]);

			// Take the FFT of the filter slice
			// IMPORTANT: subject to change in scaling factors
			
			planFFT = fftw_plan_dft_2d(SmallFilterSlice.nx, SmallFilterSlice.ny / 2, (fftw_complex*)SmallFilterSlice.GetPointer(),
				(fftw_complex*)SmallFilterSlice.GetPointer(), FFTW_FORWARD, FFTW_ESTIMATE);

			fftw_execute(planFFT);

			fftw_destroy_plan(planFFT);

    		FilterSlice.FillSubMatrix(SmallFilterSlice, ChannelIndex * nBlock, 0);
		}
	}
	catch (std::bad_alloc)
	{
		MessageLog("NdDirectionalFilterBank", "GetIrcFilterSlices", "Out of memory!");
		throw;
	}
}


//////////////////////////////////////////////////////////////////////////
//  Extract an N-D array into a list of sub-arrays
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  nDims:
//      Rank of the input array
//
//  pDims:
//      Dimension of the input array
//
//  pInData:
//      Input array
//
//  pSubbands:
//      Output arrays
//
//  levels
//      Specifying how each dimension is divided.
//
//  dim_axis
//      The axis on which the input array utilizes the Hermitian symmetry
//
//	full_dim
//		The original dimension along dim_axis

void NdDirectionalFilterBank::BoxToCell(int nDims, int pDims[], double *pInData, SurfArray pSubbands[], int levels[], int dim_axis, int full_dim)
{
    // We only work on multidimensional arrays
    assert((nDims >= 2) && (dim_axis >= 0) && (dim_axis < nDims) && (levels[dim_axis] == -1));
    
    register int n, m, j, i;
    
    // level = x ==> 2 ^ x subarrays
    levels[dim_axis] = 0;
    for (n = 0; n < nDims; n++)
    {
        levels[n] = (int)pow(2.0,  (double)levels[n]);
        // 2 ^ level must divide the corresponding dimension
        if ((pDims[n] % levels[n]) != 0)
            assert(false);
    }

    // Total number of output subbands
    int nSubbands = 1;
    for (n = 0; n < nDims; n++)
    {
        nSubbands *= levels[n];
    }

    double **pOutData, **pOutDataFixed;
    int *pSkip, *pNewDims, *pIndices;
    try
    {
        // pointers to each output subband array
        pOutData = new double * [nSubbands];
        pOutDataFixed = new double * [nSubbands];
        pSkip = new int [nDims];

        // dimension of the smaller output arrays
        pNewDims = new int [nDims];
        for (n = 0; n < nDims; n++)
        {
            pNewDims[n] = pDims[n] / levels[n];
        }   

        pIndices = new int [nDims];

        int nSubbandPoints = 2;
        for (n = 0; n < nDims; n++)
            nSubbandPoints *= pNewDims[n];

        // Allocate subband arrays
        for (n = 0; n < nSubbands; n++)
        {
            pOutData[n] = pOutDataFixed[n] = new double [nSubbandPoints];
        }
    }
    catch (std::bad_alloc)
    {
    	MessageLog("NdDirectionalFilterBank", "BoxToCell", "Out of memory!");
        throw;
    }

    // calculate skip distance    
    pSkip[nDims - 1] = 1;
    for (n = nDims - 2; n >= 0; n--)
        pSkip[n] = pSkip[n + 1] * levels[n + 1];

    // number of 2-D slices
    int nSlices = 1;
    for (n = 0; n <= nDims - 3; n++)
        nSlices *= pDims[n];

    int nDimLast = levels[nDims - 1];
    int nDimSecondLast0 = levels[nDims - 2];
    int nDimSecondLast1 = pNewDims[nDims - 2];
    
    // Initialize pIndices;
    memset(pIndices, 0, nDims * sizeof(int));
    double **ptrStartPosition;
    int nCopy = 2 * pNewDims[nDims - 1];

	double *pInFixed = pInData;
    
    for (i = 0; i < nSlices; i++)
    {
        ptrStartPosition = pOutData;
        for (j = 0; j <= nDims - 3; j++)
            ptrStartPosition += (pIndices[j] / pNewDims[j]) * pSkip[j];
        
        for (j = 0; j < nDimSecondLast0; j++)
        {
            for (m = 0; m < nDimSecondLast1; m++)
            {
                for (n = 0; n < nDimLast; n++)
                {
                    memcpy(ptrStartPosition[n], pInData, nCopy * sizeof(double));
                    ptrStartPosition[n] += nCopy;
                    pInData += nCopy;
                }
            }
            ptrStartPosition += pSkip[nDims - 2];
        }

        // Update pIndices
        for (j = nDims - 3; j >= 0; j--)
        {
            if ( (++pIndices[j]) < pDims[j])
                break;
            else
                pIndices[j] = 0;
        }
    }

	pNewDims[dim_axis] = full_dim;

    // Create SurfArray-typed subbands
    for (n = 0; n < nSubbands; n++)
    {
        pSubbands[n].AllocateSpace(nDims, pNewDims);
		pSubbands[n].SetIsRealValued(false);
        pSubbands[n].RestoreHermianSymmetryAxis(pOutDataFixed[n], dim_axis);
        // clear some memory
        delete [] pOutDataFixed[n];
    }

    // restore levels
    for (n = 0; n < nDims; n++)
        levels[n] = (int)(log((double)levels[n]) / log(2.0));
    levels[dim_axis] = -1;
    
    // release temporary memory space
    delete [] pIndices;
    delete [] pNewDims;
    delete [] pOutData;
    delete [] pOutDataFixed;
    delete [] pSkip;
}


//////////////////////////////////////////////////////////////////////////
//  Combine a list of small arrays into a larger one
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  pSubbands
//      A list of smaller array of the same dimensions
//
//  dim_axis
//      The new dimension along which the big array will utilize the Hermitian symmetry
//
//  levels
//      Specifying how each dimension is divided
//
//  RETUEN
//      A pointer to the allocated memory space for the big array

double *NdDirectionalFilterBank::CellToBox(SurfArray pSubbands[], int dim_axis, int levels[])
{
    int nDims = pSubbands[0].GetRank();
    assert((nDims >= 2) && (dim_axis >= 0) && (dim_axis < nDims) && (levels[dim_axis] == -1));

    register int n, m, j, i;

    // level = x ==> 2 ^ x sub-arrays
    levels[dim_axis] = 0;
    for (n = 0; n < nDims; n++)
    {
        levels[n] = (int)pow(2.0,  (double)levels[n]);
    }

    // Total number of output subbands
    int nSubbands = 1;
    for (n = 0; n < nDims; n++)
    {
        nSubbands *= levels[n];
    }

    double **pSubbandData, **pSubbandDataFixed, *pDataBox, *pData;
    int *pSkip, *pDims, *pLargeDims, *pIndices;
    try
    {
        // pointers to each output subband array
        pSubbandData = new double * [nSubbands];
        pSubbandDataFixed = new double * [nSubbands];
        pSkip = new int [nDims];
        pDims = new int [nDims];
        pLargeDims = new int [nDims];

        // dimension of the smaller output arrays
        pSubbands[0].GetDims(pDims);
        pDims[dim_axis] = pDims[dim_axis] / 2 + 1;

        int nTotalPoints = 2;
        for (n = 0; n < nDims; n++)
        {
            pLargeDims[n] = pDims[n] * levels[n];
            nTotalPoints *= pLargeDims[n];
        }

        pIndices = new int [nDims];

        // Allocate subband arrays
        for (n = 0; n < nSubbands; n++)
        {
            // NOTE: memory space allocated here
            pSubbandData[n] = pSubbandDataFixed[n] = pSubbands[n].NewHermitianSymmetryAxis(dim_axis);
        }

        pData = pDataBox = new double [nTotalPoints];
    }
    catch (std::bad_alloc)
    {
        MessageLog("NdDirectionalFilterBank", "BoxToCell", "Out of memory!");
        throw;
    }

    // calculate skip distance    
    pSkip[nDims - 1] = 1;
    for (n = nDims - 2; n >= 0; n--)
        pSkip[n] = pSkip[n + 1] * levels[n + 1];

    // number of 2-D slices
    int nSlices = 1;
    for (n = 0; n <= nDims - 3; n++)
        nSlices *= pLargeDims[n];

    int nDimLast = levels[nDims - 1];
    int nDimSecondLast0 = levels[nDims - 2];
    int nDimSecondLast1 = pDims[nDims - 2];

    // Initialize pIndices;
    memset(pIndices, 0, nDims * sizeof(int));
    double **ptrStartPosition;
    int nCopy = 2 * pDims[nDims - 1];

    for (i = 0; i < nSlices; i++)
    {
        ptrStartPosition = pSubbandData;
        for (j = 0; j <= nDims - 3; j++)
            ptrStartPosition += (pIndices[j] / pDims[j]) * pSkip[j];

        for (j = 0; j < nDimSecondLast0; j++)
        {
            for (m = 0; m < nDimSecondLast1; m++)
            {
                for (n = 0; n < nDimLast; n++)
                {
                    memcpy(pData, ptrStartPosition[n], nCopy * sizeof(double));
                    ptrStartPosition[n] += nCopy;
                    pData += nCopy;
                }
            }
            ptrStartPosition += pSkip[nDims - 2];
        }

        // Update pIndices
        for (j = nDims - 3; j >= 0; j--)
        {
            if ( (++pIndices[j]) < pLargeDims[j])
                break;
            else
                pIndices[j] = 0;
        }
    }

    // restore levels
    for (n = 0; n < nDims; n++)
        levels[n] = (int)(log((double)levels[n]) / log(2.0));
    levels[dim_axis] = -1;

    // release temporary memory space
    for (n = 0; n < nSubbands; n++)
        delete [] pSubbandDataFixed[n];

    delete [] pIndices;
    delete [] pLargeDims;
    delete [] pDims;
    delete [] pSubbandData;
    delete [] pSubbandDataFixed;
    delete [] pSkip;

    return pDataBox;
}


//////////////////////////////////////////////////////////////////////////
//	Filtering followed by downsampling
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nDims
//		Number of dimensions of the input data
//
//	pDims
//		Dimensions of the input data
//
//	pData
//		pointer to the input data
//
//	K
//		The major dimension
//
//	M
//		The minor dimension
//
//	Level
//		The level of multiscale pyramid