/*
 * Copyright 1993 by the University of Pennsylvania
 *
 * Permission to use, copy, and distribute for non-commercial purposes,
 * is hereby granted without fee, providing that the above copyright
 * notice appear in all copies and that both the copyright notice and this
 * permission notice appear in supporting documentation.
 *
 * The software may be modified for your own purposes, but modified versions
 * may not be distributed.
 *
 * This software is provided "as is" without any expressed or implied warranty.
 */

#include <stdio.h>
#include <string.h>
#include <memory.h>
#include <math.h>
#include <fcntl.h>
#include "global.h"
#include "mpeg.h"

double get_snr();

/*---------------------- globals ------------------------------------*/
extern int mb_type;       /* Macro block type 1: intra, 2:pred or Bi-Dir */
/****************************************************************************/
/*                     Subroutine to code a frame using                     */
/*                     predictive picture mode                              */
/****************************************************************************/
FRAME_STAT*
  do_pred(frame,recon_frame,quant_scale, dctnum, height, width)
FRAME *frame;
FRAME **recon_frame;
int quant_scale, dctnum, height, width; 
{
  MB  mblock;
  short *rec_mb,**dec_mb,*mc_mblock;
  int i,j,k;
  int n_hor_mb,n_vert_mb;
  int voff,hoff,offset1d;
  int *mvx,*mvy;
  int intra_prev = 0;
  FRAME *prev_frame;
  int prev_mvx = 0;
  int prev_mvy = 0;
  FRAME_STAT  *FS ;
/*************************************************************/
/*-------------- Initialize & Memory allocations ------------*/
/*************************************************************/
/*----- Compute number of macroblocks we need to code  ----*/
  n_vert_mb = height/MBSIZE;
  n_hor_mb  = width/MBSIZE;
/*---------------------- stat initialize --------------------*/
  FS = (FRAME_STAT *)malloc(sizeof(FRAME_STAT));
  FS->block_mode = (char *)malloc(sizeof(char)*n_vert_mb*n_hor_mb);
  FS->bit_num  = (short *)malloc(sizeof(short)*n_vert_mb*n_hor_mb);
  FS->bit_numl = NULL;
  if(dctnum < 64) 
    FS->bit_numl = (short *)malloc(sizeof(short)*n_vert_mb*n_hor_mb);
  FS->frame_mode = 1;
  FS->width  = width;
  FS->height = height;
/*-----------------------------------------------------------*/
  prev_frame = (FRAME *)make_frame(width,height);
  copy_frame(prev_frame,recon_frame[0]);

/*------------------- Motion vectors -------------------------*/
  mvx = (int *)calloc(n_vert_mb*n_hor_mb,sizeof(int));
  mvy = (int *)calloc(n_vert_mb*n_hor_mb,sizeof(int));
    
/*-------------- macroblock data -----------------------------*/
  mblock.data = make_mb();

/*----- low and high priority reconstructed macroblocks ------*/
  dec_mb = make_mb();

/*********************************************************************/
/*             Begin predictive picture mode coding                  */
/*********************************************************************/
#ifdef DEBUG
  fprintf(stderr,"PREDICTION MODE\n");
#endif
/*------------ Insert picture (frame) header ------------------------*/  
  put_picture_code(1);        /* 66 bits */
/*----------------- Get MV for MB's of current frame ----------------*/
  if (Xon)
    if(working_update(1,"Prediction mode : Making Motion Vector",
		       height/MBSIZE))  goto interrupt_jump_label;
  if (construct_mv(frame->Y,prev_frame->Y,mvx,mvy,width,height))
    return NULL;
         
#ifdef DEBUG
  fprintf(stderr,"ENCODING ");
#endif
  if (Xon)
    if(working_update(1,"Prediction mode : Encoding",height/MBSIZE))
      goto interrupt_jump_label;
/*--------------------------- Loop start ----------------------------*/
  for (i=0; i< n_vert_mb; i++) { /* slice increment */
    intra_prev = 0;
    voff = i*16;
    prev_mvx = 0; prev_mvy = 0;
/*---------------------- Insert slice header -------------------------*/
    put_slice_code(i+1,quant_scale);
/*--------------------------------------------------------------------*/
#ifdef DEBUG
    fprintf(stderr,".");
#endif
    if (Xon)
      if(working_update(1, NULL, 0)) goto interrupt_jump_label;
/*--------------------------------------------------------------------*/
    for (j=0; j < n_hor_mb; j++) {    /* block increment */
      mb_type = 2;                    /* set no-I MB */
      offset1d = i*n_hor_mb+j;
      hoff = j*16;  
      free(mblock.data[0]);
      mblock.data[0] = get_mb(frame,voff,hoff);      /* get MB */
      rec_mb         = get_mb(prev_frame,voff,hoff); /* get prev MB w/o MC */
      mc_mblock      = get_mb(prev_frame,            /* get prev MB w/  MC */
                       voff+mvy[offset1d],hoff+mvx[offset1d]);
/*--------------- decide use MC or not (MC = (0,0)) ----------------*/
      mblock.mc = mc_mb(mblock.data[0], rec_mb, mc_mblock);
      if((mvx[offset1d] || mvy[offset1d]) && mblock.mc) { /* MC */
	free(rec_mb);
	rec_mb = mc_mblock;
	mblock.mv[0] = mvy[offset1d] - prev_mvy;
	mblock.mv[1] = mvx[offset1d] - prev_mvx;
	mblock.mv[0] += big(mblock.mv[0]);
	mblock.mv[1] += big(mblock.mv[1]);
	prev_mvx = mvx[offset1d];
	prev_mvy = mvy[offset1d];
      }
      else { /* No MC */
	free(mc_mblock);
	prev_mvx = prev_mvy = 0;
	mblock.mc = 0;
      }
/*--------------- decide MB intra coded or not --------------*/
      mblock.intra = int_noint(mblock.data[0],rec_mb);
      if(mblock.intra) {  /* MB intra coding */
        mb_type = 1;
	code_i_mb(mblock.data, dec_mb, intra_prev, dctnum, quant_scale);
	mblock.mc = 0;
	mblock.type = 0;  /* Intra */
	prev_mvx = prev_mvy = 0;
	intra_prev = 1;
        FS->block_mode[offset1d] = 0;  /* I mode */
      }
      else {          /* MB Fwd Prediction coding */
	subtract_mb(mblock.data[0],rec_mb);
	code_i_mb(mblock.data, dec_mb, 0, dctnum, quant_scale);
	mblock.type = 1; /* Fwd-P */

	add_mb(dec_mb[1],rec_mb);
        add_mb(dec_mb[0],rec_mb);

        FS->block_mode[offset1d] = 1 + ((1-mblock.mc) * 10);
                /* P mode  1: w/o MC  11: w/ MC */
	intra_prev = 0;
      }

      range_mb(dec_mb[0],0,255);
      range_mb(dec_mb[1],0,255);

      for (k=0;k<NUM_P+1;k++) 
	put_mb(recon_frame[k],dec_mb[k],voff,hoff);
/*------------------------- Perform VLC ---------------------------*/
      FS->bit_num[offset1d]  =   /* perform VLC on HP code */
          (short)do_vlc(mblock,0, dctnum,((j==0)?1:0), 1, n_hor_mb);
      if(dctnum < 64)
      FS->bit_numl[offset1d] =   /* perform VLC on LP code */
          (short)do_vlc(mblock,1, dctnum,((j==0)?1:0), 1, n_hor_mb);
       
      free(rec_mb); 
    } /* end block */
  }   /* end slice */
/* compute SNR for original and reconstructed frame */
   FS->SNR = get_snr(frame,recon_frame,width,height);
/*--------------------*/
/*     FREE's         */
/*--------------------*/
  if(0) {
    interrupt_jump_label:
    free(FS->block_mode);
    free(FS->bit_num);
    free(FS->block_mode);
    if(FS->bit_numl != NULL) free(FS->bit_numl);
    free(FS);
    FS = NULL;
  }
  free(mvx);
  free(mvy);
  free_mb(dec_mb);
  free_mb(mblock.data);
  free_frame(prev_frame);
  
#ifdef DEBUG
  fprintf(stderr," Complete!\n");
#endif
  return(FS);
}
/***********************************************************************/
/* Decide whether to use MC or not (0: MC, 1: No MC)                    */
/***********************************************************************/
int
  mc_mb(or_mb,rec_mb,disp_mb)
