* \brief
 *
 * Step 2: Target Bit Calculation (Before Encoding the Current Frame)
 *
 ***************************************************************************
 */

void Target_Bit_Calculation(RCQ2* rc,			/*<--> rc struct */
							int frametype		/*<--  I or P frame type */
							)
{
	int temp,temp1;/*,temp2,temp3;*/
	int Nrp=0,Nri=0;/*number of I,P,B*/

		Nri=rc->Nr/rc->max_key_interval;
		Nrp=rc->Nr-Nri;

/*	if(frametype==0&&(rc->Nr==rc->max_key_interval)) */
		/*whq,2002.12.22*/
	if(frametype==0)
		temp1=(int)((float)rc->Rr/(float)(1+(3*Nrp/8)));
		/*設(shè)置I，P之間的位數(shù)關(guān)系，可以在此直接改變*/
	else 
		temp1=rc->Rr/(Nrp);
	
	rc->T=(temp1*95+rc->data[frametype].S*5)/100;

	temp=(2*rc->Bs-rc->B)*1000/(rc->B+rc->Bs);
	rc->T=(int)(rc->T*temp)/1000;	/* increase if less than half */
						/* decrease if more than half, don't change if half */
	if (rc->B+rc->T > rc->Bs*9/10)
		rc->T = (int)max(rc->Rs/30,9*rc->Bs/10-rc->B);	/* to avoid overflow*/
	else if (rc->B-rc->Rp+rc->T < rc->Bs/10)
		rc->T = (int)(rc->Rp-rc->B+rc->Bs/10);		/* to avoid underflow*/

}

/*!
 ****************************************************************************
 * \brief
 *
 * Step 3: Quantisation Level Calculation  (Before Encoding the Current Frame)
 *
 ***************************************************************************
 */

void Quantisation_Calculation(FILE *pfile,		/*<--  file pointer,output stat*/
							  RCQ2* rc,			/*<--> rc struct */
							  int frametype		/*<--  I or P frame type */
							  )
{	

	int delt=0;
	rc->T= max(rc->Rp/3+rc->data[frametype].Hp, rc->T);

#ifdef stat_output
	fprintf(pfile,"before  rc->T=%d \n",rc->T );
#endif

	delt=(int)rc->data[frametype].X1*(int)rc->data[frametype].X1+4*(int)rc->data[frametype].X2*(rc->T-rc->data[frametype].Hp)/rc->data[frametype].Ec;

	if ((rc->data[frametype].X2==0.0) || (delt<0))/* fall back to one order model */
	{
		rc->data[frametype].Qc = (int)(rc->data[frametype].X11*(double)rc->data[frametype].Ec/(double)(rc->T-rc->data[frametype].Hp));		/* fall back 1st order mode */
	}
	else /* quadratic rate distortion model  */
	{
		rc->data[frametype].Qc=(int)((2*rc->data[frametype].X2)/((int)sqrt((double)delt)-(int)rc->data[frametype].X1));
	} 
	if(frametype==0) /* clip */
	{
		rc->data[frametype].Qc =(int) min(ceil(rc->data[frametype].Ql*1.25), min(rc->data[frametype].Qc, 31)); /* clipping   22*/
		rc->data[frametype].Qc =(int) max(ceil(rc->data[frametype].Ql*0.75), max(rc->data[frametype].Qc, 1));  /* clipping   8 */

	}
	else if(frametype==1)
	{
		rc->data[frametype].Qc =(int) min(ceil(rc->data[frametype].Ql*1.25), min(rc->data[frametype].Qc, 31)); /* clipping  28*/
		rc->data[frametype].Qc =(int) max(ceil(rc->data[frametype].Ql*0.75), max(rc->data[frametype].Qc, 1));  /* clipping  8*/
	}

#ifdef stat_output	
	fprintf(pfile," %dth frametype=%d ,rc->Qc=%d \n",rc->N_Total,frametype,rc->data[frametype].Qc );
#endif

}

