/*
***********************************************************************
* COPYRIGHT AND WARRANTY INFORMATION
*
* Copyright 2001, International Telecommunications Union, Geneva
*
* DISCLAIMER OF WARRANTY
*
* These software programs are available to the user without any
* license fee or royalty on an "as is" basis. The ITU disclaims
* any and all warranties, whether express, implied, or
* statutory, including any implied warranties of merchantability
* or of fitness for a particular purpose.  In no event shall the
* contributor or the ITU be liable for any incidental, punitive, or
* consequential damages of any kind whatsoever arising from the
* use of these programs.
*
* This disclaimer of warranty extends to the user of these programs
* and user's customers, employees, agents, transferees, successors,
* and assigns.
*
* The ITU does not represent or warrant that the programs furnished
* hereunder are free of infringement of any third-party patents.
* Commercial implementations of ITU-T Recommendations, including
* shareware, may be subject to royalty fees to patent holders.
* Information regarding the ITU-T patent policy is available from
* the ITU Web site at http://www.itu.int.
*
* THIS IS NOT A GRANT OF PATENT RIGHTS - SEE THE ITU-T PATENT POLICY.
************************************************************************
*/

/*!
 ***************************************************************************
 * \file ratectl.c
 *
 * \brief
 *    Rate Control algorithm
 *
 * \major contributors:
 *  Siwei Ma <swma@jdl.ac.cn>
 *  Zhengguo LI<ezgli@lit.a-star.edu.sg>
 *
 * \date
 *   16 Jan. 2003
 **************************************************************************
 */
#include <stdio.h>
#include <math.h>
#include "global.h"
#include "ratectl.h"


const double THETA=1.3636;//LIZG


/* rate control variables */
int Xp, Xb;
static int R,T_field;
static int Np, Nb, bits_topfield, Q;
long T,T1;

double Wp,Wb; //LIZG, JVT019r1, 25/10/2002
int TotalPFrame;//LIZG, JVT019r1, 25/10/2002
int DuantQp; //LIZG, JVT019r1
int PDuantQp;//LIZG,JVT019r1
FILE *BitRate;//LIZG,JVT019r1
double DeltaP;//LIZG, JVT019r1



// Initiate rate control parameters
void rc_init_seq()
{

 //LIZG, JVT019r1, our proposed adaptive rate control


	 int i;
	 Xp=0;
	 Xb=0;
	 
	 bit_rate=input->bit_rate;
	 frame_rate = (float)(img->framerate *(input->successive_Bframe + 1)) / (float) (input->jumpd + 1);
	 PreviousBit_Rate=bit_rate;

/*compute the total number of MBs in a frame*/

	 img->Frame_Total_Number_MB=img->height*img->width/256;
	 if(input->basicunit>img->Frame_Total_Number_MB)
		 input->basicunit=img->Frame_Total_Number_MB;
	 if(input->basicunit<img->Frame_Total_Number_MB)
		 TotalNumberofBasicUnit=img->Frame_Total_Number_MB/input->basicunit;

	
	 if(input->basicunit<img->Frame_Total_Number_MB)
		 MINVALUE=8.0;
	 else
		 MINVALUE=4.0;


/*LIZG, JVT019r1, initialize the parameters of fluid flow traffic model*/

	 BufferSize=bit_rate*2;
	 CurrentBufferFullness=BufferSize/8;
	 GOPTargetBufferLevel=CurrentBufferFullness;

      
/*LIZG, JVT019r1, initialize the previous window size*/
	 m_windowSize=0;
	 MADm_windowSize=0;
	 img->NumberofCodedBFrame=0;
	 img->NumberofCodedPFrame=0;
	 img->NumberofGOP=0;
 /*LIZG, JVT019r1, remaining # of bits in GOP */
	 R = 0;
/*set control parameter for buffer regulation */
	 if(input->successive_Bframe>0)
	 {
		 GAMMAP=0.25;
		 BETAP=0.9;
	 }
	 else
	 {
		 GAMMAP=0.75;
		 BETAP=0.5;
	 }
  
 /*LIZG, JVT019r1, quadratic rate-distortion model*/
	 PPreHeader=0;
	
	 Pm_X1=bit_rate*1.0;
	 Pm_X2=0.0;
 /* linear prediction model for P picture*/
	 PMADPictureC1=1.0;
	 PMADPictureC2=0.0;

	 for(i=0;i<20;i++)
	 {
		 Pm_rgQp[i]=0;
		 Pm_rgRp[i]=0.0;
		 PPictureMAD[i]=0.0;
	 }
	 PPictureMAD[20]=0.0;

	
//Define the largest variation of quantization parameters

	 
	 PDuantQp=2;
	 
	
	
/*basic unit layer rate control*/
	  PAveHeaderBits1=0;
	  PAveHeaderBits3=0;	
	  if(TotalNumberofBasicUnit>=9)
		  DDquant=1;
	  else
		  DDquant=2;

	  MBPerRow=input->img_width/16;

/*adaptive field/frame coding*/
	  img->FieldControl=0;
	  img->UpdateTimes=0;
  
}

