chenimagedct.cpp

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

		RECORD_SUCCESS( ChenImageDCT_blockDCTQuantize(nBlockSize, nImageW, nImageH, pSrcF, pTable, pDstQuantF), bSuccess );

		for (int i = 0; i < nPixels; i++)
		{
			pDstQuant[i] = ChenImage_RoundFloat(pDstQuantF[i]);
		}

		delete [] pSrcF;
		delete [] pDstQuantF;

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

		if (pDstQuantF)
			delete [] pDstQuantF;

		return BAD_RETURN;
	}

#endif

}

int ChenImageDCT_blockDCTQuantize(const int nBlockSize, const int nImageW, const int nImageH, const float* pSrc, const short* pTable, float* pDstQuant)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Initialize forward DCT structure
	IppiDCTFwdSpec_32f* pDCTSpec;
	IppiSize roiSize = {nBlockSize, nBlockSize};
	IppHintAlgorithm algHint = ippAlgHintAccurate;
	RECORD_IPP_SUCCESS( ippiDCTFwdInitAlloc_32f(&pDCTSpec, roiSize, algHint), bSuccess );
	int nDCTBuffSize;
	RECORD_IPP_SUCCESS( ippiDCTFwdGetBufSize_32f(pDCTSpec, &nDCTBuffSize), bSuccess );
	uchar* pDCTBuffer = new uchar [nDCTBuffSize];
	float* pBlock  = new float [nBlockSize*nBlockSize];
	float* pBlockDCT = new float [nBlockSize*nBlockSize];
	float* pTableF = new float [nBlockSize*nBlockSize];

	// Convert table to float format
	int nSrcStep = nBlockSize*sizeof(short);
	int nDstStep = nBlockSize*sizeof(float);
	RECORD_IPP_SUCCESS( ippiConvert_16s32f_C1R((const Ipp16s*) pTable, nSrcStep, (Ipp32f*)pTableF, nDstStep, roiSize), bSuccess );

	// Loop over blocks
	int nBlockRows = nImageH / nBlockSize;
	int nBlockCols = nImageW / nBlockSize;
	for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
		for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {

			// Copy from source image into block
			for (int nRow = 0; nRow < nBlockSize; nRow++) {
				int nRowReal = nBlockRow*nBlockSize + nRow;
				for (int nCol = 0; nCol < nBlockSize; nCol++) {
					int nColReal = nBlockCol*nBlockSize + nCol;
					pBlock[nRow*nBlockSize+nCol] = pSrc[nRowReal*nImageW+nColReal];
				}
			}

			// Perform blockwise DCT
			nSrcStep = nBlockSize*sizeof(float);
			nDstStep = nBlockSize*sizeof(float);
			RECORD_IPP_SUCCESS( ippiDCTFwd_32f_C1R((const Ipp32f*)pBlock, nSrcStep, (Ipp32f*)pBlockDCT, nDstStep, pDCTSpec, (Ipp8u*)pDCTBuffer), bSuccess );

			// Quantize
			nSrcStep = nBlockSize*sizeof(float);
			nDstStep = nBlockSize*sizeof(float);
			RECORD_IPP_SUCCESS( ippiDiv_32f_C1IR((const Ipp32f*)pTableF, nSrcStep, (Ipp32f*)pBlockDCT, nDstStep, roiSize), bSuccess );

			// Copy from block into quantized destination image
			for (int nRow = 0; nRow < nBlockSize; nRow++) {
				int nRowReal = nBlockRow*nBlockSize + nRow;
				for (int nCol = 0; nCol < nBlockSize; nCol++) {
					int nColReal = nBlockCol*nBlockSize + nCol;
					pDstQuant[nRowReal*nImageW+nColReal] = pBlockDCT[nRow*nBlockSize+nCol];
				}
			}

		} // end nBlockCol
	} // end nBlockRow

	delete [] pDCTBuffer;
	delete [] pBlock;
	delete [] pBlockDCT;
	delete [] pTableF;
	RECORD_IPP_SUCCESS( ippiDCTFwdFree_32f(pDCTSpec), bSuccess );
	return bSuccess ? GOOD_RETURN : BAD_RETURN;


