/*****************************************************************************
 *                                                                          
 * This software module was originally developed by 
 *
 * J. Ignacio Ronda (UPM-GTI / ACTS-MoMuSys)
 *                                                      
 * and edited by                                        
 *                                                      
 * Martina Eckert       (UPM-GTI / ACTS-MoMuSys)             
 * Fernando Jaureguizar (UPM-GTI / ACTS-MoMuSys)
 * Jordi Ribas          (Sharp Labs of America): 
 * Tihao Chiang         (Sarnoff Corporation): 
 *
 * in the course of development of the MPEG-4 Video (ISO/IEC 14496-2) standard.
 * This software module is an implementation of a part of one or more MPEG-4 
 * Video (ISO/IEC 14496-2) tools as specified by the MPEG-4 Video (ISO/IEC 
 * 14496-2) standard. 
 *
 * ISO/IEC gives users of the MPEG-4 Video (ISO/IEC 14496-2) standard free 
 * license to this software module or modifications thereof for use in hardware
 * or software products claiming conformance to the MPEG-4 Video (ISO/IEC 
 * 14496-2) standard. 
 *
 * Those intending to use this software module in hardware or software products
 * are advised that its use may infringe existing patents. The original 
 * developer of this software module and his/her company, the subsequent 
 * editors and their companies, and ISO/IEC have no liability for use of this 
 * software module or modifications thereof in an implementation. Copyright is 
 * not released for non MPEG-4 Video (ISO/IEC 14496-2) Standard conforming 
 * products. 
 *
 * ACTS-MoMuSys partners retain full right to use the code for his/her own 
 * purpose, assign or donate the code to a third party and to inhibit third 
 * parties from using the code for non MPEG-4 Video (ISO/IEC 14496-2) Standard
 * conforming products. This copyright notice must be included in all copies or
 * derivative works. 
 *
 * Copyright (c) 1997
 *
 *****************************************************************************/

/***********************************************************HeaderBegin*******
 *                                                                         
 * File:        rc_q2.c
 *
 * Author:      J. Ignacio Ronda, UPM-GTI
 * Created:     18-06-97
 *                                                                         
 * Description: Q2 rate control functions
 *
 * Flags:       -D_RC_DEBUG_  -  RC debugging   
 *
 * Modified:
 *      16-10-97  Fernando Jaureguizar: Decreased the numFrame in
 *                RCQ2_ExcludeIFrame() instead of rcQ2_init().
 *      12.11.97  Martina Eckert: Headers, comments, cleaning
 *      18.11.97  M.Wollborn: include unistd only for non-PC
 *      27.11.97  M. Eckert: Changes for independent frame type control
 *      19.02.98  Jordi Ribas: added RCQ2_MB_ routines for VM8 macroblock rate control
 *
 ***********************************************************HeaderEnd*********/

/************************    INCLUDE FILES    ********************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#if !defined(WIN32)
#  include <unistd.h>
#endif

#include <ctype.h>

#include "rc.h"
#include "momusys.h"
#include "mom_vop.h"
#include "vm_config.h"
#include "mom_vol.h"
#include "mom_structs.h"
#include "vm_enc_defs.h"


#include <errno.h>


/* General configuration: */
/* static Int rc_gen_buff_size;  Global buffers */
/* static Int rc_vol_buff_size[MAX_NUM_VOS][MAX_NUM_VOLS];  Per-VOL buffers */

/* VOLs history and configuration, and model parameters: */
static RC_HIST   rc_hist[MAX_NUM_VOS][MAX_NUM_VOLS][3]; /* History of each VOL
							   for B,I,P frames */
static RCQ2_DATA RCQ2_config[MAX_NUM_VOS][MAX_NUM_VOLS];/* Q2 data of each VOL*/

/* --------------------------------- */
/* SHARP_start                       */

static RCQ2_MB_DATA   Q2_mb; /* data for Q2 macroblock rate control */

/* SHARP_end                        */
/* -------------------------------  */

/*static Int  RCQ2_SepMotShapeText;*/
/* inserted by Mistubishi & SAMSUNG AIT */
Int QPfirst[MAX_NUM_VOS][MAX_NUM_VOLS];
/* inserted by Mistubishi & SAMSUNG AIT */

/***********************************************************HeaderBegin*******
 * LEVEL 1 
 ***********************************************************HeaderEnd*********/

/***********************************************************CommentBegin******
 *
 * -- RCQ2_init --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      General initialization
 *
 * Arguments in :       
 *      VOConfig *list    -  list of VO
 *
 * Modified : 
 *      27.11.97  M. Eckert: Changes for independent frame type control
 *      05.12.97  M. Eckert: Rename of function
 *
 ***********************************************************CommentEnd********/

void RCQ2_init(
	       VOConfig *list
	       )
{
  VOConfig  *vo_config;
  VolConfig *vol_config;
  Int qp_first[3];

  vo_config = list;

  while (vo_config != NULL)
    {
    Int vo_id = GetVOConfigId(vo_config);
    vol_config = GetVOConfigLayers(vo_config);

    /*    RCQ2_SepMotShapeText     = GetVolConfigSepMotShapeText(vol_config);*/

    while (vol_config != NULL)
       {
       Int vol_id = GetVolConfigId(vol_config);

#ifdef _RC_DEBUG_
       fprintf(stdout, "RC: ----->RCQ2_init_global: vol_id= %d\n", vol_id);
       fprintf(stdout, "RC: StartFrame= %d\n",
                            GetVolConfigStartFrame(vol_config));
       fprintf(stdout, "RC: EndFrame= %d\n",
                            GetVolConfigEndFrame(vol_config));
       fprintf(stdout, "RC: FrameSkip= %d\n",
                            GetVolConfigFrameSkip(vol_config));
       fprintf(stdout, "RC: Bitrate= %d\n",
                            GetVolConfigBitrate(vol_config));
       fprintf(stdout, "RC: FrameRate= %f\n",
                            GetVolConfigFrameRate(vol_config));
#endif /* _RC_DEBUG_ */

       qp_first[0] = GetVolConfigIntraQuantizer(vol_config);
       qp_first[1] = GetVolConfigQuantizer(vol_config);
       qp_first[2] = GetVolConfigBQuantizer(vol_config);

       RCQ2_init_vol (vo_id, vol_id, vol_config, qp_first);

       vol_config = GetVolConfigNext(vol_config);

       }
    vo_config = GetVOConfigNext(vo_config);
    }
}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_Free --
 *
 * Author : 
 *      M. Eckert, UPM-GTI     
 *
 * Created :            
 *      19-02-98
 *
 * Purpose :
 *      Free initialized memory
 *
 * Arguments in :       
 *      VOConfig *list    -  list of VO
 *
 * Modified : 
 *
 ***********************************************************CommentEnd********/