short *or_mb;
short *rec_mb;
short *disp_mb;
{
  float x=0,y=0;
  int i;

  for(i=0;i<MBSIZE_2;i++) {
    y += abs(or_mb[i] - disp_mb[i]);
    x += abs(or_mb[i] - rec_mb[i]);
  }
  x = x/MBSIZE_2 ; y = y/MBSIZE_2;
  if(x > 1) {
    if(((x < 3) && (x > 2*y)) || ((x >= 3) && (x > 1.1*y))){
      return(1); /* MC (1)*/
    }
  }
  return(0); /* No MC (0)*/
}
/**********************************************************************/
/* Decide whether to perform I or P Coding ( 0: Intra, 1: Fwd Pred )  */
/**********************************************************************/
int
  int_noint(or_mb,rec_mb)
short *or_mb;
short *rec_mb;
{
  int i;
  float var=0,varor=0,meanor=0;
  for(i=0;i<MBSIZE_2;i++) {
    var += ((rec_mb[i] - or_mb[i])*(rec_mb[i] - or_mb[i]));
    varor += or_mb[i]*or_mb[i];
    meanor += or_mb[i];
  }
  var = var/256;
  varor = varor/256 - meanor/256*meanor/256;

  if((varor < var) && (var > 64))
    return(1); /* Intra (1)*/
  else
    return(0); /* Non-Intra (0) */
}
/****************************************************************************/
/*                     Subroutine to code a frame using                     */
/*                     intraframe picture mode                              */
/****************************************************************************/
FRAME_STAT*
  do_intra(frame,recon_frame, quant_scale, dctnum, height, width)