#else


	float* pBlock    = new float [nBlockSize*nBlockSize];
	float* pBlockDCT = new float [nBlockSize*nBlockSize];
	float* pTableF   = new float [nBlockSize*nBlockSize];

	// Convert table to float format
	for (int i = 0; i < nBlockSize*nBlockSize; i++)
	{
		pTableF[i] = (float)(pTable[i]);
	}


	// Loop over blocks
	int nBlockRows = nImageH / nBlockSize;
	int nBlockCols = nImageW / nBlockSize;
	for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
		for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {

			// Copy from source image into block
			for (int nRow = 0; nRow < nBlockSize; nRow++) {
				int nRowReal = nBlockRow*nBlockSize + nRow;
				for (int nCol = 0; nCol < nBlockSize; nCol++) {
					int nColReal = nBlockCol*nBlockSize + nCol;
					pBlock[nRow*nBlockSize+nCol] = pSrc[nRowReal*nImageW+nColReal];
				}
			}

			// Perform blockwise DCT
			ChenImageDCT_DCT8x8Fwd_32f(nBlockSize, pBlock, pBlockDCT);

			// Quantize
			for (int i = 0; i < nBlockSize*nBlockSize; i++)
			{
				pBlockDCT[i] = pBlockDCT[i]/pTableF[i];
			}

			// Copy from block into quantized destination image
			for (int nRow = 0; nRow < nBlockSize; nRow++) {
				int nRowReal = nBlockRow*nBlockSize + nRow;
				for (int nCol = 0; nCol < nBlockSize; nCol++) {
					int nColReal = nBlockCol*nBlockSize + nCol;
					pDstQuant[nRowReal*nImageW+nColReal] = pBlockDCT[nRow*nBlockSize+nCol];
				}
			}

		} // end nBlockCol
	} // end nBlockRow

	delete [] pBlock;
	delete [] pBlockDCT;
	delete [] pTableF;

	return GOOD_RETURN;


#endif

}

