surfaceletfilterbank.cpp

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

    Pyr_Subbands[PyrLevel].GetInPlaceBackwardFFT();
	OutLowpassAray = Pyr_Subbands[PyrLevel];
	
	// Free memory space
	Pyr_Subbands[PyrLevel].Reset();

	// NDFB decomposition at each pyramid scale
	for (iLevel = 0; iLevel < PyrLevel; iLevel++)
	{
		// Prepare the output array
		OutHighpassArrays[iLevel] = new SurfArray * [nDims]; 
		
		// Get decomposition levels
		ReconstructionInfo.GetNdfbLevel(iLevel, NDFB_levels);

		// NDFB decomposition			
		NDFB_fb.GetDecomposition(Pyr_Subbands[iLevel], OutHighpassArrays[iLevel], false, NDFB_levels, ReconstructionInfo.Dec_filter0, ReconstructionInfo.Dec_center0, ReconstructionInfo.Dec_filter1, ReconstructionInfo.Dec_center1, 
			ReconstructionInfo.mSize, ReconstructionInfo.beta, ReconstructionInfo.lambda);

		// Free memory space
		Pyr_Subbands[iLevel].Reset();
	}

	delete [] NDFB_levels;
	delete [] Pyr_Subbands;

	return true;
}


//////////////////////////////////////////////////////////////////////////
//	Surfacelet filter bank reconstruction
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	InHighpassArrays
//		The input highpass directional subbands
//
//	InLowpassArray
//		The lowpass subband
//
//	OutArray
//		The reconstructed array
//
//	ReconstructionInfo
//		Filter bank parameters

void SurfaceletFilterBank::GetReconstruction(SurfArray ***InHighpassArrays, SurfArray &InLowpassArray, 
	SurfArray &OutArray, SurfaceletRecInfo &ReconstructionInfo)
{
	int nDims = InHighpassArrays[0][0][0].GetRank();
	assert(nDims >= 2);

	SurfArray *Pyr_Subbands;
	int iLevel, *NDFB_levels;
	int PyrLevel = ReconstructionInfo.PyrLevel;

	try
	{
		Pyr_Subbands = new SurfArray [PyrLevel + 1];
		NDFB_levels = new int [nDims * nDims];

	}
	catch(std::bad_alloc)
	{
		MessageLog("SurfaceletFilterBank", "GetDecomposition", "Out of memory!");
		throw;
	}

	// NDFB reconstruction at each pyramid scale
	for (iLevel = 0; iLevel < PyrLevel; iLevel++)
	{
		// Get decomposition levels
		ReconstructionInfo.GetNdfbLevel(iLevel, NDFB_levels);

		NDFB_fb.GetReconstruction(InHighpassArrays[iLevel], Pyr_Subbands[iLevel], true, NDFB_levels,
			ReconstructionInfo.Rec_filter0, ReconstructionInfo.Rec_center0, ReconstructionInfo.Rec_filter1,
			ReconstructionInfo.Rec_center1);
		
	}

	// Copy the lowpass subband
	Pyr_Subbands[PyrLevel] = InLowpassArray;

	// Multidimensional multiscale pyramid decomposition
	pyr_fb.GetReconstruction(Pyr_Subbands, OutArray, PyrLevel, false, ReconstructionInfo.pyr_mode, ReconstructionInfo.func_type);
	
	delete [] NDFB_levels;
	delete [] Pyr_Subbands;

}


//////////////////////////////////////////////////////////////////////////
//	Check the validity of the decomposition parameters
//////////////////////////////////////////////////////////////////////////
//
//  PARAMETERS:
//
//  InArray
//      The input multidimensional array.
//
//  ReconstructionInfo
//      Decomposition and reconstruction parameters for the surfacelet filter bank
//
//  Return Value
//      False if the parameters are not yet supported by the current implementation.
//
//  NOTE:
//