void RCQ2_Free(
	       VOConfig *vo_list
	       )
{
  VolConfig *vol_list;

  while (vo_list != NULL)
  {
    Int vo_id = GetVOConfigId(vo_list);
    vol_list = GetVOConfigLayers(vo_list);
    while (vol_list != NULL)
    {
      Int vol_id = GetVolConfigId(vol_list);
      RCQ2_DATA *rcd = &RCQ2_config[vo_id][vol_id];

      free(rcd->X1);
      free(rcd->X2);
      free(rcd->QPfirst);

      rch_free(rc_hist[vo_id][vol_id]); /* OJO! ok? */

      vol_list = GetVolConfigNext(vol_list);
    }
    vo_list = GetVOConfigNext(vo_list);
    
  }
}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_init_vol --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Initialization of Q2 rate control
 *
 * Arguments in :       
 *      Int      vo_id        -  VO Id
 *      Int      vol_id       -  VOL Id
 *      Int      start        -  First input frame to process
 *      Int      end          -  Last input frame to process
 *      Int      intra_period -  # of frames between I frames +1
 *      Int      M            -  # of B frames between P frames +1
 *      Int      frame_inc    -  Frame temporal subsampling factor
 *      Float    frame_rate   -  Input frame rate (frames per second)
 *      Int      bit_rate     -  Individual bit-rate (bits per second) 
 *      Int      qp_first     -  QP to use the first time RC is employed  
 *
 * Modified : 
 *      27.11.97  M. Eckert: Changes for independent frame type control
 *
 *
 ***********************************************************CommentEnd********/

void RCQ2_init_vol(
                   Int   vo_id,
                   Int   vol_id,
		   VolConfig *vol_config,
                   Int   *qp_first
                   )
{
  rcQ2_init(&RCQ2_config[vo_id][vol_id], vol_config, qp_first);
  rch_init(rc_hist[vo_id][vol_id], MAX_SLIDING_WINDOW);
}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_QuantAdjust --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Quantizer obtainment: Called after mad is known.
 *
 * Arguments in :       
 *      Int      vo_id        -  VO Id
 *      Int      vol_id       -  VOL Id
 *      Double   mad          -  Pred. error mad
 *      UInt     num_pels_vop -  # pels per VOP
 *      Int      frame_type   -  I=0, P=1, B=2
 *
 * Modified : 
 *      27.11.97  M. Eckert: Changes for independent frame type control
 *
 *
 ***********************************************************CommentEnd********/

Int RCQ2_QuantAdjust(
                     Int    vo_id,
                     Int    vol_id,
                     Double mad,       /* Pred. error mad */
                     UInt   num_pels_vop,
		     Int    frame_type
                     )
{
  Int qp, bits_P, bits_frame=0;
  RC_HIST   *rch=&rc_hist[vo_id][vol_id][frame_type]; /* Histories for frame type */
  RCQ2_DATA *rcd=&RCQ2_config[vo_id][vol_id];


/* Computation of QP */
  if (rch_first(rch))
    qp = rcQ2_get_QPfirst(rcd, frame_type);
  else
  {
    /* RC-BIP: Calculate bits_frame (target0) depending on frame type */

    bits_P = (Int)( rcd->R / ( rcd->alpha_I * rcd->num_I_Frame
			+ rcd->alpha_B * rcd->num_B_Frame
			+ rcd->num_P_Frame));
    
    switch (frame_type)    
      {
      case I_VOP:
	printf("Coding I_frame \n");
	bits_frame = (Int)( rcd->alpha_I * bits_P);
	break;
      case B_VOP:
	printf("Coding B_frame \n");
	bits_frame = (Int)( rcd->alpha_B * bits_P);
	break;
      case P_VOP:
	printf("Coding P_frame \n");
	bits_frame = bits_P;
	break;
      default:
	error_exit("RCQ2_QuantAdjust: Error: frame_type? \n");
      }
    printf("RC-BIP ------> a_I = %f, a_B = %f, #I = %d, #B = %d, #P = %d,  bits_P = %d \n",
	   rcd->alpha_I, rcd->alpha_B, rcd->num_I_Frame, rcd->num_B_Frame, rcd->num_P_Frame, bits_P );
    
#ifdef _RC_DEBUG_
    fprintf(stdout, "RCQ2_QuantAdjust: R= %f  numFrame= %d  target0= %f\n",
	    (double)(rcd->R), rcd->numFrame, (double)rcd->R/rcd->numFrame); 
#endif


    qp=RCQ2_compute_quant(rcd->bit_rate,
	     /* RC-BIP */ bits_frame, /* rcd->R/rcd->numFrame,*/
	     /* RC-BIP */ rch_get_n(rch)>1 ? rch_get_last_bits_vop(rch) : 0,
			  rcd->VBV_size,
			  rcd->VBV_fullness,
	     /* RC-BIP */ rcd->X1[frame_type],
	     /* RC-BIP */ rcd->X2[frame_type],
			  mad,
	     /* RC-BIP */ rch_get_last_qp(rch)
			  );

  }
  
  /* Storage of result */
  rch_store_before(rch, (Int) num_pels_vop, mad, qp);
  
  return qp;
}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_FrameSkipping --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Returns 1 if next frame should be skipped, 0 otherwise.
 *      Called before deciding what VOL to code next.
 *
 * Arguments in :       
 *      Int      vo_id        -  VO Id
 *      Int      vol_id       -  VOL Id
 *
 * Modified : 
 *
 *
 ***********************************************************CommentEnd********/

