ratectl.c

h.264官方测试软件

    else
    {
      /*compute the MAD of the current basic unit*/
      if((input->MbInterlace)&&(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->PicInterlace==ADAPTIVE_CODING)||(input->MbInterlace))\
        &&(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;
  } 
}

/*! 
 *************************************************************************************
 * \brief
 *    Model Estimator
 *
 *************************************************************************************
*/
void RCModelEstimator (int n_windowSize)
{
  int n_realSize = n_windowSize;
  int i;
  double oneSampleQ = 0;
  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;
  }
}


/*! 
 *************************************************************************************
 * \brief
 *    Compute Frame MAD
 *
 *************************************************************************************
*/
double ComputeFrameMAD()
{
  double TotalMAD;
  int i;
  TotalMAD=0.0;
  for(i=0;i<img->Frame_Total_Number_MB;i++)
    TotalMAD +=img->MADofMB[i];
  TotalMAD /=img->Frame_Total_Number_MB;
  return TotalMAD;
}


/*! 
 *************************************************************************************
 * \brief
 *    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)
  {
    //assert (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->PicInterlace==ADAPTIVE_CODING)||(input->MbInterlace))\
        &&(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);
  }
}


/*! 
 *************************************************************************************
 * \brief
 *    MAD mode estimator
 *
 *************************************************************************************
*/
void MADModelEstimator (int n_windowSize)
{
  int n_realSize = n_windowSize;
  int i;
  double oneSampleQ = 0;
  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;
  }
}


/*! 
 *************************************************************************************
 * \brief
 *    map QP to Qstep
 *
 *************************************************************************************
*/
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;
}


/*! 
 *************************************************************************************
 * \brief
 *    map Qstep to QP
 *
 *************************************************************************************
*/
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);
}