ratectl.c

h.264标准和jm

									   PAverageQP=(int)(1.0*TotalFrameQP/TotalNumberofBasicUnit+0.5);
									   FrameQPBuffer=PAverageQP;
									   FrameAveHeaderBits=PAveHeaderBits2;
								   }
								   else
								   {
									   PAverageQP=(int)(1.0*TotalFrameQP/TotalNumberofBasicUnit+0.5);
									   FieldQPBuffer=PAverageQP;
									   FieldAveHeaderBits=PAveHeaderBits2;
								   }
							   }
						   }
					   }
					   if(GOPOverdue==TRUE)
						   Pm_Qp=PAveFrameQP;
					   else
						   Pm_Qp=m_Qc;
					   return m_Qc;
				   }
				   else
				   {
		         /*predict the MAD of current picture*/
					   if(((input->InterlaceCodingOption==1)||(input->InterlaceCodingOption==3))\
						   &&(img->FieldControl==1))
					   {
						   CurrentFrameMAD=MADPictureC1*FCBUPFMAD[TotalNumberofBasicUnit-NumberofBasicUnit]+MADPictureC2;
						   TotalBUMAD=0;
						   for(i=TotalNumberofBasicUnit-1; i>=(TotalNumberofBasicUnit-NumberofBasicUnit);i--)
						   {
							   CurrentBUMAD=MADPictureC1*FCBUPFMAD[i]+MADPictureC2;
							   TotalBUMAD +=CurrentBUMAD*CurrentBUMAD;
						   }
					   }
					   else
					   {
						   CurrentFrameMAD=MADPictureC1*BUPFMAD[TotalNumberofBasicUnit-NumberofBasicUnit]+MADPictureC2;
						   TotalBUMAD=0;
						   for(i=TotalNumberofBasicUnit-1; i>=(TotalNumberofBasicUnit-NumberofBasicUnit);i--)
						   {
							   CurrentBUMAD=MADPictureC1*BUPFMAD[i]+MADPictureC2;
							   TotalBUMAD +=CurrentBUMAD*CurrentBUMAD;
						   }
					   }
					      
			       /*compute the total number of bits for the current basic unit*/
					   m_Bits =(int)(T*CurrentFrameMAD*CurrentFrameMAD/TotalBUMAD);
					   /*compute the number of texture bits*/
					   m_Bits -=PAveHeaderBits2;
					   
					   m_Bits=MAX(m_Bits,(int)(bit_rate/(MINVALUE*frame_rate*TotalNumberofBasicUnit)));
				
					   dtmp = CurrentFrameMAD * m_X1 * CurrentFrameMAD * m_X1 \
						   + 4 * m_X2 * CurrentFrameMAD * m_Bits;
					   if ((m_X2 == 0.0) || (dtmp < 0) || ((sqrt (dtmp) - m_X1 * CurrentFrameMAD) <= 0.0))	// fall back 1st order mode
						   m_Qstep = (float)(m_X1 * CurrentFrameMAD / (double) m_Bits);
					   else // 2nd order mode
						   m_Qstep = (float) ((2 * m_X2 * CurrentFrameMAD) / (sqrt (dtmp) - m_X1 * CurrentFrameMAD));

					   m_Qc=Qstep2QP(m_Qstep);
					   m_Qc = MIN(m_Qp+DDquant,  m_Qc);	// control variation

						 if(input->basicunit>=MBPerRow)
						   m_Qc = MIN(PAveFrameQP+6, m_Qc);
					   else
						   m_Qc = MIN(PAveFrameQP+3, m_Qc);

					   m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
					   m_Qc = MAX(m_Qp-DDquant, m_Qc);	// control variation 
					   if(input->basicunit>=MBPerRow)
						   m_Qc = MAX(PAveFrameQP-6, m_Qc);
					   else
						   m_Qc = MAX(PAveFrameQP-3, m_Qc);

					   m_Qc = MAX(RC_MIN_QUANT, m_Qc);
					   TotalFrameQP +=m_Qc;
					   Pm_Qp=m_Qc;
					   NumberofBasicUnit--;
					   if((NumberofBasicUnit==0)&&(img->type==P_SLICE))
					   {
						   if((!topfield)||(img->FieldControl==0))
						   {
							   /*frame coding or field coding*/
							   if((input->InterlaceCodingOption==0)||(input->InterlaceCodingOption==2))
							   {
								   PAverageQP=(int)(1.0*TotalFrameQP/TotalNumberofBasicUnit+0.5);
								   if (img->NumberofPPicture == (input->intra_period - 2))
							            QPLastPFrame = PAverageQP;

								   img->TotalQpforPPicture +=PAverageQP;
								   PreviousQp1=PreviousQp2;
								   PreviousQp2=PAverageQP; 
								   PAveFrameQP=PAverageQP;
								   PAveHeaderBits3=PAveHeaderBits2;
							   }
							   else if((input->InterlaceCodingOption==1)\
								   ||(input->InterlaceCodingOption==3))
							   {
								   if(img->FieldControl==0)
								   {
									   PAverageQP=(int)(1.0*TotalFrameQP/TotalNumberofBasicUnit+0.5);
									   FrameQPBuffer=PAverageQP;
									   FrameAveHeaderBits=PAveHeaderBits2;
								   }
								   else
								   {
									   PAverageQP=(int)(1.0*TotalFrameQP/TotalNumberofBasicUnit+0.5);
									   FieldQPBuffer=PAverageQP;
									   FieldAveHeaderBits=PAveHeaderBits2;
								   }
							   }
						   }
					   }
					   return m_Qc;
				   }
			   }
			}
		}	
	}
}