Int RCQ2_FrameSkipping(
                       Int vo_id,
                       Int vol_id
                       )
{
   RCQ2_DATA *rcd=&RCQ2_config[vo_id][vol_id];

#ifdef _RC_DEBUG_
   Char b_file_name[15];

   sprintf(b_file_name,"b_script_%d_%d",vo_id,vol_id);
   
   WRITE_INT(b_file_name, rcd->VBV_fullness);  
#endif


   if (rcd->VBV_fullness
       > SKIP_MARGIN*rcd->VBV_size)
      {
#ifdef _RC_DEBUG_
      fprintf(stdout,
         "RC: RCQ2_FrameSkipping: VBV fullness (vo= %d vol= %d): %d of %d\n",
         vo_id, vol_id, (int)rcd->VBV_fullness, (int)rcd->VBV_size);
      fflush(stdout);
#endif
      return 1;
      }
   else
      return 0;
}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_UpdateAfterSkipping --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Updates control algorithm parameters after deciding that next frame
 *      should be skipped.
 *
 * Arguments in :       
 *      Int      vo_id        -  VO Id
 *      Int      vol_id       -  VOL Id
 *      Int      frame_type   -  I=0, P=1, B=2
 *
 * Modified : 
 *      27.11.97  M. Eckert: Changes for independent frame type control
 *
 *
 ***********************************************************CommentEnd********/

void RCQ2_UpdateAfterSkipping(
                              Int vo_id, 
                              Int vol_id,
			      Int frame_type
                              )
{
   RCQ2_DATA *rcd=&RCQ2_config[vo_id][vol_id];

   rcd->numFrame --;

   if (frame_type == I_VOP) rcd->num_I_Frame --;
   else if (frame_type == B_VOP) rcd->num_B_Frame --;
   else if (frame_type == P_VOP) rcd->num_P_Frame --;
   else error_exit("ERROR RCQ2_UpdateAfterSkipping: Unknown frame type ! \n");

      if (frame_type == P_VOP)
	 rcd->VBV_fullness -= (Int) rcd->Bpp;
      else
	 if (frame_type == I_VOP)
	    rcd->VBV_fullness -= (Int) (rcd->alpha_I * rcd->Bpp);
	 else if (frame_type == B_VOP)
	    rcd->VBV_fullness -= (Int) (rcd->alpha_B * rcd->Bpp);

}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_ExcludeIFrame --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Updates control algorithm parameters after encoding the first I frame.
 *
 * Arguments in :       
 *      Int      vo_id        -  VO Id
 *      Int      vol_id       -  VOL Id
 *      Int      vopbits      -  Total number of bits of I frame
 *
 * Modified : 
 *     16-10-97  Fernando Jaureguizar: Decreased the numFrame in
 *               RCQ2_ExcludeIFrame() instead of rcQ2_init().
 *     27.11.97  M. Eckert: Changes for independent frame type control
 *
 ***********************************************************CommentEnd********/

void RCQ2_ExcludeIFrame(
                        Int vo_id,
                        Int vol_id,
                        Int vopbits  /* Total number of bits for I frame */
                        )
{
   RCQ2_DATA *rcd=&RCQ2_config[vo_id][vol_id];

   rcd->R -= vopbits;
   rcd->numFrame--;
   rcd->num_I_Frame--; /* RC-BIP */
   rcd->Bpp = rcd->R/rcd->numFrame;
#ifdef _RC_DEBUG_
   fprintf(stdout,"RCQ2_ExcludeIFrame: rcd->Bpp=%d (R=%d numFrame=%d)\n",
           (int)rcd->Bpp, (int)rcd->R, (int)rcd->numFrame);
#endif
   if (rcd->Bpp <= 0)
      {
      fprintf(stdout,
            "RCQ2_ExcludeIFrame: All bits spent in the I_VOP (aborting).\n");
      exit(1);
      }
   
}

/***********************************************************CommentBegin******
 *
 * -- rc_update_Q2_model --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Obtainment of Q2 model
 *
 * Arguments in :       
 *      RC_HIST     *rch     -  VO Id
 *
 * Arguments in/out :       
 *      Double      *X1      -  Old and computed parameter
 *      Double      *X2      -  Old and computed parameter
 *
 * Modified : 
 *
 *
 ***********************************************************CommentEnd********/

void rc_update_Q2_model(
                        RC_HIST *rch,
                        Double  *X1,   /* <-> Old and computed params. */
                        Double  *X2
                        )
{
Double RE[MAX_SLIDING_WINDOW];     /* R/E data */
Double invQP[MAX_SLIDING_WINDOW];  /* Inverse QP   */
Double invQP2[MAX_SLIDING_WINDOW]; /* Inverse QP^2 */

Int    n;                          /* Number of data to employ */
Double last_mad;
Double plast_mad=0;
Int    i, count;

Int *bits_text=rch_get_bits_text(rch);
Double *mad_text = rch_get_mad_text(rch);
Int *qp = rch_get_qp(rch);



#ifdef _RC_DEBUG_
   fprintf(stdout, "RC: RC_Q2: --> rc_update_Q2_model\n");
   fflush(stdout);
   rch_print(rch, stdout, "rc_update_Q2_model");
#endif

   if (rch_get_n(rch)<START_MODEL_UPDATING)
     return;

   last_mad =rch_get_last_mad_text(rch);
   if (rch_get_n(rch)>1)
      plast_mad=rch_get_plast_mad_text(rch);

   /* Computation of number of data to employ */
   n=MIN(rch_get_n(rch),MAX_SLIDING_WINDOW);

#ifdef _RC_DEBUG_
   fprintf(stdout, "rc_update_Q2_model: n= %d\n", n);
#endif
   if (rch_get_n(rch)>1)
     n = (plast_mad/last_mad > 1.0) ? (Int)((last_mad/plast_mad)*n+ 1)
                             : (Int)((plast_mad/last_mad)*n + 1);
   else
     n = 1;

#ifdef _RC_DEBUG_
   fprintf(stdout, "rc_update_Q2_model: plastmad= %f  lastmad= %f  n= %d\n",  plast_mad, last_mad, n);
#endif

   n=MIN(n,MIN(rch_get_n(rch),MAX_SLIDING_WINDOW));

/* Arrays initialization */
   for (i=rch_dec_mod(rch,rch_get_ptr(rch)),count=0; count<n; count++,i=rch_dec_mod(rch,i))
   {
     /* As the VM says, model only takes into account texture bits, although it is
        used to estimate the total number of bits !!! */
     RE[count]=(Double)bits_text[i]/mad_text[i]; 

     invQP[count] = 1./(Double)qp[i];
     invQP2[count] = invQP[count] * invQP[count];

#ifdef _RC_DEBUG_log_sgicc_dbg

     fprintf(stdout, "rc_update_Q2_model: %d: RE= %f  IQ= %f  IQ2= %f\n", 
                          count, RE[count], invQP[count], invQP2[count]);
#endif
   }

   /* LS computation */
   rc_2ls_comp_q2(RE, invQP, n, X1, X2);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RC: RC_Q2: <-- rc_2ls_comp_q2: X1= %f  X2= %f \n", *X1, *X2);
#endif


   if (n >= START_OUTLIER_REMOVING)
      {
      /* Outlier removing (arrays reinitialization) */
      rc_2ls_remove_out(RE, invQP, invQP2, &n, X1, X2);

      /* LS recomputation */
      rc_2ls_comp_q2(RE, invQP, n, X1, X2);
      }
   else
      printf("No he entrado en rc_2ls_remove_out: n= %d\n", n);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RC: RC_Q2: <-- rc_update_Q2_model: X1= %f  X2= %f  \n", *X1, *X2);
#endif
}

