src/gdcmmpeg2/src/mpeg2enc/ratectl.c

   1 /* ratectl.c, bitrate control routines (linear quantization only currently) */
   2
   3 /* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */
   4
   5 /*
   6  * Disclaimer of Warranty
   7  *
   8  * These software programs are available to the user without any license fee or
   9  * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
  10  * any and all warranties, whether express, implied, or statuary, including any
  11  * implied warranties or merchantability or of fitness for a particular
  12  * purpose.  In no event shall the copyright-holder be liable for any
  13  * incidental, punitive, or consequential damages of any kind whatsoever
  14  * arising from the use of these programs.
  15  *
  16  * This disclaimer of warranty extends to the user of these programs and user's
  17  * customers, employees, agents, transferees, successors, and assigns.
  18  *
  19  * The MPEG Software Simulation Group does not represent or warrant that the
  20  * programs furnished hereunder are free of infringement of any third-party
  21  * patents.
  22  *
  23  * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
  24  * are subject to royalty fees to patent holders.  Many of these patents are
  25  * general enough such that they are unavoidable regardless of implementation
  26  * design.
  27  *
  28  */
  29
  30 #include <stdio.h>
  31 #include <math.h>
  32
  33 #include "config.h"
  34 #include "global.h"
  35
  36 /* private prototypes */
  37 static void calc_actj _ANSI_ARGS_((unsigned char *frame));
  38 static double var_sblk _ANSI_ARGS_((unsigned char *p, int lx));
  39
  40 /* rate control variables */
  41 int Xi, Xp, Xb, r, d0i, d0p, d0b;
  42 double avg_act;
  43 static int R, T, d;
  44 static double actsum;
  45 static int Np, Nb, S, Q;
  46 static int prev_mquant;
  47
  48 void rc_init_seq()
  49 {
  50   /* reaction parameter (constant) */
  51   if (r==0)  r = (int)floor(2.0*bit_rate/frame_rate + 0.5);
  52
  53   /* average activity */
  54   if (avg_act==0.0)  avg_act = 400.0;
  55
  56   /* remaining # of bits in GOP */
  57   R = 0;
  58
  59   /* global complexity measure */
  60   if (Xi==0) Xi = (int)floor(160.0*bit_rate/115.0 + 0.5);
  61   if (Xp==0) Xp = (int)floor( 60.0*bit_rate/115.0 + 0.5);
  62   if (Xb==0) Xb = (int)floor( 42.0*bit_rate/115.0 + 0.5);
  63
  64   /* virtual buffer fullness */
  65   if (d0i==0) d0i = (int)floor(10.0*r/31.0 + 0.5);
  66   if (d0p==0) d0p = (int)floor(10.0*r/31.0 + 0.5);
  67   if (d0b==0) d0b = (int)floor(1.4*10.0*r/31.0 + 0.5);
  68 /*
  69   if (d0i==0) d0i = (int)floor(10.0*r/(qscale_tab[0] ? 56.0 : 31.0) + 0.5);
  70   if (d0p==0) d0p = (int)floor(10.0*r/(qscale_tab[1] ? 56.0 : 31.0) + 0.5);
  71   if (d0b==0) d0b = (int)floor(1.4*10.0*r/(qscale_tab[2] ? 56.0 : 31.0) + 0.5);
  72 */
  73
  74   fprintf(statfile,"\nrate control: sequence initialization\n");
  75   fprintf(statfile,
  76     " initial global complexity measures (I,P,B): Xi=%d, Xp=%d, Xb=%d\n",
  77     Xi, Xp, Xb);
  78   fprintf(statfile," reaction parameter: r=%d\n", r);
  79   fprintf(statfile,
  80     " initial virtual buffer fullness (I,P,B): d0i=%d, d0p=%d, d0b=%d\n",
  81     d0i, d0p, d0b);
  82   fprintf(statfile," initial average activity: avg_act=%.1f\n", avg_act);
  83 }
  84
  85 void rc_init_GOP(np,nb)
  86 int np,nb;
  87 {
  88   R += (int) floor((1 + np + nb) * bit_rate / frame_rate + 0.5);
  89   Np = fieldpic ? 2*np+1 : np;
  90   Nb = fieldpic ? 2*nb : nb;
  91
  92   fprintf(statfile,"\nrate control: new group of pictures (GOP)\n");
  93   fprintf(statfile," target number of bits for GOP: R=%d\n",R);
  94   fprintf(statfile," number of P pictures in GOP: Np=%d\n",Np);
  95   fprintf(statfile," number of B pictures in GOP: Nb=%d\n",Nb);
  96 }
  97
  98 /* Note: we need to substitute K for the 1.4 and 1.0 constants -- this can
  99    be modified to fit image content */
 100
 101 /* Step 1: compute target bits for current picture being coded */
 102 void rc_init_pict(frame)
 103 unsigned char *frame;
 104 {
 105   double Tmin;
 106
 107   switch (pict_type)
 108   {
 109   case I_TYPE:
 110     T = (int) floor(R/(1.0+Np*Xp/(Xi*1.0)+Nb*Xb/(Xi*1.4)) + 0.5);
 111     d = d0i;
 112     break;
 113   case P_TYPE:
 114     T = (int) floor(R/(Np+Nb*1.0*Xb/(1.4*Xp)) + 0.5);
 115     d = d0p;
 116     break;
 117   case B_TYPE:
 118     T = (int) floor(R/(Nb+Np*1.4*Xp/(1.0*Xb)) + 0.5);
 119     d = d0b;
 120     break;
 121   }
 122
 123   Tmin = (int) floor(bit_rate/(8.0*frame_rate) + 0.5);
 124
 125   if (T<Tmin)
 126     T = (int)Tmin;
 127
 128   S = bitcount();
 129   Q = 0;
 130
 131   calc_actj(frame);
 132   actsum = 0.0;
 133
 134   fprintf(statfile,"\nrate control: start of picture\n");
 135   fprintf(statfile," target number of bits: T=%d\n",T);
 136 }
 137
 138 static void calc_actj(frame)
 139 unsigned char *frame;
 140 {
 141   int i,j,k;
 142   unsigned char *p;
 143   double actj,var;
 144
 145   k = 0;
 146
 147   for (j=0; j<height2; j+=16)
 148     for (i=0; i<width; i+=16)
 149     {
 150       p = frame + ((pict_struct==BOTTOM_FIELD)?width:0) + i + width2*j;
 151
 152       /* take minimum spatial activity measure of luminance blocks */
 153
 154       actj = var_sblk(p,width2);
 155       var = var_sblk(p+8,width2);
 156       if (var<actj) actj = var;
 157       var = var_sblk(p+8*width2,width2);
 158       if (var<actj) actj = var;
 159       var = var_sblk(p+8*width2+8,width2);
 160       if (var<actj) actj = var;
 161
 162       if (!fieldpic && !prog_seq)
 163       {
 164         /* field */
 165         var = var_sblk(p,width<<1);
 166         if (var<actj) actj = var;
 167         var = var_sblk(p+8,width<<1);
 168         if (var<actj) actj = var;
 169         var = var_sblk(p+width,width<<1);
 170         if (var<actj) actj = var;
 171         var = var_sblk(p+width+8,width<<1);
 172         if (var<actj) actj = var;
 173       }
 174
 175       actj+= 1.0;
 176
 177       mbinfo[k++].act = actj;
 178     }
 179 }
 180
 181 void rc_update_pict()
 182 {
 183   double X;
 184
 185   S = bitcount() - S; /* total # of bits in picture */
 186   R-= S; /* remaining # of bits in GOP */
 187   X = (int) floor(S*((0.5*(double)Q)/(mb_width*mb_height2)) + 0.5);
 188   d+= S - T;
 189   avg_act = actsum/(mb_width*mb_height2);
 190
 191   switch (pict_type)
 192   {
 193   case I_TYPE:
 194     Xi = (int)X;
 195     d0i = d;
 196     break;
 197   case P_TYPE:
 198     Xp = (int)X;
 199     d0p = d;
 200     Np--;
 201     break;
 202   case B_TYPE:
 203     Xb = (int)X;
 204     d0b = d;
 205     Nb--;
 206     break;
 207   }
 208
 209   fprintf(statfile,"\nrate control: end of picture\n");
 210   fprintf(statfile," actual number of bits: S=%d\n",S);
 211   fprintf(statfile," average quantization parameter Q=%.1f\n",
 212     (double)Q/(mb_width*mb_height2));
 213   fprintf(statfile," remaining number of bits in GOP: R=%d\n",R);
 214   fprintf(statfile,
 215     " global complexity measures (I,P,B): Xi=%d, Xp=%d, Xb=%d\n",
 216     Xi, Xp, Xb);
 217   fprintf(statfile,
 218     " virtual buffer fullness (I,P,B): d0i=%d, d0p=%d, d0b=%d\n",
 219     d0i, d0p, d0b);
 220   fprintf(statfile," remaining number of P pictures in GOP: Np=%d\n",Np);
 221   fprintf(statfile," remaining number of B pictures in GOP: Nb=%d\n",Nb);
 222   fprintf(statfile," average activity: avg_act=%.1f\n", avg_act);
 223 }
 224
 225 /* compute initial quantization stepsize (at the beginning of picture) */
 226 int rc_start_mb()
 227 {
 228   int mquant;
 229
 230   if (q_scale_type)
 231   {
 232     mquant = (int) floor(2.0*d*31.0/r + 0.5);
 233
 234     /* clip mquant to legal (linear) range */
 235     if (mquant<1)
 236       mquant = 1;
 237     if (mquant>112)
 238       mquant = 112;
 239
 240     /* map to legal quantization level */
 241     mquant = non_linear_mquant_table[map_non_linear_mquant[mquant]];
 242   }
 243   else
 244   {
 245     mquant = (int) floor(d*31.0/r + 0.5);
 246     mquant <<= 1;
 247
 248     /* clip mquant to legal (linear) range */
 249     if (mquant<2)
 250       mquant = 2;
 251     if (mquant>62)
 252       mquant = 62;
 253
 254     prev_mquant = mquant;
 255   }
 256
 257 /*
 258   fprintf(statfile,"rc_start_mb:\n");
 259   fprintf(statfile,"mquant=%d\n",mquant);
 260 */
 261
 262   return mquant;
 263 }
 264
 265 /* Step 2: measure virtual buffer - estimated buffer discrepancy */
 266 int rc_calc_mquant(j)
 267 int j;
 268 {
 269   int mquant;
 270   double dj, Qj, actj, N_actj;
 271
 272   /* measure virtual buffer discrepancy from uniform distribution model */
 273   dj = d + (bitcount()-S) - j*(T/(mb_width*mb_height2));
 274
 275   /* scale against dynamic range of mquant and the bits/picture count */
 276   Qj = dj*31.0/r;
 277 /*Qj = dj*(q_scale_type ? 56.0 : 31.0)/r;  */
 278
 279   actj = mbinfo[j].act;
 280   actsum+= actj;
 281
 282   /* compute normalized activity */
 283   N_actj = (2.0*actj+avg_act)/(actj+2.0*avg_act);
 284
 285   if (q_scale_type)
 286   {
 287     /* modulate mquant with combined buffer and local activity measures */
 288     mquant = (int) floor(2.0*Qj*N_actj + 0.5);
 289
 290     /* clip mquant to legal (linear) range */
 291     if (mquant<1)
 292       mquant = 1;
 293     if (mquant>112)
 294       mquant = 112;
 295
 296     /* map to legal quantization level */
 297     mquant = non_linear_mquant_table[map_non_linear_mquant[mquant]];
 298   }
 299   else
 300   {
 301     /* modulate mquant with combined buffer and local activity measures */
 302     mquant = (int) floor(Qj*N_actj + 0.5);
 303     mquant <<= 1;
 304
 305     /* clip mquant to legal (linear) range */
 306     if (mquant<2)
 307       mquant = 2;
 308     if (mquant>62)
 309       mquant = 62;
 310
 311     /* ignore small changes in mquant */
 312     if (mquant>=8 && (mquant-prev_mquant)>=-4 && (mquant-prev_mquant)<=4)
 313       mquant = prev_mquant;
 314
 315     prev_mquant = mquant;
 316   }
 317
 318   Q+= mquant; /* for calculation of average mquant */
 319
 320 /*
 321   fprintf(statfile,"rc_calc_mquant(%d): ",j);
 322   fprintf(statfile,"bitcount=%d, dj=%f, Qj=%f, actj=%f, N_actj=%f, mquant=%d\n",
 323     bitcount(),dj,Qj,actj,N_actj,mquant);
 324 */
 325
 326   return mquant;
 327 }
 328
 329 /* compute variance of 8x8 block */
 330 static double var_sblk(p,lx)
 331 unsigned char *p;
 332 int lx;
 333 {
 334   int i, j;
 335   unsigned int v, s, s2;
 336
 337   s = s2 = 0;
 338
 339   for (j=0; j<8; j++)
 340   {
 341     for (i=0; i<8; i++)
 342     {
 343       v = *p++;
 344       s+= v;
 345       s2+= v*v;
 346     }
 347     p+= lx - 8;
 348   }
 349
 350   return s2/64.0 - (s/64.0)*(s/64.0);
 351 }
 352
 353 /* VBV calculations
 354  *
 355  * generates warnings if underflow or overflow occurs
 356  */
 357
 358 /* vbv_end_of_picture
 359  *
 360  * - has to be called directly after writing picture_data()
 361  * - needed for accurate VBV buffer overflow calculation
 362  * - assumes there is no byte stuffing prior to the next start code
 363  */
 364
 365 static int bitcnt_EOP;
 366
 367 void vbv_end_of_picture()
 368 {
 369   bitcnt_EOP = bitcount();
 370   bitcnt_EOP = (bitcnt_EOP + 7) & ~7; /* account for bit stuffing */
 371 }
 372
 373 /* calc_vbv_delay
 374  *
 375  * has to be called directly after writing the picture start code, the
 376  * reference point for vbv_delay
 377  */
 378
 379 void calc_vbv_delay()
 380 {
 381   double picture_delay;
 382   static double next_ip_delay; /* due to frame reordering delay */
 383   static double decoding_time;
 384
 385   /* number of 1/90000 s ticks until next picture is to be decoded */
 386   if (pict_type == B_TYPE)
 387   {
 388     if (prog_seq)
 389     {
 390       if (!repeatfirst)
 391         picture_delay = 90000.0/frame_rate; /* 1 frame */
 392       else
 393       {
 394         if (!topfirst)
 395           picture_delay = 90000.0*2.0/frame_rate; /* 2 frames */
 396         else
 397           picture_delay = 90000.0*3.0/frame_rate; /* 3 frames */
 398       }
 399     }
 400     else
 401     {
 402       /* interlaced */
 403       if (fieldpic)
 404         picture_delay = 90000.0/(2.0*frame_rate); /* 1 field */
 405       else
 406       {
 407         if (!repeatfirst)
 408           picture_delay = 90000.0*2.0/(2.0*frame_rate); /* 2 flds */
 409         else
 410           picture_delay = 90000.0*3.0/(2.0*frame_rate); /* 3 flds */
 411       }
 412     }
 413   }
 414   else
 415   {
 416     /* I or P picture */
 417     if (fieldpic)
 418     {
 419       if(topfirst==(pict_struct==TOP_FIELD))
 420       {
 421         /* first field */
 422         picture_delay = 90000.0/(2.0*frame_rate);
 423       }
 424       else
 425       {
 426         /* second field */
 427         /* take frame reordering delay into account */
 428         picture_delay = next_ip_delay - 90000.0/(2.0*frame_rate);
 429       }
 430     }
 431     else
 432     {
 433       /* frame picture */
 434       /* take frame reordering delay into account*/
 435       picture_delay = next_ip_delay;
 436     }
 437
 438     if (!fieldpic || topfirst!=(pict_struct==TOP_FIELD))
 439     {
 440       /* frame picture or second field */
 441       if (prog_seq)
 442       {
 443         if (!repeatfirst)
 444           next_ip_delay = 90000.0/frame_rate;
 445         else
 446         {
 447           if (!topfirst)
 448             next_ip_delay = 90000.0*2.0/frame_rate;
 449           else
 450             next_ip_delay = 90000.0*3.0/frame_rate;
 451         }
 452       }
 453       else
 454       {
 455         if (fieldpic)
 456           next_ip_delay = 90000.0/(2.0*frame_rate);
 457         else
 458         {
 459           if (!repeatfirst)
 460             next_ip_delay = 90000.0*2.0/(2.0*frame_rate);
 461           else
 462             next_ip_delay = 90000.0*3.0/(2.0*frame_rate);
 463         }
 464       }
 465     }
 466   }
 467
 468   if (decoding_time==0.0)
 469   {
 470     /* first call of calc_vbv_delay */
 471     /* we start with a 7/8 filled VBV buffer (12.5% back-off) */
 472     picture_delay = ((vbv_buffer_size*16384*7)/8)*90000.0/bit_rate;
 473     if (fieldpic)
 474       next_ip_delay = (int)(90000.0/frame_rate+0.5);
 475   }
 476
 477   /* VBV checks */
 478
 479   /* check for underflow (previous picture) */
 480   if (!low_delay && (decoding_time < bitcnt_EOP*90000.0/bit_rate))
 481   {
 482     /* picture not completely in buffer at intended decoding time */
 483     if (!quiet)
 484       fprintf(stderr,"vbv_delay underflow! (decoding_time=%.1f, t_EOP=%.1f\n)",
 485         decoding_time, bitcnt_EOP*90000.0/bit_rate);
 486   }
 487
 488   /* when to decode current frame */
 489   decoding_time += picture_delay;
 490
 491   /* warning: bitcount() may overflow (e.g. after 9 min. at 8 Mbit/s */
 492   vbv_delay = (int)(decoding_time - bitcount()*90000.0/bit_rate);
 493
 494   /* check for overflow (current picture) */
 495   if ((decoding_time - bitcnt_EOP*90000.0/bit_rate)
 496       > (vbv_buffer_size*16384)*90000.0/bit_rate)
 497   {
 498     if (!quiet)
 499       fprintf(stderr,"vbv_delay overflow!\n");
 500   }
 501
 502   fprintf(statfile,
 503     "\nvbv_delay=%d (bitcount=%d, decoding_time=%.2f, bitcnt_EOP=%d)\n",
 504     vbv_delay,bitcount(),decoding_time,bitcnt_EOP);
 505
 506   if (vbv_delay<0)
 507   {
 508     if (!quiet)
 509       fprintf(stderr,"vbv_delay underflow: %d\n",vbv_delay);
 510     vbv_delay = 0;
 511   }
 512
 513   if (vbv_delay>65535)
 514   {
 515     if (!quiet)
 516       fprintf(stderr,"vbv_delay overflow: %d\n",vbv_delay);
 517     vbv_delay = 65535;
 518   }
 519 }