int ChenImageDCT_blockInvDCTDequantize(const int nBlockSize, const int nImageW, const int nImageH, const short* pSrc, const short* pTable, short* pDst)
{

#ifdef _USEIPP

	bool bSuccess = true;
	int nPixels = nImageW * nImageH;
	float* pSrcF = new float[nPixels];
	float* pDstF = new float[nPixels];

	RECORD_IPP_SUCCESS( ippsConvert_16s32f((const Ipp16s*)pSrc, (Ipp32f*)pSrcF, nPixels), bSuccess );

	RECORD_SUCCESS( ChenImageDCT_blockInvDCTDequantize(nBlockSize, nImageW, nImageH, pSrcF, pTable, pDstF), bSuccess );

	IppRoundMode rndMode = ippRndNear;
	int nScaleFactor = 0;
	RECORD_IPP_SUCCESS( ippsConvert_32f16s_Sfs((const Ipp32f*)pDstF, (Ipp16s*)pDst, nPixels, rndMode, nScaleFactor), bSuccess );

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


#else


	bool bSuccess = true;
	int nPixels = nImageW * nImageH;
	float* pSrcF = new float[nPixels];
	float* pDstF = new float[nPixels];

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

		RECORD_SUCCESS( ChenImageDCT_blockInvDCTDequantize(nBlockSize, nImageW, nImageH, pSrcF, pTable, pDstF), bSuccess );

		for (int i = 0; i < nPixels; i++)
		{
			pDst[i] = ChenImage_RoundFloat(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 ChenImageDCT_blockInvDCTDequantize(const int nBlockSize, const int nImageW, const int nImageH, const float* pSrc, const short* pTable, float* pDst)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Initialize forward DCT structure
	IppiDCTInvSpec_32f* pDCTSpec;
	IppiSize roiSize = {nBlockSize, nBlockSize};
	IppHintAlgorithm algHint = ippAlgHintAccurate;
	RECORD_IPP_SUCCESS( ippiDCTInvInitAlloc_32f(&pDCTSpec, roiSize, algHint), bSuccess );
	int nDCTBuffSize;
	RECORD_IPP_SUCCESS( ippiDCTInvGetBufSize_32f(pDCTSpec, &nDCTBuffSize), bSuccess );
	uchar* pDCTBuffer = new uchar [nDCTBuffSize];
	float* pBlock  = new float [nBlockSize*nBlockSize];
	float* pBlockDCT = new float [nBlockSize*nBlockSize];
	float* pTableF = new float [nBlockSize*nBlockSize];

	// Convert table to float format
	int nSrcStep = nBlockSize*sizeof(short);
	int nDstStep = nBlockSize*sizeof(float);
	RECORD_IPP_SUCCESS( ippiConvert_16s32f_C1R((const Ipp16s*) pTable, nSrcStep, (Ipp32f*)pTableF, nDstStep, roiSize), bSuccess );

	// Loop over blocks
	int nBlockRows = nImageH / nBlockSize;
	int nBlockCols = nImageW / nBlockSize;
	for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
		for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {

			// Copy from source image into block
			for (int nRow = 0; nRow < nBlockSize; nRow++) {
				int nRowReal = nBlockRow*nBlockSize + nRow;
				for (int nCol = 0; nCol < nBlockSize; nCol++) {
					int nColReal = nBlockCol*nBlockSize + nCol;
					pBlockDCT[nRow*nBlockSize+nCol] = pSrc[nRowReal*nImageW+nColReal];
				}
			}

			// Dequantize
			nSrcStep = nBlockSize*sizeof(float);
			nDstStep = nBlockSize*sizeof(float);
			RECORD_IPP_SUCCESS( ippiMul_32f_C1IR((const Ipp32f*)pTableF, nSrcStep, (Ipp32f*)pBlockDCT, nDstStep, roiSize), bSuccess );

			// Perform blockwise DCT
			nSrcStep = nBlockSize*sizeof(float);
			nDstStep = nBlockSize*sizeof(float);
			RECORD_IPP_SUCCESS( ippiDCTInv_32f_C1R((const Ipp32f*)pBlockDCT, nSrcStep, (Ipp32f*)pBlock, nDstStep, pDCTSpec, (Ipp8u*)pDCTBuffer), bSuccess );
			
			// Copy from block into output
			for (int nRow = 0; nRow < nBlockSize; nRow++) {
				int nRowReal = nBlockRow*nBlockSize + nRow;
				for (int nCol = 0; nCol < nBlockSize; nCol++) {
					int nColReal = nBlockCol*nBlockSize + nCol;
					pDst[nRowReal*nImageW+nColReal] = pBlock[nRow*nBlockSize+nCol];
				}
			}

		} // end nBlockCol
	} // end nBlockRow

	delete [] pDCTBuffer;
	delete [] pBlock;
	delete [] pBlockDCT;
	delete [] pTableF;
	RECORD_IPP_SUCCESS( ippiDCTInvFree_32f(pDCTSpec), bSuccess );
	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#else

	bool bSuccess = true;

	float* pBlock    = new float [nBlockSize*nBlockSize];
	float* pBlockDCT = new float [nBlockSize*nBlockSize];

	if ((pBlock != NULL) && (pBlockDCT != NULL))
	{
		// Loop over blocks
		int nBlockRows = nImageH / nBlockSize;
		int nBlockCols = nImageW / nBlockSize;
		for (int nBlockRow = 0; nBlockRow < nBlockRows; nBlockRow++) {
			for (int nBlockCol = 0; nBlockCol < nBlockCols; nBlockCol++) {

				// Copy from source image into block
				for (int nRow = 0; nRow < nBlockSize; nRow++) {
					int nRowReal = nBlockRow*nBlockSize + nRow;
					for (int nCol = 0; nCol < nBlockSize; nCol++) {
						int nColReal = nBlockCol*nBlockSize + nCol;
						pBlockDCT[nRow*nBlockSize+nCol] = pSrc[nRowReal*nImageW+nColReal];
					}
				}

				// Dequantize
				for (int i = 0; i < nBlockSize*nBlockSize; i++)
				{
					pBlockDCT[i] = pTable[i]*pBlockDCT[i];
				}

				// Perform blockwise DCT
				ChenImageDCT_DCT8x8Inv_32f(nBlockSize, pBlockDCT, pBlock);
				
				// Copy from block into output
				for (int nRow = 0; nRow < nBlockSize; nRow++) {
					int nRowReal = nBlockRow*nBlockSize + nRow;
					for (int nCol = 0; nCol < nBlockSize; nCol++) {
						int nColReal = nBlockCol*nBlockSize + nCol;
						pDst[nRowReal*nImageW+nColReal] = pBlock[nRow*nBlockSize+nCol];
					}
				}

			} // end nBlockCol
		} // end nBlockRow

		delete [] pBlock;
		delete [] pBlockDCT;

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

		if (pBlockDCT)
			delete [] pBlockDCT;

		return BAD_RETURN;
	}

#endif

}

int ChenImageDCT_blockDCTQuantizeRange(const int nBlockSize, const short* pTable, short* pRange)
{

#ifdef _USEIPP

	bool bSuccess = true;

	// Convert table to float precision
	int nLen = nBlockSize * nBlockSize;
	float* pA = new float[nLen];
	float* pTableF = new float[nLen];
	RECORD_IPP_SUCCESS( ippsConvert_16s32f((const Ipp16s*) pTable, (Ipp32f*) pTableF, nLen), bSuccess );

	// Form |A| matrix
	RECORD_IPP_SUCCESS( ippsSet_32f((Ipp32f) (1.0/sqrt(8.0)), (Ipp32f*) pA, nLen), bSuccess );
	for (int nI = 1; nI <= nBlockSize-1; nI++) {
		for (int nK = 0; nK <= nBlockSize-1; nK++) {
			pA[nI*nBlockSize + nK] = (float)( 0.25*cos(PI*(2*nK+1)*nI/16.0) );
		}
	}
	RECORD_IPP_SUCCESS( ippsAbs_32f_I((Ipp32f*) pA, nLen), bSuccess );

	// Form ones matrix
	float* pOnes = new float[nLen];
	RECORD_IPP_SUCCESS( ippsSet_32f(1.0, (Ipp32f*) pOnes, nLen), bSuccess );

	// Calculate |A| * Ones(8,8) * |A|^T
	int nSrcWidth = nBlockSize;
	int nSrcHeight = nBlockSize;
	int nSrcStride2 = sizeof(Ipp32f);
	int nSrcStride1 = nBlockSize*sizeof(Ipp32f);
	float* pATranspose = new float[nLen];
	RECORD_IPP_SUCCESS( ippmTranspose_m_32f((const Ipp32f*) pA, nSrcStride1, nSrcStride2, nSrcWidth, nSrcHeight, 
							(Ipp32f*) pATranspose, nSrcStride1, nSrcStride2), bSuccess );
	float* pInterimProduct = new float[nLen];
	RECORD_IPP_SUCCESS( ippmMul_mm_32f((const Ipp32f*) pOnes, nSrcStride1, nSrcStride2, nSrcWidth, nSrcHeight, 
						(const Ipp32f*) pATranspose, nSrcStride1, nSrcStride2, nSrcWidth, nSrcHeight,
						pInterimProduct, nSrcStride1, nSrcStride2), bSuccess );
	float* pRangeF = new float[nLen];
	RECORD_IPP_SUCCESS( ippmMul_mm_32f((const Ipp32f*) pA, nSrcStride1, nSrcStride2, nSrcWidth, nSrcHeight, 
						(const Ipp32f*) pInterimProduct, nSrcStride1, nSrcStride2, nSrcWidth, nSrcHeight,
						pRangeF, nSrcStride1, nSrcStride2), bSuccess );

	// Scale and round
	RECORD_IPP_SUCCESS( ippsMulC_32f_I(255.0, (Ipp32f*) pRangeF, nLen), bSuccess );
	RECORD_IPP_SUCCESS( ippsMulC_32f_I(0.5, (Ipp32f*) pRangeF, nLen), bSuccess );
	RECORD_IPP_SUCCESS( ippsDiv_32f_I((const Ipp32f*) pTableF, (Ipp32f*)pRangeF, nLen), bSuccess );
	pRangeF[0] = (float)( 2*pRangeF[0] + 1.01 );
	IppRoundMode rndMode = ippRndNear;
	int nScaleFactor = 0;
	RECORD_IPP_SUCCESS( ippsConvert_32f16s_Sfs((const Ipp32f*) pRangeF, (Ipp16s*) pRange, nLen, rndMode, nScaleFactor), bSuccess );

	// Clean up
	delete [] pA;
	delete [] pOnes;
	delete [] pATranspose;
	delete [] pInterimProduct;
	delete [] pRangeF;
	delete [] pTableF;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;


#else

	bool bSuccess = true;

	// Convert table to float precision
	int nLen = nBlockSize * nBlockSize;
	float* pA = new float[nLen];

	// Form |A| matrix
	float fV = (float)(1.0/sqrt(8.0));
	for (int i = 0; i < nLen; i++)
	{
		pA[i] = fV;
	}

	for (int nI = 1; nI <= nBlockSize-1; nI++) {
		for (int nK = 0; nK <= nBlockSize-1; nK++) {
			pA[nI*nBlockSize + nK] = (float)( fabs(0.25*cos(PI*(2*nK+1)*nI/16.0)) );
		}
	}

	// Form ones matrix
	float* pOnes = new float[nLen];
	for (int i = 0; i < nLen; i++)
	{
		pOnes[i] = 1.0;
	}

	// Calculate |A| * Ones(8,8) * |A|^T
	Engine* pEngine = ChenImageMatlab_getEngine();


	//Send the input image and parameters to MATLAB
	ChenImageMatlab_putImageInWorkspace(nBlockSize, nBlockSize, 1, pA, "A");

	// Carry out image interpolation
	engEvalString(pEngine, "pRangeF = A * ones(8,8) * A'");

	float* pRangeF = new float[nLen];

	//Get the result image
	ChenImageMatlab_getImageFromWorkspace(nBlockSize, nBlockSize, 1, "pRangeF", pRangeF);

	// Scale and round
	for (int i = 0; i < nLen; i++)
	{
		pRangeF[i] = (float)(255.0*pRangeF[i]);
		pRangeF[i] = (float)(  0.5*pRangeF[i]);
		pRangeF[i] = pRangeF[i]/pTable[i];
	}

	pRangeF[0] = (float)( 2*pRangeF[0] + 1.01 );
	for (int i = 0; i < nLen; i++)
	{
		pRange[i] = ChenImage_RoundFloat(pRangeF[i]);
	}

	// Clean up
	delete [] pA;
	delete [] pOnes;
	delete [] pRangeF;

	return bSuccess ? GOOD_RETURN : BAD_RETURN;

#endif

}