/***********************************************************CommentBegin******
 *
 * -- RCQ2_Update2OrderModel --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Updates model after encoding a frame.
 *
 * Arguments in :       
 *      Int      vo_id        -  VO Id
 *      Int      vol_id       -  VOL Id
 *      Int      vopbits      -  Total number of bits
 *      Int      vopheader    -  Total minus texture
 *      Double   dist         -  Distortion
 *      Int      frame_type   -  I=0, P=1, B=2
 *
 * Modified : 
 *
 *      19.02.98  Jordi Ribas: added call to update VM8 macroblock rate control
 *
 ***********************************************************CommentEnd********/

void RCQ2_Update2OrderModel(
                            Int vo_id,
                            Int vol_id,
                            Int vopbits,   /* Total number of bits */
                            Int vopheader, /* Total minus texture */
                            Double dist,
			    Int frame_type
                            )
{
   RCQ2_DATA *rcd=&RCQ2_config[vo_id][vol_id];
   RC_HIST   *rch=&rc_hist[vo_id][vol_id][frame_type];

   /* ------------------------------------------ */
   /* SHARP_start                                */
   if (Q2_mb.active_vop == Q2MB_ACTIVE)
     {
       /* If the Q2 VM8 macroblock rate control was used for the last VOP
          should only do this */

       rch->total_dist      += dist;
       rch->total_bits_text += vopbits-vopheader;
       rch->total_bits      += vopbits;
       rch->total_frames++;
       rch->prev= 0;

       Q2_mb.bits_prev= vopbits;
       rcQ2_update(rcd, rch, vopbits, frame_type);
       fprintf(stdout,"RCQ2_Update2OrderModel: Buffer fullness: %d \n\n",
               (Int) rcd->VBV_fullness);

       Q2_mb.active_vop= Q2MB_INACTIVE;
     }
   /* SHARP_end                                */
   /* ------------------------------------------ */
   else
     { 
       rch_store_after(rch, vopbits-vopheader, vopbits, dist);
       
       rcQ2_update(rcd, rch, vopbits, frame_type);
       
       rc_update_Q2_model(rch, &rcd->X1[frame_type], &rcd->X2[frame_type]);
     }
}

/* ------------------------------------------------------------------------- */
/* SHARP_start                                                               */

/***********************************************************CommentBegin******
 *
 * -- RCQ2_MB_active --
 *
 * Author : 
 *      Jordi Ribas, Sharp Labs of America, jordi@sharplabs.com
 *
 *      The macroblock-layer rate control is based on the technology 
 *      jointly developed by Sarnoff Corporation and Sharp Labs 
 *      of America.
 *
 * Created :            
 *      19-02-98
 *
 * Purpose :
 *      Can be used to turn the Q2 macroblock rate control active or inactive.
 *      Initializes the A1, A2, A3, parameters in the Q2_mb structure.
 *      (The Q2_mb structure is static in rc_q2.c)
 *      It also returns current state: 1 if active and 0 if not active
 *
 * Arguments in :  
 *      Int     turn_active    -  0   : turn inactive
 *                                1   : turn active
 *                                else: any other number does not change state 
 *
 * Modified : 
 *
 ***********************************************************CommentEnd********/

Int RCQ2_MB_active(Int turn_active)
{
  
  if (turn_active == 0) 
    {
      Q2_mb.active=0;
#ifdef _RC_DEBUG_
      fprintf(stdout, "RCQ2_MB_active: macroblock rate control is off \n");
#endif
    }
  else 
    if (turn_active == 1) 
      {
        Q2_mb.active=     Q2MB_ACTIVE;
        Q2_mb.active_vop= Q2MB_INACTIVE;  /* this is only activated at VOP level, in Q2MB_init */
        Q2_mb.A1=         100.0;
        Q2_mb.A1_prev=    100.0;
        Q2_mb.A1_curr=      0.0;
        Q2_mb.A2=         400.0;
        Q2_mb.A3=           0.0;
        Q2_mb.T_sms=        0.0;
        Q2_mb.bits_prev=    0.0;
#ifdef _RC_DEBUG_
        fprintf(stdout, "RCQ2_MB_active: macroblock rate control is on \n");
#endif
      }
  
  return (Q2_mb.active);
}