//update the parameters of quadratic R-D model
void updateRCModel ()
{

	int n_windowSize;
	int i;
	double error[20], std = 0.0, threshold;
	int m_Nc;
	Boolean MADModelFlag = FALSE;
   
	if(img->type==P_SLICE)
	{
		/*frame layer rate control*/
		if(img->BasicUnit==img->Frame_Total_Number_MB)
		{
			CurrentFrameMAD=ComputeFrameMAD();
			m_Nc=img->NumberofCodedPFrame;
		}
		/*basic unit layer rate control*/
		else
		{
			/*compute the MAD of the current basic unit*/
			if((input->InterlaceCodingOption==3)&&(img->FieldControl==0))
				CurrentFrameMAD=img->TotalMADBasicUnit/img->BasicUnit/2;
			else
				CurrentFrameMAD=img->TotalMADBasicUnit/img->BasicUnit;

				
			img->TotalMADBasicUnit=0;
		    			
			/* compute the average number of header bits*/
			
				CodedBasicUnit=TotalNumberofBasicUnit-NumberofBasicUnit;
				if(CodedBasicUnit>0)
				{
					PAveHeaderBits1=(int)(1.0*(PAveHeaderBits1*(CodedBasicUnit-1)+\
						+img->NumberofBasicUnitHeaderBits)/CodedBasicUnit+0.5);
					if(PAveHeaderBits3==0)
						PAveHeaderBits2=PAveHeaderBits1;
					else
						PAveHeaderBits2=(int)(1.0*(PAveHeaderBits1*CodedBasicUnit+\
						+PAveHeaderBits3*NumberofBasicUnit)/TotalNumberofBasicUnit+0.5);
				}
					/*update the record of MADs for reference*/
			    if(((input->InterlaceCodingOption==1)||(input->InterlaceCodingOption==3))\
				  &&(img->FieldControl==1))
     				FCBUCFMAD[TotalNumberofBasicUnit-1-NumberofBasicUnit]=CurrentFrameMAD;
	    		else
		    		BUCFMAD[TotalNumberofBasicUnit-1-NumberofBasicUnit]=CurrentFrameMAD;

				if(NumberofBasicUnit!=0)
					m_Nc=img->NumberofCodedPFrame*TotalNumberofBasicUnit+CodedBasicUnit;
				else
					m_Nc=(img->NumberofCodedPFrame-1)*TotalNumberofBasicUnit+CodedBasicUnit;
			
		}
		
		
	 
		if(m_Nc>1)
			MADModelFlag=TRUE;
		PPreHeader=img->NumberofHeaderBits;
		for (i = 19; i > 0; i--) {// update the history
			Pm_rgQp[i] = Pm_rgQp[i - 1];
			m_rgQp[i]=Pm_rgQp[i];
		    Pm_rgRp[i] = Pm_rgRp[i - 1];
			m_rgRp[i]=Pm_rgRp[i];
		}
		Pm_rgQp[0] = QP2Qstep(m_Qc); //*1.0/CurrentFrameMAD;
		/*frame layer rate control*/
		if(img->BasicUnit==img->Frame_Total_Number_MB)
			Pm_rgRp[0] = img->NumberofTextureBits*1.0/CurrentFrameMAD;
		/*basic unit layer rate control*/
		else
			Pm_rgRp[0]=img->NumberofBasicUnitTextureBits*1.0/CurrentFrameMAD;

		m_rgQp[0]=Pm_rgQp[0];
		m_rgRp[0]=Pm_rgRp[0];
		m_X1=Pm_X1;
    	m_X2=Pm_X2;
	

/*compute the size of window*/


		n_windowSize = (CurrentFrameMAD>PreviousFrameMAD)?(int)(PreviousFrameMAD/CurrentFrameMAD*20)\
			:(int)(CurrentFrameMAD/PreviousFrameMAD*20);
		n_windowSize=MAX(n_windowSize, 1);
		n_windowSize=MIN(n_windowSize,m_Nc);
		n_windowSize=MIN(n_windowSize,m_windowSize+1);
		n_windowSize=MIN(n_windowSize,20);

			/*update the previous window size*/
	m_windowSize=n_windowSize;
	


	for (i = 0; i < 20; i++) {
		m_rgRejected[i] = FALSE;
	}

	// initial RD model estimator
	RCModelEstimator (n_windowSize);
 
	n_windowSize = m_windowSize;
	// remove outlier	
	
	for (i = 0; i < (int) n_windowSize; i++) {
		error[i] = m_X1 / m_rgQp[i] + m_X2 / (m_rgQp[i] * m_rgQp[i]) - m_rgRp[i];
		std += error[i] * error[i];	
	}
	threshold = (n_windowSize == 2) ? 0 : sqrt (std / n_windowSize);
	for (i = 0; i < (int) n_windowSize; i++) {
		if (fabs(error[i]) > threshold)
			m_rgRejected[i] = TRUE;
	}
    // always include the last data point
	m_rgRejected[0] = FALSE;

	// second RD model estimator
	RCModelEstimator (n_windowSize);

	if(MADModelFlag)
		updateMADModel();
	else if(img->type==P_SLICE)
		PPictureMAD[0]=CurrentFrameMAD;
	}	
}

