chenimagewavelet.cpp

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

{
	bool bSuccess = true;

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

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_wavelet53MultiLevelForward(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, "5/3", pDst), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_wavelet53MultiLevelInverse(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, "5/3", pDst), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

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

	RECORD_SUCCESS( ChenImageWavelet_waveletForward(nWidth, nHeight, pSrc, "5/3", pDst), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

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

	RECORD_SUCCESS( ChenImageWavelet_waveletInverse(nWidth, nHeight, pSrc, "5/3", pDst), bSuccess );
	
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletRateAllocate(const int nWidth, const int nHeight, const short** pSrc, const int nLevels, const float fRate, const float fConditionalReduction, float pStepSizes[4], float pRates[4])
{
	bool bSuccess = true;

	// Get band dimensions
	int nBands = 3*nLevels + 1;
	int* pWidths = new int[nBands];
	int* pHeights = new int[nBands];
	RECORD_SUCCESS( ChenImageWavelet_waveletMultiLevelAlloc(nWidth, nHeight, nLevels, pWidths, pHeights), bSuccess );

	// Calculate statistics for each subband
	float* pFilterGain = new float[nBands];
	float* pEpsilonSqr = new float[nBands];
	float* pSigmaSqr = new float[nBands];
	float* pEta = new float[nBands];
	float* pCondReduc = new float[nBands];
	float* pWeightedVar = new float[nBands];
	short* pMin = new short[nBands];
	short* pMax = new short[nBands];
	float fDenomProd = 1.0;
	for (int nBand = 0; nBand < nBands; nBand++) {
		// Filter gain
		int nLevel = MAX(0,(nBand-1)/3);
		pFilterGain[nBand] = (float)1.0/(float)pow(nLevel+1.0,2.0);

		// Standard deviation
		float fMean, fStdDev;
		RECORD_SUCCESS( ChenImage_imageMeanStdDev(pWidths[nBand], pHeights[nBand], pSrc[nBand], &fMean, &fStdDev), bSuccess );
		pSigmaSqr[nBand] = fStdDev * fStdDev;

		// Distortion epsilon
		pEpsilonSqr[nBand] = 1.0;
		RECORD_SUCCESS( ChenImage_imageDistortionEpsilonSqr(pWidths[nBand], pHeights[nBand], pSrc[nBand], pEpsilonSqr+nBand), bSuccess );

		// Fraction of total coefficients
		pEta[nBand] = (float)(pWidths[nBand]*pHeights[nBand]) / (float)(nWidth*nHeight);

		// Conditional-entropy-induced rate reduction
		pCondReduc[nBand] = fConditionalReduction;

		// Weighted variance
		pWeightedVar[nBand] = pFilterGain[nBand] * pEpsilonSqr[nBand] * pSigmaSqr[nBand] * pCondReduc[nBand];
		fDenomProd *= pow(pWeightedVar[nBand], pEta[nBand]);

		// Min and max
		int nSrcStep = pWidths[nBand] * sizeof(short);
		IppiSize roiSize = {pWidths[nBand], pHeights[nBand]};
		RECORD_IPP_SUCCESS( ippiMinMax_16s_C1R((const Ipp16s*)pSrc[nBand], nSrcStep, roiSize, pMin+nBand, pMax+nBand), bSuccess );
	}

	// Calculate subband rates (Tonomura)
	bool* pPositiveRate = new bool[nBands];
	for (int nBand = 0; nBand < nBands; nBand++) {
		pRates[nBand] = (float) MAX(0, fRate + 0.5 * log(pWeightedVar[nBand] / fDenomProd) / log(2.0));
		pPositiveRate[nBand] = pRates[nBand] > 0.0;
	}

	// Constrain rates to be nonnegative (Bradley)
	int nMaxIterations = 100;
	for (int nIteration = 1; nIteration <= nMaxIterations; nIteration++) {
		float fS = 0.0;
		float fDenomProd = 1.0;
		for (int nBand = 0; nBand < nBands; nBand++) {
			fS += (float) (pPositiveRate[nBand] ? pEta[nBand] : 0.0);
			fDenomProd *= (float) (pPositiveRate[nBand] ? pow(pWeightedVar[nBand], pEta[nBand]) : 1.0);
		}
		fDenomProd = pow(fDenomProd, (float)1.0/fS);
		bool bAllNonnegative = true;
		for (int nBand = 0; nBand < nBands; nBand++) {
			pRates[nBand] = (float) (fRate/fS + 0.5 * log(pWeightedVar[nBand] / fDenomProd) / log(2.0));
			pPositiveRate[nBand] = pRates[nBand] > 0.0;
			bAllNonnegative = bAllNonnegative && pRates[nBand] >= 0.0;
			pRates[nBand] = MAX(0, pRates[nBand]);
		}
		if (bAllNonnegative) break;
	}

	// Choose uniform scalar quantizer step sizes
	for (int nBand = 0; nBand < nBands; nBand++) {
		float fBins = (float) MAX(1.0, floor( pow((float)2.0, pRates[nBand]) ) );
		float fGamma = 2.5;
		float fRange = MIN( 2*fGamma*sqrt(pSigmaSqr[nBand]), pMax[nBand]-pMin[nBand] );
		pStepSizes[nBand] = (float) MAX(1.0, fRange / fBins);
	}
	
	delete [] pWidths;
	delete [] pHeights;
	delete [] pFilterGain;
	delete [] pEpsilonSqr;
	delete [] pSigmaSqr;
	delete [] pEta;
	delete [] pCondReduc;
	delete [] pWeightedVar;
	delete [] pMin;
	delete [] pMax;
	delete [] pPositiveRate;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

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

	// Get band dimensions
	int nBands = 3*nLevels + 1;
	int* pWidths = new int[nBands];
	int* pHeights = new int[nBands];
	RECORD_SUCCESS( ChenImageWavelet_waveletMultiLevelAlloc(nWidth, nHeight, nLevels, pWidths, pHeights), bSuccess );

	// Quantize subbands
	for (int nBand = 0; nBand < nBands; nBand++) {
		// Divide unquantized subband by step size, in floating point
		float* pSubband = new float[pWidths[nBand] * pHeights[nBand]];
		int nSrcStep = pWidths[nBand] * sizeof(short);
		int nDstStep = pWidths[nBand] * sizeof(float);
		IppiSize roiSize = {pWidths[nBand], pHeights[nBand]};
		RECORD_IPP_SUCCESS( ippiConvert_16s32f_C1R((const Ipp16s*)pSrc[nBand], nSrcStep, (Ipp32f*)pSubband, nDstStep, roiSize), bSuccess );
		int nSrcDstStep = pWidths[nBand] * sizeof(float);
		RECORD_IPP_SUCCESS( ippiDivC_32f_C1IR(pStepSizes[nBand], (Ipp32f*)pSubband, nSrcDstStep, roiSize), bSuccess );

		// Convert to short output
		nSrcStep = pWidths[nBand] * sizeof(float);
		nDstStep = pWidths[nBand] * sizeof(short);
		IppRoundMode roundMode = ippRndNear;
		RECORD_IPP_SUCCESS( ippiConvert_32f16s_C1R((const Ipp32f*)pSubband, nSrcStep, (Ipp16s*)pDst[nBand], nDstStep, roiSize, roundMode), bSuccess );
		delete [] pSubband;
	}

	delete [] pWidths;
	delete [] pHeights;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

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

	// Get band dimensions
	int nBands = 3*nLevels + 1;
	int* pWidths = new int[nBands];
	int* pHeights = new int[nBands];
	RECORD_SUCCESS( ChenImageWavelet_waveletMultiLevelAlloc(nWidth, nHeight, nLevels, pWidths, pHeights), bSuccess );

	for (int nBand = 0; nBand < nBands; nBand++) {
		// Multiply by step size, in floating point
		float* pSubband = new float[pWidths[nBand] * pHeights[nBand]];
		int nSrcStep = pWidths[nBand] * sizeof(short);
		int nDstStep = pWidths[nBand] * sizeof(float);
		IppiSize roiSize = {pWidths[nBand], pHeights[nBand]};
		RECORD_IPP_SUCCESS( ippiConvert_16s32f_C1R((const Ipp16s*)pSrc[nBand], nSrcStep, (Ipp32f*)pSubband, nDstStep, roiSize), bSuccess );
		int nSrcDstStep = pWidths[nBand] * sizeof(float);
		RECORD_IPP_SUCCESS( ippiMulC_32f_C1IR(pStepSizes[nBand], (Ipp32f*)pSubband, nSrcDstStep, roiSize), bSuccess );

		// Convert to short output
		nSrcStep = pWidths[nBand] * sizeof(float);
		nDstStep = pWidths[nBand] * sizeof(short);
		IppRoundMode roundMode = ippRndNear;
		RECORD_IPP_SUCCESS( ippiConvert_32f16s_C1R((const Ipp32f*)pSubband, nSrcStep, (Ipp16s*)pDst[nBand], nDstStep, roiSize, roundMode), bSuccess );
		delete [] pSubband;
	}

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

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

	for (int nRow = 0; nRow < nHeight; nRow++) {
		for (int nCol = 0; nCol < nWidth; nCol++) {
			int nPix = nRow*nWidth + nCol;
			int nBand, nCoeff;

			// Subband 0
			if ((nRow % 2) == 0 && (nCol % 2) == 0) {
				nBand = 0;
			}
			// Subband 1
			else if ((nRow % 2) == 0 && (nCol % 2) == 1) {
				nBand = 1;
			}
			// Subband 2
			else if ((nRow % 2) == 1 && (nCol % 2) == 0) {
				nBand = 2;
			}
			// Subband 3
			else {
				nBand = 3;
			}
			nCoeff = (nRow/2)*nWidth/2 + nCol/2;

			pDst[pDstMap[nPix]] = pSrc[nBand][nCoeff];

		} // end nCol
	} // end nRow

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

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

	for (int nRow = 0; nRow < nHeight/2; nRow++) {
		for (int nCol = 0; nCol < nWidth/2; nCol++) {
			int nCoeff = nRow*nWidth/2 + nCol;

			// Subband 0
			int nPix = 2*nRow*nWidth + 2*nCol;
			pDst[0][nCoeff] = pSrc[pSrcMap[nPix]];

			// Subband 1
			nPix = 2*nRow*nWidth + 2*nCol + 1;
			pDst[1][nCoeff] = pSrc[pSrcMap[nPix]];

			// Subband 2
			nPix = (2*nRow + 1)*nWidth + 2*nCol;
			pDst[2][nCoeff] = pSrc[pSrcMap[nPix]];

			// Subband 3
			nPix = (2*nRow + 1)*nWidth + 2*nCol + 1;
			pDst[3][nCoeff] = pSrc[pSrcMap[nPix]];
		}
	}

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletMultiLevelInterleavePositions(const int nWidth, const int nHeight, const int nLevels, int* pInterleavedToStacked, int* pStackedToInterleaved, int* pInterleavedToSubband)
{
	bool bSuccess = true;

	// Generate initial position map
	int* pInterleavedPositions = new int[nWidth * nHeight]; // positions in interleaved image
	for (int nPos = 0; nPos < nWidth * nHeight; nPos++) {
		pInterleavedPositions[nPos] = nPos;
	}

	// Loop over levels
	int nLevelWidth = nWidth;
	int nLevelHeight = nHeight;
	int nSubbandWidth = nWidth/2;
	int nSubbandHeight = nHeight/2;
	for (int nLevel = nLevels-1; nLevel >= 0; nLevel--) {
		// Loop over interleaved positions
		for (int nInterleavedRow = 0; nInterleavedRow < nLevelHeight; nInterleavedRow++) {
			for (int nInterleavedCol = 0; nInterleavedCol < nLevelWidth; nInterleavedCol++) {
				int nInterleaved = nInterleavedRow*nLevelWidth + nInterleavedCol;
				int nBand, nBandStart;
				// Subband 0
				if ((nInterleavedRow % 2) == 0 && (nInterleavedCol % 2) == 0) {
					nBand = 3*nLevel;
					nBandStart = 0;
				}
				// Subband 1
				else if ((nInterleavedRow % 2) == 0 && (nInterleavedCol % 2) == 1) {
					nBand = 3*nLevel + 1;
					nBandStart = nSubbandWidth;
				}
				// Subband 2
				else if ((nInterleavedRow % 2) == 1 && (nInterleavedCol % 2) == 0) {
					nBand = 3*nLevel + 2;
					nBandStart = nSubbandHeight*nWidth;
				}
				// Subband 3
				else {
					nBand = 3*nLevel + 3;
					nBandStart = nSubbandHeight*nWidth + nSubbandWidth;
				}
				int nBandRow = nInterleavedRow/2;
				int nBandCol = nInterleavedCol/2;
				int nStacked = nBandStart + nBandRow*nWidth + nBandCol;

				pInterleavedToStacked[pInterleavedPositions[nInterleaved]] = nStacked;
				pInterleavedToSubband[pInterleavedPositions[nInterleaved]] = nBand;
				pStackedToInterleaved[nStacked] = pInterleavedPositions[nInterleaved];

			} // end nInterleavedCol
		} // end nInterleavedRow

		// Generate new position map by downsampling current position map
		int* pInterleavedPositionsDown = new int[nSubbandWidth * nSubbandHeight];
		for (int nRow = 0; nRow < nSubbandHeight; nRow++) {
			for (int nCol = 0; nCol < nSubbandWidth; nCol++) {
				int nPixDst = nRow*nSubbandWidth + nCol;
				int nPixSrc = 2*nRow*nLevelWidth + 2*nCol;
				pInterleavedPositionsDown[nPixDst] = pInterleavedPositions[nPixSrc];
			} // end nCol
		} // end nRow
		delete [] pInterleavedPositions;
		pInterleavedPositions = pInterleavedPositionsDown;

		// Prepare for next level
		nSubbandWidth /= 2;
		nSubbandHeight /= 2;
		nLevelWidth /= 2;
		nLevelHeight /= 2;
	} // end nLevel
	delete [] pInterleavedPositions;

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImageWavelet_waveletMultiLevelInterleave(const int nWidth, const int nHeight, const short** pSrc, const int nLevels, short* pDst)
{
	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 );

	short* pImgStacked = new short[nWidth * nHeight];
	RECORD_SUCCESS( ChenImageWavelet_waveletStack(nWidth, nHeight, pSrc, nLevels, pImgStacked), bSuccess );
	for (int nRow = 0; nRow < nHeight; nRow++) {
		for (int nCol = 0; nCol < nWidth; nCol++) {
			int nStacked = nRow*nWidth + nCol;
			pDst[pStackedToInterleaved[nStacked]] = pImgStacked[nStacked];
		} // end nCol
	} // end nRow

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

int ChenImageWavelet_waveletMultiLevelDeinterleave(const int nWidth, const int nHeight, const short* pSrc, const int nLevels, short** pDst)
{
	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 );

	short* pImgStacked = new short[nWidth * nHeight];
	for (int nRow = 0; nRow < nHeight; nRow++) {
		for (int nCol = 0; nCol < nWidth; nCol++) {
			int nStacked = nRow*nWidth + nCol;
			pImgStacked[nStacked] = pSrc[pStackedToInterleaved[nStacked]];
		} // end nCol
	} // end nRow
	ChenImageWavelet_waveletDestack(nWidth, nHeight, pImgStacked, nLevels, pDst);

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