// Initiate one GOP
void rc_init_GOP(int np, int nb)
{
	int OverBits;
	int OverDuantQp;
	int AllocatedBits;
	int GOPDquant;
/*check if the last GOP over uses its budget. If yes, the initial QP of the I frame in 
 the coming  GOP will be increased.*/
/*recent change*/
	OverBits=-R;

 /*compute the total number of bits for the current GOP*/ 
	AllocatedBits = (int) floor((1 + np + nb) * bit_rate / frame_rate + 0.5);
	R +=AllocatedBits;
	Np  = np;
	Nb  = nb;
/*recent change*/

	OverDuantQp=(int)(8*OverBits/AllocatedBits+0.5);


	GOPOverdue=FALSE;
  
/*field coding*/
	img->IFLAG=1;

/*LIZG, JVT019r1, Compute the starting quantization parameter for each GOP.
	The quantization parameter is adaptive to the available channel bandwidth and the 
	GOP length*/

	
	TotalPFrame=np;
	img->NumberofGOP++;
	if(img->NumberofGOP==1)
	{
		InitialQp=input->initialQP;
		PreviousQp2=InitialQp;
	}
	else
	{
/*adaptive field/frame coding*/
		if((input->InterlaceCodingOption==1)\
			||(input->InterlaceCodingOption==3))
		{
			if(img->FieldFrame==1)
				img->TotalQpforPPicture +=FrameQPBuffer;
			else
				img->TotalQpforPPicture +=FieldQPBuffer;
		}
/*compute the average QP of P frames in the previous GOP*/

		if(img->SwitchFrameRate>30)
			PAverageQp=(int)(1.0*img->TotalQpforPPicture/img->NumberofPPicture+0.5);//recent change -1
		else
			PAverageQp=(int)(1.0*img->TotalQpforPPicture/img->NumberofPPicture+0.5)-1;//recent change -1
		GOPDquant = (int)(0.5+1.0*(np+nb+1)/15);
		//PAverageQp -=(int)(0.5+1.0*(np+nb+1)/15);
		if(GOPDquant>2)
			GOPDquant=2;
		PAverageQp-=GOPDquant;
		PAverageQp +=OverDuantQp;
		InitialQp=PAverageQp;
		Pm_Qp=PAverageQp;
		PAveFrameQP=PAverageQp;
		PreviousQp1=PreviousQp2;
		PreviousQp2=InitialQp;
		
	}

	img->TotalQpforPPicture=0;
	img->NumberofPPicture=0;
	NumberofBFrames=0;

 
}

