chenimagewavelet.cpp

LDPC码的实现,包括编码器和解码器,使用了DCT.

// ChenImageWavelet.h

#include <ChenImage.h>
#include <ChenImageWavelet.h>
#include <ChenImageMatlab.h>
#include <ippi.h>
#include <ippj.h>

int ChenImageWavelet_showWavelet(const int nWidth, const int nHeight, const short** pSubbands, const int nLevels, const char* pName, const bool bNewFigure, const bool bScaleDC, const bool bScaleAC)
{
	bool bSuccess = true;
	short* pImage = new short[nWidth * nHeight];
	int nDivisor = ROUND(pow(2.0,nLevels));
	int nSubbandWidth = nWidth/nDivisor;
	int nSubbandHeight = nHeight/nDivisor;
	
	// Stack image
	RECORD_SUCCESS( ChenImageWavelet_waveletStack(nWidth, nHeight, pSubbands, nLevels, pImage), bSuccess );

	// Copy low-low
	short* pLowLow = new short[nSubbandWidth * nSubbandHeight];
	int nSrcStep = nWidth * sizeof(short);
	int nDstStep = nSubbandWidth * sizeof(short);
	int nScaleFactor = 0;
	IppiSize roiSize = {nSubbandWidth, nSubbandHeight};
	RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pImage, nSrcStep, (Ipp16s*)pLowLow, nDstStep, roiSize), bSuccess );

	// For AC bands, optionally scale them to full display range
	if (bScaleAC) {
		// Set low-low to zeros
		RECORD_IPP_SUCCESS( ippiSet_16s_C1R(0, (Ipp16s*)pImage, nDstStep, roiSize), bSuccess );

		// Scale
		float fMean, fStdDev;
		RECORD_SUCCESS( ChenImage_imageMeanStdDev(nWidth, nHeight, pImage, &fMean, &fStdDev), bSuccess );
		short nScale = ROUND( 255.0/(3*fStdDev) );
		nSrcStep = nWidth * sizeof(short);
		roiSize.width = nWidth;
		roiSize.height = nHeight;
		nScaleFactor = 0;
		RECORD_IPP_SUCCESS( ippiMulC_16s_C1IRSfs(nScale, (Ipp16s*)pImage, nSrcStep, roiSize, nScaleFactor), bSuccess );

		// Copy back low-low
		nSrcStep = nSubbandWidth * sizeof(short);
		nDstStep = nWidth * sizeof(short);
		roiSize.width = nSubbandWidth;
		roiSize.height = nSubbandHeight;
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pLowLow, nSrcStep, (Ipp16s*)pImage, nDstStep, roiSize), bSuccess );
	}
	
	// Offset AC bands by 128
	nDstStep = nWidth * sizeof(short);
	roiSize.width = nWidth;
	roiSize.height = nHeight;
	nScaleFactor = 0;
	RECORD_IPP_SUCCESS( ippiAddC_16s_C1IRSfs(128, (Ipp16s*)pImage, nDstStep, roiSize, nScaleFactor), bSuccess );

	// For DC band, undo shift of 128
	roiSize.width = nSubbandWidth;
	roiSize.height = nSubbandHeight;
	RECORD_IPP_SUCCESS( ippiSubC_16s_C1IRSfs(128, (Ipp16s*)pImage, nDstStep, roiSize, nScaleFactor), bSuccess );

	// For DC band, optionally scale it to full display range
	if (bScaleDC) {
		float fMean, fStdDev;
		RECORD_SUCCESS( ChenImage_imageMeanStdDev(nSubbandWidth, nSubbandHeight, pLowLow, &fMean, &fStdDev), bSuccess );

		nSrcStep = nWidth * sizeof(short);
		nScaleFactor = 0;
		roiSize.width = nSubbandWidth;
		roiSize.height = nSubbandHeight;
		RECORD_IPP_SUCCESS( ippiSubC_16s_C1IRSfs(ROUND(fMean), (Ipp16s*)pImage, nSrcStep, roiSize, nScaleFactor), bSuccess ); 

		short nScale = ROUND( 255.0/(3*fStdDev) );
		RECORD_IPP_SUCCESS( ippiMulC_16s_C1IRSfs(nScale, (Ipp16s*)pImage, nSrcStep, roiSize, nScaleFactor), bSuccess );

		RECORD_IPP_SUCCESS( ippiAddC_16s_C1IRSfs(128, (Ipp16s*)pImage, nSrcStep, roiSize, nScaleFactor), bSuccess );
	}

	RECORD_SUCCESS( ChenImageMatlab_showImage(nWidth, nHeight, 1, pImage, pName, bNewFigure), bSuccess );
	delete [] pImage;
	delete [] pLowLow;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletStack(const int nWidth, const int nHeight, const short** pSubbands, const int nLevels, short* pStacked)
{
	bool bSuccess = true;

	// Fill image with zeros
	int nDstStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	ippiSet_16s_C1R(0, pStacked, nDstStep, roiSize);

	// Fill image with subbands
	int nDivisor = ROUND(pow(2.0,nLevels));
	int nSubbandWidth = nWidth/nDivisor;
	int nSubbandHeight = nHeight/nDivisor;
	for (int nLevel = 0; nLevel < nLevels; nLevel++) {
		// Select subbands
		const short* pLowLow = pSubbands[3*nLevel];
		const short* pLowHigh = pSubbands[3*nLevel+1];
		const short* pHighLow = pSubbands[3*nLevel+2];
		const short* pHighHigh = pSubbands[3*nLevel+3];

		// Copy subbands into output image
		int nSrcStep = nSubbandWidth * sizeof(short);
		int nDstStep = nWidth * sizeof(short);
		IppiSize roiSize = {nSubbandWidth, nSubbandHeight};
		if (nLevel == 0) {
			RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pLowLow, nSrcStep, (Ipp16s*)pStacked, nDstStep, roiSize), bSuccess );
		}
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pLowHigh, nSrcStep, (Ipp16s*)(pStacked + nSubbandWidth), nDstStep, roiSize), bSuccess );
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pHighLow, nSrcStep, (Ipp16s*)(pStacked + nSubbandHeight*nWidth), nDstStep, roiSize), bSuccess );
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pHighHigh, nSrcStep, (Ipp16s*)(pStacked + nSubbandHeight*nWidth + nSubbandWidth), nDstStep, roiSize), bSuccess );

		// Prepare for next level
		// ChenImage_showImage(nWidth, nHeight, 1, pStacked, "stacked");
		nSubbandWidth *= 2;
		nSubbandHeight *= 2;
	} // end nLevel

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletDestack(const int nWidth, const int nHeight, const short* pStacked, const int nLevels, short** pSubbands)
{
	bool bSuccess = true;

	int nDivisor = ROUND(pow(2.0,nLevels));
	int nSubbandWidth = nWidth/nDivisor;
	int nSubbandHeight = nHeight/nDivisor;
	for (int nLevel = 0; nLevel < nLevels; nLevel++) {
		// Select subbands
		short* pLowLow = pSubbands[3*nLevel];
		short* pLowHigh = pSubbands[3*nLevel+1];
		short* pHighLow = pSubbands[3*nLevel+2];
		short* pHighHigh = pSubbands[3*nLevel+3];

		// Copy subbands into output image
		int nSrcStep = nWidth * sizeof(short);
		int nDstStep = nSubbandWidth * sizeof(short);
		IppiSize roiSize = {nSubbandWidth, nSubbandHeight};
		if (nLevel == 0) {
			RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pStacked, nSrcStep, (Ipp16s*)pLowLow, nDstStep, roiSize), bSuccess );
		}
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pStacked + nSubbandWidth), nSrcStep, (Ipp16s*)pLowHigh, nDstStep, roiSize), bSuccess );
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pStacked + nSubbandHeight*nWidth), nSrcStep, (Ipp16s*)pHighLow, nDstStep, roiSize), bSuccess );
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pStacked + nSubbandHeight*nWidth + nSubbandWidth), nSrcStep, (Ipp16s*)pHighHigh, nDstStep, roiSize), bSuccess );

		// Prepare for next level
		// ChenImage_showImage(nWidth, nHeight, 1, pStacked, "stacked");
		nSubbandWidth *= 2;
		nSubbandHeight *= 2;
	} // end nLevel
	
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletStackPositions(const int nWidth, const int nHeight, const int nLevels, int* pStackedToSubband)
{
	bool bSuccess = true;

	int* pInterleavedToStacked = new int[nWidth * nHeight];
	int* pStackedToInterleaved = new int[nWidth * nHeight];
	int* pInterleavedToSubband = new int[nWidth * nHeight];
	RECORD_SUCCESS( ChenImageWavelet_waveletMultiLevelInterleavePositions(nWidth, nHeight, nLevels, pInterleavedToStacked, pStackedToInterleaved, pInterleavedToSubband), bSuccess );

	for (int nPix = 0; nPix < nWidth * nHeight; nPix++) {
		pStackedToSubband[nPix] = pInterleavedToSubband[pStackedToInterleaved[nPix]];
	}

	delete [] pInterleavedToStacked;
	delete [] pStackedToInterleaved;
	delete [] pInterleavedToSubband;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletMultiLevelAlloc(const int nWidth, const int nHeight, const int nLevels, int* pWidths, int* pHeights)
{
	bool bSuccess = true;

	int nDivisor = ROUND(pow(2.0,nLevels));
	for (int nLevel = 0; nLevel < nLevels; nLevel++) {
		if (nLevel == 0) {
			pWidths[0] = nWidth/nDivisor;
			pHeights[0] = nHeight/nDivisor;
		}
		for (int nBand = 3*nLevel+1; nBand <= 3*nLevel+3; nBand++) {
			pWidths[nBand] = nWidth/nDivisor;
			pHeights[nBand] = nHeight/nDivisor;
		} // end nBand
		nDivisor /= 2;
	} // end nLevel

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletMultiLevelForward(const int nWidth, const int nHeight, const short* pSrc, const int nLevels, const char* pFilter, short** pDst)
{
	bool bSuccess = true;

	// Perform multi-level decomposition
	int nSrcWidth = nWidth, nSrcHeight = nHeight;
	int nDstWidth = nWidth/2, nDstHeight = nHeight/2;
	const short* pCurrSrc = pSrc;
	short* pCurrDst[4];
	for (int nLevel = nLevels; nLevel > 0; nLevel--) {
		// Allocate space
		int nDstPixels = nDstWidth * nDstHeight;
		for (int nBand = 0; nBand < 4; nBand++) {
			pCurrDst[nBand] = new short[nDstPixels];
		}

		// Perform single-level decomposition
		if (strcmp("9/7", pFilter) == 0) {
			RECORD_SUCCESS( ChenImageWavelet_wavelet97Forward(nSrcWidth, nSrcHeight, pCurrSrc, pCurrDst), bSuccess );
		}
		else {
			RECORD_SUCCESS( ChenImageWavelet_wavelet53Forward(nSrcWidth, nSrcHeight, pCurrSrc, pCurrDst), bSuccess );
		}

		// Copy subbands into real output
		short* pLowLow = pDst[3*(nLevel-1)];
		short* pLowHigh = pDst[3*(nLevel-1)+1];
		short* pHighLow = pDst[3*(nLevel-1)+2];
		short* pHighHigh = pDst[3*(nLevel-1)+3];
		int nSrcStep = nDstWidth * sizeof(short);
		int nDstStep = nDstWidth * sizeof(short);
		IppiSize roiSize = {nDstWidth, nDstHeight};
		if (nLevel == 1) {
			RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pCurrDst[0], nSrcStep, (Ipp16s*)pLowLow, nDstStep, roiSize), bSuccess );
		}
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pCurrDst[1], nSrcStep, (Ipp16s*)pLowHigh, nDstStep, roiSize), bSuccess );
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pCurrDst[2], nSrcStep, (Ipp16s*)pHighLow, nDstStep, roiSize), bSuccess );
		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pCurrDst[3], nSrcStep, (Ipp16s*)pHighHigh, nDstStep, roiSize), bSuccess );
		// ChenImage_showWavelet(nSrcWidth, nSrcHeight, (const short**)pCurrDst);
		
		// Prepare for next level
		pCurrSrc = pCurrDst[0];
		nSrcWidth /= 2;
		nSrcHeight /= 2;
		nDstWidth /= 2;
		nDstHeight /= 2;
		
		// Clean up space, but don't delete low-low band
		for (int nBand = 1; nBand < 4; nBand++) {
			delete [] pCurrDst[nBand];
		}
	} // end nLevel
	delete [] pCurrDst[0];

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_wavelet97MultiLevelForward(const int nWidth, const int nHeight, const short* pSrc, const int nLevels, short** pDst)
{
	bool bSuccess = true;

	RECORD_SUCCESS( ChenImageWavelet_waveletMultiLevelForward(nWidth, nHeight, pSrc, nLevels, "9/7", pDst), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletMultiLevelInverse(const int nWidth, const int nHeight, const short** pSrc, const int nLevels, const char* pFilter, short* pDst)
{
	bool bSuccess = true;

	// Perform multi-level synthesis
	int nDivisor = ROUND(pow(2.0,nLevels));
	int nSrcWidth = nWidth/nDivisor, nSrcHeight = nHeight/nDivisor;
	int nDstWidth = nSrcWidth*2, nDstHeight = nSrcHeight*2;
	short* pCurrDst = (short*)pSrc[0];
	const short* pCurrSrc[4];
	for (int nLevel = 0; nLevel < nLevels; nLevel++) {
		// Allocate space
		int nDstPixels = nDstWidth * nDstHeight;
		pCurrSrc[0] = pCurrDst;
		for (int nBand = 1; nBand < 4; nBand++) {
			pCurrSrc[nBand] = pSrc[3*nLevel+nBand];
		}
		pCurrDst = new short[nDstPixels];
		// ChenImage_showWavelet(nDstWidth, nDstHeight, pCurrSrc);

		// Perform single-level synthesis
		if (strcmp("9/7", pFilter) == 0) {
			RECORD_SUCCESS( ChenImageWavelet_wavelet97Inverse(nDstWidth, nDstHeight, (const short**)pCurrSrc, pCurrDst), bSuccess );
		}
		else {
			RECORD_SUCCESS( ChenImageWavelet_wavelet53Inverse(nDstWidth, nDstHeight, (const short**)pCurrSrc, pCurrDst), bSuccess );
		}
		// ChenImage_showImage(nDstWidth, nDstHeight, 1, pCurrDst);

		// Clean up space
		if (nLevel > 0) {
			delete [] pCurrSrc[0];
		}

		// Prepare for next level
		nSrcWidth *= 2;
		nSrcHeight *= 2;
		nDstWidth *= 2;
		nDstHeight *= 2;
	} // end nLevel

	// Copy final coarse image into real output
	int nSrcStep = nWidth * sizeof(short);
	int nDstStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	ippiCopy_16s_C1R((const Ipp16s*)pCurrDst, nSrcStep, (Ipp16s*)pDst, nDstStep, roiSize);
	delete [] pCurrDst;

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}


int ChenImageWavelet_wavelet97MultiLevelInverse(const int nWidth, const int nHeight, const short** pSrc, const int nLevels, short* pDst)
{
	bool bSuccess = true;

	RECORD_SUCCESS( ChenImageWavelet_waveletMultiLevelInverse(nWidth, nHeight, pSrc, nLevels, "9/7", pDst), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletForward(const int nWidth, const int nHeight, const short* pSrc, const char* pFilter, short* pDst[4])
{
	bool bSuccess = true;

	// Allocate space
	IppiRect waveletTileRect = {0, 0, nWidth, nHeight};
	int nBufferSize;
	if (strcmp("9/7", pFilter) == 0) {
		RECORD_IPP_SUCCESS( ippiWTGetBufSize_D97_JPEG2K_16s_C1R(&waveletTileRect, &nBufferSize), bSuccess );
	}
	else {
		RECORD_IPP_SUCCESS( ippiWTGetBufSize_B53_JPEG2K_16s_C1R(&waveletTileRect, &nBufferSize), bSuccess );
	}
	int nSrcStep = nWidth * sizeof(short);
	int pDstStep[4];
	for (int nBand = 0; nBand < 4; nBand++) {
		pDstStep[nBand] = waveletTileRect.width/2 * sizeof(short);
	}
	unsigned char* pBuffer = new unsigned char[nBufferSize];

	// Perform actual DWT
	if (strcmp("9/7", pFilter) == 0) {
		RECORD_IPP_SUCCESS( ippiWTFwd_D97_JPEG2K_16s_C1R(pSrc, nSrcStep, &waveletTileRect, pDst, pDstStep, pBuffer), bSuccess );
	}
	else {
		RECORD_IPP_SUCCESS( ippiWTFwd_B53_JPEG2K_16s_C1R(pSrc, nSrcStep, &waveletTileRect, pDst, pDstStep, pBuffer), bSuccess );
	}
	delete [] pBuffer;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_wavelet97Forward(const int nWidth, const int nHeight, const short* pSrc, short* pDst[4])
{
	bool bSuccess = true;

	RECORD_SUCCESS( ChenImageWavelet_waveletForward(nWidth, nHeight, pSrc, "9/7", pDst), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletInverse(const int nWidth, const int nHeight, const short** pSrc, const char* pFilter, short* pDst)
{
	bool bSuccess = true;

	// Allocate space
	IppiRect waveletTileRect = {0, 0, nWidth, nHeight};
	int nBufferSize;
	if (strcmp("9/7", pFilter) == 0) {
		RECORD_IPP_SUCCESS( ippiWTGetBufSize_D97_JPEG2K_16s_C1R(&waveletTileRect, &nBufferSize), bSuccess );
	}
	else {
		RECORD_IPP_SUCCESS( ippiWTGetBufSize_B53_JPEG2K_16s_C1R(&waveletTileRect, &nBufferSize), bSuccess );
	}
	int nDstStep = nWidth * sizeof(short);
	int pSrcStep[4];
	for (int nBand = 0; nBand < 4; nBand++) {
		pSrcStep[nBand] = waveletTileRect.width/2 * sizeof(short);
	}
	unsigned char* pBuffer = new unsigned char[nBufferSize];

	// Perform actual IDWT
	if (strcmp("9/7", pFilter) == 0) {
		RECORD_IPP_SUCCESS( ippiWTInv_D97_JPEG2K_16s_C1R(pSrc, pSrcStep, pDst, nDstStep, &waveletTileRect, pBuffer), bSuccess );
	}
	else {
		RECORD_IPP_SUCCESS( ippiWTInv_B53_JPEG2K_16s_C1R(pSrc, pSrcStep, pDst, nDstStep, &waveletTileRect, pBuffer), bSuccess );
	}
	delete [] pBuffer;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_wavelet97Inverse(const int nWidth, const int nHeight, const short** pSrc, short* pDst)