chenimage.cpp

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

// ChenImage.cpp

#include <ChenImage.h>
#include <ChenImageMatlab.h>

#ifdef _USEIPP
#include <ippi.h>
#include <ippcv.h>
#include <ipps.h>
#include <ippdc.h>
#endif


short ChenImage_RoundFloat(float fV)
{
	if (fV < 0.0)
	{
		return (short)(fV - 0.5);
	}
	else
	{
		return (short)(fV + 0.5);
	}
}

int ChenImage_imageMeanStdDev(const int nWidth, const int nHeight, const short* pSrc, float* pMean, float* pStdDev)
{
#ifdef _USEIPP

	bool bSuccess = true;

	float* pSrcF = new float[nWidth * nHeight];
	int nSrcStep = nWidth * sizeof(short);
	int nDstStep = nWidth * sizeof(float);
	IppiSize roiSize = {nWidth, nHeight};
	RECORD_IPP_SUCCESS( ippiConvert_16s32f_C1R((const Ipp16s*)pSrc, nSrcStep, (Ipp32f*)pSrcF, nDstStep, roiSize), bSuccess );
	nSrcStep = nWidth * sizeof(float);
	double dMean, dStdDev;
	RECORD_IPP_SUCCESS( ippiMean_StdDev_32f_C1R((const Ipp32f*)pSrcF, nSrcStep, roiSize, &dMean, &dStdDev), bSuccess );
	*pMean = (float)dMean;
	*pStdDev = (float)dStdDev;

	delete [] pSrcF;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else // use C code

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

		int nNumPixels = nWidth*nHeight;

		double dMean (0); 
		for (i = 0; i < nNumPixels; i++)
		{
			dMean += pSrc[i];
		}
		*pMean  = (float) dMean/nNumPixels;

		double dStdDev (0);
		for (i = 0; i < nNumPixels; i++)
		{
			double dV = (pSrc[i] - *pMean); 
			dStdDev += dV*dV; 
		}
		// (nNumPixels-1) is not used for the baseline match with IPP mode
		*pStdDev = (float)(sqrt(dStdDev/(nNumPixels))); 

		return GOOD_RETURN;
	}
	else
	{
		return BAD_RETURN; 
	}

#endif

}

int ChenImage_imageBlockMeanStdDev(const int nWidth, const int nHeight, const short* pSrc, const int nBlockSize, float* pMean, float* pStdDev)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Allocate space
	short* pBlock = new short[nBlockSize * nBlockSize];

	// Loop over blocks
	int nBlockCols = nWidth / nBlockSize;
	int nBlockRows = nHeight / nBlockSize;
	for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
		for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {
			
			// Copy block from source
			int nSrcStep = nWidth * sizeof(short);
			int nDstStep = nBlockSize * sizeof(short);
			IppiSize roiSize = {nBlockSize, nBlockSize};
			int nRow = nBlockRow*nBlockSize;
			int nCol = nBlockCol*nBlockSize;
			RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pSrc + nRow*nWidth + nCol), nSrcStep, (Ipp16s*)pBlock, nDstStep, roiSize), bSuccess );

			// Calculate mean and standard deviation of block
			int nBlock = nBlockRow*nBlockCols + nBlockCol;
			RECORD_SUCCESS( ChenImage_imageMeanStdDev(nBlockSize, nBlockSize, pBlock, pMean+nBlock, pStdDev+nBlock), bSuccess );

		}// end nBlockCol
	} // end nBlockRow

	delete [] pBlock;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else

	bool bSuccess = true;

	// Allocate space
	short* pBlock = new short [nBlockSize * nBlockSize];

	if (pBlock != NULL)
	{
		// Loop over blocks
		int nBlockCols = nWidth / nBlockSize;
		int nBlockRows = nHeight / nBlockSize;
		for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
			for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {
				
				// Copy block from source
				int nRow = nBlockRow*nBlockSize;
				int nCol = nBlockCol*nBlockSize;

				short *pSrcHere = ((short*)pSrc + nRow*nWidth + nCol);
				short *pDstHere = pBlock;
				size_t nCount = nBlockSize*sizeof(short);

				for (int i = 0; i < nBlockSize; i++)
				{
					memcpy(pDstHere, pSrcHere, nCount);
					pSrcHere += nWidth;
					pDstHere += nBlockSize;
				}

				// Calculate mean and standard deviation of block
				int nBlock = nBlockRow*nBlockCols + nBlockCol;
				RECORD_SUCCESS( ChenImage_imageMeanStdDev(nBlockSize, nBlockSize, pBlock, pMean+nBlock, pStdDev+nBlock), bSuccess );

			}// end nBlockCol
		} // end nBlockRow

		delete [] pBlock;

		return bSuccess ? GOOD_RETURN : BAD_RETURN;
	}
	else
	{
		return BAD_RETURN;
	}