bool SurfaceletFilterBank::CheckParameters(SurfArray &InArray, SurfaceletRecInfo &ReconstructionInfo)
{
	int nDims = InArray.GetRank();

	if (nDims < 2)
	{
		MessageLog("SurfaceletFilterBank", "CheckParameters", "The input array must be multidimensional.");
		return false;
	}

	if (ReconstructionInfo.PyrLevel < 1)
	{
		MessageLog("SurfaceletFilterBank", "CheckParameters", "The levels of pyramid decomposition must be at least one.");
		return false;
	}

	bool IsValid = true;

	int *pDims, *pLevels;

	pDims = new int [nDims];
	pLevels = new int [nDims * nDims];

	InArray.GetDims(pDims);

	int i, j, downsample_factor, iLevel, Level, subband_nx, max_shiftingfactor, filter_ny, max_filter_ny;
	
	for (iLevel = 0; iLevel < ReconstructionInfo.PyrLevel; iLevel++)
	{
		// Get the NDFB decomposition levels
		ReconstructionInfo.GetNdfbLevel(iLevel, pLevels);

        // Check conditions for the hourglass filter bank
        for (i = 0; i < nDims; i++)
        {
            if (pDims[i] < 2 * ReconstructionInfo.mSize - 1)
            {
                MessageLog("SurfaceletFilterBank", "CheckParameters", "The dimensions of the array is smaller than 2 * mSize - 1. Try reducing mSize.");
                IsValid = false;
                break;
            }
        }

        if (!IsValid) break;

		for (i = 0; i < nDims; i++)
		{
			if (pLevels[i * nDims + i] != -1)
			{
				MessageLog("SurfaceletFilterBank", "CheckParameters", "The diagonal elements of the NDFB level matrices must be -1.");
				IsValid = false;
				break;
			}

			for (j = 0; j < nDims; j++)
			{
				if (i == j) continue;

				if ((Level = pLevels[i * nDims + j]) < 0)
				{
					MessageLog("SurfaceletFilterBank", "CheckParameters", "The off-diagonal elements of the NDFB level matrices must be nonnegative.");
					IsValid = false;
					break;
				}

				if (Level == 0) continue;

				Level = (int)pow(2.0, Level - 1);

				max_shiftingfactor = Level - 1;
				max_filter_ny = -1;

				filter_ny = ReconstructionInfo.Dec_filter0.ny + (ReconstructionInfo.Dec_filter0.nx - 1) * max_shiftingfactor;
				if (filter_ny > max_filter_ny)
					max_filter_ny = filter_ny;

				filter_ny = ReconstructionInfo.Dec_filter1.ny + (ReconstructionInfo.Dec_filter1.nx - 1) * max_shiftingfactor;
				if (filter_ny > max_filter_ny)
					max_filter_ny = filter_ny;

				filter_ny = ReconstructionInfo.Rec_filter0.ny + (ReconstructionInfo.Rec_filter0.nx - 1) * max_shiftingfactor;
				if (filter_ny > max_filter_ny)
					max_filter_ny = filter_ny;

				filter_ny = ReconstructionInfo.Rec_filter1.ny + (ReconstructionInfo.Rec_filter1.nx - 1) * max_shiftingfactor;
				if (filter_ny > max_filter_ny)
					max_filter_ny = filter_ny;

				if (pDims[i] < max_filter_ny)
				{
					MessageLog("SurfaceletFilterBank", "CheckParameters", "The resampled (sheared) filter size if larger than the downsampled array size.");
					IsValid = false;
					break;
				}


				if (pDims[j] % (Level * 2))
				{
					MessageLog("SurfaceletFilterBank", "CheckParameters", "The downsampling factor of the multidimensional pyramid must be able to divide the dimensions of the input array");
					IsValid = false;
					break;
				}
				
				subband_nx = pDims[j] / Level;

				if ((subband_nx < ReconstructionInfo.Dec_filter0.nx) || (subband_nx < ReconstructionInfo.Dec_filter1.nx)
					|| (subband_nx < ReconstructionInfo.Rec_filter0.nx) || (subband_nx < ReconstructionInfo.Rec_filter1.nx))
				{
					MessageLog("SurfaceletFilterBank", "CheckParameters", "The checkerboard filter size is larger than the downsampled subband size. Try reducing the level of decompositions or bo.");
					IsValid = false;
					break;
				}

			}

			if (!IsValid) break;
		}


		// Update the dimensions according to the downsample operations in the pyramid decomposition
		switch (ReconstructionInfo.pyr_mode)
		{
		case DOWNSAMPLE_1:
			
			downsample_factor = (iLevel == 0)? 1 : 2;
			break;

		case DOWNSAMPLE_15:
			
			downsample_factor = (iLevel == 0)? 3 : 2;
			break;

		case DOWNSAMPLE_2:

			downsample_factor = 2;
			break;

		default:
			assert(false);
		}

		for (i = 0; i < nDims; i++)
		{
			if (pDims[i] % downsample_factor)
			{
				MessageLog("SurfaceletFilterBank", "CheckParameters", "The downsampling factor of the multidimensional pyramid must be able to divide the dimensions of the input array");
				IsValid = false;
				break;
			}
			else
				pDims[i] /= downsample_factor;

			if ((iLevel == 0) && (ReconstructionInfo.pyr_mode == DOWNSAMPLE_15))
				pDims[i] *= 2;
		}

		if (!IsValid) break;

	}
	
	delete [] pLevels;
	delete [] pDims;

	return IsValid;

}


//	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.