FRAME *frame;
FRAME **recon_frame;
int quant_scale, dctnum, height, width;
{
  int i,j,k;
  int n_hor_mb,n_vert_mb,voff,hoff,offset1d;
  MB mblock;
  short **recon_mb = NULL;
  FRAME_STAT  *FS;
/*************************************************************/
/*-------------- Initialize & Memory allocations ------------*/
/*************************************************************/
/*----- Compute number of macroblocks we need to code  ----*/
  n_vert_mb = height/MBSIZE;
  n_hor_mb  = width/MBSIZE;
/*---------------------- stat initialize --------------------*/
  FS = (FRAME_STAT *)malloc(sizeof(FRAME_STAT));
  FS->block_mode = (char *)malloc(sizeof(char)*n_vert_mb*n_hor_mb);
  FS->bit_num  = (short *)malloc(sizeof(short)*n_vert_mb*n_hor_mb);
  FS->bit_numl = NULL;
  if(dctnum < 64)
    FS->bit_numl = (short *)malloc(sizeof(short)*n_vert_mb*n_hor_mb);
  FS->frame_mode = 0;
  FS->width = width;
  FS->height = height;

/*-------------- mblock data ----------------------------------*/
  mblock.data = make_mb();
  recon_mb = make_mb();

/*********************************************************************/
/*             Begin predictive picture mode coding                  */
/*********************************************************************/
#ifdef DEBUG
  fprintf(stderr,"INTRAFRAME MODE\nENCODING ");
#endif
  if (Xon)
    if(working_update(1,"Intra mode : Encoding",height/MBSIZE))
      goto interrupt_jump_label;
/*------------------------------------------------------------*/
  put_group_code(); /* I frame is always beginning of group */
/*---------------- Set macroblock type to INTRA -------------*/
  mb_type = 1;
  mblock.intra = 1; /* intra mode */
  mblock.type  = 1; /* intra mode */
  mblock.mc    = 0; /* no MC (Obviously) */
/*------------------------------------------------------------*/
  put_picture_code(0);      /* 62 bits */
/*--------------------------- Loop start ----------------------------*/
  for (i=0; i<n_vert_mb; i++) { /* slice increment */
#ifdef DEBUG
    fprintf(stderr,".");
#endif
    if (Xon)    
      if(working_update(1,NULL,0)) goto interrupt_jump_label;
    voff = i*16;
/*---------------------- Insert slice header -------------------------*/
    put_slice_code(i+1,quant_scale);
    for (j=0; j< n_hor_mb; j++) { /* block increment */
      hoff = j*16;
      offset1d = i*n_hor_mb+j;
      free(mblock.data[0]); 
      mblock.data[0] = get_mb(frame,voff,hoff);
      code_i_mb(mblock.data, recon_mb, j, dctnum,quant_scale); 

      range_mb(recon_mb[0],0,255);
      range_mb(recon_mb[1],0,255);

      for (k=0;k<NUM_P+1;k++) 
	put_mb(recon_frame[k], recon_mb[k], voff, hoff);
      FS->block_mode[offset1d] = 0; /* I mode MB */
/*------------------------- Perform VLC ---------------------------*/
      FS->bit_num[offset1d] =   /* perform VLC on HP code */
              (short)do_vlc(mblock,0, dctnum,((j==0)?1:0), 0, n_hor_mb);
      if(dctnum < 64)
        FS->bit_numl[offset1d] =   /* perform VLC on LP code */
              (short)do_vlc(mblock,1, dctnum,((j==0)?1:0), 0, n_hor_mb);
      
    }
  }
/* compute SNR for original and reconstructed frame */
  FS->SNR = get_snr(frame,recon_frame,width,height);
/*--------------------*/
/*     FREE's         */
/*--------------------*/
  if(0) {
  interrupt_jump_label:
    free(FS->block_mode);
    free(FS->bit_num);
    free(FS->block_mode);
    if(FS->bit_numl != NULL) free(FS->bit_numl);
    free(FS);
    FS = NULL;
  }

#ifdef DEBUG
  fprintf(stderr," Complete!\n");
#endif
  free_mb(recon_mb);
  free_mb(mblock.data);
  return(FS);
}
/****************************************************************************/
/* Take macroblock as input, performs DCT, variable length coding, 
   outputs coded data and returns reconstructed macroblock */
