Merge CV, ver=5.1; ops/methods from HV, and interface from CV where possible
[goodguy/history.git] / cinelerra-5.1 / mpeg2enc / ratectl.c
diff --git a/cinelerra-5.1/mpeg2enc/ratectl.c b/cinelerra-5.1/mpeg2enc/ratectl.c
new file mode 100644 (file)
index 0000000..7683843
--- /dev/null
@@ -0,0 +1,639 @@
+/* ratectl.c, bitrate control routines (linear quantization only currently) */
+
+/* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */
+
+/*
+ * Disclaimer of Warranty
+ *
+ * These software programs are available to the user without any license fee or
+ * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
+ * any and all warranties, whether express, implied, or statuary, including any
+ * implied warranties or merchantability or of fitness for a particular
+ * purpose.  In no event shall the copyright-holder 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 MPEG Software Simulation Group does not represent or warrant that the
+ * programs furnished hereunder are free of infringement of any third-party
+ * patents.
+ *
+ * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
+ * are subject to royalty fees to patent holders.  Many of these patents are
+ * general enough such that they are unavoidable regardless of implementation
+ * design.
+ *
+ */
+
+#include <stdio.h>
+#include <limits.h>
+#include <math.h>
+#include <pthread.h>
+
+#include "config.h"
+#include "global.h"
+#include "fastintfns.h"
+
+/* rate control variables */
+/*
+ * static double R, T, d;
+ * static double actsum;
+ * static int Np, Nb;
+ * static double S, Q;
+ * static int prev_mquant;
+ * static double bitcnt_EOP;
+ * static double next_ip_delay; // due to frame reordering delay
+ * static double decoding_time;
+ * static int Xi, Xp, Xb, r, d0i, d0p, d0b;
+ * static double avg_act;
+ */
+
+void ratectl_init_seq(ratectl_t *ratectl)
+{
+       pthread_mutexattr_t mutex_attr;
+       pthread_mutexattr_init(&mutex_attr);
+       pthread_mutex_init(&(ratectl->ratectl_lock), &mutex_attr);
+
+       ratectl->avg_KI = 2.5;   /* TODO: These values empirically determined            */
+       ratectl->avg_KB = 10.0;   /* for MPEG-1, may need tuning for MPEG-2   */
+       ratectl->avg_KP = 10.0;
+
+       ratectl->bits_per_mb = (double)bit_rate / (mb_per_pict);
+       /* reaction parameter (constant) decreased to increase response
+          rate as encoder is currently tending to under/over-shoot... in
+          rate TODO: Reaction parameter is *same* for every frame type
+          despite different weightings...  */
+
+       if (ratectl->r == 0)  
+               ratectl->r = (int)floor(2.0 * bit_rate / frame_rate + 0.5);
+
+       ratectl->Ki = 1.2;  /* EXPERIMENT: ADJUST activities for I MB's */
+       ratectl->Kb = 1.4;
+       ratectl->Kp = 1.1;
+
+/* average activity */
+       if (ratectl->avg_act == 0.0)  ratectl->avg_act = 400.0;
+
+/* remaining # of bits in GOP */
+       ratectl->R = 0;
+       ratectl->IR = 0;
+
+       /* Heuristic: In constant bit-rate streams we assume buffering
+          will allow us to pre-load some (probably small) fraction of the
+          buffers size worth of following data if earlier data was
+          undershot its bit-rate allocation
+        
+       */
+
+       ratectl->CarryR = 0;
+       ratectl->CarryRLim = video_buffer_size / 3;
+       /* global complexity (Chi! not X!) measure of different frame types */
+       /* These are just some sort-of sensible initial values for start-up */
+
+       ratectl->Xi = 1500*mb_per_pict;   /* Empirically derived values */
+       ratectl->Xp = 550*mb_per_pict;
+       ratectl->Xb = 170*mb_per_pict;
+       ratectl->d0i = -1;                              /* Force initial Quant prediction */
+       ratectl->d0pb = -1;
+
+       ratectl->current_quant = 1;
+}
+
+void ratectl_init_GOP(ratectl_t *ratectl, int np, int nb)
+{
+       double per_gop_bits = 
+               (double)(1 + np + nb) * (double)bit_rate / frame_rate;
+
+       /* A.Stevens Aug 2000: at last I've found the wretched
+          rate-control overshoot bug...  Simply "topping up" R here means
+          that we can accumulate an indefinately large pool of bits
+          "saved" from previous low-activity frames.  This is of
+          course nonsense.
+
+          In CBR we can only accumulate as much as our buffer allows, after that
+          the eventual system stream will have to be padded.  The automatic padding
+          will make this calculation fairly reasonable but since that's based on 
+          estimates we again impose our rough and ready heuristic that we can't
+          accumulate more than approximately half  a video buffer full.
+
+          In VBR we actually do nothing different.  Here the bitrate is
+          simply a ceiling rate which we expect to undershoot a lot as
+          our quantisation floor cuts in. We specify a great big buffer
+          and simply don't pad when we undershoot.
+
+          However, we don't want to carry over absurd undershoots as when it
+          does get hectic we'll breach our maximum.
+               
+          TODO: For CBR we should do a proper video buffer model and use
+          it to make bit allocation decisions.
+
+       */
+
+       if( ratectl->R > 0  )
+       {
+               /* We replacing running estimate of undershoot with
+                  *exact* value and use that for calculating how much we
+                  may "carry over"
+               */
+               ratectl->gop_undershoot = intmin( video_buffer_size/2, (int)ratectl->R );
+
+               ratectl->R = ratectl->gop_undershoot + per_gop_bits;            
+       }
+       else 
+       {
+               /* Overshoots are easy - we have to make up the bits */
+               ratectl->R +=  per_gop_bits;
+               ratectl->gop_undershoot = 0;
+       }
+       ratectl->IR = ratectl->R;
+       ratectl->Np = fieldpic ? 2 * np + 1 : np;
+       ratectl->Nb = fieldpic ? 2 * nb : nb;
+}
+
+static int scale_quant(pict_data_s *picture, double quant )
+{
+       int iquant;
+       
+       if (picture->q_scale_type  )
+       {
+               iquant = (int) floor(quant+0.5);
+
+               /* clip mquant to legal (linear) range */
+               if (iquant<1)
+                       iquant = 1;
+               if (iquant>112)
+                       iquant = 112;
+
+               iquant =
+                       non_linear_mquant_table_hv[map_non_linear_mquant_hv[iquant]];
+       }
+       else
+       {
+               /* clip mquant to legal (linear) range */
+               iquant = (int)floor(quant+0.5);
+               if (iquant<2)
+                       iquant = 2;
+               if (iquant>62)
+                       iquant = 62;
+               iquant = (iquant/2)*2; // Must be *even*
+       }
+       return iquant;
+}
+
+
+
+
+/* compute variance of 8x8 block */
+static double var_sblk(p, lx)
+unsigned char *p;
+int lx;
+{
+       int j;
+       register unsigned int v, s, s2;
+
+       s = s2 = 0;
+
+       for (j=0; j<8; j++)
+       {
+               v = p[0];   s += v;    s2 += v * v;
+               v = p[1];   s += v;    s2 += v * v;
+               v = p[2];   s += v;    s2 += v * v;
+               v = p[3];   s += v;    s2 += v * v;
+               v = p[4];   s += v;    s2 += v * v;
+               v = p[5];   s += v;    s2 += v * v;
+               v = p[6];   s += v;    s2 += v * v;
+               v = p[7];   s += v;    s2 += v * v;
+               p += lx;
+       }
+
+       return (double)s2 / 64.0 - ((double)s / 64.0) * ((double)s / 64.0);
+}
+
+
+static double calc_actj(pict_data_s *picture)
+{
+       int i,j,k,l;
+       double actj,sum;
+       uint16_t *i_q_mat;
+       int actsum;
+       sum = 0.0;
+       k = 0;
+
+       for (j=0; j<height2; j+=16)
+               for (i=0; i<width; i+=16)
+               {
+                       /* A.Stevens Jul 2000 Luminance variance *has* to be a rotten measure
+                          of how active a block in terms of bits needed to code a lossless DCT.
+                          E.g. a half-white half-black block has a maximal variance but 
+                          pretty small DCT coefficients.
+
+                          So.... we use the absolute sum of DCT coefficients as our
+                          variance measure.  
+                       */
+                       if( picture->mbinfo[k].mb_type  & MB_INTRA )
+                       {
+                               i_q_mat = i_intra_q;
+                               /* EXPERIMENT: See what happens if we compensate for
+                                the wholly disproprotionate weight of the DC
+                                coefficients.  Shold produce more sensible results...  */
+                               actsum =  -80*COEFFSUM_SCALE;
+                       }
+                       else
+                       {
+                               i_q_mat = i_inter_q;
+                               actsum = 0;
+                       }
+
+                       /* It takes some bits to code even an entirely zero block...
+                          It also makes a lot of calculations a lot better conditioned
+                          if it can be guaranteed that activity is always distinctly
+                          non-zero.
+                        */
+
+
+                       for( l = 0; l < 6; ++l )
+                               actsum += 
+                                       (*pquant_weight_coeff_sum)
+                                           ( cur_picture.mbinfo[k].dctblocks[l], i_q_mat ) ;
+                       actj = (double)actsum / (double)COEFFSUM_SCALE;
+                       if( actj < 12.0 )
+                               actj = 12.0;
+
+                       picture->mbinfo[k].act = (double)actj;
+                       sum += (double)actj;
+                       ++k;
+               }
+       return sum;
+}
+
+/* Note: we need to substitute K for the 1.4 and 1.0 constants -- this can
+   be modified to fit image content */
+
+/* Step 1: compute target bits for current picture being coded */
+void ratectl_init_pict(ratectl_t *ratectl, pict_data_s *picture)
+{
+       double avg_K;
+       double target_Q;
+       double current_Q;
+       double Si, Sp, Sb;
+       /* TODO: A.Stevens  Nov 2000 - This modification needs testing visually.
+
+          Weird.  The original code used the average activity of the
+          *previous* frame as the basis for quantisation calculations for
+          rather than the activity in the *current* frame.  That *has* to
+          be a bad idea..., surely, here we try to be smarter by using the
+          current values and keeping track of how much of the frames
+          activitity has been covered as we go along.
+
+          We also guesstimate the relationship between  (sum
+          of DCT coefficients) and actual quantisation weighted activty.
+          We use this to try to predict the activity of each frame.
+       */
+
+       ratectl->actsum =  calc_actj(picture );
+       ratectl->avg_act = (double)ratectl->actsum/(double)(mb_per_pict);
+       ratectl->sum_avg_act += ratectl->avg_act;
+       ratectl->actcovered = 0.0;
+
+       /* Allocate target bits for frame based on frames numbers in GOP
+          weighted by global complexity estimates and B-frame scale factor
+          T = (Nx * Xx/Kx) / Sigma_j (Nj * Xj / Kj)
+       */
+       ratectl->min_q = ratectl->min_d = INT_MAX;
+       ratectl->max_q = ratectl->max_d = INT_MIN;
+       switch (picture->pict_type)
+       {
+       case I_TYPE:
+               
+               /* There is little reason to rely on the *last* I-frame
+                  as they're not closely related.  The slow correction of
+                  K should be enough to fine-tune...
+               */
+
+               ratectl->d = ratectl->d0i;
+               avg_K = ratectl->avg_KI;
+               Si = (ratectl->Xi + 3.0*avg_K*ratectl->actsum)/4.0;
+               ratectl->T = ratectl->R/(1.0+ratectl->Np*ratectl->Xp*ratectl->Ki/(Si*ratectl->Kp)+ratectl->Nb*ratectl->Xb*ratectl->Ki/(Si*ratectl->Kb));
+
+               break;
+       case P_TYPE:
+               ratectl->d = ratectl->d0pb;
+               avg_K = ratectl->avg_KP;
+               Sp = (ratectl->Xp + avg_K*ratectl->actsum) / 2.0;
+               ratectl->T =  ratectl->R/(ratectl->Np+ratectl->Nb*ratectl->Kp*ratectl->Xb/(ratectl->Kb*Sp)) + 0.5;
+               break;
+       case B_TYPE:
+               ratectl->d = ratectl->d0pb;            // I and P frame share ratectl virtual buffer
+               avg_K = ratectl->avg_KB;
+               Sb = ratectl->Xb /* + avg_K * ratectl->actsum) / 2.0 */;
+               ratectl->T =  ratectl->R/(ratectl->Nb+ratectl->Np*ratectl->Kb*ratectl->Xp/(ratectl->Kp*Sb));
+               break;
+       }
+
+       /* Undershot bits have been "returned" via R */
+       if( ratectl->d < 0 )
+               ratectl->d = 0;
+
+       /* We don't let the target volume get absurdly low as it makes some
+          of the prediction maths ill-condtioned.  At these levels quantisation
+          is always minimum anyway
+       */
+       if( ratectl->T < 4000.0 )
+       {
+               ratectl->T = 4000.0;
+       }
+       target_Q = scale_quant(picture, 
+                                                  avg_K * ratectl->avg_act *(mb_per_pict) / ratectl->T);
+       current_Q = scale_quant(picture,62.0*ratectl->d / ratectl->r);
+#ifdef DEBUG
+       if( !quiet )
+       {
+               /* printf( "AA=%3.4f T=%6.0f K=%.1f ",avg_act, (double)T, avg_K  ); */
+               printf( "AA=%3.4f SA==%3.4f ",avg_act, sum_avg_act  ); 
+       }
+#endif 
+       
+       if ( current_Q < 3 && target_Q > 12 )
+       {
+               /* We're undershooting and a serious surge in the data_flow
+                  due to lagging adjustment is possible...
+               */
+               ratectl->d = (int) (target_Q * ratectl->r / 62.0);
+       }
+
+       ratectl->S = bitcount();
+       ratectl->frame_start = bitcount();
+//     ratectl->current_quant = ratectl->d * 62.0 / ratectl->r;
+       if(ratectl->current_quant < 1) ratectl->current_quant = 1;
+       if(ratectl->current_quant > 100) ratectl->current_quant = 100;
+}
+
+/* compute initial quantization stepsize (at the beginning of picture) */
+int ratectl_start_mb(ratectl_t *ratectl, pict_data_s *picture)
+{
+       double Qj;
+       int mquant;
+       
+       if(fixed_mquant) 
+               Qj = fixed_mquant;
+       else
+               Qj = ratectl->current_quant;
+//             Qj = ratectl->d * 62.0 / ratectl->r;
+
+       mquant = scale_quant( picture, Qj);
+       mquant = intmax(mquant, quant_floor);
+
+       return mquant;
+}
+
+void ratectl_update_pict(ratectl_t *ratectl, pict_data_s *picture)
+{
+       double X;
+       double K;
+       int64_t AP,PP;          /* Actual and padded picture bit counts */
+       int    i;
+       int    Qsum;
+       int frame_overshoot;
+       double avg_bitrate;
+       int last_size;
+       double new_weight;
+       double old_weight;
+       
+       if(fixed_mquant) return;
+       
+       AP = bitcount() - ratectl->S;
+       frame_overshoot = (int)AP-(int)ratectl->T;
+
+       /* For the virtual buffers for quantisation feedback it is the
+          actual under/overshoot that counts, not what's left after padding
+       */
+       ratectl->d += frame_overshoot;
+       
+       /* If the cummulative undershoot is getting too large (as
+          a rough and ready heuristic we use 1/2 video buffer size)
+          we start padding the stream.  Or, in the case of VBR,
+          we pretend we're padding but don't actually write anything!
+
+        */
+
+       if( ratectl->gop_undershoot-frame_overshoot > video_buffer_size/2 )
+       {
+               int padding_bytes = 
+                       ((ratectl->gop_undershoot - frame_overshoot) - video_buffer_size/2)/8;
+               if( quant_floor != 0 )  /* VBR case pretend to pad */
+               {
+                       PP = AP + padding_bytes;
+               }
+               else
+               {
+//                     printf( "PAD" );
+//                     alignbits();
+                       for( i = 0; i < padding_bytes/2; ++i )
+                       {
+//                             putbits(0, 16);
+                       }
+                       PP = bitcount() - ratectl->S;                   /* total # of bits in picture */
+               }
+               frame_overshoot = (int)PP - (int)ratectl->T;
+       }
+       else
+               PP = AP;
+       
+       /* Estimate cummulative undershoot within this gop. 
+          This is only an estimate because T includes an allocation
+          from earlier undershoots causing multiple counting.
+          Padding and an exact calculation each gop prevent the error
+          in the estimate growing too excessive...
+          */
+       ratectl->gop_undershoot -= frame_overshoot;
+       ratectl->gop_undershoot = ratectl->gop_undershoot > 0 ? ratectl->gop_undershoot : 0;
+       ratectl->R -= PP;                                               /* remaining # of bits in GOP */
+
+       Qsum = 0;
+       for( i = 0; i < mb_per_pict; ++i )
+       {
+               Qsum += picture->mbinfo[i].mquant;
+       }
+
+
+       ratectl->AQ = (double)Qsum/(double)mb_per_pict;
+       /* TODO: The X are used as relative activity measures... so why
+          bother dividing by 2?  
+          Note we have to be careful to measure the actual data not the
+          padding too!
+       */
+       ratectl->SQ += ratectl->AQ;
+       X = (double)AP*(ratectl->AQ/2.0);
+       
+       K = X / ratectl->actsum;
+#ifdef DEBUG
+       if( !quiet )
+       {
+               printf( "AQ=%.1f SQ=%.2f",  AQ,SQ);
+       }
+#endif
+       /* Bits that never got used in the past can't be resurrected
+          now...  We use an average of past (positive) virtual buffer
+          fullness in the event of an under-shoot as otherwise we will
+          tend to over-shoot heavily when activity picks up.
+        
+          TODO: We should really use our estimate K[IBP] of
+          bit_usage*activity / quantisation ratio to set a "sensible"
+          initial d to achieve a reasonable initial quantisation. Rather
+          than have to cut in a huge (lagging correction).
+
+          Alternatively, simply requantising with mean buffer if there is
+          a big buffer swing would work nicely...
+
+       */
+       
+       /* EXPERIMENT: Xi are used as a guesstimate of likely *future* frame
+          activities based on the past.  Thus we don't want anomalous outliers
+          due to scene changes swinging things too much.  Introduce moving averages
+          for the Xi...
+          TODO: The averaging constants should be adjust to suit relative frame
+          frequencies...
+       */
+       switch (picture->pict_type)
+       {
+       case I_TYPE:
+               ratectl->avg_KI = (K + ratectl->avg_KI * K_AVG_WINDOW_I) / (K_AVG_WINDOW_I+1.0) ;
+               ratectl->d0i = ratectl->d;
+               ratectl->Xi = (X + 3.0 * ratectl->Xi) / 4.0;
+               break;
+       case P_TYPE:
+               ratectl->avg_KP = (K + ratectl->avg_KP * K_AVG_WINDOW_P) / (K_AVG_WINDOW_P+1.0) ;
+               ratectl->d0pb = ratectl->d;
+               ratectl->Xp = (X + ratectl->Xp * 12.0) / 13.0;
+               ratectl->Np--;
+               break;
+       case B_TYPE:
+               ratectl->avg_KB = (K + ratectl->avg_KB * K_AVG_WINDOW_B) / (K_AVG_WINDOW_B + 1.0) ;
+               ratectl->d0pb = ratectl->d;
+               ratectl->Xb = (X + ratectl->Xb * 24.0) / 25.0;
+               ratectl->Nb--;
+               break;
+       }
+
+       ratectl->frame_end = bitcount();
+
+       last_size = ratectl->frame_end - ratectl->frame_start;
+       avg_bitrate = (double)last_size * frame_rate;
+       switch(picture->pict_type)
+       {
+               case I_TYPE:
+                       new_weight = avg_bitrate / bit_rate * 1 / N;
+                       old_weight = (double)(N - 1) / N;
+                       break;
+
+               default:
+               case P_TYPE:
+                       new_weight = avg_bitrate / bit_rate * (N - 1) / N;
+                       old_weight = (double)1 / N;
+                       break;
+       }
+       ratectl->current_quant *= (old_weight + new_weight);
+
+/*
+ * printf("ratectl_update_pict %f %f\n", 
+ *     avg_bitrate,
+ *     ratectl->current_quant);
+ */
+
+}
+
+/* Step 2: measure virtual buffer - estimated buffer discrepancy */
+int ratectl_calc_mquant(ratectl_t *ratectl, pict_data_s *picture, int j)
+{
+       int mquant;
+       double dj, Qj, actj, N_actj; 
+
+//     pthread_mutex_lock(&(ratectl->ratectl_lock));
+       /* A.Stevens 2000 : we measure how much *information* (total activity)
+          has been covered and aim to release bits in proportion.  Indeed,
+          complex blocks get an disproprortionate boost of allocated bits.
+          To avoid visible "ringing" effects...
+
+       */
+
+       actj = picture->mbinfo[j].act;
+       /* Guesstimate a virtual buffer fullness based on
+          bits used vs. bits in proportion to activity encoded
+       */
+
+
+       dj = ((double)ratectl->d) + 
+               ((double)(bitcount() - ratectl->S) - ratectl->actcovered * ((double)ratectl->T) / ratectl->actsum);
+
+
+
+       /* scale against dynamic range of mquant and the bits/picture count.
+          quant_floor != 0.0 is the VBR case where we set a bitrate as a (high)
+          maximum and then put a floor on quantisation to achieve a reasonable
+          overall size.
+          Not that this *is* baseline quantisation.  Not adjust for local activity.
+          Otherwise we end up blurring active macroblocks. Silly in a VBR context.
+       */
+
+       Qj = dj * 62.0 / ratectl->r;
+
+//printf("ratectl_calc_mquant %f\n", Qj);
+       if(fixed_mquant)
+               Qj = fixed_mquant;
+       else
+               Qj = ratectl->current_quant;
+
+       Qj = (Qj > quant_floor) ? Qj : quant_floor;
+       /*  Heuristic: Decrease quantisation for blocks with lots of
+               sizeable coefficients.  We assume we just get a mess if
+               a complex texture's coefficients get chopped...
+       */
+               
+       N_actj =  actj < ratectl->avg_act ? 
+               1.0 : 
+               (actj + act_boost * ratectl->avg_act) /
+               (act_boost * actj + ratectl->avg_act);
+   
+       mquant = scale_quant(picture, Qj * N_actj);
+
+
+       /* Update activity covered */
+
+       ratectl->actcovered += actj;
+//     pthread_mutex_unlock(&(ratectl->ratectl_lock));
+
+       return mquant;
+}
+
+/* VBV calculations
+ *
+ * generates warnings if underflow or overflow occurs
+ */
+
+/* vbv_end_of_picture
+ *
+ * - has to be called directly after writing picture_data()
+ * - needed for accurate VBV buffer overflow calculation
+ * - assumes there is no byte stuffing prior to the next start code
+ */
+
+void vbv_end_of_picture()
+{
+}
+
+/* calc_vbv_delay
+ *
+ * has to be called directly after writing the picture start code, the
+ * reference point for vbv_delay
+ */
+
+void calc_vbv_delay()
+{
+}
+
+void stop_ratectl(ratectl_t *ratectl)
+{
+       pthread_mutex_destroy(&(ratectl->ratectl_lock));
+}