#endif

}

int ChenImage_imageMSE(const int nWidth, const int nHeight, const short* pSrc1, const short* pSrc2, float* pMSE)
{
#ifdef _USEIPP

	bool bSuccess = true;

	// Calculate difference image
	short* pDiffImg = new short[nWidth * nHeight];
	int nSrcStep = nWidth * sizeof(short);
	int nDstStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	int nScaleFactor = 0;
	ippiSub_16s_C1RSfs((const Ipp16s*)pSrc1, nSrcStep, (const Ipp16s*)pSrc2, nSrcStep, (Ipp16s*)pDiffImg, nDstStep, roiSize, nScaleFactor);

	// Calculate MSE
	float fMean, fStdDev;
	RECORD_SUCCESS( ChenImage_imageMeanStdDev(nWidth, nHeight, pDiffImg, &fMean, &fStdDev), bSuccess );
	*pMSE = fStdDev * fStdDev;

	delete [] pDiffImg;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else

	bool bSuccess = true;

	// Calculate difference image

	int nNumPixels = (nWidth * nHeight);

	short* pDiffImg = new short [nNumPixels];

	if (pDiffImg != NULL)
	{
		for (int i = 0; i < nNumPixels; i++)
		{
			pDiffImg[i] = (pSrc2[i] - pSrc1[i]);
		}

		// Calculate MSE
		float fMean (0), fStdDev (0);
		RECORD_SUCCESS( ChenImage_imageMeanStdDev(nWidth, nHeight, pDiffImg, &fMean, &fStdDev), bSuccess );
		*pMSE = fStdDev * fStdDev;

		delete [] pDiffImg;

		return GOOD_RETURN;
	}
	else
	{
		return BAD_RETURN;
	}


#endif

}

