chenimage.cpp

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

	delete [] pLevels;
	delete [] pHist;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

	#else

	// Calculate histogram
	int nLevels = 510;
	int nShift  = 255;

	float* pHist = new float [nLevels];

	if (pHist != NULL)
	{
		int i (0);

		// Calculate the histogram of the input image block

		memset(pHist, 0, nLevels*sizeof(float));

		int nNumPixels = nWidth*nHeight;

		for (i = 0; i < nNumPixels; i++)
		{
			pHist[(pSrc[i]+nShift)] += 1.0;
		}

		// Calculate PMF
		float fHistSum (0.0);
		for (i = 0; i < nLevels; i++)
		{
			fHistSum += pHist[i];
		}

		for (i = 0; i < nLevels; i++)
		{
			pPMF[i] = pHist[i]/fHistSum;
		}

		delete [] pHist;

		return GOOD_RETURN;
	}
	else
	{
		return BAD_RETURN;
	}

#endif

}

int ChenImage_imageEntropy(const int nWidth, const int nHeight, const short* pSrc, float* pEntropy)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Calculate PMF
	int nLevels = 256;
	float* pPMF = new float[nLevels];
	RECORD_SUCCESS( ChenImage_imagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	for (int nLevel = 0; nLevel < nLevels; nLevel++) {
		pPMF[nLevel] = (float) MAX(EPS, pPMF[nLevel]);
	}
	
	// Calculate entropy
	float* pLnPMF = new float[nLevels];
	RECORD_IPP_SUCCESS( ippsLn_32f((const Ipp32f*)pPMF, (Ipp32f*)pLnPMF, nLevels), bSuccess );
	RECORD_IPP_SUCCESS( ippsMul_32f_I((const Ipp32f*)pPMF, (Ipp32f*)pLnPMF, nLevels), bSuccess );
	IppHintAlgorithm hint = ippAlgHintAccurate;
	RECORD_IPP_SUCCESS( ippsSum_32f((const Ipp32f*)pLnPMF, nLevels, pEntropy, hint), bSuccess );
	*pEntropy /= -(float)log(2.0);
	
	delete [] pLnPMF;
	delete [] pPMF;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else

	bool bSuccess = true;

	// Calculate PMF
	int nLevels = 256;
	float* pPMF   = new float [nLevels];
	float* pLnPMF = new float [nLevels];

	if ((pPMF != NULL) && (pLnPMF != NULL))
	{
		RECORD_SUCCESS( ChenImage_imagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	
		*pEntropy = 0;
		for (int nLevel = 0; nLevel < nLevels; nLevel++) 
		{
			pPMF[nLevel]   = (float) MAX(EPS, pPMF[nLevel]);
			pLnPMF[nLevel] = (float)(log(pPMF[nLevel]));
			*pEntropy += pLnPMF[nLevel];
		}
		*pEntropy /= -(float)log(2.0);
		
		delete [] pLnPMF;
		delete [] pPMF;

		return GOOD_RETURN;
	}
	else
	{
		if (pPMF)
			delete [] pPMF;

		if (pLnPMF)
			delete [] pLnPMF;

		return BAD_RETURN;
	}

#endif
}

int ChenImage_diffImageEntropy(const int nWidth, const int nHeight, const short* pSrc, float* pEntropy)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Calculate PMF
	int nLevels = 511;
	float* pPMF = new float[nLevels];
	RECORD_SUCCESS( ChenImage_diffImagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	for (int nLevel = 0; nLevel < nLevels; nLevel++) {
		pPMF[nLevel] = (float)MAX(EPS, pPMF[nLevel]);
	}
	
	// Calculate entropy
	float* pLnPMF = new float[nLevels];
	RECORD_IPP_SUCCESS( ippsLn_32f((const Ipp32f*)pPMF, (Ipp32f*)pLnPMF, nLevels), bSuccess );
	RECORD_IPP_SUCCESS( ippsMul_32f_I((const Ipp32f*)pPMF, (Ipp32f*)pLnPMF, nLevels), bSuccess );
	IppHintAlgorithm hint = ippAlgHintAccurate;
	RECORD_IPP_SUCCESS( ippsSum_32f((const Ipp32f*)pLnPMF, nLevels, pEntropy, hint), bSuccess );
	*pEntropy /= -(float)log(2.0);
	float entropy = *pEntropy;
	
	delete [] pLnPMF;
	delete [] pPMF;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;


#else

	bool bSuccess = true;

	// Calculate PMF
	int nLevels = 511;
	float* pPMF   = new float[nLevels];
	float* pLnPMF = new float[nLevels];

	if ((pPMF != NULL) && (pLnPMF != NULL))
	{
		RECORD_SUCCESS( ChenImage_diffImagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	
		for (int nLevel = 0; nLevel < nLevels; nLevel++) {
			pPMF[nLevel] = (float)MAX(EPS, pPMF[nLevel]);
		}

		*pEntropy = 0;
		for (int nLevel = 0; nLevel < nLevels; nLevel++) 
		{
			pPMF[nLevel]   = (float) MAX(EPS, pPMF[nLevel]);
			pLnPMF[nLevel] = (float)(log(pPMF[nLevel]));
			*pEntropy += pLnPMF[nLevel];
		}
		*pEntropy /= -(float)log(2.0);
		
		delete [] pLnPMF;
		delete [] pPMF;

		return GOOD_RETURN;
	}
	else
	{
		if (pPMF)
			delete [] pPMF;

		if (pLnPMF)
			delete [] pLnPMF;

		return BAD_RETURN;
	}

#endif

}

int ChenImage_imageBlockMatch1D(const int nWidth, const int nHeight, const int nBlockSize, const int nMaxShift, const short* pSrc1, const short* pSrc2, short* pShifted, short* pResidual, short* pShifts)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Initialize
	int nBlockRows = nHeight / nBlockSize;
	int nBlockCols = nWidth / nBlockSize;
	short* pBlock1 = new short[nBlockSize * nBlockSize];
	short* pBlock2 = new short[nBlockSize * nBlockSize];
	short* pBlockRes = new short[nBlockSize * nBlockSize];
	short* pBlockResAbs = new short[nBlockSize * nBlockSize];

	// Copy src 2 into output by default
	int nSrcStep = nWidth * sizeof(short);
	int nDstStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pSrc2, nSrcStep, (Ipp16s*)pShifted, nDstStep, roiSize), bSuccess );

	// Iterate across all blocks
	double dMinResAbsSum;
	int nSrc1StartCol, nSrc1StartRow, nSrc2StartCol, nSrc2StartRow;
	roiSize.width = nBlockSize; roiSize.height = nBlockSize;
	for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
		for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {
			int nBlock = nBlockRow*nBlockCols + nBlockCol;
			dMinResAbsSum = LARGE;
			nSrc1StartCol = nBlockSize * nBlockCol;
			nSrc1StartRow = nBlockSize * nBlockRow;

			// Copy block from X
			nSrcStep = nWidth * sizeof(short);
			nDstStep = nBlockSize * sizeof(short);
			RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pSrc1 + nSrc1StartRow*nWidth + nSrc1StartCol), nSrcStep, (Ipp16s*)pBlock1, nDstStep, roiSize), bSuccess );

			// Iterate across all local horizontal
			nSrc2StartRow = nSrc1StartRow;
			for (nSrc2StartCol = MAX(0,nSrc1StartCol-nMaxShift); nSrc2StartCol <= MIN(nSrc1StartCol+nMaxShift,nWidth-nBlockSize); nSrc2StartCol++) {

				// Copy block from Y
				nSrcStep = nWidth * sizeof(short);
				nDstStep = nBlockSize * sizeof(short);
				RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pSrc2 + nSrc2StartRow*nWidth + nSrc2StartCol), nSrcStep, (Ipp16s*)pBlock2, nDstStep, roiSize), bSuccess );

				// Create residual block
				nSrcStep = nBlockSize * sizeof(short);
				nDstStep = nBlockSize * sizeof(short);
				int nScaleFactor = 0;
				RECORD_IPP_SUCCESS( ippiSub_16s_C1RSfs((const Ipp16s*)pBlock2, nSrcStep, (const Ipp16s*)pBlock1, nSrcStep, (Ipp16s*)pBlockRes, nDstStep, roiSize, nScaleFactor), bSuccess );

				// Calculate sum of absolute residuals
				double dResAbsSum = LARGE;
				nSrcStep = nBlockSize * sizeof(short);
				nDstStep = nBlockSize * sizeof(short);
				RECORD_IPP_SUCCESS( ippiAbs_16s_C1R((const Ipp16s*)pBlockRes, nSrcStep, (Ipp16s*)pBlockResAbs, nDstStep, roiSize), bSuccess );
				RECORD_IPP_SUCCESS( ippiSum_16s_C1R((const Ipp16s*)pBlockResAbs, nSrcStep, roiSize, &dResAbsSum), bSuccess );

				// Compare with previous min sum of absolute residuals
				if (dResAbsSum < dMinResAbsSum) {
					dMinResAbsSum = dResAbsSum;
					// Copy shifted block into output image
					nSrcStep = nBlockSize * sizeof(short);
					nDstStep = nWidth * sizeof(short);
					RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pBlock2, nSrcStep, (Ipp16s*)(pShifted + nSrc1StartRow*nWidth + nSrc1StartCol), nDstStep, roiSize), bSuccess );

					// Copy block residual into output residual image
					nSrcStep = nBlockSize * sizeof(short);
					nDstStep = nWidth * sizeof(short);
					RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pBlockRes, nSrcStep, (Ipp16s*)(pResidual + nSrc1StartRow*nWidth + nSrc1StartCol), nDstStep, roiSize), bSuccess );

					// Record optimal shift value
					if (pShifts != NULL) {
						pShifts[nBlock] = nSrc2StartCol - nSrc1StartCol;
					}
				}

			} // end nSrc2StartCol
		
		} // end nBlockCol
	} // end nBlockRow
	delete [] pBlock1;
	delete [] pBlock2;
	delete [] pBlockRes;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;


