hourglassfilterbank.cpp

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

						
		}

	}

	delete [] pHourglassFilters;
	delete [] pFilterValue;
	delete [] pIndices;

    return;

}


//////////////////////////////////////////////////////////////////////////
//	Filtering operation for the reconstruction step.
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nDims
//		Dimension of the input signal, e.g., 2-D, 3-D, etc.
// 
//	pDims
//		Actual array dimensions, e.g., 320 * 240 * 256, etc.
//
//	pDataIn
//		Pointers to the N input data.
//      Note: Input data are in the Fourier domain, and the last dimension is
//      almost cut in half. See SurfArray.h for more details.
// 
//	pDataOut
//		Pointer to the output array.
// 
//	pHourglass
//		Pointer to the hourglass filter coefficients at different 2-D planes
//
//	Explanation: (tricky!)
//
//	See DecompositionFilter for more details.

void HourglassFilterBank::ReconstructionFiltering(int nDims, int* pDims, 
        double *pDataIn[], double *pDataOut, SurfMatrix *pHourglass[])
{
	int *pIndices = NULL;
	double *pFilterValue = NULL;
	double **pHourglassFilters = NULL;

	try
	{
		pIndices = new int [nDims];
		pFilterValue = new double [nDims];
		pHourglassFilters = new double * [nDims * (nDims - 1)];
	}
	catch (std::bad_alloc)
	{
		MessageLog("HourglassFilterBank", "DecompositionFiltering", "Out of memory!");
		if (pIndices) delete [] pIndices;
		if (pFilterValue) delete [] pFilterValue;
	}
    register int m, k; // These two have the priority to get the register
    register int i; 
    register int j;
    int iSlice, nSlices, nDimLast, nDimSecondLast;

    for (k = 0; k < nDims * (nDims - 1); k++)
        pHourglassFilters[k] = pHourglass[k]->GetPointer();

    // Initialize array indices
    for (k = 0; k < nDims; k++)
        pIndices[k] = 0;

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

    // Number of complex numbers along the last dimension.
    // Note: we utilize the conjugate symmetry in the Fourier transform of real signals,
    // so the data size is nearly reduced by half.
    nDimLast = pDims[nDims - 1] / 2 + 1;
    nDimSecondLast = pDims[nDims - 2];

    double denominator, FilterValue, SumReal, SumImag;
    int idxHourglass;

    // Start the computation ...
    for (iSlice = 0; iSlice < nSlices; iSlice++)
    {
        for (j = 0; j < nDimSecondLast; j++)
        {
            pIndices[nDims - 2] = j;

            for (i = 0; i < nDimLast; i++)
            {
                pIndices[nDims - 1] = i;

                denominator = 0.0;

                idxHourglass = 0;

                // Fill in pFilterValue;
                for (k = nDims - 1; k >= 0; k--)
                {
                    FilterValue = 1.0;

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

                        FilterValue *= *(pHourglassFilters[idxHourglass++] 
                            + pIndices[k] + pIndices[m] * pDims[k]);

                    }

                    pFilterValue[nDims - 1 - k] = FilterValue;

                    denominator += FilterValue * FilterValue;
                }

                denominator = sqrt(denominator * double(nDims));

                for (k = 0; k < nDims; k++)
                    pFilterValue[k] /= denominator;

                SumReal = SumImag = 0.0;

                for (k = 0; k < nDims; k++)
                {
                    // Input and Output data are complex-valued, while the filter is real-valued.
                    SumReal += *(pDataIn[k]++) * (FilterValue = pFilterValue[k]);
                    SumImag += *(pDataIn[k]++) * (FilterValue);
                }

                *(pDataOut++) = SumReal;
                *(pDataOut++) = SumImag;

            }

        }

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

        }

    }

	delete [] pHourglassFilters;
	delete [] pFilterValue;
	delete [] pIndices;

    return;

}


//////////////////////////////////////////////////////////////////////////
//	Get the hourglass filters along all possible 2-D planes
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nDims
//		Dimension of the input signal, e.g., 2-D, 3-D, etc.
// 
//	pDims
//		Actual array dimensions, e.g., 320 * 240 * 256, etc.
//
//	pMappingMatrix
//		A matrix containing the diamond mapping kernel
//
//	pHourglassFilters
//		The resulting hourglass filters along different 2-D planes
//
//	lambda
//		See paper for more details