void rc_init_pict(int fieldpic,int topfield,int targetcomputation)
{
	int i;
/*compute the total number of basic units in a frame*/

	if(input->InterlaceCodingOption==3)
		TotalNumberofBasicUnit=img->Frame_Total_Number_MB/img->BasicUnit;
	img->NumberofCodedMacroBlocks=0;


/*LIZG, JVT019r1, Normally, the bandwith for the VBR case is estimated by 
a congestion control algorithm. A bandwidth curve can be predefined if we 
only want to test the proposed algorithm*/
	if(input->channel_type==1)
	{
		if(img->NumberofCodedPFrame==58)
			bit_rate *=1.5;
		else if(img->NumberofCodedPFrame==59)
			PreviousBit_Rate=bit_rate;
	}

	
	
/*LIZG, JVT019r1, predefine a target buffer level for each P frame*/
	if((fieldpic||topfield)&&targetcomputation)
	{
		switch (img->type)
		{
			case INTER_IMG:
/*Since the available bandwidth may vary at any time, the total number of 
bits is updated frame  by frame*/
				if(PreviousBit_Rate!=bit_rate)
					R +=(int) floor((bit_rate-PreviousBit_Rate)*(Np+Nb)/frame_rate+0.5);
							
/* predefine the  target buffer level for each picture.
frame layer rate control*/
				if(img->BasicUnit==img->Frame_Total_Number_MB)
				{
					if(img->NumberofPPicture==1)
					{
						TargetBufferLevel=CurrentBufferFullness;
						DeltaP=(CurrentBufferFullness-GOPTargetBufferLevel)/(TotalPFrame-1);
						TargetBufferLevel -=DeltaP;
					}
					else if(img->NumberofPPicture>1)
						TargetBufferLevel -=DeltaP;
				}

/*basic unit layer rate control*/
				else
				{
					if(img->NumberofCodedPFrame>0)
					{
						/*adaptive frame/filed coding*/
						if(((input->InterlaceCodingOption==1)||(input->InterlaceCodingOption==3))\
							&&(img->FieldControl==1))
						{
							for(i=0;i<TotalNumberofBasicUnit;i++)
								FCBUPFMAD[i]=FCBUCFMAD[i];
						}
						else
						{
							for(i=0;i<TotalNumberofBasicUnit;i++)
								BUPFMAD[i]=BUCFMAD[i];
						}
							
					}

					if(img->NumberofGOP==1)
					{
						if(img->NumberofPPicture==1)
						{
							TargetBufferLevel=CurrentBufferFullness;
							DeltaP=(CurrentBufferFullness-GOPTargetBufferLevel)/(TotalPFrame-1);
							TargetBufferLevel -=DeltaP;
						}
						else if(img->NumberofPPicture>1)
							TargetBufferLevel -=DeltaP;
					}
					else if(img->NumberofGOP>1)
					{
						if(img->NumberofPPicture==0)
						{
							TargetBufferLevel=CurrentBufferFullness;
							DeltaP=(CurrentBufferFullness-GOPTargetBufferLevel)/TotalPFrame;
							TargetBufferLevel -=DeltaP;
						}
						else if(img->NumberofPPicture>0)
							TargetBufferLevel -=DeltaP;
					}
//					BitRate=fopen("c:\\mybitrate101.dat","a+");
//                  fprintf(BitRate,"%5.2f\n",TargetBufferLevel);
//                  fclose(BitRate);

				}
					
//compute the average weight of P frame
				if(input->successive_Bframe>0)
				{
					if((img->NumberofCodedPFrame<8)&&(img->NumberofCodedPFrame>1))
						AWp=Wp*(img->NumberofCodedPFrame-1)/img->NumberofCodedPFrame+\
							AWp/img->NumberofCodedPFrame;
					else if(img->NumberofCodedPFrame>1)
						AWp=Wp/8+7*AWp/8;
//compute the target buffer level
					TargetBufferLevel +=(AWp*(input->successive_Bframe+1)*bit_rate\
						/(frame_rate*(AWp+AWb*input->successive_Bframe))-bit_rate/frame_rate);
				}
				break;


			   case B_IMG:
/* update the total number of bits if the bandwidth is changed*/
				   if(PreviousBit_Rate!=bit_rate)
					   R +=(int) floor((bit_rate-PreviousBit_Rate)*(Np+Nb)/frame_rate+0.5);
				   	if((img->NumberofCodedPFrame==1)&&(img->NumberofCodedBFrame==1))
					{
						AWp=Wp;
						AWb=Wb;
					}
					else if(img->NumberofCodedBFrame>1)
					{
//compute the average weight of B frame
						if(img->NumberofCodedBFrame<8)
							AWb=Wb*(img->NumberofCodedBFrame-1)/img->NumberofCodedBFrame+\
								AWb/img->NumberofCodedBFrame;
						else
							AWb=Wb/8+7*AWb/8;
					}

				    break;
		}

		 /*Compute the target bit for each frame*/
		if(img->type==INTER_IMG)
		{
/*frame layer rate control*/
			if(img->BasicUnit==img->Frame_Total_Number_MB)
			{
				if(img->NumberofCodedPFrame>0)
				{
					T = (int) floor(Wp*R/(Np*Wp+Nb*Wb) + 0.5);
					T1 = (int) floor(bit_rate/frame_rate-GAMMAP*(CurrentBufferFullness-TargetBufferLevel)+0.5);
					T = (int)(floor(BETAP*T+(1-BETAP)*T1+0.5));
				}
			 }
 /*basic unit layer rate control*/
			else
			{
				if((img->NumberofGOP==1)&&(img->NumberofCodedPFrame>0))
				{
					T = (int) floor(Wp*R/(Np*Wp+Nb*Wb) + 0.5);
					T1 = (int) floor(bit_rate/frame_rate-GAMMAP*(CurrentBufferFullness-TargetBufferLevel)+0.5);
					T = (int)(floor(BETAP*T+(1-BETAP)*T1+0.5));
				}
				else if(img->NumberofGOP>1)
				{
					T = (int) floor(Wp*R/(Np*Wp+Nb*Wb) + 0.5);
					T1 = (int) floor(bit_rate/frame_rate-GAMMAP*(CurrentBufferFullness-TargetBufferLevel)+0.5);
					T = (int)(floor(BETAP*T+(1-BETAP)*T1+0.5));
				}
			}
			if((topfield)||(fieldpic&&((input->InterlaceCodingOption==1)\
				||(input->InterlaceCodingOption==3))))
				T_field=T;
		}
	}

	if(fieldpic||topfield)
	{
		/*frame layer rate control*/
		img->NumberofHeaderBits=0;
		img->NumberofTextureBits=0;


		/*basic unit layer rate control*/
		if(img->BasicUnit<img->Frame_Total_Number_MB)
		{
			TotalFrameQP=0;
			img->NumberofBasicUnitHeaderBits=0;
			img->NumberofBasicUnitTextureBits=0;
			img->TotalMADBasicUnit=0;
			if(img->FieldControl==0)
				NumberofBasicUnit=TotalNumberofBasicUnit;
			else
				NumberofBasicUnit=TotalNumberofBasicUnit/2;
		}
	}
		
	if((img->type==INTER_IMG)&&(img->BasicUnit<img->Frame_Total_Number_MB)\
		&&(img->FieldControl==1))
	{
	/*top filed at basic unit layer rate control*/
		if(topfield)
		{
			bits_topfield=0;
			T=(long)(T_field*0.6);
		}
  /*bottom filed at basic unit layer rate control*/
		else
		{
			T=T_field-bits_topfield;
			img->NumberofBasicUnitHeaderBits=0;
			img->NumberofBasicUnitTextureBits=0;
			img->TotalMADBasicUnit=0;
			NumberofBasicUnit=TotalNumberofBasicUnit/2;
		}
	}
}



