decoderquad.cpp

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

            pTxSeq[k] = pTxSeqRead[nCode][k];
		}
		for (int k = 0; k < nCodeEdges; k++) {
            pIR[k] = pIrRead[nCode][k];
		}
        int nSubcodeLength = pNumEdgesRead[nCode];
        
        //------------------------------------------------------------
		// Construct syndrome
		//------------------------------------------------------------
        int nCurrIndex = pTxSeq[0], nPrevIndex;
        pSyndrome[0] = pParity[nCurrIndex];
        for (int k = 1; k < nSubcodeLength; k++) {
            nPrevIndex = nCurrIndex;
            nCurrIndex = pTxSeq[k%nNumInc] + (k/nNumInc)*nTotalNumInc;
            pSyndrome[k] = (double) (((int) (pParity[nCurrIndex] + pParity[nPrevIndex])) % 2);
        }
                
        //------------------------------------------------------------
		// Initialize LDPC memory
		//------------------------------------------------------------
		memset(pDecodedBits, 0, nCodeLength * sizeof(double));
		memset(pLLRMemory, 0, nCodeEdges * sizeof(double));
		memset(pLLROutput, 0, nCodeLength * sizeof(double));
        
        //------------------------------------------------------------
		// Initialize motion field
		//------------------------------------------------------------
		// Motion oracle
		if (pOption->bMotionOracle) {
			RECORD_SUCCESS( pixelQuads_initializeMotionOracleQuad(pOracleShiftsX, pOracleShiftsY, pMotionBlock), bSuccess );
		}
		// Motion learning
		else if (pOption->bMotionEstimate) {
			RECORD_SUCCESS( pixelQuads_initializeMotionLearningQuad(pMotionBlock), bSuccess );
		}
		// No motion compensation
		else { 
			RECORD_SUCCESS( pixelQuads_initializeMotionNoCompensationQuad(pMotionBlock), bSuccess );
		}

		// Interpolate motion probabilities from block to pixel
		RECORD_SUCCESS( pixelQuads_interpolateMotionQuad(pOption->nInterpolateMethod, pMotionBlock, pMotionPixel), bSuccess );

		//------------------------------------------------------------
		// Loop over iterations of EM
		//------------------------------------------------------------        
		sprintf(pTextBuffer, "code = %d: ", nCode); SCREEN_AND_LOG(pTextBuffer, pOption->pLogFile);

		float fSuccessRateMax = 0;
		for (int nIteration = 0; nIteration < MAX_ITERATIONS; nIteration++) {

			//------------------------------------------------------------
			// Use LDPC pixel and bit decoders
			//------------------------------------------------------------
			// Joint bitplane decoder
            RECORD_SUCCESS( pixelQuads_pixelDecoderQuad(nCodeLength, pSideInfo, pChannel, pMotionPixel, pQuantCoeffRange, pLLROutput, pLLRSideInfo, pProbArray), bSuccess );
            
			// Message-passing bit decoder
            RECORD_SUCCESS( pixelQuads_ldpcDecoderQuad(pIR, pJC, nSubcodeLength, nCodeLength, nCodeEdges, pSyndrome, pLLRSideInfo, pLLROutput, pLLRMemory, pCheckLLR, pCheckLLRMag, pRowTotal), bSuccess );  
            
			//------------------------------------------------------------
			// Use MAP decoder
			//------------------------------------------------------------
            int l = 0;
			for (int k = 0; k < nCodeLength; k++) {
                if (pDecodedBits[k] == ((pLLRSideInfo[k] + pLLROutput[k] < 0) ? 1 : 0))
                    l++;
                else
                    pDecodedBits[k] = ((pLLRSideInfo[k] + pLLROutput[k] < 0) ? 1 : 0);       
			}
            
			//------------------------------------------------------------
			// Check row totals and estimate partial decoding errors
			//------------------------------------------------------------
			for (int k = 0; k < nSubcodeLength; k++) {
                pRowTotal[k] = pSyndrome[k];
			}
			for (int k = 0; k < nCodeLength; k++) {
				for (int l = pJC[k]; l < pJC[k+1]; l++) {
                    pRowTotal[pIR[l]] += pDecodedBits[k];
				}
			}
               
			int nErrors = 0;
			for (int k = 0; k < nSubcodeLength; k++) {
                if (((int) pRowTotal[k] % 2) != 0)
					nErrors++;
                else if (k == nSubcodeLength-1 && nErrors == 0)
                {
					bConverged = true;
					break;
					// all syndrome checks satisfied
                }
			}

			//------------------------------------------------------------
			// Decide whether or not to terminate decoding
			//------------------------------------------------------------
			float fSuccessRate = ((float)(nSubcodeLength-nErrors))/((float)nSubcodeLength);
			fSuccessRateMax = MAX( fSuccessRateMax, fSuccessRate );
			if (bConverged || (nIteration == MAX_ITERATIONS-1)) {
				sprintf(pTextBuffer, "iter. = %d, success = %.4f, ", nIteration+1, fSuccessRateMax); SCREEN_AND_LOG(pTextBuffer, pOption->pLogFile);
			
				time(&nTimeEnd);
				int nTimeElapsed = (int)(nTimeEnd-nTimeStart);
				sprintf(pTextBuffer, "time (sec) = %d \n", nTimeElapsed); SCREEN_AND_LOG(pTextBuffer, pOption->pLogFile);

				*pDecodedRate = pRateRead[nCode];
				break;
			}
            
			//------------------------------------------------------------
			// Blockwise disparity estimate followed by interpolation
			//------------------------------------------------------------
			if (pOption->bMotionEstimate && !pOption->bMotionOracle) {
				// Estimate blockwise motion probabilities
				RECORD_SUCCESS( pixelQuads_motionEstimatorQuad(pLLROutput, pSideInfo, pChannel, pQuantCoeffRange, pThetaArray, pMotionBlock), bSuccess );
				
				// Interpolate motion probabilities from block to pixel
				RECORD_SUCCESS( pixelQuads_interpolateMotionQuad(pOption->nInterpolateMethod, pMotionBlock, pMotionPixel), bSuccess );
			}
			
        } // end nIteration

		if (bConverged) break;

    } // end nCode

	//------------------------------------------------------------
	// Clean up memory
	//------------------------------------------------------------
	delete [] pIR;
    delete [] pJC;
	delete [] pTxSeq;
	delete [] pSyndrome;
	
    delete [] pLLRMemory;
	delete [] pCheckLLR;
	delete [] pCheckLLRMag;
	delete [] pRowTotal;
	delete [] pLLROutput;
	delete [] pLLRSideInfo;

	delete [] pChannel;
    delete [] pProbArray;
    delete [] pThetaArray;

	for (int nCode = 0; nCode < nNumCodes; nCode++) {
		delete [] pTxSeqRead[nCode];
		delete [] pIrRead[nCode];
	}
	delete [] pNumIncRead;
	delete [] pTxSeqRead;
	delete [] pIrRead;
	delete [] pNumEdgesRead;
	delete [] pRateRead;
	
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int pixelQuads_generateSideInfoQuad(const short* pImgY, const short* pScaledQArray, short* pSideInfoQuant)
{
	bool bSuccess = true;

	// Initialize side info
	memset(pSideInfoQuant, 0, PIXELS_QUAD*NUM_SHIFTS_2D*sizeof(short));

	// Calculate shifted versions
	short* pImgYShift = new short[PIXELS_QUAD];
	short* pTransQuantYShift = new short[PIXELS_QUAD];
	for (int nShiftY = 0; nShiftY < NUM_SHIFTS; nShiftY++) {
		for (int nShiftX = 0; nShiftX < NUM_SHIFTS; nShiftX++) {
			// Create shifted image
			RECORD_SUCCESS( ChenImage_shiftImage(WIDTH_QUAD, HEIGHT_QUAD, pImgY, -(nShiftX-MAX_SHIFT), -(nShiftY-MAX_SHIFT), pImgYShift), bSuccess );	

			// Perform 8-by-8 block DCT and quantized
			RECORD_SUCCESS( ChenImageDCT_block88DCTQuantize(WIDTH_QUAD, HEIGHT_QUAD, pImgYShift, pScaledQArray, pTransQuantYShift), bSuccess );

			// Copy quantized DCT coefficients into output
			int nShift = nShiftY*NUM_SHIFTS + nShiftX;
			int nShiftImageStart = nShift*PIXELS_QUAD;
			for (int nRow = 0; nRow < HEIGHT_QUAD; nRow++) {
				for (int nCol = 0; nCol < WIDTH_QUAD; nCol++) {
					pSideInfoQuant[nShiftImageStart + nRow*WIDTH_QUAD + nCol] = pTransQuantYShift[nRow*WIDTH_QUAD + nCol];
				}
			}
		} // end nShiftX
	} // end nShiftY
	delete [] pImgYShift;
	delete [] pTransQuantYShift;

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int pixelQuads_motionCompensateReconstructQuad(const short* pImgY, const short* pDecodedTransQuant, const double* pProbDisp, const short* pScaledQArray, short* pImgRecover)
{
	bool bSuccess = true;

	// Find max shift probability in each block
	short* pMaxProbIndexX = new short[BLOCKS_QUAD];
	short* pMaxProbIndexY = new short[BLOCKS_QUAD];
	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;
			
			double dMaxProb = 0.0;
			for (int nShiftY = -MAX_SHIFT; nShiftY <= MAX_SHIFT; nShiftY++) {
				double dProbSum = 0.0;
				for (int nShiftX = -MAX_SHIFT; nShiftX <= MAX_SHIFT; nShiftX++) {
					int nShift = (nShiftY+MAX_SHIFT)*NUM_SHIFTS + nShiftX+MAX_SHIFT;
					dProbSum += pProbDisp[nBlock*NUM_SHIFTS_2D + nShift];
				} // end nShiftX
				if (dProbSum > dMaxProb) {
					dMaxProb = dProbSum;
					pMaxProbIndexY[nBlock] = (short)nShiftY;
				}
			} // end nShiftY

			dMaxProb = 0.0;
			for (int nShiftX = -MAX_SHIFT; nShiftX <= MAX_SHIFT; nShiftX++) {
				double dProbSum = 0.0;
				for (int nShiftY = -MAX_SHIFT; nShiftY <= MAX_SHIFT; nShiftY++) {
					int nShift = (nShiftY+MAX_SHIFT)*NUM_SHIFTS + nShiftX+MAX_SHIFT;
					dProbSum += pProbDisp[nBlock*NUM_SHIFTS_2D + nShift];
				} // end nShiftY
				if (dProbSum > dMaxProb) {
					dMaxProb = dProbSum;
					pMaxProbIndexX[nBlock] = (short)nShiftX;
				}
			} // end nShiftX

		} // end nBlockCol
	} // end nBlockRow

	// Use disparity-compensated reconstruction
	short* pImgYComp = new short[PIXELS_QUAD];
	for (int nPix = 0; nPix < PIXELS_QUAD; nPix++) {
		pImgYComp[nPix] = pImgY[nPix];
	}
	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;

			// Copy disparity-compensated block
			int nShiftX = pMaxProbIndexX[nBlock];
			int nShiftY = pMaxProbIndexY[nBlock];
			for (int nRow = 0; nRow < BLOCK_SIZE; nRow++) {
				int nRowReal = nBlockRow*BLOCK_SIZE + nRow;
				int nRowRealShift = nRowReal + nShiftY;
				if ((nRowRealShift < 0) || (nRowRealShift >= HEIGHT_QUAD)) continue;
				for (int nCol = 0; nCol < BLOCK_SIZE; nCol++) {
					int nColReal = nBlockCol*BLOCK_SIZE + nCol;
					int nColRealShift = nColReal + nShiftX;
					if ((nColRealShift < 0) || (nColRealShift >= WIDTH_QUAD)) continue;
					pImgYComp[nRowReal*WIDTH_QUAD + nColReal] = pImgY[nRowRealShift*WIDTH_QUAD + nColRealShift];
				} // end nCol
			} // end nRow

		} // end nBlockCol
	} // end nBlockRow
	// ChenImage_showImage(WIDTH_QUAD, HEIGHT_QUAD, 1, pImgYComp, "y-comp");

	// Transform and quantize
	short* pTransYComp = new short[PIXELS_QUAD];
	RECORD_SUCCESS( ChenImageDCT_block88DCTQuantize(WIDTH_QUAD, HEIGHT_QUAD, pImgYComp, QArrayUnity, pTransYComp), bSuccess );
	short* pTransQuantYComp = new short[PIXELS_QUAD];
	RECORD_SUCCESS( ChenImageDCT_block88DCTQuantize(WIDTH_QUAD, HEIGHT_QUAD, pImgYComp, pScaledQArray, pTransQuantYComp), bSuccess );

	// Decide final reconstruction to use
	float* pDecodedTransF = new float[PIXELS_QUAD];
	for (int nBlockRow = 0; nBlockRow < BLOCK_ROWS_QUAD; nBlockRow++) {
		for (int nBlockCol = 0; nBlockCol < BLOCK_COLS_QUAD; nBlockCol++) {

			// Use nearest bin
			for (int nRow = 0; nRow < BLOCK_SIZE; nRow++) {
				int nRowReal = nBlockRow*BLOCK_SIZE + nRow;
				for (int nCol = 0; nCol < BLOCK_SIZE; nCol++) {
					int nColReal = nBlockCol*BLOCK_SIZE + nCol;
					short nValueQuant = pDecodedTransQuant[nRowReal*WIDTH_QUAD + nColReal];
					short nValueQuantComp = pTransQuantYComp[nRowReal*WIDTH_QUAD + nColReal];
					short nValueComp = pTransYComp[nRowReal*WIDTH_QUAD + nColReal];
					if (nValueQuant < nValueQuantComp)
						pDecodedTransF[nRowReal*WIDTH_QUAD + nColReal] = (float)( (nValueQuant + 0.5) * pScaledQArray[nRow*BLOCK_SIZE + nCol] );
					else if (nValueQuant > nValueQuantComp)
						pDecodedTransF[nRowReal*WIDTH_QUAD + nColReal] = (float)( (nValueQuant - 0.5) * pScaledQArray[nRow*BLOCK_SIZE + nCol] );
					else
						pDecodedTransF[nRowReal*WIDTH_QUAD + nColReal] = nValueComp;
				} // end nCol
			} // end nRow

		} // end nBlockCol
	} // end nBlockRow

	// Dequantize and inverse-transform
	float* pImgRecoverF = new float[PIXELS_QUAD];
	RECORD_SUCCESS( ChenImageDCT_block88InvDCTDequantize(WIDTH_QUAD, HEIGHT_QUAD, pDecodedTransF, QArrayUnity, pImgRecoverF), bSuccess );

	// Round 
	for (int nPix = 0; nPix < PIXELS_QUAD; nPix++) {
		pImgRecover[nPix] = (short) ROUND( pImgRecoverF[nPix] );
	}

	// Clean up
	delete [] pMaxProbIndexX;
	delete [] pMaxProbIndexY;
	delete [] pImgYComp;
	delete [] pTransQuantYComp;
	delete [] pTransYComp;
	delete [] pDecodedTransF;
	delete [] pImgRecoverF;
	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}

int pixelQuads_findMaxMotionProbabilityQuad(const double* pProbDispBlock, short* pDecodedShiftsX, short* pDecodedShiftsY, float* pDecodedShiftProb)
{
	bool bSuccess = true;

	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;

			double dMaxProb = 0.0;
			for (int nShiftY = 0; nShiftY < NUM_SHIFTS; nShiftY++) {
				for (int nShiftX = 0; nShiftX < NUM_SHIFTS; nShiftX++) {
					int nShift = nShiftY*NUM_SHIFTS + nShiftX;
					
					if (pProbDispBlock[nBlock*NUM_SHIFTS_2D + nShift] > dMaxProb) {
						dMaxProb = pProbDispBlock[nBlock*NUM_SHIFTS_2D + nShift];
						pDecodedShiftsX[nBlock] = (short)(nShiftX - MAX_SHIFT);
						pDecodedShiftsY[nBlock] = (short)(nShiftY - MAX_SHIFT);
						pDecodedShiftProb[nBlock] = (float)dMaxProb;
					}

				} // end nShiftX
			} // end nShiftY

		} // end nBlockCol
	} // end nBlockRow

	return bSuccess ? GOOD_RETURN : BAD_RETURN;
}