void HourglassFilterBank::GetHourglassFilters(int nDims, int *pDims, SurfMatrix *pMappingMatrix, 
            SurfMatrix *pHourglassFilters[], double lambda)
{
    register int k, m;
	register int i, j;
	SurfArray *pDiamond = NULL;
	SurfMatrix *pHourglass;
	int Dims2D[2];
	double *pMapping, *pDiamondFilter;
	int padding;

	assert(pMappingMatrix->nx == pMappingMatrix->ny);

	int DiamondSize = pMappingMatrix->nx;
	
	for (i = 0; i < nDims * (nDims - 1); i++)
		pHourglassFilters[i] = NULL;

	try
	{
		int idx = 0;
		for (i = nDims - 1; i >= 0; i--)
			for (j = nDims - 1; j >= 0; j--)
			{
				if (i == j) continue;
				
				pDiamond = new SurfArray;
				Dims2D[0] = pDims[j];
				Dims2D[1] = pDims[i];
				pDiamond->AllocateSpace(2, Dims2D);

				// Import values
				assert((Dims2D[0] >= DiamondSize) && (Dims2D[1] >= DiamondSize));
				assert(Dims2D[1] % 2 == 0);

				pDiamondFilter = pDiamond->GetPointer(padding);
				pMapping = pMappingMatrix->GetPointer();
				
				for (k = 0; k < DiamondSize; k++)
				{
					for (m = 0; m < DiamondSize; m++)
						*(pDiamondFilter++) = *(pMapping++);
					memset(pDiamondFilter, 0, (Dims2D[1] - DiamondSize) * sizeof(double));
					pDiamondFilter += (Dims2D[1] - DiamondSize + padding);

				}
				memset(pDiamondFilter, 0, (Dims2D[1] + padding) * (Dims2D[0] - DiamondSize) * sizeof(double));
				// Get its Fourier transform
				pDiamond->GetInPlaceForwardFFT();
				pDiamond->GetMagnitudeResponse();

				pHourglassFilters[idx++] = pHourglass = new SurfMatrix;
				pHourglass->AllocateSpace(Dims2D[0], Dims2D[1]);
				pDiamond->ExportRealValues(pHourglass->GetPointer());

				// Delete the diamond filter
				delete pDiamond;
				pDiamond = NULL;

				// Get the hourglass filter by shifting the diamond filter
				pHourglass->CircularShift(0, Dims2D[1] / 2);
                
                // Raise the matrix values to the power of lambda
                (*pHourglass) ^ lambda;
                
	
			}
	}
	catch (std::bad_alloc)
	{
		MessageLog("HourglassFilterBank", "GetHourglassFilters", "Out of memory!");
		for (i = 0; i < nDims * (nDims - 1); i++)
			if (pHourglassFilters[i]) delete pHourglassFilters[i];

		if (pDiamond) delete pDiamond;
		
		throw;
	}
	
}


//////////////////////////////////////////////////////////////////////////
//	Diamond mapping used in the hourglass filter bank
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nOrder
//		The order of the diamond mapping kernel
//
//	beta
//		For Kaiser window
//
//  pMappingMatrix
//		A matrix containing the resulting diamond mapping kernel

void HourglassFilterBank::GetDiamondMapping(int nOrder, double beta, SurfMatrix* pMappingMatrix)
{
	// nOrder must be an odd number greater than or equal to 3
	assert((nOrder >= 3) && ((nOrder / 2) * 2) != nOrder);
	
	double *pMapping = pMappingMatrix->GetPointer();

	double *pWin = new double [nOrder];
	int i, center;
	
    GetKaiserWindow(nOrder, beta, pWin);
	// Modulate the Kaiser window
	center = (nOrder - 1) / 2;
	for (i = 0; i < nOrder; i++)
	{
		pWin[i] *= Sinc((i - center) / 2.0);
	}

	int x1, x2, idx;
	double val;
	for (x1 = - (nOrder - 1); x1 <= nOrder - 1; x1++)
		for (x2 = -(nOrder - 1); x2 <= nOrder - 1; x2++)
		{
			idx = x1 + x2 + center;
			if ((idx < 0) || (idx >= nOrder))
				val = 0.0;
			else
				val = pWin[idx];

			idx = x2 - x1 + center;
			if ((idx < 0) || (idx >= nOrder))
				val = 0.0;
			else
				val *= pWin[idx];

			idx = (2 * nOrder - 1) * (x2 + nOrder - 1) + (x1 + nOrder - 1);
			pMapping[idx] = val;
		}
	
	delete [] pWin;

	return;
}


//////////////////////////////////////////////////////////////////////////
//	Get the Kaiser smoothing window
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	nOrder
//		The order of the diamond mapping
//
//	beta
//		For Kaiser window
//
//	pWin
//		Pointer to the Kaiser window

void HourglassFilterBank::GetKaiserWindow(int nOrder, double beta, double* pWin)
{
	// nOrder should be an odd integer greater than or equal to 3.
	// We will not check these conditions here.

	double denominator = CalculateI0(beta);
	int i, center;
	double val;


	center = (nOrder - 1) / 2;
	for (i = 0; i <= center; i++)
	{
		val = 2 * i / (double)(nOrder - 1) - 1;
		pWin[i] = CalculateI0(beta * sqrt(1 - val * val)) / denominator;
	}

	// Utilize the symmetry
	for (i = 1; i <= center; i++)
	{
		pWin[center + i] = pWin[center - i];
	}
}


//////////////////////////////////////////////////////////////////////////
//	Calculate the 0th-order modified Bessel function I_0
//////////////////////////////////////////////////////////////////////////

double HourglassFilterBank::CalculateI0(double x)
{
	// y = sum_{k=0}^\infty ((x/2)^k / k!)^2
	
	double y = 1.0, z = 1.0, t;
	int k = 1;
	
	x = fabs(x) / 2.0;

	do 
	{
		z *= x / (k++);
		t = z * z;
		y += t;
	} while(t > EPSILON);

	return y;
}


//////////////////////////////////////////////////////////////////////////
//	y = sinc(x) = sin(\pi x) / (\pi x)
//////////////////////////////////////////////////////////////////////////

double HourglassFilterBank::Sinc(double x)
{
	double y;
	if (fabs(x) <= EPSILON)
	{
		y = 1.0;
	}
	else
	{
		y = sin(PI * x) / (PI * x);
	}

	return y;
}


//	This software is provided "as-is", without any express or implied
//	warranty. In no event will the authors be held liable for any 
//	damages arising from the use of this software.