/*Boolean skipThisFrame ()
{
	if (m_B > (int) ((RC_SKIP_MARGIN / 100.0) * m_Bs)) {	// buffer full!
		m_skipNextFrame = TRUE;						// set the status
		m_Nr--;										// skip one frame
		m_B -= m_Rp;								// decrease current buffer level
	} else
		m_skipNextFrame = FALSE;
	return m_skipNextFrame;
}
*/

void RCModelEstimator (int n_windowSize)
{
	int n_realSize = n_windowSize;
	int i;
	double oneSampleQ;
	double a00 = 0.0, a01 = 0.0, a10 = 0.0, a11 = 0.0, b0 = 0.0, b1 = 0.0;
	double MatrixValue;
	Boolean estimateX2 = FALSE;

	for (i = 0; i < n_windowSize; i++) {// find the number of samples which are not rejected
		if (m_rgRejected[i])
			n_realSize--;
	}

	// default RD model estimation results

	m_X1 = m_X2 = 0.0;

	for (i = 0; i < n_windowSize; i++)	{
		if (!m_rgRejected[i])
			oneSampleQ = m_rgQp[i];
	}
	for (i = 0; i < n_windowSize; i++)	{// if all non-rejected Q are the same, take 1st order model
		if ((m_rgQp[i] != oneSampleQ) && !m_rgRejected[i])
			estimateX2 = TRUE;
		if (!m_rgRejected[i])
			m_X1 += (m_rgQp[i] * m_rgRp[i]) / n_realSize;
	}

	// take 2nd order model to estimate X1 and X2
	if ((n_realSize >= 1) && estimateX2) {
			for (i = 0; i < n_windowSize; i++) {
			if (!m_rgRejected[i]) {
				a00 = a00 + 1.0;
				a01 += 1.0 / m_rgQp[i];
				a10 = a01;
				a11 += 1.0 / (m_rgQp[i] * m_rgQp[i]);
				b0 += m_rgQp[i] * m_rgRp[i];
				b1 += m_rgRp[i];
			}
		}
		// solve the equation of AX = B
			MatrixValue=a00*a11-a01*a10;
			if(fabs(MatrixValue)>0.000001)
			{
				m_X1=(b0*a11-b1*a01)/MatrixValue;
				m_X2=(b1*a00-b0*a10)/MatrixValue;
			}
			else
			{
				m_X1=b0/a00;
				m_X2=0.0;
			}
	
	}
	if(img->type==P_SLICE)
	{
		Pm_X1=m_X1;
		Pm_X2=m_X2;
	}
}

