www.pudn.com > lencod.rar > ratectl.c


/*
***********************************************************************
* COPYRIGHT AND WARRANTY INFORMATION
*
* Copyright 2003, Advanced Audio Video Coding Standard, Part II
*
* DISCLAIMER OF WARRANTY
*
* These software programs are available to the users without any
* license fee or royalty on an "as is" basis. The AVS 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 contributors or
* the AVS be liable for any incidental, punitive, or consequential
* damages of any kind whatsoever arising from the use of this program.
*
* This disclaimer of warranty extends to the user of this program
* and user's customers, employees, agents, transferees, successors,
* and assigns.
*
* The AVS does not represent or warrant that the program furnished
* hereunder are free of infringement of any third-party patents.
* Commercial implementations of AVS, including shareware, may be
* subject to royalty fees to patent holders. Information regarding
* the AVS patent policy is available from the AVS Web site at
* http://www.avs.org.cn
*
* THIS IS NOT A GRANT OF PATENT RIGHTS - SEE THE AVS 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"
#include "header.h"

const double THETA=1.3636;
const int Switch=0;

int Iprev_bits=0;
int Pprev_bits=0;

/* rate control variables */
int Xp, Xb;
static int R,T_field;
static int Np, Nb, bits_topfield, Q;
long T,T1;
//HRD consideration
long UpperBound1, UpperBound2, LowerBound;
double InitialDelayOffset;
const double OMEGA=0.9;

double Wp,Wb;
int TotalPFrame;
int DuantQp;
int PDuantQp;
FILE *BitRate;
double DeltaP;