/***********************************************************CommentBegin******
 *
 * -- RCQ2_MB_init --
 *
 * Author : 
 *      Jordi Ribas, Sharp Labs of America, jordi@sharplabs.com
 *
 * Created :            
 *      19-02-98
 *
 * Purpose :
 *      Initializes the data structure used for Q2 macroblock rate control
 *      at the beginning of each VOP.  The data structure is Q2_mb and is 
 *      static in rc_q2.c. There are four steps: 
 *
 *      1) check whether VOP is rectangular and P (currently only this is supported), 
 *      2) compute the target number of bits for the VOP (as in Q2 rate control), 
 *      3) computation of the macroblock perceptual weights, 
 *      3) computation of the MADs for the VOP macroblocks.
 *
 * Arguments in :  
 *      Vop     *error_vop      -  pointer to the VOP's prediction error
 *      Int     vo_id           -  VO Id
 *      Int     vol_id          -  VOL Id
 *      Int     MB_h            -  number of macroblocks in VOP along horizontal 
 *      Int     MB_v            -  number of macroblocks in VOP along vertical
 * 
 * Returns: the value of QP for the VOP layer and first macroblock
 *
 * Modified : 
 *
 ***********************************************************CommentEnd********/

Int RCQ2_MB_init(Vop *error_vop, Int vo_id, Int vol_id, Int MB_h, Int MB_v)
{
  RCQ2_DATA *rcd=&RCQ2_config[vo_id][vol_id];
  Int j,i, m,n;
  SInt  *curr_in=0,  *tmp_in;
  Float *curr_fin=0, *tmp_fin;
  Int    pos, row_pos, type, line;
  Double mad=0.0, mean= 0.0;

  /* variables for computing target bits */
  Double target, a, b;
  Double bit_rate=   rcd->bit_rate;             /* Bit rate (bits per second) */
  Int    bits_frame= (rcd->R/rcd->numFrame);    /* Target bits per frame when no buffer */
  Int    bits_prev;                             /* Total bits previous frame */
  Int    buff_size=  rcd->VBV_size;             /* Buffer size */
  Int    buff_occup= rcd->VBV_fullness;         /* Buffer occupancy */                      

  /* STEP 1: check proper mode for macroblock rate control is in place.
     Currently, macroblock rate control is only implemented for P and rectangular VOPs */
  if (error_vop->prediction_type != P_VOP)
    {
      /* return to default: VM5 rate control mode and buffer */
      Q2_mb.active=      Q2MB_INACTIVE;
      Q2_mb.active_vop=  Q2MB_INACTIVE; 

      rcd->VBV_size = rcd->bit_rate/2;        
      rcd->VBV_fullness += (rcd->bit_rate/4 - rcd->bit_rate/16);  

#ifdef _RC_DEBUG_
        fprintf(stdout, "RCQ2_MB_init: Q2 VM8 macroblock rate control off. \n");
        fprintf(stdout, "              Currently only P VOPs are supported. \n");
        fprintf(stdout, "              Q2 VM5 rate control with default buffer used instead. \n");
#endif
      return (0);
    }
  
  if (error_vop->arbitrary_shape != RECTANGULAR)
    {
      /* return to default: VM5 rate control mode and buffer */
      Q2_mb.active=      Q2MB_INACTIVE;
      Q2_mb.active_vop=  Q2MB_INACTIVE; 

      rcd->VBV_size = rcd->bit_rate/2;        
      rcd->VBV_fullness += (rcd->bit_rate/4 - rcd->bit_rate/16);  

#ifdef _RC_DEBUG_
        fprintf(stdout, "RCQ2_MB_init: Q2 VM8 macroblock rate control off. \n");
        fprintf(stdout, "              Currently only rectangular VOPs are supported. \n");
        fprintf(stdout, "              Q2 VM5 rate control with default buffer used instead. \n");
#endif
      return (0);  
    }
 
  Q2_mb.active_vop=  Q2MB_ACTIVE; 

  /* STEP 2: Computing the Target number of bits for the VOP. 
     Done equivalently as in Q2 rate control */
  bits_prev= (Int) Q2_mb.bits_prev;
  target= (Double) bits_frame;
  target= floor(target);
 
  if (bits_prev > 0)
    target = target*(1.-PAST_PERCENT)+bits_prev*PAST_PERCENT;
  target = floor(target);

  a = buff_occup;
  b = buff_size-buff_occup;
  target = (a+2.0*b)*target/(2.0*a+b);
  target = floor(target);

  target = MAX(bit_rate/30, target);
  target = floor(target);

  if (buff_occup+target > (1.-SAFETY_MARGIN)*buff_size)
    target = MAX(bit_rate/30, (1.-SAFETY_MARGIN)*buff_size-buff_occup);
  else if (buff_occup-bits_frame+target < SAFETY_MARGIN*buff_size)
    target = SAFETY_MARGIN*buff_size-buff_occup+bits_frame;
  target = floor(target);
  
  Q2_mb.T_bits= (Float) target;
  Q2_mb.T_text=  Q2_mb.T_bits - Q2_mb.T_sms;
  
  if (Q2_mb.T_bits < 0) Q2_mb.T_bits= 0;
  if (Q2_mb.T_text < 0) Q2_mb.T_text= 0;
#ifdef _RC_DEBUG_
  fprintf(stdout,"RCQ2_MB_init: Target for macroblock rate control: %.1f - %.1f = %.1f \n", 
          Q2_mb.T_bits, Q2_mb.T_sms, Q2_mb.T_text);
#endif
  /* reset value of bits for shape, motion and syntax */
  Q2_mb.T_sms = 0;
  
  if (Q2_mb.T_text > Q2MB_MODE_THR*((Float) MB_SIZE*MB_SIZE*MB_v*MB_h))  Q2_mb.mode= 1;
  else Q2_mb.mode= 2;

#ifdef _RC_DEBUG_
  fprintf(stdout,"RCQ2_MB_init: T_text= %.1f, Threshold: %.1f, mode used = %d", 
          Q2_mb.T_text, Q2MB_MODE_THR*((Float) MB_SIZE*MB_SIZE*MB_v*MB_h),  Q2_mb.mode);
#endif

  /* STEP 3: Memory allocation and computation of macroblocks' perceptual weights */
  Q2_mb.weight= (Float **)  calloc(MB_v,sizeof(Float *));
  for (j=0; j < MB_v; j++) 
    {
      Q2_mb.weight[j]= (Float *) calloc(MB_h,sizeof(Float));
      for (i=0; i < MB_h; i++)
        Q2_mb.weight[j][i]= 1.0;        /* weights are chosen to be 1 by default */
    }
  
  /* STEP 4: Memory allocation and computation of mad for each macroblock */
  Q2_mb.mad= (Float **) calloc(MB_v,sizeof(Float *));
  
  type= GetImageType(error_vop->y_chan);
  switch(type)
    {
    case SHORT_TYPE:
      curr_in = GetImageIData(error_vop->y_chan); 
      break;
    case FLOAT_TYPE:
      curr_fin = GetImageFData(error_vop->y_chan);
      break;
    }

  pos = 0;             /* pixel position at the top-left corner of (j,i) macroblock */
  row_pos= 0;          /* pixel position at the top-left corner of (j,0) macroblock */ 
  line= MB_h*MB_SIZE;  /* number of pixels in a line                                */
  Q2_mb.sum_mad= 0;
  Q2_mb.energy_S= 0;

#ifdef _RC_DEBUG_
  fprintf(stdout,"\nRCQ2_MB_init: MAD values for the macroblocks:");
#endif

  for (j=0; j < MB_v; j++) 
    {
#ifdef _RC_DEBUG_
      printf("\n");
#endif
      Q2_mb.mad[j]= (Float *) calloc(MB_h,sizeof(Float));  
      pos= row_pos;
      for (i=0; i < MB_h; i++)
        {
          mean= 0;
          mad= 0;
          switch (type)
            {
            case SHORT_TYPE:
              tmp_in= curr_in + pos;
              for (m=0; m < MB_SIZE; m++)
                {
                  tmp_in= curr_in + pos + m*line;
                  for (n=0; n < MB_SIZE; n++)
                    {
                      mean += *tmp_in;
                      tmp_in++;
                    }
                }
              mean /= (Double) MB_SIZE*MB_SIZE;

              tmp_in= curr_in + pos;
              for (m=0; m < MB_SIZE; m++)
                {
                  tmp_in= curr_in + pos + m*line;
                  for (n=0; n < MB_SIZE; n++)
                    {
                      mad += fabs(mean - (Double) *tmp_in);
                      tmp_in++;
                    }
                }
              break;
            case FLOAT_TYPE:
              tmp_fin = curr_fin + pos;
              for (m=0; m < MB_SIZE; m++)
                {
                  tmp_fin= curr_fin + pos + m*line;
                  for (n=0; n < MB_SIZE; n++)
                    {
                      mean += *tmp_fin;
                      tmp_fin++;
                    }
                }
              mean /= (Double) MB_SIZE*MB_SIZE;
              
              tmp_fin = curr_fin + pos;
              for (m=0; m < MB_SIZE; m++)
                {
                  tmp_fin= curr_fin + pos + m*line;
                  for (n=0; n < MB_SIZE; n++)
                    {
                      mad += fabs(mean - (Double) *tmp_fin);
                      tmp_fin++;
                    }
                }
              break;
            default: 
              break;
            }
          
          mad /= (Double) MB_SIZE*MB_SIZE;
          pos += MB_SIZE;
          Q2_mb.mad[j][i]= (Float) mad;

#ifdef _RC_DEBUG_
          fprintf(stdout," %4.1f", Q2_mb.mad[j][i]);
#endif
          Q2_mb.sum_mad   +=  Q2_mb.mad[j][i];
          Q2_mb.energy_S  +=  Q2_mb.mad[j][i]*Q2_mb.weight[j][i];
        }
      row_pos += MB_h*MB_SIZE*MB_SIZE;
    }

  /* finally we also initialize the fields MB_h, MB_v, and MB_skip */
  Q2_mb.MB_h=    MB_h;
  Q2_mb.MB_v=    MB_v;
  Q2_mb.MB_skip= 0;

#ifdef _RC_DEBUG_
  fprintf(stdout,"\nRCQ2_MB_init: sum_mad: %.2f, energy_S: %.2f", Q2_mb.sum_mad, Q2_mb.energy_S);
  fprintf(stdout,"\nRCQ2_MB_init: A1: %.1f = A1_prev: %.1f ", Q2_mb.A1, Q2_mb.A1_prev);
  fprintf(stdout,"\nRCQ2_MB_init: A2: %.1f, A3: %.1f \n", Q2_mb.A2, Q2_mb.A3);
#endif

  /* returns the QP quantizer value for the first macroblock */
  return (RCQ2_MB_compute_QP(0,0));

}