void
  code_i_mb(mblock, recon_mb, first, dctnum, quant_scale)
short **mblock;
short **recon_mb;
int first;
int dctnum;
int quant_scale;
{
  int j,k;
  int btype;  /* represents MB type Y:0, Cb:1, Cr:2 */
  int fflag; /* takes 1 if MB is first slice and first block of Y or Cr or Cb 
                otherwise, takes 1 */

  /* type used to determine if macroblock is Y,Cb or Cr */
  for( btype = 0, k = 0; k < 6; k++ ) {
    fflag = ((first == 0) && ((k == 0) || (k > 3)));
    if ( k > 3 ) btype++; /*increment btype for Cb and Cr blocks*/
    code_intra( mblock[0] + k*DCT_SIZE_2, fflag, btype, quant_scale); 
  }

  do_split(mblock,recon_mb,dctnum,first, quant_scale);

  for(j=0;j<NUM_P+1;j++)
    for(k=0;k<6;k++) 
      make_zz(mblock[j]+k*DCT_SIZE_2);
}
/****************************************************************************/
void 
  code_intra(block, fflag, type, quant_scale)
short *block;
int fflag; /* 1: first MB of each slice */
int type;  /* 0: Y, 1: Cb, 2: Cr */
int quant_scale;
{
  j_fwd_dct(block);
  if (mb_type != 1)
    var_threshold(block,quant_scale);
  if(mb_type == 1)
    quantize_dc(block,fflag,type); 
  quantize_ac(block,quant_scale);
}
/****************************************************************************/
void 
  quantize_dc(block, fflag, type)