// Initiate rate control parameters
void rc_init_seq()
{
double L1,L2,L3,bpp;
int qp;
int i;

Xp=0;
Xb=0;

bit_rate=input->bit_rate;
frame_rate =(int)((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;

MINVALUE=4.0;
/*initialize the parameters of fluid flow traffic model*/

BufferSize=bit_rate*2.56;
CurrentBufferFullness=0;
GOPTargetBufferLevel=CurrentBufferFullness;
/*HRD consideration*/
InitialDelayOffset=BufferSize*0.8;

/*initialize the previous window size*/
m_windowSize=0;
MADm_windowSize=0;
img->NumberofCodedBFrame=0;
img->NumberofCodedPFrame=0;
img->NumberofGOP=0;
/*remaining # of bits in GOP */
R = 0;
/*control parameter */
if(input->successive_Bframe>0)
{
//GAMMAP=0.25;
BETAP=0.9;
}
else
{
GAMMAP=0.5;
BETAP=0.5;
}

/*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;

RC_MAX_QUANT = 63; // clipping
RC_MIN_QUANT = 0;//clipping

/*compute thei initial QP*/
bpp = 1.0*bit_rate /(frame_rate*img->width*img->height);
if (img->width == 176)
{
L1 = 0.1;
L2 = 0.3;
L3 = 0.6;
}else if (img->width == 352)
{
L1 = 0.2;
L2 = 0.6;
L3 = 1.2;
}else
{
L1 = 0.6;
L2 = 1.4;
L3 = 2.4;
}

if (input->SeinitialQP==0)
{
if(bpp<= L1)
qp = 35;
else
if(bpp<=L2)
qp = 25;
else
if(bpp<=L3)
qp = 20;
else
qp =10;
input->SeinitialQP = qp;
}
}

// Initiate one GOP
void rc_init_GOP(int np, int nb)
{
Boolean Overum=FALSE;
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.*/

if(R<0)
Overum=TRUE;
OverBits=-R;

/*initialize the lower bound and the upper bound for the target bits of each frame, HRD consideration*/
LowerBound=(long)(R+bit_rate/frame_rate);
UpperBound1=(long)(R+InitialDelayOffset);

/*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;

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

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

/*Compute InitialQp for each GOP*/
TotalPFrame=np;
img->NumberofGOP++;
if(img->NumberofGOP==1)
{
MyInitialQp=input->SeinitialQP;
PreviousQp2=MyInitialQp-1; //recent change -0;
QPLastGOP=MyInitialQp;
}
else
{
/*adaptive field/frame coding*/
if((input->InterlaceCodingOption==2))
{
if (img->FieldFrame == 1)
{
img->TotalQpforPPicture += FrameQPBuffer;
QPLastPFrame = FrameQPBuffer;
}
else
{
img->TotalQpforPPicture += FieldQPBuffer;
QPLastPFrame = FieldQPBuffer;
}

}
/*compute the average QP of P frames in the previous GOP*/
PAverageQp=(int)(1.0*img->TotalQpforPPicture/img->NumberofPPicture+0.5);

GOPDquant=(int)(0.5+1.0*(np+nb+1)/15);
if(GOPDquant>2)
GOPDquant=2;

PAverageQp-=GOPDquant;

if (PAverageQp > (QPLastPFrame - 2))
PAverageQp--;
PAverageQp = MAX(QPLastGOP-2, PAverageQp);
PAverageQp = MIN(QPLastGOP+2, PAverageQp);
PAverageQp = MIN(RC_MAX_QUANT, PAverageQp);
PAverageQp = MAX(RC_MIN_QUANT, PAverageQp);


MyInitialQp=PAverageQp;
QPLastGOP = MyInitialQp;
Pm_Qp=PAverageQp;
PAveFrameQP=PAverageQp;
PreviousQp1=PreviousQp2;
PreviousQp2=MyInitialQp-1;
}

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*/
img->NumberofCodedMacroBlocks=0;

/*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;
}

/*predefine a target buffer level for each frame*/
if((fieldpic||topfield)&amt;&amt;targetcomputation)
{
switch (img->type)
{
case INTER_IMG:
/*Since the available bandwidth may vary at any time, the total number of
bits is updated picture by picture*/
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==2))\
&amt;&amt;(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;
}
}

if(img->NumberofCodedPFrame==1)
AWp=Wp;
if((img->NumberofCodedPFrame<8)&amt;&amt;(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 average complexity of B frames
if(input->successive_Bframe>0)
{
//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)&amt;&amt;(img->NumberofCodedBFrame==1))
{
AWp=Wp;
AWb=Wb;
}
else if(img->NumberofCodedBFrame>1)
{
//compute the average weight
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 = (long) floor(Wp*R/(Np*Wp+Nb*Wb) + 0.5);

T1 = (long) floor(bit_rate/frame_rate-GAMMAP*(CurrentBufferFullness-TargetBufferLevel)+0.5);
T1=MAX(0,T1);
T = (long)(floor(BETAP*T+(1.0-BETAP)*T1+0.5));
}
}
/*basic unit layer rate control*/
else
{
if((img->NumberofGOP==1)&amt;&amt;(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);
T1=MAX(0,T1);
T = (int)(floor(BETAP*T+(1.0-BETAP)*T1+0.5));
}
else if(img->NumberofGOP>1)
{
T = (long) floor(Wp*R/(Np*Wp+Nb*Wb) + 0.5);
T1 = (long) floor(bit_rate/frame_rate-GAMMAP*(CurrentBufferFullness-TargetBufferLevel)+0.5);
T1 = MAX(0,T1);
T = (long)(floor(BETAP*T+(1.0-BETAP)*T1+0.5));
}
}

/*reserve some bits for smoothing*/

T=(long)((1.0-0.0*input->successive_Bframe)*T);
/*HRD consideration*/
T = MAX(T, (long) LowerBound);
T = MIN(T, (long) UpperBound2);

if((topfield)||(fieldpic&amt;&amt;((input->InterlaceCodingOption==2))))
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)&amt;&amt;(img->BasicUnit<img->Frame_Total_Number_MB)\
&amt;&amt;(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;
}
}
}

//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 one picture after frame/field encoding
void rc_update_pict(int nbits)
{
R-= nbits; /* remaining # of bits in GOP */
CurrentBufferFullness += nbits - bit_rate/frame_rate;

/*update the lower bound and the upper bound for the target bits of each frame, HRD consideration*/
LowerBound +=(long)(bit_rate/frame_rate-nbits);
UpperBound1 +=(long)(bit_rate/frame_rate-nbits);
UpperBound2 = (long)(OMEGA*UpperBound1);

return;
}

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

/*update the
complexity weight of I, P, B frame*/
int Avem_Qc;
int X;

/*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)&amt;&amt;(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*/
if(((img->IFLAG==0)&amt;&amt;(img->FieldControl==1))\
||(img->FieldControl==0))
{
Xp = X;
Np--;
Wp=Xp;
Pm_Hp=img->NumberofHeaderBits;
img->NumberofCodedPFrame++;
img->NumberofPPicture++;
}
else if((img->IFLAG!=0)&amt;&amt;(img->FieldControl==1))
img->IFLAG=0;
break;
case B_IMG:
Xb = X;
Nb--;
Wb=Xb/THETA;

img->NumberofCodedBFrame++;
NumberofBFrames++;

break;
}
}

// 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;
int i;