/***********************************************************CommentBegin******
 *
 * -- RCQ2_MB_compute_QP --
 *
 * Author : 
 *      Jordi Ribas, Sharp Labs of America, jordi@sharplabs.com
 *
 * Created :            
 *      19-02-98
 *
 * Purpose :
 *      Returns a quantization parameter QP for a given macroblock
 *
 * Arguments in :
 *      Float   index_y         -  vertical index of macroblock
 *      Float   index_x         -  horizontal index of macroblock
 *
 * Modified : 
 *
 ***********************************************************CommentEnd********/

Int RCQ2_MB_compute_QP(Int index_y, Int index_x)
{
  Float mad=      Q2_mb.mad[index_y][index_x];
  Float w=        Q2_mb.weight[index_y][index_x];
  Float energy=   Q2_mb.energy_S;
  Float num_bits= Q2_mb.T_text;
  Float A1=       Q2_mb.A1;
  Float A2=       Q2_mb.A2;
  Float A3=       Q2_mb.A3;
  Int   mode=     Q2_mb.mode;

  Int QP;  
  Float Q, a, b, c;
  Float target_mb;

  if (num_bits <= 0) return(31);
  
  target_mb= num_bits*mad*w/energy;
  if (mode == 1)  
    {
      Q=  mad*sqrt(A1/target_mb);
      QP= (Int) (Q + 0.5);
    }
  else  
    {
      a= - target_mb;
      b= A3*mad;
      c= A2*mad;

      if ((b*b - 4*a*c) > 0)
        Q= (-b - sqrt(b*b - 4*a*c))/(2*a);
      else { /* VM5 fall back mode */ 
        Q = mad * Q2_mb.X11 / target_mb ; 
      }
      
      QP= (Int) (Q + 0.5);
    }

  if (QP < 1) QP=1;
  if (QP > 31) QP= 31;

  return(QP);
}


