chenimage.cpp

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

		RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)(pSrcExtra - nYShift*nWidthExtra - nXShift), nSrcStep, (Ipp16s*)pDst, nDstStep, roiSize), bSuccess );
	}

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

#else

	bool bSuccess = true;
	if (nXShift >= nWidth || nXShift <= -nWidth || nYShift >= nHeight || nYShift <= -nHeight) return 1;

	// Gray out destination image
	int nSrcStep = nWidth * sizeof(short);
	int nDstStep = nWidth * sizeof(short);
	memset(pDst, 0, nWidth*nHeight*sizeof(short));

	// Extrapolate missing regions
	int nXShiftAbs = nXShift < 0 ? -nXShift : nXShift;
	int nYShiftAbs = nYShift < 0 ? -nYShift : nYShift;
	int nWidthExtra = nWidth + nXShiftAbs;
	int nHeightExtra = nHeight + nYShiftAbs;

	short* pSrcExtra = new short[nWidthExtra * nHeightExtra];
	
	if (pSrcExtra != NULL)
	{
		RECORD_SUCCESS( ChenImage_extrapolateImage((short)nWidth, (short)nHeight, pSrc, -nXShift, -nYShift, pSrcExtra), bSuccess );

		// Copy subregion
		nSrcStep = nWidthExtra * sizeof(short);
		nDstStep = nWidth * sizeof(short);

		short *pSrcHere = NULL;
		short *pDstHere = pDst;

		if (nXShift >= 0 && nYShift >= 0) {
			pSrcHere = (short*)pSrcExtra;
		}
		else if (nXShift < 0 && nYShift >= 0) {
			pSrcHere = (short*)(pSrcExtra - nXShift);
		}
		else if (nXShift >= 0 && nYShift < 0) {
			pSrcHere = (short*)(pSrcExtra - nYShift*nWidthExtra);
		}
		else if (nXShift < 0 && nYShift < 0) {
			pSrcHere = (short*)(pSrcExtra - nYShift*nWidthExtra - nXShift);
		}

		for (int i = 0; i < nHeight; i++)
		{
			memcpy(pDstHere, pSrcHere, nWidth*sizeof(short));
			pSrcHere += nWidthExtra;
			pDstHere += nWidth;
		}

		delete [] pSrcExtra;

		return GOOD_RETURN;
	}
	else
	{
		return BAD_RETURN;
	}


#endif

}

int ChenImage_interpolateImage(const int nInterpolateMethod, const int nWidthSmall, const int nHeightSmall, const short* pSrc, const int nWidthLarge, const int nHeightLarge, short* pDst)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Perform interpolation in float precision
	float* pSrcF = new float[nWidthSmall * nHeightSmall];
	float* pDstF = new float[nWidthLarge * nHeightLarge];
	IppiSize roiSize = {nWidthSmall, nHeightSmall};
	RECORD_IPP_SUCCESS( ippiConvert_16s32f_C1R((const Ipp16s*)pSrc, nWidthSmall*sizeof(short), (Ipp32f*)pSrcF, nWidthSmall*sizeof(float), roiSize), bSuccess );
	RECORD_SUCCESS( ChenImage_interpolateImage(nInterpolateMethod, nWidthSmall, nHeightSmall, pSrcF, nWidthLarge, nHeightLarge, pDstF), bSuccess );
	roiSize.width = nWidthLarge; roiSize.height = nHeightLarge;
	RECORD_IPP_SUCCESS( ippiConvert_32f16s_C1R((const Ipp32f*)pDstF, nWidthLarge*sizeof(float), (Ipp16s*)pDst, nWidthLarge*sizeof(short), roiSize, ippRndNear), bSuccess );

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

#else

	bool bSuccess = true;

	// Perform interpolation in float precision
	float* pSrcF = new float [nWidthSmall * nHeightSmall];
	float* pDstF = new float [nWidthLarge * nHeightLarge];

	if ((pSrcF != NULL) && (pDstF != NULL))
	{
		for (int i = 0; i < nWidthSmall*nHeightSmall; i++)
		{
			pSrcF[i] = (float)(pSrc[i]);
		}

		RECORD_SUCCESS( ChenImage_interpolateImage(nInterpolateMethod, nWidthSmall, nHeightSmall, pSrcF, nWidthLarge, nHeightLarge, pDstF), bSuccess );
	
		for (int i = 0; i < nWidthLarge*nHeightLarge; i++)
		{
			pDst[i] = (short)(pDstF[i]);
		}

		delete [] pSrcF;
		delete [] pDstF;

		return bSuccess ? GOOD_RETURN : BAD_RETURN;
	}
	else
	{
		if (pSrcF)
			delete [] pSrcF;

		if (pDstF)
			delete [] pDstF;

		return BAD_RETURN;
	}