double ComputeFrameMAD()
{
	double TotalMAD;
	int i;
	TotalMAD=0.0;
//	CurrentFrameMAD=0.0;
	for(i=0;i<img->Frame_Total_Number_MB;i++)
		TotalMAD +=img->MADofMB[i];
	TotalMAD /=img->Frame_Total_Number_MB;
	return TotalMAD;
}


//update the parameters of linear prediction model
void updateMADModel ()
{
	
	int n_windowSize;
	int i;
	double error[20], std = 0.0, threshold;
	int m_Nc;
	
	if(img->NumberofCodedPFrame>0)
	{
		
		if(img->type==P_SLICE)
		{
			/*frame layer rate control*/
			if(img->BasicUnit==img->Frame_Total_Number_MB)
				m_Nc=img->NumberofCodedPFrame;
			/*basic unit layer rate control*/
			else
				m_Nc=img->NumberofCodedPFrame*TotalNumberofBasicUnit+CodedBasicUnit;
			for (i = 19; i > 0; i--) {// update the history
				PPictureMAD[i] = PPictureMAD[i - 1];
				PictureMAD[i]=PPictureMAD[i];
				ReferenceMAD[i]= ReferenceMAD[i-1];
			}
			PPictureMAD[0] = CurrentFrameMAD;
			PictureMAD[0]=PPictureMAD[0];
			if(img->BasicUnit==img->Frame_Total_Number_MB)
				ReferenceMAD[0]=PictureMAD[1];
			else
			{
				if(((input->InterlaceCodingOption==1)||(input->InterlaceCodingOption==3))\
					&&(img->FieldControl==1))
					ReferenceMAD[0]=FCBUPFMAD[TotalNumberofBasicUnit-1-NumberofBasicUnit];
				else
					ReferenceMAD[0]=BUPFMAD[TotalNumberofBasicUnit-1-NumberofBasicUnit];
			}
			MADPictureC1=PMADPictureC1;
			MADPictureC2=PMADPictureC2;
		}
		
		
		/*compute the size of window*/
		
		n_windowSize = (CurrentFrameMAD>PreviousFrameMAD)?(int)(PreviousFrameMAD/CurrentFrameMAD*20)\
			:(int)(CurrentFrameMAD/PreviousFrameMAD*20);
		n_windowSize=MIN(n_windowSize,(m_Nc-1));
		n_windowSize=MAX(n_windowSize, 1);
		n_windowSize=MIN(n_windowSize,MADm_windowSize+1);
		n_windowSize=MIN(20,n_windowSize);
		/*update the previous window size*/
		MADm_windowSize=n_windowSize;
		
		for (i = 0; i < 20; i++) {
			PictureRejected[i] = FALSE;
		}
		//update the MAD for the previous frame
		if(img->type==P_SLICE)
			PreviousFrameMAD=CurrentFrameMAD;
		
		// initial MAD model estimator
		MADModelEstimator (n_windowSize);
		
		// remove outlier	
		
		for (i = 0; i < (int) n_windowSize; i++) {
			error[i] = MADPictureC1*ReferenceMAD[i]+MADPictureC2-PictureMAD[i];
			std += error[i] * error[i];	
		}
		threshold = (n_windowSize == 2) ? 0 : sqrt (std / n_windowSize);
		for (i = 0; i < (int) n_windowSize; i++) {
			if (fabs(error[i]) > threshold)
				PictureRejected[i] = TRUE;
		}
    // always include the last data point
		PictureRejected[0] = FALSE;
		
		// second MAD model estimator
		MADModelEstimator (n_windowSize);
	}
}