/***********************************************************CommentBegin******
 *
 * --  RCQ2_MB_update --
 *
 * Author : 
 *      Jordi Ribas, Sharp Labs of America, jordi@sharplabs.com
 *
 * Created :            
 *      19-02-98
 *
 * Purpose :
 *      Updates the model parameters A1, A2, A3 used for Q2 macroblock 
 *      rate control. If it is the last macroblock of the VOP, it also
 *      deallocates memory and prepares for next VOP.
 *
 * Arguments in :  
 *
 *      Float   bits_MB         - bits used for luma and chroma DCT for current macroblock
 *      Float   bits_sms        - total bits used for motion, shape and syntax so far in the frame
 *      Float   QP              - value of QP used for encoding the macroblock
 *      Int     index_y         - vertical index of macroblock 
 *      Int     index_x         - horizontal index of macroblock
 *
 * Modified : 
 *
 *
 ***********************************************************CommentEnd********/

Void RCQ2_MB_update(Float bits_MB, Float bits_sms, Float QP, Int index_y, Int index_x)
{
  Float  mad=       Q2_mb.mad[index_y][index_x];       
  Float *A1=       &Q2_mb.A1;      
  Float *A1_curr=  &Q2_mb.A1_curr;  
  Float  A1_prev=   Q2_mb.A1_prev;  
  Float *A2=       &Q2_mb.A2;      
  Float *A3=       &Q2_mb.A3;       
  Int    mode=      Q2_mb.mode;
  
  Float MB_total=   (Float) (Q2_mb.MB_v*Q2_mb.MB_h);
  Float MB_coded=   (Float) (index_y*Q2_mb.MB_h + index_x + 1);
  Float MB_noskp=   MB_coded - (Float) Q2_mb.MB_skip;
  Float A1_mb=0;
  Int j;

  static  Int   N_ls= 0;
  Int     k;
  Double  invX[2][2], corrXY[2];
  Double  a=0, b=0, c=0, d=0, determ=0;

  if (mode == 1)  
    {
      if (bits_MB > 0)
        {
          /* Updating A1*/
          A1_mb= bits_MB*QP*QP/(mad*mad);
          *A1_curr= (*A1_curr)*(MB_noskp-1)/MB_noskp + A1_mb/MB_noskp;
        }
      else Q2_mb.MB_skip++;
      *A1= (*A1_curr)*MB_coded/MB_total + A1_prev*(MB_total-MB_coded)/MB_total;
    }   
  else  
    {
      /* COMMENT TO SARNOFF:  the estimation of A2 and A3 is here. 
         Please, check this carefully and let me know if you change anything.
         From discussions with Tihao in Stockholm and Fribourg, 
         I understand that only the linear regression is done at MB level
         (e.g., no outlier rejection) in order to save complexity.         */

      if ( bits_MB > 0 ) 
        {
          if (N_ls < Q2MB_N) 
            {
              Q2_mb.X_ls[N_ls][0]= (double) 1/QP;
              Q2_mb.X_ls[N_ls][1]= (double) 1;
              Q2_mb.Y_ls[N_ls]=    (double) (bits_MB*QP/mad);
              
              { int i;  /* VM5 */
              float sum_Y=0;
              for (i=0; i<=N_ls; i++) {
                sum_Y += Q2_mb.Y_ls[N_ls];
              }
              Q2_mb.X11 = sum_Y/(N_ls+1);
              /* fprintf(stdout,"Q2_mb.Y_ls[%2d]=%f\n",N_ls,Q2_mb.Y_ls[N_ls]); */
              }
              
              N_ls++;
            }
          else
            { /* Q2_mb.X_ls and Q2_mb.Y_ls are full with the Q2MB_N data points, 
                 so we just shift them to the left and add new data at the end */

              for (k=1; k < Q2MB_N; k++)
                {  
                  Q2_mb.X_ls[k-1][0]= Q2_mb.X_ls[k][0];
                  Q2_mb.X_ls[k-1][1]= Q2_mb.X_ls[k][1];
                  Q2_mb.Y_ls[k-1]   = Q2_mb.Y_ls[k];
                }
              /* N_ls is Q2MB_N now: the last index in array is N_ls-1 */
              Q2_mb.X_ls[N_ls-1][0]= (double)  1/QP;
              Q2_mb.X_ls[N_ls-1][1]= (double)  1;
              Q2_mb.Y_ls[N_ls-1]=    (double)  (bits_MB*QP/mad);
              
              { int i;  /* VM 5 */
              float sum_Y=0;
              for (i=0; i<Q2MB_N; i++)
                sum_Y += Q2_mb.Y_ls[N_ls];
              Q2_mb.X11 = sum_Y/(Q2MB_N);
              }
              /* fprintf(stdout,"Q2_mb.Y_ls[%2d]=%f\n",N_ls-1,Q2_mb.Y_ls[N_ls-1]); */
            }
          
          /* ----------- LS fit or Linear regression */
          corrXY[0]= 0;
          corrXY[1]= 0;
          for (k=0; k < N_ls; k++)
            {
              corrXY[0] += Q2_mb.Y_ls[k]*Q2_mb.X_ls[k][0];
              corrXY[1] += Q2_mb.Y_ls[k]*Q2_mb.X_ls[k][1];
            }
          
          /* LS part 2: Q2_mb.X_ls transposed by itself, inverted */
          /* Inverse of Autocorrelation matrix of Q2_mb.X_ls */
          for (k=0; k < N_ls; k++)
            {
              a += Q2_mb.X_ls[k][0]*Q2_mb.X_ls[k][0];
              b += Q2_mb.X_ls[k][1]*Q2_mb.X_ls[k][0];
              d += Q2_mb.X_ls[k][1]*Q2_mb.X_ls[k][1];
            }
          c= b;
          determ= a*d - c*b;
          if (fabs(determ) == 0)
            {
              /* fprintf(stdout,"\n ERROR: singular matrix. Keeping old values of A2 and A3");
               */
            }
          else
            {
              invX[0][0]=  d/determ;
              invX[1][0]= -c/determ;
              invX[0][1]= -b/determ;
              invX[1][1]=  a/determ;
              *A2= (float) (invX[0][0]*corrXY[0] + invX[0][1]*corrXY[1]);
              *A3= (float) (invX[1][0]*corrXY[0] + invX[1][1]*corrXY[1]);
            }
        }
    }                   
  
  /* Update bits left and energy left */
  Q2_mb.energy_S -= Q2_mb.mad[index_y][index_x]*Q2_mb.weight[index_y][index_x];
  Q2_mb.T_text   -= bits_MB;
  
  /* If last macroblock, must deallocate memory and change value of A1_prev */
  if (MB_coded == MB_total)
    {
      for (j=0; j <= index_y; j++) 
        {   
          free(Q2_mb.mad[j]);
          free(Q2_mb.weight[j]);
        }
      
      free(Q2_mb.mad);
      free(Q2_mb.weight);
      
      /* A1_prev for next frame is the value of A1 at the end of VOP */
      Q2_mb.A1_prev= Q2_mb.A1;  
      /* bits_sms is only used here 
         -- 55 are the typical vop syntax bits for P VOP*/
      Q2_mb.T_sms = bits_sms + 55; 

#ifdef _RC_DEBUG_
      fprintf(stdout,"\n");     /* adds one line after the QP array */
#endif                            

    }       
}