#else

	// Initialize
	int nBlockRows = nHeight / nBlockSize;
	int nBlockCols = nWidth  / nBlockSize;
	
	short* pBlock1      = new short[nBlockSize * nBlockSize];
	short* pBlock2      = new short[nBlockSize * nBlockSize];
	short* pBlockRes    = new short[nBlockSize * nBlockSize];
	short* pBlockResAbs = new short[nBlockSize * nBlockSize];

	if ((pBlock1 != NULL) && (pBlock2 != NULL) && (pBlockRes != NULL) && (pBlockResAbs != NULL))
	{
		// Copy src 2 into output by default
		memcpy(pShifted, pSrc2, nWidth*nHeight*sizeof(short));

		// Iterate across all blocks
		double dMinResAbsSum;
		int nSrc1StartCol, nSrc1StartRow, nSrc2StartCol, nSrc2StartRow;

		for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
			for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {
				int nBlock = nBlockRow*nBlockCols + nBlockCol;
				dMinResAbsSum = LARGE;
				nSrc1StartCol = nBlockSize * nBlockCol;
				nSrc1StartRow = nBlockSize * nBlockRow;

				// Copy block from X
				short* pSrcHere1 = (short*)(pSrc1) + (nSrc1StartRow*nWidth + nSrc1StartCol);
				short* pDstHere1 = pBlock1;

				for (int i = 0; i < nBlockSize; i++)
				{
					memcpy(pDstHere1, pSrcHere1, nBlockSize*sizeof(short));
					pSrcHere1 += nWidth;
					pDstHere1 += nBlockSize;
				}

				// Iterate across all local horizontal
				nSrc2StartRow = nSrc1StartRow;
				for (nSrc2StartCol = MAX(0,nSrc1StartCol-nMaxShift); nSrc2StartCol <= MIN(nSrc1StartCol+nMaxShift,nWidth-nBlockSize); nSrc2StartCol++) {

					// Copy block from Y
					short* pSrcHere2 = (short*)(pSrc2) + (nSrc2StartRow*nWidth + nSrc2StartCol);
					short* pDstHere2 = pBlock2;

					for (int i = 0; i < nBlockSize; i++)
					{
						memcpy(pDstHere2, pSrcHere2, nBlockSize*sizeof(short));
						pSrcHere2 += nWidth;
						pDstHere2 += nBlockSize;
					}

					// Create residual block and sum of absolute residuals
					double dResAbsSum (0);
					for (int k = 0; k < nBlockSize*nBlockSize; k++)
					{
						pBlockRes[k]    = pBlock1[k] - pBlock2[k];
						pBlockResAbs[k] = (short)(abs(pBlockRes[k]));
						dResAbsSum += pBlockResAbs[k] ;
					}

					// Compare with previous min sum of absolute residuals
					if (dResAbsSum < dMinResAbsSum) {
						dMinResAbsSum = dResAbsSum;

						// Copy shifted block into output image
						pSrcHere2 = pBlock2;
						pDstHere2 = (pShifted + nSrc1StartRow*nWidth + nSrc1StartCol);
						for (int i = 0; i < nBlockSize; i++)
						{
							memcpy(pDstHere2, pSrcHere2, nBlockSize*sizeof(short));
							pSrcHere2 += nBlockSize;
							pDstHere2 += nWidth;
						}

						// Copy block residual into output residual image
						pSrcHere2 = pBlockRes;
						pDstHere2 = (pResidual + nSrc1StartRow*nWidth + nSrc1StartCol);
						for (int i = 0; i < nBlockSize; i++)
						{
							memcpy(pDstHere2, pSrcHere2, nBlockSize*sizeof(short));
							pSrcHere2 += nBlockSize;
							pDstHere2 += nWidth;
						}

						// Record optimal shift value
						if (pShifts != NULL) {
							pShifts[nBlock] = (short)(nSrc2StartCol - nSrc1StartCol);
						}
					}

				} // end nSrc2StartCol
			
			} // end nBlockCol
		} // end nBlockRow
		
		delete [] pBlock1;
		delete [] pBlock2;
		delete [] pBlockRes;
		delete [] pBlockResAbs;