/*!
 ****************************************************************************
 * \brief
 *
 * Step 4: After encoding the current frame,update rate control param and 
 *         model after coding one frame
 *
 ***************************************************************************
 */
void After_encoding_current(FILE *pfile,		/*<--  file pointer,output stat*/
							RCQ2* rc,			/*<--> rc struct */
							int frametype,		/*<--  I or P frame type */
							int length,			/*<--  frame length */
							int head_vector,	/*<--  motion vector length */
							int mad				/*<--  frame mad */
							)
{	
	rc->N_Total++; 
	rc->data[frametype].Rc=length;
	rc->data[frametype].S=rc->data[frametype].Rc;
	rc->data[frametype].Hc=head_vector;
	rc->data[frametype].Hp=rc->data[frametype].Hc;	/* update the previous header and motion bit counts */
	rc->data[frametype].Ep=rc->data[frametype].Ec;
	rc->data[frametype].Ec=mad;
	rc->data[frametype].Ql= rc->data[frametype].Qc;
	rc->Nr--;	
	if(frametype==0)
	{
		rc->Rf=length;
		rc->Rr-=rc->Rf;                            /* total number of bits available for this segment */
		rc->B +=rc->Rf-rc->Rp;	
	}
	else
	{
		rc->B += rc->data[frametype].Rc - rc->Rp;	/* update buffer fullness */
		rc->Rr -= rc->data[frametype].Rc;			/* update the remaining bit counts */

	}


	rch_store_after(rc->rc_queue,frametype,rc->data[frametype].Rc,rc->data[frametype].Qc,rc->data[frametype].Hc,rc->data[frametype].Ec);
	UpdateRDModel(rc,frametype);

	
}


/*!
 ****************************************************************************
 * \brief
 *
 * update the model after coding one frame 
 *
 ***************************************************************************
 */

void UpdateRDModel (RCQ2* rc,			/*<--> rc struct */
					int frametype		/*<--  I or P frame type */
					)
{
/*	
	n: number of encoded past frames;
	x: matrix contains Q;
	y: matrix contains Q*(R-H)/E;
	Ep:	mean absolute difference for the previous frame. This is computed after motion compensation for the Y component only.
*/
	int w,n; 
	n=rch_get_size(rc->rc_queue,frametype);
	w = min(n, 20);	/* Maximum data number set to 20 */
	if(w>1){
	
		if(rc->data[frametype].Ep>rc->data[frametype].Ec)
			w = (int)ceil((double)rc->data[frametype].Ec/(double)rc->data[frametype].Ep*(double)w);		/* sliding window for scene change */

		else
			w = (int)ceil((double)rc->data[frametype].Ep/(double)rc->data[frametype].Ec*(double)w);		/* sliding window for scene change */
	}

	rc->data[frametype].Ep = rc->data[frametype].Ec;				/* update mad */
	
	rch_get_data(rc->rc_queue,frametype,w,rc->data[frametype].QpArray,rc->data[frametype].RpArray);
	Estimator(rc,frametype,&w);
	RemoveOutlier(rc,frametype,rc->data[frametype].X1,rc->data[frametype].X2,
		rc->data[frametype].QpArray,rc->data[frametype].RpArray,rc->data[frametype].Ec,&w);
	Estimator(rc,frametype,&w);
}


/*!
 ****************************************************************************
 * \brief
 *
 * compute quadratic rate distortion model param X1,X2
 *
 ***************************************************************************
 */

