ratectl.c

G729、h263、h264、MPEG4四种最流行的音频和视频标准的压缩和解压算法的源代码.rar

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

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

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

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

  // default RD model estimation results

  m_X1 = m_X2 = 0.0;

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

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

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


//update the parameters of linear prediction model
void updateMADModel ()
{
  
  int n_windowSize;
  int i;
  double error[20], std = 0.0, threshold;
  int m_Nc;
  
  if(img->NumberofCodedPFrame>0)
  {
    
    if(img->type==P_SLICE)
    {
      /*frame layer rate control*/
      if(img->BasicUnit==img->Frame_Total_Number_MB)
        m_Nc=img->NumberofCodedPFrame;
      /*basic unit layer rate control*/
      else
        m_Nc=img->NumberofCodedPFrame*TotalNumberofBasicUnit+CodedBasicUnit;
      for (i = 19; i > 0; i--) {// update the history
        PPictureMAD[i] = PPictureMAD[i - 1];
        PictureMAD[i]=PPictureMAD[i];
        ReferenceMAD[i]= ReferenceMAD[i-1];
      }
      PPictureMAD[0] = CurrentFrameMAD;
      PictureMAD[0]=PPictureMAD[0];
      if(img->BasicUnit==img->Frame_Total_Number_MB)
        ReferenceMAD[0]=PictureMAD[1];
      else
      {
        if(((input->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);
  }
}

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


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

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