void MADModelEstimator (int n_windowSize)
{
	int n_realSize = n_windowSize;
	int i;
	double oneSampleQ;
	double a00 = 0.0, a01 = 0.0, a10 = 0.0, a11 = 0.0, b0 = 0.0, b1 = 0.0;
	double MatrixValue;
	Boolean estimateX2 = FALSE;
	
	for (i = 0; i < n_windowSize; i++) {// find the number of samples which are not rejected
		if (PictureRejected[i])
			n_realSize--;
	}
	
	// default MAD model estimation results
	
	MADPictureC1 = MADPictureC2 = 0.0;
	
	for (i = 0; i < n_windowSize; i++)	{
		if (!PictureRejected[i])
			oneSampleQ = PictureMAD[i];
	}
	for (i = 0; i < n_windowSize; i++)	{// if all non-rejected MAD are the same, take 1st order model
		if ((PictureMAD[i] != oneSampleQ) && !PictureRejected[i])
			estimateX2 = TRUE;
		if (!PictureRejected[i])
			MADPictureC1 += PictureMAD[i] / (ReferenceMAD[i]*n_realSize);
	}
	
	// take 2nd order model to estimate X1 and X2
	if ((n_realSize >= 1) && estimateX2) {
		for (i = 0; i < n_windowSize; i++) {
			if (!PictureRejected[i]) {
				a00 = a00 + 1.0;
				a01 += ReferenceMAD[i];
				a10 = a01;
				a11 += ReferenceMAD[i]*ReferenceMAD[i];
				b0 += PictureMAD[i];
				b1 += PictureMAD[i]*ReferenceMAD[i];
			}
		}
		// solve the equation of AX = B
		MatrixValue=a00*a11-a01*a10;
		if(fabs(MatrixValue)>0.000001)
		{
			MADPictureC2=(b0*a11-b1*a01)/MatrixValue;
			MADPictureC1=(b1*a00-b0*a10)/MatrixValue;
		}
		else
		{
			MADPictureC1=b0/a01;
			MADPictureC2=0.0;
		}
		
	}
	if(img->type==P_SLICE)
	{
		PMADPictureC1=MADPictureC1;
		PMADPictureC2=MADPictureC2;
	}
}


double QP2Qstep( int QP )
{
	int i; 
	double Qstep;
	static const double QP2QSTEP[6] = { 0.625, 0.6875, 0.8125, 0.875, 1.0, 1.125 };
	
	Qstep = QP2QSTEP[QP % 6];
	for( i=0; i<(QP/6); i++)
		Qstep *= 2;
	
	return Qstep;
}

int Qstep2QP( double Qstep )
{
	int q_per = 0, q_rem = 0;
	
	//	assert( Qstep >= QP2Qstep(0) && Qstep <= QP2Qstep(51) );
	if( Qstep < QP2Qstep(0))
		return 0;
	else if (Qstep > QP2Qstep(51) )
		return 51;
	
	while( Qstep > QP2Qstep(5) )
	{
		Qstep /= 2;
		q_per += 1;
	}
	
	if (Qstep <= (0.625+0.6875)/2) 
	{
		Qstep = 0.625;
		q_rem = 0;
	}
	else if (Qstep <= (0.6875+0.8125)/2)
	{
		Qstep = 0.6875;
		q_rem = 1;
	}
	else if (Qstep <= (0.8125+0.875)/2)
	{
		Qstep = 0.8125;
		q_rem = 2;
	}
	else if (Qstep <= (0.875+1.0)/2)
	{
		Qstep = 0.875;
		q_rem = 3;
	}
	else if (Qstep <= (1.0+1.125)/2)
	{
		Qstep = 1.0;	
		q_rem = 4;
	}
	else 
	{
		Qstep = 1.125;
		q_rem = 5;
	}
	
	return (q_per * 6 + q_rem);
}