void Estimator(RCQ2* rc,			/*<--> rc struct */
			   int frametype,		/*<--  I or P frame type */
			   int * w				/*<--  number of points*/
			   )
{


	double a,b,c,x,y;
	int i;
	int  temp_cmp;
	BOOL cmp_result;
	
	a=b=c=x=y=0.0;
	
	for(i=0;i<*w;i++){
		rc->array_x[i][0]=1.0;
		rc->array_x[i][1]=1.0/(double)rc->data[frametype].QpArray[i];	

		rc->array_y[i]=(double)rc->data[frametype].QpArray[i]*rc->data[frametype].RpArray[i];
	}
	
	cmp_result=FALSE;
	temp_cmp=rc->data[frametype].QpArray[0];
	rc->data[frametype].X11=0.0;
	for (i=0; i<*w; i++)
	{   
		rc->data[frametype].X11 += rc->array_y[i]/(*w);/*可以最后做個(gè)除法，提到外面去*/
		
		if((!cmp_result)&&(temp_cmp!=rc->data[frametype].QpArray[i]))
			cmp_result=TRUE;
		else
			temp_cmp=rc->data[frametype].QpArray[i];
	}

	rc->data[frametype].X11=(int)(rc->data[frametype].X11+0.5) ;
	if (cmp_result) 
	{
		
		/* Dimension of the matrix */
		/* 2x1 = (2xw * wx2)**(-1)*(2xw)*(w*1) */
		
		for (i=0;i<*w;i++){
			a+=rc->array_x[i][0]*rc->array_x[i][0];
			b+=rc->array_x[i][0]*rc->array_x[i][1];
			c+=rc->array_x[i][1]*rc->array_x[i][1];
			x+=rc->array_x[i][0]*rc->array_y[i];
			y+=rc->array_x[i][1]*rc->array_y[i];
		}
		rc->data[frametype].X1=(c*x-b*y)/(a*c-b*b);
		rc->data[frametype].X2=(a*y-b*x)/(a*c-b*b);
	} 
}

/*!
 ****************************************************************************
 * \brief
 *
 * remove the bad point out of the set
 *
 ***************************************************************************
 */

void RemoveOutlier(RCQ2* rc,	/*<-- rc struct */
				   int frametype,/*<-- I or P frame type */
				   double X1,	/*<-- quadratic model param X1 */
				   double X2,	/*<-- quadratic model param X1 */
				   int Qp[],	/*<-- previous frame Qp */
				   double Rp[],	/*<-- previous frame complexity scale */
				   double Ec,	/*<-- current frame complexity */
				   int * w		/*<-->in: number of total points/out:number of usable points */
				   )
{
	
	int new_w, i;
	double threshold=0.0;
	double std=0.0;

	for (i=0; i<*w; i++) {
		std += (X1*Ec/Qp[i]+X2*Ec/(Qp[i]*Qp[i])-Rp[i]*Ec)*
			(X1*Ec/Qp[i]+X2*Ec/(Qp[i]*Qp[i])-Rp[i]*Ec);
		rc->rc_error[i] = X1*Ec/Qp[i]+X2*Ec/(Qp[i]*Qp[i])-Rp[i]*Ec;

	}
	
	/* When w is 2 the threshold should be 0, some implementation may resolve a small number due to precision problem. */
	threshold = sqrt(std/(*w));			/* Setup rejection threshold */
	
	new_w = 0;

	for (i=0; i<*w; i++){
		if(rc->rc_error[i]<0.0)
			rc->rc_error[i]=-rc->rc_error[i];
		if(rc->rc_error[i]<= threshold){
			rc->data[frametype].QpArray[new_w] = rc->data[frametype].QpArray[i];
			rc->data[frametype].RpArray[new_w] = rc->data[frametype].RpArray[i];
			new_w++;
		}
	}

	*w = new_w;
}

/*!
 ****************************************************************************
 * \brief
 *
 *  decide whether skip the frame according to the buffer fullness
 *	If the buffer has reached 80% of the buffer size, the encoder will skip the 
 *  upcoming frame to be encoded. Thus the buffer is reduced to prevent from the 
 *  buffer overflow. 
 ***************************************************************************
 */

BOOL IsSkipFrame(				/*==> return flag showing whether skip the frame */
				 RCQ2* rc       /*<-- rc struct */
				 )
{ 

	if(rc->B > 0.8 * rc->Bs){

		rc->Nr--;
		rc->N_Total++; 
		rc->B -= rc->Rp;
		return TRUE;
	}
	else
		return FALSE;
}