short *block;
int fflag; /* 1: first MB of each slice */
int type;  /* 0: Y, 1: Cb, 2: Cr */
{
  static int dct_dc_past[3];
  int level;

  /* round to nearest integer */
  level = (int)((double)block[0]/(double)intra_quant[0] + 0.5); 
  block[0] = RANGE((level - ( fflag ? 128 : dct_dc_past[type] )),-255,255);
  dct_dc_past[type] = level;
}
/****************************************************************************/
void 
  quantize_ac(block,quant_scale)
short *block;
int quant_scale;
{
  int i,j, level;
  double quant;

  j = (mb_type == 1 ? 1 : 0);



  for ( i = j ; i < DCT_SIZE_2 ; i++ ) {
    if (block[i]) {
      quant = (double)( j ? intra_quant[i] : 16.0);
      level = (int)((double)(16 * block[i])/(double)(2 * quant * quant_scale) 
		    + (1 - j)*sgn(block[i])/(double)(2*quant_scale));
      block[i] = RANGE(level,-255,255); 
    }
  }
}
/****************************************************************************/
void
  make_zz(block)
short *block;
{
  int i,j;
  short tmp[DCT_SIZE_2];
  
   for(i=0;i<DCT_SIZE;i++)
     for(j=0;j<DCT_SIZE;j++)
       tmp[scan[i][j]] = block[i*DCT_SIZE+j];

  for(i=0;i<DCT_SIZE_2;i++)
    block[i] = tmp[i] ;
}
/****************************************************************************/
void
  make_unzz(block)
short *block;
{
  int i,j;
  short tmp[DCT_SIZE_2];
  
   for(i=0;i<DCT_SIZE;i++)
     for(j=0;j<DCT_SIZE;j++)
       tmp[i*DCT_SIZE+j] = block[scan[i][j]];

  for(i=0;i<DCT_SIZE_2;i++)
    block[i] = tmp[i] ;
}
/****************************************************************************/
/* decode a quantized macroblock */
void
  decode_i(mblock, recon_mb, first, testrun,quant_scale)
short *mblock;
short *recon_mb;
int first;
int testrun;
int quant_scale;
{
  int i,l,m,n,index;
  float quant=16;
  static short y;
  static short cb;
  static short cr;
  short dct_dc_y_past,dct_dc_cb_past,dct_dc_cr_past;

  dct_dc_y_past = y;
  dct_dc_cb_past = cb;
  dct_dc_cr_past = cr;

  i = (mb_type == 1 ? 1 : 0);

  for(l=0; l<6; l++){
    for(m=0; m<DCT_SIZE; m++)
      for(n=0; n<DCT_SIZE; n++){
	index = n + m*DCT_SIZE + l*DCT_SIZE_2;
	recon_mb[index] = (2 * mblock[index] + (1-i)*sgn(mblock[index]))*
	  quant_scale*(i ? intra_quant[n + m*DCT_SIZE] : quant)/16.0;
	if ((recon_mb[index] & 1) == 0)
	  recon_mb[index] = recon_mb[index] - sgn(recon_mb[index]);
	if(recon_mb[index] > 2047) recon_mb[index] = 2047;
	if(recon_mb[index] < -2048) recon_mb[index] = -2048;
      }
    if (mb_type == 1) {
      switch(l) {
      case 0: 
	recon_mb[0] = (first == 0 ? 128*8 : dct_dc_y_past) + mblock[0]*8 ;
	dct_dc_y_past = recon_mb[0];
	break;
      case 1: case 2: case 3: 
	recon_mb[l*DCT_SIZE_2] = dct_dc_y_past + mblock[l*DCT_SIZE_2]*8;
	dct_dc_y_past = recon_mb[l*DCT_SIZE_2];
	break;
      case 4:
	recon_mb[l*DCT_SIZE_2] = (first == 0 ? 128*8 : dct_dc_cb_past)
	  + mblock[l*DCT_SIZE_2]*8;
	dct_dc_cb_past = recon_mb[l*DCT_SIZE_2];
	break;
      case 5:
	recon_mb[l*DCT_SIZE_2] = (first == 0 ? 128*8 : dct_dc_cr_past)
	  + mblock[l*DCT_SIZE_2]*8;
	dct_dc_cr_past = recon_mb[l*DCT_SIZE_2];
	break;
      }
    }
    j_rev_dct(recon_mb + l*DCT_SIZE_2);
  }
  if (testrun != 1) {
    y = dct_dc_y_past;
    cb = dct_dc_cb_past;
    cr = dct_dc_cr_past;
  }
}
/****************************************************************************/
/* reads a MB from an image as specified by voff,hoff                       */
/****************************************************************************/
short *
  get_mb(frame, voff, hoff)