// LIZG, calculate MAD for the current macroblock 
double calc_MAD()
{
	int k,l;
   	int s = 0;
	double MAD;

	for (k = 0; k < 16; k++)
		for (l = 0; l < 16; l++)
			s+= abs(diffy[k][l]);
	
	MAD=s*1.0/256;
	return MAD;
}

// update the buffer fullness and the number of remaining bits after frame/field encoding
void rc_update_pict(int nbits)
{
	R-= nbits; /* remaining # of bits in GOP */
//	BitRate=fopen("c:\\mybitrate103.dat","a+");
//   fprintf(BitRate,"%d\n",nbits);
//   fclose(BitRate);
/*update the buffer fullness*/
	CurrentBufferFullness += nbits - bit_rate/frame_rate;
//	BitRate=fopen("c:\\mybitrate102.dat","a+");
//  fprintf(BitRate,"%5.2f\n",CurrentBufferFullness);
//  fclose(BitRate);

	return;
}

// update after frame encoding
void rc_update_pict_frame(int nbits)
{

/*LIZG, our proposed adaptive rate control, update the
complexity weight of I, P, B frame*/
	int X;
	int Avem_Qc;

		
/*frame layer rate control*/
	if(img->BasicUnit==img->Frame_Total_Number_MB)
		X = (int) floor(nbits*m_Qc+ 0.5);
/*basic unit layer rate control*/
	else
	{
		if(img->type==INTER_IMG)
		{
			if(((img->IFLAG==0)&&(img->FieldControl==1))\
				||(img->FieldControl==0))
			{
				Avem_Qc=TotalFrameQP/TotalNumberofBasicUnit;
				X=(int)floor(nbits*Avem_Qc+0.5);
			}
		}
		else if(img->type==B_IMG)
			X = (int) floor(nbits*m_Qc+ 0.5);
	}

	switch (img->type)
	{
		case INTER_IMG:
 /*filed coding, need to consider the case that the bottom filed of an I frame is a P-type field*/
			if(((img->IFLAG==0)&&(img->FieldControl==1))\
				||(img->FieldControl==0))
			{
				Xp = X;
				Np--;
				Wp=Xp;
				Pm_Hp=img->NumberofHeaderBits;
				img->NumberofCodedPFrame++;
				img->NumberofPPicture++;
			}
			else if((img->IFLAG!=0)&&(img->FieldControl==1))
				img->IFLAG=0;
		    break;
	      case B_IMG:
		    Xb = X;
			Nb--;
		    Wb=Xb/THETA; 
		
			img->NumberofCodedBFrame++;

			NumberofBFrames++;

	
		    break;
	  }

}




//set coded bits for top field
void setbitscount(int nbits)
{
	bits_topfield = nbits;
}