/* SHARP_end                                  */
/* ------------------------------------------ */



/***********************************************************HeaderBegin*******
 * LEVEL 2 (NO ACCESS TO STATIC VARS)
 ***********************************************************HeaderEnd*********/

/***********************************************************CommentBegin******
 *
 * -- RCQ2_compute_quant --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Q2 quantizer computation
 *      Note: floor is used to produce the same results of original 
 *      implementation.
 *
 * Arguments in :       
 *      Double   bit_rate     -  Bit rate (bits per second)
 *      Int      bits_frame   -  Target bits per frame 
 *      Int      bits_prev    -  Total bits previous frame 
 *      Int      buff_size    -  Buffer size
 *      Int      buff_occup   -  Buffer occupancy
 *      Double   X1           -  Model parameter
 *      Double   X2           -  Model parameter
 *      Double   mad          -  Pred. error. magnitude
 *      Int      last_qp      -  Previous QP
 *
 * Modified : 
 *
 *
 ***********************************************************CommentEnd********/

Int RCQ2_compute_quant(
                       Double bit_rate,     /* Bit rate (bits per second) */
                       Int    bits_frame,   /* Target bits per frame */
                       Int    bits_prev,    /* Total bits previous frame */
                       Int    buff_size,    /* Buffer size */
                       Int    buff_occup,   /* Buffer occupancy */
                       Double X1,           /* Model parameters */
                       Double X2,
                       Double mad,          /* Pred. error. magnitude */
                       Int    last_qp       /* Previous QP */
                       )
{
Double target=(Double)bits_frame;
Double qp;
Int a, b;

target = floor(target);


#ifdef _RC_DEBUG_
   fprintf(stdout, "--> RCQ2_compute_quant\n");
   fprintf(stdout, "RC: RC_Q2: --> RCQ2_compute_quant\n");
   fprintf(stdout, "RC:        bits_frame= %d\n", bits_frame);
   fprintf(stdout, "RC:        bits_prev= %d\n", bits_prev);
   fprintf(stdout, "RC:        buff_size= %d\n", buff_size);
   fprintf(stdout, "RC:        buff_occup= %d\n", buff_occup);
   fprintf(stdout, "RC:        X1= %f\n", X1);
   fprintf(stdout, "RC:        X2= %f\n", X2);
   fprintf(stdout, "RC:        mad= %f\n", mad);
   fprintf(stdout, "RC:        last_qp= %d\n", last_qp);
   fprintf(stdout, "<-- RCQ2_compute_quant\n");
   fflush(stdout);
#endif

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant:        target0= %f\n", target);
   fflush(stdout);
#endif

/*** Low-pass filtering of the target ***/
if (bits_prev > 0)
   target = target*(1.-PAST_PERCENT)+bits_prev*PAST_PERCENT;
target = floor(target);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant:        target1= %f\n", target);
   fflush(stdout);
#endif

/*** Target tuning depending on buffer occupancy ***/
  a = buff_occup;
  b = buff_size-buff_occup;
  target = (a+2.0*b)*target/(2.0*a+b);
  target = floor(target);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant:        target2= %f\n", target);
   fflush(stdout);
#endif

/*** Lower-bounding of target ***/
target = MAX(bit_rate/30, target);
target = floor(target);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant:        target3= %f\n", target);
   fflush(stdout);
#endif

/*** Target clipping depending on buffer occupancy ***/
if (buff_occup+target > (1.-SAFETY_MARGIN)*buff_size)
   target = MAX(bit_rate/30, (1.-SAFETY_MARGIN)*buff_size-buff_occup);
else if (buff_occup-bits_frame+target < SAFETY_MARGIN*buff_size)
   target = SAFETY_MARGIN*buff_size-buff_occup+bits_frame;
target = floor(target);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant:        target4= %f\n", target);
   fflush(stdout);
#endif


if (target<=0)
   return MAX_QUANT;

if (X2==0)
   qp = X1*mad/target;
else
   {
   Double disc=X1*X1+4*X2*target/mad;
   if (disc < 0)
      {
#ifdef _RC_DEBUG_
      fprintf(stdout, "RCQ2_compute_quant: Second order model without solution: X2 set to 0\n");
      fflush(stdout);
#endif
      qp=X1*mad/target;
      }
   else
      {
      qp = (2*X2*mad)/(sqrt((X1*mad)*(X1*mad)+4*X2*mad*target)-X1*mad);
      }
   }
#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant: qp0= %f\n", qp);
#endif
qp = (Int) floor(qp+0.5);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant: qp1= %f\n", qp);
   fflush(stdout);
#endif

/*** Bounding of QP variation ***/
qp = MIN3(ceil(last_qp*1.25), qp, MAX_QUANT);
qp = MAX3(ceil(last_qp*0.75), qp, MIN_QUANT);

#ifdef _RC_DEBUG_
   fprintf(stdout, "RCQ2_compute_quant:        qp2= %f\n", qp);
   fflush(stdout);
#endif

#ifdef _RC_DEBUG_
   fprintf(stdout, "<-- RCQ2_compute_quant\n");
   fflush(stdout);
#endif


return (Int)(qp+.5);
}

/*********************************************************** End of file ***/