FRAME *frame;
int voff;
int hoff;

{
  int i,j,k,index;
  short *mblock;
  int width = frame->width;
  int cindex,yindex;

  mblock = (short *)calloc(6*DCT_SIZE_2,sizeof(short));
   
  /* Macroblock structure consists of 6 blocks:
     Y : Blocks 0-3; Cb : Block 4 ; Cr : Block 5 */

  for(k=0;k<3;k++)
    for(i=0;i<DCT_SIZE;i++)
      for(j=0;j<DCT_SIZE;j++) {
	index = j + i*DCT_SIZE +  k*DCT_SIZE_2;
	switch(k) {
	case 0 : case 1 : /* Luminance components */
	  yindex=(voff+i)*width + hoff + j + k*DCT_SIZE;
	  mblock[index] = frame->Y[yindex];
          mblock[index + 2*DCT_SIZE_2] = frame->Y[yindex + DCT_SIZE*width];
	  break;
	case 2 : /* Chrominance Components */
	  cindex = (voff/2 + i)*width/2 + hoff/2 + j;
	  if (hoff % 2 && voff % 2) {
	    mblock[index + 2*DCT_SIZE_2] = 
	      ((int)frame->Cb[cindex] + (int)frame->Cb[cindex + 1] + 
	       (int)frame->Cb[cindex + width/2] + 
	       (int)frame->Cb[cindex + 1 + width/2] + 0.5)/4.0;
	    mblock[index + 3*DCT_SIZE_2] = 
	      ((int)frame->Cr[cindex] + (int)frame->Cr[cindex + 1] + 
	       (int)frame->Cr[cindex + width/2] + 
	       (int)frame->Cr[cindex + 1 + width/2] +
	       0.5)/4.0;
	  }
	  else if (hoff % 2) {
	    mblock[index + 2*DCT_SIZE_2] = 
	      ((int)frame->Cb[cindex] + (int)frame->Cb[cindex + 1] + 0.5)/2.0; 
	    mblock[index + 3*DCT_SIZE_2] = 
	      ((int)frame->Cr[cindex] + (int)frame->Cr[cindex + 1] + 0.5)/2.0;
	  }
	  else if (voff % 2) {
	    mblock[index + 2*DCT_SIZE_2] = 
	      ((int)frame->Cb[cindex] + (int)frame->Cb[cindex + width/2] 
	       + 0.5)/2.0;
	    mblock[index + 3*DCT_SIZE_2] = 
	      ((int)frame->Cr[cindex] + (int)frame->Cr[cindex + width/2] 
	       + 0.5)/2.0;
	  }
	  else { 
	    mblock[index + 2*DCT_SIZE_2] = (int)frame->Cb[cindex];
	    mblock[index + 3*DCT_SIZE_2] = (int)frame->Cr[cindex];
	  }
	  break;
	}
      }
  return mblock;
}	
/****************************************************************************/
/* places a macroblock in an image array at position specified by hoff,voff */
/****************************************************************************/
void 
  put_mb(frame, mblock, voff,  hoff)