#endif

}

int ChenImage_interpolateImage(const int nInterpolateMethod, const int nWidthSmall, const int nHeightSmall, const float* pSrc, const int nWidthLarge, const int nHeightLarge, float* pDst)
{

#ifdef _USEIPP

	bool bSuccess = true;

	double dFactorX = nWidthLarge/nWidthSmall;
	double dFactorY = nHeightLarge/nHeightSmall;
	IppiSize srcSize = {nWidthSmall, nHeightSmall};
	IppiRect srcRoi = {0, 0, nWidthSmall, nHeightSmall};
	IppiSize dstSize = {nWidthLarge, nHeightLarge};
	RECORD_SUCCESS( ippiResize_32f_C1R((const Ipp32f*)pSrc, srcSize, nWidthSmall*sizeof(float), srcRoi, (Ipp32f*)pDst, nWidthLarge*sizeof(float), dstSize, dFactorX, dFactorY, nInterpolateMethod), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN; 

#else

	Engine* pEngine = ChenImageMatlab_getEngine();

	if (pEngine != NULL)
	{
		//Send the input image and parameters to MATLAB

		ChenImageMatlab_putImageInWorkspace(nWidthSmall, nHeightSmall, 1, pSrc, "InputImage");

		mxArray* pMxWidthLarge = mxCreateDoubleMatrix(1, 1, mxREAL);
		*mxGetPr(pMxWidthLarge) = nWidthLarge;
		engPutVariable(pEngine, "nWidthLarge",  pMxWidthLarge);

		mxArray* pMxHeightLarge = mxCreateDoubleMatrix(1, 1, mxREAL);
		*mxGetPr(pMxHeightLarge) = nHeightLarge;
		engPutVariable(pEngine, "nHeightLarge",  pMxHeightLarge);

		// Carry out image interpolation
		switch(nInterpolateMethod)
		{
			case 1:
			engEvalString(pEngine, "OutputImage = imresize(InputImage, [nHeightLarge nWidthLarge], 'nearest') ");
			break;

			case 2:
			engEvalString(pEngine, "OutputImage = imresize(InputImage, [nHeightLarge nWidthLarge], 'bilinear') ");
			break;

			case 4:
			engEvalString(pEngine, "OutputImage = imresize(InputImage, [nHeightLarge nWidthLarge], 'bicubic') ");
			break;

			default:
			engEvalString(pEngine, "OutputImage = imresize(InputImage, [nHeightLarge nWidthLarge], 'nearest') ");
			break;
		}

		//Get the result image
		ChenImageMatlab_getImageFromWorkspace(nWidthLarge, nHeightLarge, 1, "OutputImage", pDst);

		engEvalString(pEngine, "clear InputImage");
		engEvalString(pEngine, "clear OutputImage");


		return GOOD_RETURN;
	}
	else
	{
		return BAD_RETURN;
	}

#endif

}



int ChenImage_copyImage(const int nWidth, const int nHeight, const short* pSrc, short* pDst)
{
#ifdef _USEIPP

	bool bSuccess = true;

	IppiSize roiSize = {nWidth, nHeight};
	RECORD_IPP_SUCCESS( ippiCopy_16s_C1R((const Ipp16s*)pSrc, nWidth*sizeof(short), pDst, nWidth*sizeof(short), roiSize), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else

	memcpy(pDst, pSrc, nWidth*nHeight*sizeof(short));

	return GOOD_RETURN;

#endif
}

int ChenImage_laplacianPDF(const int nLength, const float fLambda, const float* pSrc, float* pDst)
{
#ifdef _USEIPP

	bool bSuccess = true;

	float* pSrcAbs = new float[nLength];
	RECORD_IPP_SUCCESS( ippsAbs_32f((const Ipp32f*)pSrc, (Ipp32f*)pSrcAbs, nLength), bSuccess );
	RECORD_IPP_SUCCESS( ippsMulC_32f((const Ipp32f*)pSrcAbs, -fLambda, pDst, nLength), bSuccess );
	RECORD_IPP_SUCCESS( ippsExp_32f_I(pDst, nLength), bSuccess );
	RECORD_IPP_SUCCESS( ippsMulC_32f_I((Ipp32f)0.5*fLambda, (Ipp32f*)pDst, nLength), bSuccess );

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

#else

	for (int i = 0; i < nLength; i++)
	{
		pDst[i] = fabs(pSrc[i]);
		pDst[i] = -fLambda*pDst[i];
		pDst[i] = (float)(exp(pDst[i]));
		pDst[i] = (float)(0.5*fLambda*pDst[i]);
	}

	return GOOD_RETURN;

#endif
}

int ChenImage_imageFitBlockLaplacianModel(const int nWidth, const int nHeight, const short* pSrc, const int nBlockSize, float* pLaplacianML, float* pEntropyML, float* pLaplacianME, float* pEntropyME)
{
#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 );

			// Fit Laplacian model for block
			int nBlock = nBlockRow*nBlockCols + nBlockCol;
			RECORD_SUCCESS( ChenImage_imageFitLaplacianModel(nBlockSize, nBlockSize, pBlock, pLaplacianML == NULL ? NULL : pLaplacianML+nBlock, pEntropyML == NULL ? NULL : pEntropyML+nBlock, pLaplacianME == NULL ? NULL : pLaplacianME+nBlock, pEntropyME == NULL ? NULL : pEntropyME+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;

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

				// Fit Laplacian model for block
				int nBlock = nBlockRow*nBlockCols + nBlockCol;
				RECORD_SUCCESS( ChenImage_imageFitLaplacianModel(nBlockSize, nBlockSize, pBlock, pLaplacianML == NULL ? NULL : pLaplacianML+nBlock, pEntropyML == NULL ? NULL : pEntropyML+nBlock, pLaplacianME == NULL ? NULL : pLaplacianME+nBlock, pEntropyME == NULL ? NULL : pEntropyME+nBlock), bSuccess );

			}// end nBlockCol
		} // end nBlockRow

		delete [] pBlock;
		return bSuccess ? GOOD_RETURN : BAD_RETURN;
	}
	else
	{
		return  BAD_RETURN;
	}

#endif


}

int ChenImage_imageFitLaplacianModel(const int nWidth, const int nHeight, const short* pSrc, float* pLaplacianML, float* pEntropyML, float* pLaplacianME, float* pEntropyME)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Calculate histogram
	int nLevels = 511;
	float* pPMF = new float[nLevels];
	RECORD_SUCCESS( ChenImage_diffImagePMF(nWidth, nHeight, pSrc, pPMF), bSuccess );
	int* pLevels = new int[nLevels];
	float fOffset = -255.0, fSlope = 1.0;
	RECORD_IPP_SUCCESS( ippsVectorRamp_32s((Ipp32s*)pLevels, nLevels, fOffset, fSlope), bSuccess );
	float* pLevelsF = new float[nLevels];
	RECORD_IPP_SUCCESS( ippsConvert_32s32f((const Ipp32s*)pLevels, (Ipp32f*)pLevelsF, nLevels), bSuccess );

	//////////////////////////////////////////////////////////////////
	// Use maximum likelihood criterion
	//////////////////////////////////////////////////////////////////
	float fMinLaplacian = (float)0.1, fMaxLaplacian = (float)5.0;
	if (pLaplacianML != NULL && pEntropyML != NULL) {
		short* pSrcAbs = new short[nWidth * nHeight];
		int nSrcStep = nWidth * sizeof(short);
		int nDstStep = nWidth * sizeof(short);
		IppiSize roiSize = {nWidth, nHeight};
		RECORD_IPP_SUCCESS( ippiAbs_16s_C1R((const Ipp16s*)pSrc, nSrcStep, (Ipp16s*)pSrcAbs, nDstStep, roiSize), bSuccess );
		double dSrcAbsSum;
		nSrcStep = nWidth * sizeof(short);
		RECORD_IPP_SUCCESS( ippiSum_16s_C1R((const Ipp16s*)pSrcAbs, nSrcStep, roiSize, &dSrcAbsSum), bSuccess );
		*pLaplacianML = nWidth * nHeight / (float)MAX(EPS,(float)dSrcAbsSum);
		*pLaplacianML = CLIP(*pLaplacianML, fMinLaplacian, fMaxLaplacian);
		delete [] pSrcAbs;
	}

	//////////////////////////////////////////////////////////////////
	// Use min entropy criterion; test range of Laplacian parameters
	//////////////////////////////////////////////////////////////////
	Engine* pEngine = ChenImageMatlab_getEngine();
	if (pLaplacianME != NULL && pEntropyME != NULL) {
		float* pQPMF = new float[nLevels];
		*pEntropyME = LARGE;
		float fMinDistToLaplacianML = LARGE;
		for (float fLaplacian = fMinLaplacian; fLaplacian <= fMaxLaplacian; fLaplacian += (float)0.1) {
			// Form q vector
			RECORD_SUCCESS( ChenImage_laplacianPDF(nLevels, fLaplacian, pLevelsF, pQPMF), bSuccess );