decoderquad.cpp

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

// decoderQuad.cpp

//------------------------------------------------------------
// David Chen*, David Varodayan, Markus Flierl, Bernd Girod
// Image, Video, and Multimedia Systems Group
// Information Systems Laboratory
// Stanford University
//
// *Contact: [e-mail]
//------------------------------------------------------------

#include "decoderQuad.h"
#include "encoderQuad.h"
#include <time.h>
#include <ChenImageDCT.h>

static char pTextBuffer[1024];

int pixelQuads_decodeFrameYUV(const CodingOption* pOption, const short** pFrame, const uchar** pAccumSyndrome, const short* pScaledQArray, const float* pLowerRateBound, const float** pLaplacian, const char* pLadderFile, const short** pOracleShiftsX, const short** pOracleShiftsY, short** pDecodedFrame, float* pRate, short** pDecodedShiftsX, short** pDecodedShiftsY, float** pDecodedShiftProb)
{
	bool bSuccess = true;

	// Decode Y component using quadrant decoder
	RECORD_SUCCESS( pixelQuads_decodeImage(pOption, pFrame[0], pAccumSyndrome, pScaledQArray, pLowerRateBound, pLaplacian, pLadderFile, pOracleShiftsX, pOracleShiftsY, pDecodedFrame[0], pRate, pDecodedShiftsX, pDecodedShiftsY, pDecodedShiftProb), bSuccess );

	// Decode U,V components using regular decoder
	sprintf(pTextBuffer, "decoding quad U \n"); SCREEN_AND_LOG(pTextBuffer, pOption->pLogFile);
	RECORD_SUCCESS( pixelQuads_decodeImageQuad(pOption, pFrame[1], pAccumSyndrome[4], pScaledQArray, pLowerRateBound[4], pLaplacian[4], pLadderFile, pOracleShiftsX[4], pOracleShiftsY[4], pDecodedFrame[1], pRate+4, pDecodedShiftsX[4], pDecodedShiftsY[4], pDecodedShiftProb[4]), bSuccess );
	sprintf(pTextBuffer, "decoding quad V \n"); SCREEN_AND_LOG(pTextBuffer, pOption->pLogFile);
	RECORD_SUCCESS( pixelQuads_decodeImageQuad(pOption, pFrame[2], pAccumSyndrome[5], pScaledQArray, pLowerRateBound[5], pLaplacian[5], pLadderFile, pOracleShiftsX[5], pOracleShiftsY[5], pDecodedFrame[2], pRate+5, pDecodedShiftsX[5], pDecodedShiftsY[5], pDecodedShiftProb[5]), bSuccess );

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int pixelQuads_decodeImage(const CodingOption* pOption, const short* pImgY, const uchar** pAccumSyndrome, const short* pScaledQArray, const float* pLowerRateBound, const float** pLaplacian, const char* pLadderFile, const short** pOracleShiftsX, const short** pOracleShiftsY, short* pDecodedImg, float* pRate, short** pDecodedShiftsX, short** pDecodedShiftsY, float** pDecodedShiftProb)
{
	bool bSuccess = true;

	// Break side info image into four quadrants
	short* pImgYQuads[4];
	short* pImgXQuads[4];
	for (int nQuad = 0; nQuad < 4; nQuad++) {
		pImgYQuads[nQuad] = new short[PIXELS_QUAD];
		pImgXQuads[nQuad] = new short[PIXELS_QUAD];
	}
	RECORD_SUCCESS( ChenImage_imageToQuadrants(WIDTH, HEIGHT, pImgY, pImgYQuads), bSuccess );

	// Decode each quadrant
	for (int nQuad = 0; nQuad < 4; nQuad++) {
		sprintf(pTextBuffer, "decoding quad %d \n", nQuad); SCREEN_AND_LOG(pTextBuffer, pOption->pLogFile);
		RECORD_SUCCESS( pixelQuads_decodeImageQuad(pOption, pImgYQuads[nQuad], pAccumSyndrome[nQuad], pScaledQArray, pLowerRateBound[nQuad], pLaplacian[nQuad], pLadderFile, ((pOracleShiftsX == NULL) ? NULL : pOracleShiftsX[nQuad]), ((pOracleShiftsY == NULL) ? NULL : pOracleShiftsY[nQuad]), pImgXQuads[nQuad], pRate+nQuad, pDecodedShiftsX[nQuad], pDecodedShiftsY[nQuad], pDecodedShiftProb[nQuad]), bSuccess );
	}

	// Synthesize entire image from decoded quadrants
	RECORD_SUCCESS( ChenImage_quadrantsToImage(WIDTH, HEIGHT, (const short**)pImgXQuads, pDecodedImg), bSuccess );

	// Clean up
	for (int nQuad = 0; nQuad < 4; nQuad++) {
		delete [] pImgYQuads[nQuad];
		delete [] pImgXQuads[nQuad];
	}
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int pixelQuads_decodeImageQuad(const CodingOption* pOption, const short* pImgY, const unsigned char* pAccumSyndrome, const short* pScaledQArray, const float fLowerRateBound, const float* pLaplacian, const char* pLadderFile, const short* pOracleShiftsX, const short* pOracleShiftsY, short* pDecodedImg, float* pRate, short* pDecodedShiftsX, short* pDecodedShiftsY, float* pDecodedShiftProb)
{
	bool bSuccess = true;

	// Find quantized coefficient range
	short* pQCoeffRange = new short[NUM_COEFFS];
	RECORD_SUCCESS( ChenImageDCT_blockDCTQuantizeRange(BLOCK_SIZE, pScaledQArray, pQCoeffRange), bSuccess );

	// Generate side info
	short* pSideInfoQuant = new short[PIXELS_QUAD * NUM_SHIFTS_2D];
	RECORD_SUCCESS( pixelQuads_generateSideInfoQuad(pImgY, pScaledQArray, pSideInfoQuant), bSuccess );

	// Apply block-based decoder
	double* pDecodedBitsBlockD = new double[BITS_QUAD];
	double* pProbDispBlock = new double[BLOCKS_QUAD * NUM_SHIFTS_2D]; // Resultant disparity PMF
	double* pProbDispPixel = new double[PIXELS_QUAD * NUM_SHIFTS_2D];
	RECORD_SUCCESS( pixelQuads_decodeBitsQuad(pOption, pSideInfoQuant, pAccumSyndrome, fLowerRateBound, pLaplacian, pQCoeffRange, pOracleShiftsX, pOracleShiftsY, pDecodedBitsBlockD, pRate, pProbDispBlock, pProbDispPixel), bSuccess );
	unsigned char* pDecodedBitsBlock = new unsigned char[BITS_QUAD];
	for (int nBit = 0; nBit < BITS_QUAD; nBit++) {
		pDecodedBitsBlock[nBit] = (unsigned char)(pDecodedBitsBlockD[nBit] + 0.5);
	}

	// Convert bitstream back to transform coefficients
	short* pDecodedTransQuantBlock = new short[PIXELS_QUAD];
	RECORD_SUCCESS( pixelQuads_bitstreamToSignedCoefficientsQuad(pDecodedBitsBlock, pDecodedTransQuantBlock), bSuccess );

	// Reconstruct using motion-compensated nearest-neighbor binning + inverse transform
	if (pOption->nReconstructMethod == RECON_MOTION) {
		RECORD_SUCCESS( pixelQuads_motionCompensateReconstructQuad(pImgY, pDecodedTransQuantBlock, pProbDispBlock, pScaledQArray, pDecodedImg), bSuccess );
	}
	// Reconstruct using centroid of quantization interval + inverse transform
	else {
		RECORD_SUCCESS( ChenImageDCT_block88InvDCTDequantize(WIDTH_QUAD, HEIGHT_QUAD, pDecodedTransQuantBlock, pScaledQArray, pDecodedImg), bSuccess );
	}

	// Find max-probability elements in motion fields
	RECORD_SUCCESS( pixelQuads_findMaxMotionProbabilityQuad(pProbDispBlock, pDecodedShiftsX, pDecodedShiftsY, pDecodedShiftProb), bSuccess );
	
	// Clean up
	delete [] pQCoeffRange;
	delete [] pSideInfoQuant;
	delete [] pDecodedBitsBlockD;
	delete [] pProbDispBlock;
	delete [] pProbDispPixel;
	delete [] pDecodedBitsBlock;
	delete [] pDecodedTransQuantBlock;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int pixelQuads_pixelDecoderQuad(const int nCodeLength, const short* pSideInfo, const double* pChannel, const double* pMotionPixel, const short* pQuantCoeffRange, const double* pLLROutput, double* pLLRSideInfo, double* pProbArray)
{
    double pProbSum[2];    
    for (int nPixel = 0; nPixel < nCodeLength/BITPLANES; nPixel++) // loop over pixels
    {
        // blend
        int nCol = nPixel % WIDTH_QUAD;
        int nRow = nPixel / WIDTH_QUAD;
        int nCoeffCol = nPixel % BLOCK_SIZE;
        int nCoeffRow = (nPixel/WIDTH_QUAD) % BLOCK_SIZE;
        int nCoeffNum = BLOCK_SIZE*nCoeffRow + nCoeffCol;
        
		//------------------------------------------------------------
		// Motion probabilities to gray-level probabilities
		//------------------------------------------------------------
		memset(pProbArray, 0, LEVELS*sizeof(double));

		for (int nShiftY = 0; nShiftY < NUM_SHIFTS; nShiftY++) {
			for (int nShiftX = 0; nShiftX < NUM_SHIFTS; nShiftX++) {
				int nShift = nShiftY*NUM_SHIFTS + nShiftX;
				if( pMotionPixel[NUM_SHIFTS_2D*nPixel + nShift] > 0 )
				{
					int nSideInfoStart = nShift*PIXELS_QUAD;
					int nSideInfo = pSideInfo[nSideInfoStart + nRow*WIDTH_QUAD + nCol];
	                
					// DC coefficient
					if (nCoeffNum == 0)
					{
						for (int nLevel = 0; nLevel < LEVELS; nLevel++) {
							if (nLevel <= pQuantCoeffRange[nCoeffNum]) {
								pProbArray[nLevel] += pMotionPixel[NUM_SHIFTS_2D*nPixel + nShift] * pChannel[(2*LEVELS-1)*nCoeffNum + (LEVELS-1) - nSideInfo + nLevel];
							}
						} // end nLevel
					}
					// AC coefficient
					else
					{
						for (int nLevel = -LEVELS/2; nLevel < LEVELS/2; nLevel++) {
							if (abs(nLevel) <= pQuantCoeffRange[nCoeffNum]) {
								pProbArray[(nLevel+LEVELS) % LEVELS] += pMotionPixel[NUM_SHIFTS_2D*nPixel + nShift] * pChannel[(2*LEVELS-1)*nCoeffNum + (LEVELS-1) - nSideInfo + nLevel]; 
							}
						} // end nLevel
					}
				}  
			} // end nShiftX
		} // end nShiftY
        
		//------------------------------------------------------------
		// Gray-level probabilities to bit probabilities
		//------------------------------------------------------------
        for (int nBit = 0; nBit < BITPLANES; nBit++)
        {
            double p0 = min(max((1 + tanh(pLLROutput[nPixel*BITPLANES + nBit]/2))/2, 0.0001), 0.9999);
			for (int nLevel = 0; nLevel < LEVELS; nLevel++) {
				if (nLevel & ((LEVELS/2)>>nBit)) {
                    pProbArray[nLevel] *= 1-p0;
				}
				else {
                    pProbArray[nLevel] *= p0;
				}
			} // end nLevel
        } // end nBit
        
        for (int nBit = 0; nBit < BITPLANES; nBit++)
        {
            pProbSum[0] = 0;
            pProbSum[1] = 0;
            double p0 = MIN(MAX((1 + tanh(pLLROutput[nPixel*BITPLANES + nBit]/2))/2, 0.0001), 0.9999);
			for (int nLevel = 0; nLevel < LEVELS; nLevel++) {
				if (nLevel & ((LEVELS/2)>>nBit)) {
                    pProbSum[1] += MAX(EPS, pProbArray[nLevel]/(1-p0));
				}
				else {
                    pProbSum[0] += pProbArray[nLevel]/p0;
				}
			} // end nLevel
            pLLRSideInfo[nPixel*BITPLANES + nBit] = log(pProbSum[0]/pProbSum[1]);
        } // end nBit
        
    } // end nPixel
	return GOOD_RETURN;
}


int pixelQuads_ldpcDecoderQuad(const int* pIR, const int* pJC, const int nSubcodeLength, const int nCodeLength, const int nCodeEdges, const double* pSyndrome, const double* pLLRSideInfo, double* pLLROutput, double* pLLRMemory, double* pCheckLLR, double* pCheckLLRMag, double* pRowTotal)
{
    //------------------------------------------------------------
	// Step 1: compute check-to-variable messages
	//------------------------------------------------------------
	for (int k = 0; k < nCodeLength; k++) {
		for (int l = pJC[k]; l < pJC[k+1]; l++) {
            pLLRMemory[l] += pLLRSideInfo[k];
		}
	}
    
    for (int k = 0; k < nCodeEdges; k++) {
        pCheckLLR[k] = ((pLLRMemory[k] < 0) ? -1 : 1);
        pCheckLLRMag[k] = fabs(pLLRMemory[k]);
    }
    
	for (int k = 0; k < nSubcodeLength; k++) {
        pRowTotal[k] = ((pSyndrome[k] == 1) ? -1 : 1);
	}
	for (int k = 0; k < nCodeEdges; k++) {
        pRowTotal[pIR[k]] *= pCheckLLR[k];
	}
	for (int k = 0; k < nCodeEdges; k++) {
        pCheckLLR[k] = pCheckLLR[k] * pRowTotal[pIR[k]];
        //sign of check-to-variable messages
	}
        
	for (int k = 0; k < nCodeEdges; k++) {
        pCheckLLRMag[k] = -log( tanh( max(pCheckLLRMag[k], 0.000000001)/2 ) );
	}
	memset(pRowTotal, 0, nSubcodeLength * sizeof(double));
	for (int k = 0; k < nCodeEdges; k++) {
        pRowTotal[pIR[k]] += pCheckLLRMag[k];
	}
	for (int k = 0; k < nCodeEdges; k++) {
        pCheckLLRMag[k] = -log( tanh( max(pRowTotal[pIR[k]] - pCheckLLRMag[k], 0.000000001)/2 ) );
        //magnitude of check-to-variable messages
	}
        
	for (int k = 0; k < nCodeEdges; k++) {
        pCheckLLR[k] = pCheckLLR[k] * pCheckLLRMag[k];
        //check-to-variable messages
	}
       
    //------------------------------------------------------------
	// Step 2: compute variable-to-check messages
	//------------------------------------------------------------
    for (int k = 0; k < nCodeLength; k++) {
        pLLROutput[k] = 0.0;
		for (int l = pJC[k]; l < pJC[k+1]; l++) {
            pLLROutput[k] += pCheckLLR[l];
		}
    }
        
	for (int k = 0; k < nCodeLength; k++) {
		for (int l = pJC[k]; l < pJC[k+1]; l++) {
            pLLRMemory[l] = pLLROutput[k] - pCheckLLR[l];
            //variable-to-check messages
		}
	}

	return GOOD_RETURN;
}

int pixelQuads_motionEstimatorQuad(const double* pLLROutput, const short* pSideInfo, const double* pChannel, const short* pQuantCoeffRange, double* pThetaArray, double* pMotionBlock)
{
    double* pThetaNew = new double[BITPLANES];
	for (int nBlockRow = 0; nBlockRow < BLOCK_ROWS_QUAD; nBlockRow++) {
        for (int nBlockCol=0; nBlockCol < BLOCK_COLS_QUAD; nBlockCol++) {    
            int nBlock = nBlockRow*BLOCK_COLS_QUAD + nBlockCol;

			for (int l = 0; l < LEVELS*BLOCK_SIZE*BLOCK_SIZE; l++) {
                pThetaArray[l] = 1.0;
			}
            
			//------------------------------------------------------------
			// Bit probabilities to gray-level probabilities
			//------------------------------------------------------------
			for (int nRow = 0; nRow < BLOCK_SIZE; nRow++) {
                for (int nCol = 0; nCol < BLOCK_SIZE; nCol++) {
					for (int nBit = 0; nBit < BITPLANES; nBit++) {
                        pThetaNew[nBit] = (1 + tanh(pLLROutput[((nBlockRow*BLOCK_SIZE + nRow)*WIDTH_QUAD + nBlockCol*BLOCK_SIZE + nCol)*BITPLANES + nBit]/2))/2;
					}

					for (int nLevel = 0; nLevel < LEVELS; nLevel++) {
						for (int nBit = 0; nBit < BITPLANES; nBit++) {
							if (nLevel & ((LEVELS/2)>>nBit)) {
                                pThetaArray[(nRow*BLOCK_SIZE + nCol)*LEVELS + nLevel] *= 1-pThetaNew[nBit];
							}
							else {
                                pThetaArray[(nRow*BLOCK_SIZE + nCol)*LEVELS + nLevel] *= pThetaNew[nBit];
							}
						} // end nBit
					} // end nLevel
                } // end nCol
			} // nRow
         
			//------------------------------------------------------------
			// Gray-level probabilities to disparity probabilities
			//------------------------------------------------------------
            double dProbabilitySum = 0;
			for (int nShiftY = -MAX_SHIFT; nShiftY <= MAX_SHIFT; nShiftY++) {
				for (int nShiftX = -MAX_SHIFT; nShiftX <= MAX_SHIFT; nShiftX++) {
					int nShift = (nShiftY+MAX_SHIFT)*NUM_SHIFTS + nShiftX+MAX_SHIFT;
					if ( pMotionBlock[nBlock*NUM_SHIFTS_2D + nShift] != 0 )
					{
						for (int nRow = 0; nRow < BLOCK_SIZE; nRow++) {
							for (int nCol = 0; nCol  <BLOCK_SIZE; nCol++) {

								double dProbability = 0;
								int nSideInfoStart = nShift*PIXELS_QUAD;
								int nSideInfo = pSideInfo[nSideInfoStart + (nBlockRow*BLOCK_SIZE + nRow)*WIDTH_QUAD + nBlockCol*BLOCK_SIZE + nCol];
	                            
								// DC coefficient
								if (nRow == 0 && nCol == 0)
								{
									for (int nLevel = 0; nLevel < LEVELS; nLevel++) {
										if (nLevel <= pQuantCoeffRange[nRow*BLOCK_SIZE + nCol]) {
											dProbability += pThetaArray[(nRow*BLOCK_SIZE + nCol)*LEVELS + nLevel] * pChannel[(2*LEVELS-1)*(nRow*BLOCK_SIZE + nCol) + (LEVELS-1) - nSideInfo + nLevel];
										}
									} // end nLevel
								}
								// AC coefficient
								else
								{
									for (int nLevel = -LEVELS/2; nLevel < LEVELS/2; nLevel++) {
										if (abs(nLevel) <= pQuantCoeffRange[nRow*BLOCK_SIZE + nCol]) {
											dProbability += pThetaArray[(nRow*BLOCK_SIZE + nCol)*LEVELS + ((nLevel + LEVELS) % LEVELS)] * pChannel[(2*LEVELS-1)*(nRow*BLOCK_SIZE + nCol) + (LEVELS-1) - nSideInfo + nLevel];
										}
									}
								}
	                            
								pMotionBlock[nBlock*NUM_SHIFTS_2D + nShift] *= dProbability;

							} // end nCol
						} // end nRow
	                    
						dProbabilitySum += pMotionBlock[nBlock*NUM_SHIFTS_2D + nShift];
					}
				} // end nShiftX
			} // end nShiftY
            
			//------------------------------------------------------------
			// Normalize disparity probabilities
			//------------------------------------------------------------