int ChenImage_imagePSNR(const int nWidth, const int nHeight, const short* pSrc1, const short* pSrc2, float* pPSNR)
{
	bool bSuccess = true;
	
	float fMSE;
	RECORD_SUCCESS( ChenImage_imageMSE(nWidth, nHeight, pSrc1, pSrc2, &fMSE), bSuccess );
	*pPSNR = (float) (10.0*log10(255.0*255.0 / fMSE));

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int ChenImage_imageDistortionEpsilonSqr(const int nWidth, const int nHeight, const short* pSrc, float* pEpsilonSqr)
{	
#ifdef _USEIPP

	bool bSuccess = true;

	// Get min, max
	short nMin, nMax;
	int nSrcStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	RECORD_IPP_SUCCESS( ippiMinMax_16s_C1R((const Ipp16s*)pSrc, nSrcStep, roiSize, &nMin, &nMax), bSuccess );

	// Calculate mean and standard deviation
	float fStdDev, fMean;
	RECORD_SUCCESS( ChenImage_imageMeanStdDev(nWidth, nHeight, pSrc, &fMean, &fStdDev), bSuccess );

	// Calculate PMF
	// Difference image
	float* pPMF;
	int nLevels, nZeroCutoff;
	if (nMin < 0) { 
		nZeroCutoff = 255;
		nLevels = 511;
		pPMF = new float[nLevels];
		RECORD_SUCCESS( ChenImage_diffImagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	}
	// Regular image
	else {
		nZeroCutoff = ROUND(fMean)-1;
		nLevels = 256;
		pPMF = new float[nLevels];
		RECORD_SUCCESS( ChenImage_imagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	}

	// Calculate integral of cube root of PMF (beyond zero cutoff point)
	float* pPMFCubrt = new float[nLevels];
	RECORD_IPP_SUCCESS( ippsCubrt_32f((const Ipp32f*)pPMF, (Ipp32f*)pPMFCubrt, nLevels), bSuccess );
	float fPMFCubrtSum;
	IppHintAlgorithm hint = ippAlgHintAccurate;
	RECORD_IPP_SUCCESS( ippsSum_32f((const Ipp32f*)pPMFCubrt+nZeroCutoff, nLevels-nZeroCutoff, &fPMFCubrtSum, hint), bSuccess );

	// Calculate epsilon-squared
	*pEpsilonSqr = (float) (2.0/3.0 * pow(fPMFCubrtSum,3) / (fStdDev*fStdDev));

	delete [] pPMF;
	delete [] pPMFCubrt;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else

	return BAD_RETURN;

#endif

}

int ChenImage_imagePMF(const int nWidth, const int nHeight, const short* pSrc, float* pPMF)
{
#ifdef _USEIPP

	bool bSuccess = true;

	// Calculate histogram
	int nLevels = 256;
	int nSrcStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	int* pLevels = new int[nLevels];
	float fOffset = 0.0, fSlope = 1.0;
	RECORD_IPP_SUCCESS( ippsVectorRamp_32s((Ipp32s*)pLevels, nLevels, fOffset, fSlope), bSuccess );
	int* pHist = new int[nLevels];
	RECORD_IPP_SUCCESS( ippsZero_32s((Ipp32s*)pHist, nLevels), bSuccess );
	RECORD_IPP_SUCCESS( ippiHistogramRange_16s_C1R((const Ipp16s*)pSrc, nSrcStep, roiSize, (Ipp32s*)pHist, (const Ipp32s*)pLevels, nLevels), bSuccess );

	// Calculate PMF
	int nScaleFactor = 0, nHistSum;
	RECORD_IPP_SUCCESS( ippsSum_32s_Sfs((const Ipp32s*)pHist, nLevels, &nHistSum, nScaleFactor), bSuccess );
	RECORD_IPP_SUCCESS( ippsConvert_32s32f((const Ipp32s*)pHist, (Ipp32f*)pPMF, nLevels), bSuccess );
	RECORD_IPP_SUCCESS( ippsDivC_32f_I((float)nHistSum, (Ipp32f*)pPMF, nLevels), bSuccess );

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


#else

	int nLevels = 256;
	float* pHist = new float[nLevels];

	if (pHist != NULL)
	{
		int i;

		// Calculate the histogram
		memset(pHist, 0, nLevels*sizeof(float));

		int nNumPixels = nWidth*nHeight;

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

		// Calculate PMF
		float fHistSum (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_diffImagePMF(const int nWidth, const int nHeight, const short* pSrc, float* pPMF)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Calculate histogram
	int nLevels = 511;
	int nSrcStep = nWidth * sizeof(short);
	IppiSize roiSize = {nWidth, nHeight};
	int* pLevels = new int[nLevels];
	float fOffset = -255.0, fSlope = 1.0;
	RECORD_IPP_SUCCESS( ippsVectorRamp_32s((Ipp32s*)pLevels, nLevels, fOffset, fSlope), bSuccess );
	int* pHist = new int[nLevels];
	RECORD_IPP_SUCCESS( ippsZero_32s((Ipp32s*)pHist, nLevels), bSuccess );
	RECORD_IPP_SUCCESS( ippiHistogramRange_16s_C1R((const Ipp16s*)pSrc, nSrcStep, roiSize, (Ipp32s*)pHist, (const Ipp32s*)pLevels, nLevels), bSuccess );

	// Calculate PMF
	int nScaleFactor = 0, nHistSum;
	RECORD_IPP_SUCCESS( ippsSum_32s_Sfs((const Ipp32s*)pHist, nLevels, &nHistSum, nScaleFactor), bSuccess );
	RECORD_IPP_SUCCESS( ippsConvert_32s32f((const Ipp32s*)pHist, (Ipp32f*)pPMF, nLevels), bSuccess );
	RECORD_IPP_SUCCESS( ippsDivC_32f_I((float)nHistSum, (Ipp32f*)pPMF, nLevels), bSuccess );