FRAME *frame;
short *mblock;
int voff;
int hoff;

{
  int k,l,m,index;
  int width = frame->width;
  
  for(m=0;m<3;m++)
    for(k=0;k<DCT_SIZE;k++)
      for(l=0;l<DCT_SIZE;l++){
	index = m*DCT_SIZE_2 + k*DCT_SIZE + l;
	switch (m) {
	case 0: case 1:
	  frame->Y[(k + voff)*width + l + hoff + m*DCT_SIZE] = mblock[index];
	  frame->Y[(k + voff + DCT_SIZE)*width + l + hoff + m*DCT_SIZE] = 
	    mblock[index + 2*DCT_SIZE_2];
	  break;
	case 2:
	  frame->Cb[(voff/2 + k)*width/2 + hoff/2 + l] = 
	    mblock[index + 2*DCT_SIZE_2];
	  frame->Cr[(voff/2 + k)*width/2 + hoff/2 + l] = 
	    mblock[index + 3*DCT_SIZE_2];
	  break;
	}
      }
}

do_split(orig_mb,recon_mb,dctnum,fflag, quant_scale)
short **orig_mb,**recon_mb;
int dctnum,fflag;
int quant_scale;
{
  int j,k,l,m,pos;

  for(j=1;j<NUM_P+1;j++) {
    switch (j) { /* Priority specific operations */
    case 1:
      for(k=0;k<6;k++)
	for(l=0;l < DCT_SIZE;l++)
	  for(m = 0; m < DCT_SIZE; m++) {
	    pos = m + l*DCT_SIZE + k*DCT_SIZE_2;
	    if (scan[l][m] < dctnum) 
	      orig_mb[j][pos] = orig_mb[0][pos];
	    else 
	      orig_mb[j][pos] = 0;
	  }
      decode_i(orig_mb[j], recon_mb[j], fflag, 1, quant_scale);
      break;
    case 2:
      for(k=0;k<6;k++)
	for(l=0;l < DCT_SIZE;l++)
	  for(m = 0; m < DCT_SIZE; m++) {
	    pos = m + l*DCT_SIZE + k*DCT_SIZE_2;
	    if (scan[l][m] < dctnum) 
	      orig_mb[j][pos] = 0;
	    else
	      orig_mb[j][pos] = orig_mb[0][pos];
	  }
      decode_i(orig_mb[j], recon_mb[j], 1, 1, quant_scale);
      break;
    }
  }
  decode_i(orig_mb[0],recon_mb[0],fflag,0, quant_scale); 
}


var_threshold(block,quant_scale)
short *block;
int quant_scale;
{
  int i,j;
  int T = quant_scale*2;
  int Tmax = quant_scale*3;
  
  for (i=0;i<8;i++)
    for (j=0;j<8;j++)
      if (abs(block[scan[i][j]]) < T) {
	block[scan[i][j]] = 0;
	T = ((T+1 < Tmax) ? T+1 : Tmax);
      }
      else 
	T = quant_scale*2;
}




/****************************************************************************/
/* Compute SNR for original and constructed                                 */
/****************************************************************************/
double
  get_snr(frame1,frame2,width,height)
byte **frame1,**frame2;
int width, height;
{
  int i,j;
  double noise=0.0;
  double signal=0.0;

  return 0;

  for(i=1;i<=height;i++)
    for(j=1;j<=2*width;j++) {
      noise += ((frame1[i][j]-frame2[i][j])*(frame1[i][j]-frame2[i][j]));
      signal += frame1[i][j]*frame1[i][j];
    }

  return (10*log10(signal/noise));

}