//compute a  quantization parameter for each frame
int updateQuantizationParameter(int topfield)
{
	double dtmp;
	int m_Bits;
	int BFrameNumber;
	int StepSize;
	int PAverageQP;
	int SumofBasicUnit;

	
/*frame layer rate control*/
	if(img->BasicUnit==img->Frame_Total_Number_MB)
	{

/*fixed quantization parameter is used to coded I frame, the first P frame and the first B frame
	the quantization parameter is adjusted according the available channel bandwidth and 
	the type of vide*/	
/*top field. The quantization parameter of top field is the same as that of the bottom field*/
		if((topfield)||(img->FieldControl==0))
		{
			if(img->type==INTRA_IMG)
			{
				m_Qc=InitialQp;
				return m_Qc;
			}
			else if(img->type==B_IMG)
			{
				if(input->successive_Bframe==1)
				{
					BFrameNumber=(NumberofBFrames+1)%input->successive_Bframe;
					if(BFrameNumber==0)
						BFrameNumber=input->successive_Bframe;
/*adaptive field/frame coding*/
					else if(BFrameNumber==1)
					{
						if((input->InterlaceCodingOption==1)\
							||(input->InterlaceCodingOption==3))
						{
							if(img->FieldControl==0)
							{
																		
/*previous choice is frame coding*/
								if(img->FieldFrame==1)
								{
									PreviousQp1=PreviousQp2;
									PreviousQp2=FrameQPBuffer;
								}

							
/*previous choice is field coding*/
								else
								{
									PreviousQp1=PreviousQp2;
									PreviousQp2=FieldQPBuffer;
								}
							}
						}
					}
					if(PreviousQp1==PreviousQp2)
						m_Qc=PreviousQp1+2;
					else
						m_Qc=(PreviousQp1+PreviousQp2)/2+1;
					m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
					m_Qc = MAX(RC_MIN_QUANT, m_Qc);//clipping
				}
				else
				{
					BFrameNumber=(NumberofBFrames+1)%input->successive_Bframe;
					if(BFrameNumber==0)
						BFrameNumber=input->successive_Bframe;
/*adaptive field/frame coding*/
					else if(BFrameNumber==1)
					{
						if((input->InterlaceCodingOption==1)\
							||(input->InterlaceCodingOption==3))
						{
							if(img->FieldControl==0)
							{
/*previous choice is frame coding*/
								if(img->FieldFrame==1)
								{
									PreviousQp1=PreviousQp2;
									PreviousQp2=FrameQPBuffer;
								}
							
				
/*previous choice is field coding*/
								else
								{
									PreviousQp1=PreviousQp2;
									PreviousQp2=FieldQPBuffer;
								}
							}
						}
					}
					if((PreviousQp2-PreviousQp1)<=(-2*input->successive_Bframe-3))
						StepSize=-3;
					else  if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe-2))
						StepSize=-2;
					else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe-1))
						StepSize=-1;
					else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe))
						StepSize=0;
					else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe+1))
						StepSize=1;
					else
						StepSize=2;

					m_Qc=PreviousQp1+StepSize;
					m_Qc +=MIN(2*(BFrameNumber-1),MAX(-2*(BFrameNumber-1), \
							(BFrameNumber-1)*(PreviousQp2-PreviousQp1)/(input->successive_Bframe-1)));
					m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
					m_Qc = MAX(RC_MIN_QUANT, m_Qc);//clipping
				}
				return m_Qc;
			}
		
		else if((img->type==INTER_IMG)&&(img->NumberofPPicture==0))
		{
			
			if(img->NumberofGOP==1)
				m_Qc=InitialQp; 
			else 
				m_Qc=PAverageQp;
    		if(img->FieldControl==0)
			{
				if(input->InterlaceCodingOption==0)
				{
					img->TotalQpforPPicture +=m_Qc;
					PreviousQp1=PreviousQp2;
					PreviousQp2=m_Qc;
					Pm_Qp=m_Qc;
				}
				/*adaptive field/frame coding*/
				else
					FrameQPBuffer=m_Qc;
			}
				
			return m_Qc;
			
		}
		else
		{
			/*adaptive field/frame coding*/
			if(((input->InterlaceCodingOption==1)\
				||(input->InterlaceCodingOption==3))\
				&&(img->FieldControl==0))
			{
				/*previous choice is frame coding*/
				if(img->FieldFrame==1)
				{
					img->TotalQpforPPicture +=FrameQPBuffer;
					Pm_Qp=FrameQPBuffer;
				}
				/*previous choice is field coding*/
				else
				{
					img->TotalQpforPPicture +=FieldQPBuffer;
					Pm_Qp=FieldQPBuffer;
				}
			}
			m_X1=Pm_X1;
			m_X2=Pm_X2;
			m_Hp=PPreHeader;
			m_Qp=Pm_Qp;
			DuantQp=PDuantQp;
			MADPictureC1=PMADPictureC1;
			MADPictureC2=PMADPictureC2;
			PreviousPictureMAD=PPictureMAD[0];
                
/*predict the MAD of current frame*/
			CurrentFrameMAD=MADPictureC1*PreviousPictureMAD+MADPictureC2;
/*compute the number of bits for the texture*/			
		
		if(T<0)
		{
			m_Qc=m_Qp+DuantQp;
			m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
		}
	   else
		{
		   m_Bits =T-m_Hp;
		   m_Bits = MAX(m_Bits, (int)(bit_rate/(MINVALUE*frame_rate)));
		   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_Qc = (int) (m_X1 * CurrentFrameMAD / (double) m_Bits);
	       else // 2nd order mode
			   m_Qc = (int) ((2 * m_X2 * CurrentFrameMAD) / (sqrt (dtmp) - m_X1 * CurrentFrameMAD));
		   m_Qc = MIN(m_Qp+DuantQp,  m_Qc);	// control variation
		   m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
		   m_Qc = MAX(m_Qp-DuantQp, m_Qc);	// control variation
		   m_Qc = MAX(RC_MIN_QUANT, m_Qc);
	   }
	    
		if(img->FieldControl==0)
		{
			/*frame coding*/
			if(input->InterlaceCodingOption==0)
			{
				img->TotalQpforPPicture +=m_Qc;
				PreviousQp1=PreviousQp2;
				PreviousQp2=m_Qc;
				Pm_Qp=m_Qc;
			}
			/*adaptive field/frame coding*/
			else
				FrameQPBuffer=m_Qc;
		}
							
		
		return m_Qc;
		}
	 }
	 /*bottom field*/
	 else
	 {
		 if((img->type==INTER_IMG)&&(img->IFLAG==0))
		 {
			 /*field coding*/
			 if(input->InterlaceCodingOption==2)
			 {
				 img->TotalQpforPPicture +=m_Qc;
				 PreviousQp1=PreviousQp2+1; 
				 PreviousQp2=m_Qc;//
				 Pm_Qp=m_Qc;
			 }
			 /*adaptive field/frame coding*/
			 else
				 FieldQPBuffer=m_Qc;
			 
		 }
		 return m_Qc;
	 }
	}
	/*basic unit layer rate control*/
	else
	{
		/*top filed of I frame*/

		if(img->type==INTRA_IMG)
		{
			m_Qc=InitialQp;
			return m_Qc;
		}
		/*bottom field of I frame*/
		else if((img->type==INTER_IMG)&&(img->IFLAG==1)&&(img->FieldControl==1))
		{
			m_Qc=InitialQp;
			return m_Qc;
		}
		else if(img->type==B_IMG)
		{
			/*top filed of B frame*/
			if((topfield)||(img->FieldControl==0))
			{
				if(input->successive_Bframe==1)
				{
					BFrameNumber=(NumberofBFrames+1)%input->successive_Bframe;
				if(BFrameNumber==0)
					BFrameNumber=input->successive_Bframe;
				/*adaptive field/frame coding*/
				else if(BFrameNumber==1)
				{
					if((input->InterlaceCodingOption==1)\
						||(input->InterlaceCodingOption==3))
					{
						if(img->FieldControl==0)
						{
											
				/*previous choice is frame coding*/
							if(img->FieldFrame==1)
							{
								PreviousQp1=PreviousQp2;
								PreviousQp2=FrameQPBuffer;
							}
				/*previous choice is field coding*/
							else
							{
								PreviousQp1=PreviousQp2;
								PreviousQp2=FieldQPBuffer;
							}
						}
					}
				}
					if(PreviousQp1==PreviousQp2)
						m_Qc=PreviousQp1+2;
					else
						m_Qc=(PreviousQp1+PreviousQp2)/2+1;
					m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
					m_Qc = MAX(RC_MIN_QUANT, m_Qc);//clipping
				}
				else
				{
					BFrameNumber=(NumberofBFrames+1)%input->successive_Bframe;
				if(BFrameNumber==0)
					BFrameNumber=input->successive_Bframe;
				/*adaptive field/frame coding*/
				else if(BFrameNumber==1)
				{
					if((input->InterlaceCodingOption==1)\
						||(input->InterlaceCodingOption==3))
					{
						if(img->FieldControl==0)
						{
					/*previous choice is frame coding*/
							if(img->FieldFrame==1)
							{
								PreviousQp1=PreviousQp2;
								PreviousQp2=FrameQPBuffer;
							}
					/*previous choice is field coding*/
							else
							{
								PreviousQp1=PreviousQp2;
								PreviousQp2=FieldQPBuffer;
							}
						}
					
					}
				}
				
				  if((PreviousQp2-PreviousQp1)<=(-2*input->successive_Bframe-3))
						StepSize=-3;
				  else 	if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe-2))
						StepSize=-2;
					else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe-1))
						StepSize=-1;
					else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe))
						StepSize=0;
					else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe+1))
						StepSize=1;
					else
						StepSize=2;
					m_Qc=PreviousQp1+StepSize;
					m_Qc +=MIN(2*(BFrameNumber-1),MAX(-2*(BFrameNumber-1), \
						(BFrameNumber-1)*(PreviousQp2-PreviousQp1)/(input->successive_Bframe-1)));
					m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
					m_Qc = MAX(RC_MIN_QUANT, m_Qc);//clipping
				}
				return m_Qc;
			}
			/*bottom field of B frame*/
			else
				return m_Qc;
		}
		else if(img->type==INTER_IMG)
		{
			if((img->NumberofGOP==1)&&(img->NumberofPPicture==0))
			{
				if((img->FieldControl==0)||((img->FieldControl==1)\
					&&(img->IFLAG==0)))
				{
				/*top field of the first P frame*/
					m_Qc=InitialQp;
					img->NumberofBasicUnitHeaderBits=0;
					img->NumberofBasicUnitTextureBits=0;
					NumberofBasicUnit--;
				/*bottom field of the first P frame*/
					if((!topfield)&&(NumberofBasicUnit==0))
					{
						/*frame coding or field coding*/
						if((input->InterlaceCodingOption==0)||(input->InterlaceCodingOption==2))
						{
							img->TotalQpforPPicture +=m_Qc;
							PreviousQp1=PreviousQp2;
							PreviousQp2=m_Qc;
							PAveFrameQP=m_Qc;
							PAveHeaderBits3=PAveHeaderBits2;
						}
						/*adaptive frame/field coding*/
						else if((input->InterlaceCodingOption==1)\
							||(input->InterlaceCodingOption==3))
						{
							if(img->FieldControl==0)
							{
								FrameQPBuffer=m_Qc;
								FrameAveHeaderBits=PAveHeaderBits2;
							}
							else
							{
								FieldQPBuffer=m_Qc;
								FieldAveHeaderBits=PAveHeaderBits2;
							}
						}
					}
					Pm_Qp=m_Qc;
					TotalFrameQP +=m_Qc;
					return m_Qc;
				}
			}
			else
			{
				m_X1=Pm_X1;
				m_X2=Pm_X2;
				m_Hp=PPreHeader;
				m_Qp=Pm_Qp;
				DuantQp=PDuantQp;
				MADPictureC1=PMADPictureC1;
				MADPictureC2=PMADPictureC2;

				if(img->FieldControl==0)
					SumofBasicUnit=TotalNumberofBasicUnit;
				else
					SumofBasicUnit=TotalNumberofBasicUnit/2;


				/*the average QP of the previous frame is used to coded the first basic unit of the current frame or field*/
				if(NumberofBasicUnit==SumofBasicUnit)
				{

					/*adaptive field/frame coding*/
					if(((input->InterlaceCodingOption==1)\
						||(input->InterlaceCodingOption==3))\
						&&(img->FieldControl==0))
					{
						
				/*previous choice is frame coding*/
						if(img->FieldFrame==1)
						{
							if(img->NumberofPPicture>0)
								img->TotalQpforPPicture +=FrameQPBuffer;
							PAveFrameQP=FrameQPBuffer;
							PAveHeaderBits3=FrameAveHeaderBits;
						}
				
				/*previous choice is field coding*/
						else
						{
							if(img->NumberofPPicture>0)
								img->TotalQpforPPicture +=FieldQPBuffer;
							PAveFrameQP=FieldQPBuffer;
							PAveHeaderBits3=FieldAveHeaderBits;
						}
					}
						
			
					if(T<=0)
					{
						m_Qc=PAveFrameQP+2;
						if(m_Qc>RC_MAX_QUANT)
							m_Qc=RC_MAX_QUANT;
						if(topfield||(img->FieldControl==0))
							GOPOverdue=TRUE;
					}
					else
						m_Qc=PAveFrameQP;
					TotalFrameQP +=m_Qc;
					NumberofBasicUnit--;
					Pm_Qp=PAveFrameQP;
					return m_Qc;
				}
               else
			   {

    	/*compute the number of remaining bits*/			
				   T -=img->NumberofBasicUnitHeaderBits+img->NumberofBasicUnitTextureBits;
				   img->NumberofBasicUnitHeaderBits=0;
				   img->NumberofBasicUnitTextureBits=0;
				   if(T<0)
				   {
					   if(GOPOverdue==TRUE)
						   m_Qc=m_Qp+2;
					   else 
						   m_Qc=m_Qp+DDquant;
					   m_Qc = MIN(m_Qc, RC_MAX_QUANT);	// clipping
					   if(input->basicunit>MBPerRow)
						   m_Qc = MIN(m_Qc, PAveFrameQP+6);
					   else
						   m_Qc = MIN(m_Qc, PAveFrameQP+3);

					   TotalFrameQP +=m_Qc;
					   NumberofBasicUnit--;
					   if(NumberofBasicUnit==0)
					   {
						   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);
								   img->TotalQpforPPicture +=PAverageQP;
								   if(GOPOverdue==TRUE)
								   {
									   PreviousQp1=PreviousQp2+1;
									   PreviousQp2=PAverageQP;
								   }
								   else
								   {
									  
									   if((img->NumberofPPicture==0)&&(img->NumberofGOP>1))
									   {
										   PreviousQp1=PreviousQp2;
										   PreviousQp2=PAverageQP;
									   }
									   else if(img->NumberofPPicture>0)
									   {
										    PreviousQp1=PreviousQp2+1;
									   /*LIZG 10/2/2003,because of the above step, the visual quality may not smooth enough*/
						    				PreviousQp2=PAverageQP;

									   }
								   }
								   PAveFrameQP=PAverageQP;
								   PAveHeaderBits3=PAveHeaderBits2;
							   }
							   /*adaptive field/frame coding*/
							   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;
								   }

							   }
						   }
					   }
					   if(GOPOverdue==TRUE)
						   Pm_Qp=PAveFrameQP;
					   else
						   Pm_Qp=m_Qc;
					   return m_Qc;
				   }
				   else
				   {
					   
		   /*LIZG JVT50V2,  predict the MAD of current picture*/
					   if(((input->InterlaceCodingOption==1)||(input->InterlaceCodingOption==3))\
						   &&(img->FieldControl==1))
						   CurrentFrameMAD=MADPictureC1*FCBUPFMAD[TotalNumberofBasicUnit-NumberofBasicUnit]+MADPictureC2;
					   else
						   CurrentFrameMAD=MADPictureC1*BUPFMAD[TotalNumberofBasicUnit-NumberofBasicUnit]+MADPictureC2;
					      
       /*compute the total number of bits for the current basic unit*/
					   m_Bits =T/NumberofBasicUnit;


					   /*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_Qc = (int) (m_X1 * CurrentFrameMAD / (double) m_Bits);
					   else // 2nd order mode
						   m_Qc = (int) ((2 * m_X2 * CurrentFrameMAD) / (sqrt (dtmp) - m_X1 * CurrentFrameMAD));
					   m_Qc = MIN(m_Qp+DDquant,  m_Qc);	// control variation +1
					   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==INTER_IMG))
					   {
						   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);
								   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==INTER_IMG)
	{
/*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 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);
			}
			/*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;
			/*update the record of MADs*/
			if(((input->InterlaceCodingOption==1)||(input->InterlaceCodingOption==3))\
				&&(img->FieldControl==1))
				FCBUCFMAD[TotalNumberofBasicUnit-1-NumberofBasicUnit]=CurrentFrameMAD;
			else
				BUCFMAD[TotalNumberofBasicUnit-1-NumberofBasicUnit]=CurrentFrameMAD;
			img->TotalMADBasicUnit=0;
			
			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] = m_Qc;
		/*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);

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


	for (i = 0; i < 20; i++) {
		m_rgRejected[i] = FALSE;
	}
//update the MAD for the previous frame
	if(img->type==INTER_IMG)
		PreviousFrameMAD=CurrentFrameMAD;

	// initial RD model estimator
	RCModelEstimator (n_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==INTER_IMG)
		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==INTER_IMG)
	{
		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
/*LIZG, JVT019r1*/
void updateMADModel ()
{

	int n_windowSize;
	int i;
	double error[20], std = 0.0, threshold;
	int m_Nc;

	if(img->NumberofCodedPFrame>0)
	{
   
	if(img->type==INTER_IMG)
	{
		/*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);

	/*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==INTER_IMG)
		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);
	}
}
/*LIZG, JVT50V2*/
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==INTER_IMG)
	{
		PMADPictureC1=MADPictureC1;
		PMADPictureC2=MADPictureC2;
	}
}