/*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*/
if((topfield)||(img->FieldControl==0))
{
if(img->type==INTRA_IMG)
{
m_Qc=MyInitialQp;
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==2))
{
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==2))
{
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)&amt;&amt;(img->NumberofPPicture==0))
{
m_Qc=MyInitialQp;

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==2))\
&amt;&amt;(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 picture*/
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_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+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)&amt;&amt;(img->IFLAG==0))
{
/*field coding*/
if(input->InterlaceCodingOption==1)
{
img->TotalQpforPPicture +=m_Qc;
PreviousQp1=PreviousQp2+1;
PreviousQp2=m_Qc;//+0 Recent change 13/1/2003
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=MyInitialQp;
return m_Qc;
}
/*bottom field of I frame*/
else if((img->type==INTER_IMG)&amt;&amt;(img->IFLAG==1)&amt;&amt;(img->FieldControl==1))
{
m_Qc=MyInitialQp;
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==2))
{
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==2))
{
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;//0
else if((PreviousQp2-PreviousQp1)==(-2*input->successive_Bframe+1))
StepSize=1;//1
else
StepSize=2;//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)&amt;&amt;(img->NumberofPPicture==0))
{
if((img->FieldControl==0)||((img->FieldControl==1)\
&amt;&amt;(img->IFLAG==0)))
{
/*top field of the first P frame*/
m_Qc=MyInitialQp;
img->NumberofBasicUnitHeaderBits=0;
img->NumberofBasicUnitTextureBits=0;
NumberofBasicUnit--;
/*bottom field of the first P frame*/
if((!topfield)&amt;&amt;(NumberofBasicUnit==0))
{
/*frame coding or field coding*/
if((input->InterlaceCodingOption==0)||(input->InterlaceCodingOption==1))
{
img->TotalQpforPPicture +=m_Qc;
PreviousQp1=PreviousQp2;
PreviousQp2=m_Qc;
PAveFrameQP=m_Qc;
PAveHeaderBits3=PAveHeaderBits2;
}
/*adaptive frame/field coding*/
else if((input->InterlaceCodingOption==2))
{
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==2))\
&amt;&amt;(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*/
TotalBasicUnitBits=img->NumberofBasicUnitHeaderBits+img->NumberofBasicUnitTextureBits;
T -=TotalBasicUnitBits;
img->NumberofBasicUnitHeaderBits=0;
img->NumberofBasicUnitTextureBits=0;
if(T<0)
{
if(GOPOverdue==TRUE)
m_Qc=m_Qp+2;
else
m_Qc=m_Qp+DDquant;//2
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==1))
{
PAverageQP=(int)(1.0*TotalFrameQP/TotalNumberofBasicUnit+0.5);
if (img->NumberofPPicture == (input->intra_period - 2))
QPLastPFrame = PAverageQP;

img->TotalQpforPPicture +=PAverageQP;
if(GOPOverdue==TRUE)
{
PreviousQp1=PreviousQp2+1;
PreviousQp2=PAverageQP;
}
else
{
if((img->NumberofPPicture==0)&amt;&amt;(img->NumberofGOP>1))
{
PreviousQp1=PreviousQp2;
PreviousQp2=PAverageQP;
}
else if(img->NumberofPPicture>0)
{
PreviousQp1=PreviousQp2+1;
PreviousQp2=PAverageQP;
}
}
PAveFrameQP=PAverageQP;
PAveHeaderBits3=PAveHeaderBits2;
}
/*adaptive field/frame coding*/
else if((input->InterlaceCodingOption==2))
{
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
{
/*predict the MAD of current picture*/
if(((input->InterlaceCodingOption==2))\
&amt;&amt;(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)&amt;&amt;(img->type==INTER_IMG))
{
if((!topfield)||(img->FieldControl==0))
{
/*frame coding or field coding*/
if((input->InterlaceCodingOption==0)||(input->InterlaceCodingOption==1))
{
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==2))
{
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 MAD of the current basic unit*/
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==2))\
&amt;&amt;(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==INTER_IMG)
PPictureMAD[0]=CurrentFrameMAD;
}
}

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) &amt;&amt; !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) &amt;&amt; 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
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==2))\
&amt;&amt;(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==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);
}
}

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) &amt;&amt; !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) &amt;&amt; 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;
}
}


double QP2Qstep( int QP )
{
int i;
double Qstep;
static const double QP2QSTEP[8] = { 1.0, 1.0905, 1.189, 1.297, 1.414, 1.542, 1.682, 1.834};

Qstep = QP2QSTEP[QP > 8];
for( i=0; i<(QP/8); i++)
Qstep *= 2;

return Qstep;
}

int Qstep2QP( double Qstep )
{
int q_per = 0, q_rem = 0;

// assert( Qstep >= QP2Qstep(0) &amt;&amt; Qstep <= QP2Qstep(51) );
if( Qstep < QP2Qstep(0))
return 0;
else if (Qstep > QP2Qstep(63) )
return 63;

while( Qstep > QP2Qstep(7) )
{
Qstep /= 2;
q_per += 1;
}

if (Qstep <= (1.0+1.0905)/2)
{
Qstep = 1.0;
q_rem = 0;
}
else if (Qstep <= (1.0905+1.1895)/2)
{
Qstep = 1.0905;
q_rem = 1;
}
else if (Qstep <= (1.189+1.297)/2)
{
Qstep = 1.189;
q_rem = 2;
}
else if (Qstep <= (1.297+1.414)/2)
{
Qstep = 1.297;
q_rem = 3;
}
else if (Qstep <= (1.414+1.542)/2)
{
Qstep = 1.414;
q_rem = 4;
}
else if (Qstep <= (1.542+1.682)/2)
{
Qstep = 1.414;
q_rem = 5;
}
else if (Qstep <= (1.682+1.834)/2)
{
Qstep = 1.682;
q_rem = 6;
}
else
{
Qstep = 1.834;
q_rem = 7;
}

return (q_per * 8 + q_rem);
}