/*****************************************************************************
 *
 * This software module was originally developed by
 *
 *   J. Ignacio Ronda (UPM-GTI / ACTS-MoMuSys)
 *
 * and edited by
 *
 *   Angel Pacheco (UPM-GTI / ACTS-MoMuSys)
 *   Robert Danielsen (Telenor / ACTS-MoMuSys)
 *   Fernando Jaureguizar (UPM-GTI / ACTS-MoMuSys)
 *   Ulrike Pestel (TUH / ACTS-MoMuSys)
 *   Minhua Zhou (HHI / ACTS-MoMuSys)
 *   Ferran Marques (UPC / ACTS-MoMuSys)
 *   Bob Eifrig (NextLevel Systems)
 *
 * 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:        mot_est.c
 *
 * Author:      J. Ignacio Ronda / Angel Pacheco, UPM-GTI
 *
 * Created:     02.02.96
 *
 * Description: Motion Estimation  of a VOP
 *
 * Flags:       -D_NO_ESTIMATION_ : dummy zero motion vectors and INTRA
 *                                  modes are produced directly
 *
 * Modified:
 *      21.04.96 Robert Danielsen: Reformatted. New headers.
 *      02.07.96 Fernando Jaureguizar: Reformatted.
 *      22.10.96 Angel Pacheco:
 *               - new define DEFAULT_EDGE
 *               - modifications to cope with the new MB-based padding scheme
 *      28.11.96 Ferran Marques: Integration of changes from (22.10.96)
 *      05.02.97 Angel Pacheco: deletion of the valid_pel_* arrays (VM5.1)
 *               having the valid positions (padded pixels) in which the
 *               search for the motion vectors must be done (as VM4.x said).
 *      14.05.97 Fernando Jaureguizar: new commentaries about new enlarged
 *               f_code range according to VM7.0
 *               Modification of B_MotionEstimation_after_padding() call to
 *               cope with enlarged f_code range.
 *      16.06.97 Angel Pacheco: unified the TRANSPARENT modes.
 *      17.07.97 Ulrike Pestel: for scalability
 *      01.08.97 Fernando Jaureguizar: changed the wrong name advanced by the
 *               correct one enable_8x8_mv.
 *      13.11.97 Fernando Jaureguizar: formating and new header.
 *      11.12.97 Bob Eifrig:  added code for interlaced video
 *               (field motion estimation)
 *      09.03.98 Fernando Jaureguizar: New formating. Deleted not used
 *               variables and parameters.
 *
 ***********************************************************HeaderEnd*********/

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

#include "momusys.h"
#include "mom_structs.h"
#include "vm_common_defs.h"
#include "mom_access.h"
#include "mom_image.h"
#include "mom_vop.h"
#include "mot_util.h"
#include "mot_padding.h"
#include "mot_est.h"
#include "io_generic.h"
#include "mot_comp.h"

#define DEFAULT_EDGE 16

/***********************************************************CommentBegin******
 *
 * -- MotionEstimation -- Estimates the motion vectors
 *
 * Author :
 *      UPM-GTI - J. Ignacio Ronda / Angel Pacheco
 *
 * Created :
 *      02.02.96
 *
 * Purpose :
 *      Estimates the motion vectors in advanced/not_advanced unrestricted
 *      mode. The output are four images containing x/y components of
 *      MV's (one per 8x8 block), modes (one value per MB), and an the
 *      8x8-subsampled alpha plane
 *
 * Arguments in :
 *      Vop     *curr_vop,      current Vop (for luminance)
 *      Vop     *curr_rec_vop,  reconstructed current  Vop (for shape)
 *      Vop     *prev_rec_vop,  reference Vop (reconstructed)(1/2 pixel)
 *      Vop     *prev_vop,      reference Vop (original) (pixel)
 *      Int     enable_8x8_mv,  8x8 MV (=1) or only 16x16 MV (=0)
 *      Int     edge,           restricted(==0)/unrestricted(==edge) mode
 *      Int     f_code,         MV search range 1/2 pel: 1=32,2=64,...,7=2048
 *
 * Arguments in/out :
 *
 *
 * Arguments out :
 *      Image   **mot_x,        horizontal MV coordinates (Float)
 *      Image   **mot_y,        vertical MV coordinates   (Float)
 *      Image   **mode,         modes for each MB         (SInt)
 *      Image   **alpha_sub     subsampled alpha plane (blocks units) (SInt)
 *
 * Return values :
 *      none
 *
 * Side effects :
 *
 *
 * Description :
 *      1)  memory is allocated for these images.
 *      2)  mot_x/mot_y have 4 identical vector for a MB coded inter 16
 *      3)  mot_x/mot_y have 4 identical zero vectors for a MB coded intra
 *      4)  if _NO_ESTIMATION_ is used, the output is the following:
 *              - mot_x     :    all MV's are 0.0
 *              - mot_y     :    all MV's are 0.0
 *              - mode      :    all modes are MB_INTRA (IGNORING THE SHAPE)
 *              - alpha_sub :    the right alpha_sub
 *
 *      Based on: CodeOneOrTwo (TMN5, coder.c)
 *
 * See also :
 *
 *
 * Modified :
 *      05.02.97 Angel Pacheco: deletion of the valid_pel_* arrays (VM5.1)
 *               having the valid positions (padded pixels) in which the
 *               search for the motion vectors must be done (as VM4.x said).
 *      01.08.97 Fernando Jaureguiza: canged the wrong name advanced by the
 *               correct one enable_8x8_mv.
 *      05.08.97 Minhua Zhou: modified to have MC carried out on the rec. VOPs
 *
 ***********************************************************CommentEnd********/

#include <stdio.h>
#include <string.h>
#include <fcntl.h>

Void
MotionEstimation (
   Vop    *curr_vop,     /* <-- current Vop (for luminance)                  */
   Vop    *curr_rec_vop, /* <-- reconstructed current  Vop (for shape)       */
   Vop    *prev_rec_vop, /* <-- reference Vop (reconstructed)(1/2 pixel)     */
   Vop    *prev_vop,     /* <-- reference Vop (original) (pixel)             */
   Int    enable_8x8_mv, /* <-- 8x8 MV (=1) or only 16x16 MV (=0)            */
   Int    edge,          /* <-- restricted(==0)/unrestricted(==edge) mode    */
   Int    f_code,        /* <-- MV search range 1/2 pel: 1=32,2=64,...,7=2048*/
   Image  **mot_x,       /* --> horizontal MV coordinates                    */
   Image  **mot_y,       /* --> vertical MV coordinates                      */
   Image  **mode,        /* --> modes for each MB                            */
   Image  **alpha_sub    /* --> subsampled alpha plane (blocks units)        */
   )
{
  Image   *pr;             /* Reference image                   */
  Image   *pr_rec;         /* Reference image (reconstructed)   */
  Image   *alpha;          /* Alpha image                       */
  Image   *pi;             /* Interp.ref.image Y                */

  Image   *mode16;
  Image   *mv16_w;
  Image   *mv16_h;
  Image   *mv8_w;
  Image   *mv8_h;

  Image   *mosaic=NULL;
  Image   *alphamosaic=NULL;
  Image   *mosaic_edge=NULL;
  Image   *alphamosaic_edge=NULL;
  Int     edgemos=0;
  Int     offset_mos_x0=0;
  Int     offset_mos_y0=0;
  Int     no_of_sprite_points=0;
  Int     warpingaccuracy=0;

  SInt    *prev_ipol,      /* Interp.ref.image Y (subimage)            */
          *prev_orig;      /* Ref. original image with edge (subimage) */

  Int     vop_width,
          vop_height;

  Int     br_x;
  Int     br_y;
  Int     br_height;
  Int     br_width;
  Int     mv_h, mv_w;
  Char    line[100];
  Image   *mvF_w;
  Image   *mvF_h;
  int pic, type, ref, rad;

#ifdef _NO_ESTIMATION_
  /*
   * allocate memory for mv's and modes data
   *
   */

  br_height=curr_rec_vop->height;
  br_width=curr_rec_vop->width;
  mv_h=br_height/MB_SIZE;
  mv_w=br_width/MB_SIZE;

  *mode=AllocImage(mv_w,mv_h,SHORT_TYPE);
  SetConstantImage (*mode,(Float)MBM_INTRA);
  edge = 16;

  /*
  **   mv16 have 2x2 repeted the mv value to share the functions of
  **   mv prediction between CodeVopVotion and MotionEstimation
  */

  *mot_x=AllocImage(mv_w*2,mv_h*2,FLOAT_TYPE);
  *mot_y=AllocImage(mv_w*2,mv_h*2,FLOAT_TYPE);
  SetConstantImage (*mot_x,+0.0);
  SetConstantImage (*mot_y,+0.0);

  *alpha_sub=AllocImage(mv_w*2,mv_h*2,SHORT_TYPE);
  subsamp_alpha((SInt*)GetImageData(GetVopA(curr_rec_vop)), br_width,br_height,
                 1 /* block-based */, (SInt*)GetImageData(*alpha_sub));


  return;

#endif

  if (GetVopSpriteUsage(curr_vop)==STATIC_SPRITE)
    {
    fprintf(stdout,
      "\t\t\t- no motion estimation for UPDATE PIECES (STATIC SPRITE usage)\n");
    fprintf(stdout,
      "\t\t\t     but Inter/Intra decision is taken using Zero vectors\n");
    f_code=0;
    }

  /* CODE */

  /*GETBBR*/
  br_y=curr_rec_vop->ver_spat_ref;
  br_x=curr_rec_vop->hor_spat_ref;
  br_height=curr_rec_vop->height;
  br_width=curr_rec_vop->width;

  /*
  ** WE SUPPOSE prev_rec_vop & prev_vop HAVE EQUAL SIZE AND POSITIONS
  */

  vop_width=prev_rec_vop->width-2*edge;
  vop_height=prev_rec_vop->height-2*edge;

  /*
  ** Get images
  */

  alpha= GetVopA(curr_rec_vop);
  /* pr = GetVopY(prev_vop); */
  pr = GetVopY(prev_rec_vop);
  /* changed by Minhua Zhou 05.08.97 */
  pr_rec = GetVopY(prev_rec_vop);

  /*
   * Create a alpha_edge image
   *
   */

  prev_orig = (SInt*)((SInt*)GetImageData(pr) + (vop_width+2*edge)*edge+edge);

  pi = AllocImage (2*(vop_width+2*edge),2*(vop_height+2*edge),SHORT_TYPE);
  InterpolateImage(pr_rec, pi,GetVopRoundingType(curr_vop));

  prev_ipol=(SInt*)((SInt*)GetImageData(pi)+2*(vop_width+2*edge)*2*edge+2*edge);

  if (GetVopSpriteUsage (curr_vop) == ONLINE_SPRITE)
    {
    edgemos = 16;
    offset_mos_x0 = - GetVopSpriteLeftEdge(curr_vop) + edgemos;
    offset_mos_y0 = - GetVopSpriteTopEdge(curr_vop) + edgemos;
    mosaic = GetVopY(GetVopSprite(curr_vop));
    alphamosaic = GetVopA(GetVopSprite(curr_vop));

    mosaic_edge = AllocImage (GetVopSpriteHdim(curr_vop)+2*edgemos,
                              GetVopSpriteVdim(curr_vop)+2*edgemos,SHORT_TYPE);
    alphamosaic_edge = AllocImage (GetVopSpriteHdim(curr_vop)+2*edgemos,
                              GetVopSpriteVdim(curr_vop)+2*edgemos,SHORT_TYPE);

    MakeImageEdge(mosaic,edgemos,mosaic_edge);
    MakeImageEdge(alphamosaic,edgemos,alphamosaic_edge);
    RepetitiveImagePadding(alphamosaic_edge,mosaic_edge);

    no_of_sprite_points = GetVopNoOfSpritePoints(curr_vop);
    warpingaccuracy = GetVopWarpingAccuracy(curr_vop);
    }

  if (GetVopSpriteUsage (curr_vop) == GMC_SPRITE)
    {
    edgemos = 16;
    offset_mos_x0 = - prev_rec_vop->hor_spat_ref;
    offset_mos_y0 = - prev_rec_vop->ver_spat_ref;

    mosaic_edge = GetVopY(prev_rec_vop);

    no_of_sprite_points = GetVopNoOfSpritePoints(curr_vop);
    warpingaccuracy = GetVopWarpingAccuracy(curr_vop);
    }

  /*
   * allocate memory for mv's and modes data
   *
   */

  mv_h=br_height/MB_SIZE;
  mv_w=br_width/MB_SIZE;

  mode16=AllocImage (mv_w,mv_h,SHORT_TYPE);
  SetConstantImage (mode16,(Float)MBM_INTRA);

  /*
  **   mv16 have 2x2 repeted the mv value to share the functions of
  **   mv prediction between CodeVopVotion and MotionEstimation
  */

  mv16_w=AllocImage (mv_w*2,mv_h*2,FLOAT_TYPE);
  mv16_h=AllocImage (mv_w*2,mv_h*2,FLOAT_TYPE);
  mv8_w =AllocImage (mv_w*2,mv_h*2,FLOAT_TYPE);
  mv8_h =AllocImage (mv_w*2,mv_h*2,FLOAT_TYPE);
  mvF_w =AllocImage (mv_w*2,mv_h*2,FLOAT_TYPE);
  mvF_h =AllocImage (mv_w*2,mv_h*2,FLOAT_TYPE);
  SetConstantImage (mvF_h,+0.0);
  SetConstantImage (mvF_w,+0.0);
  SetConstantImage (mv16_h,+0.0);
  SetConstantImage (mv16_w,+0.0);
  SetConstantImage (mv8_h,+0.0);
  SetConstantImage (mv8_w,+0.0);

  /*
   * 8x8 subsampled alpha plane
   */

  *alpha_sub=AllocImage (mv_w*2,mv_h*2,SHORT_TYPE);

  /* Compute motion vectors and modes for each MB
  ** (TMN5 functions)
  */

  switch (curr_vop->mvfileusage)
    {
    case 1:
      (*curr_vop->mvlinenop)++;
      if (fgets(line, sizeof line, curr_vop->mvfile) == NULL)
        {
        fprintf(stderr, "EOF on %s at %d\n",
                curr_vop->mvfilename, *curr_vop->mvlinenop);
        exit(1);
        }
      if ((sscanf(line + 4, "%d %d %d %d", &pic, &type, &ref, &rad) != 4) ||
          (pic != curr_vop->frame) ||
          (type != curr_vop->prediction_type) ||
          (ref != prev_vop->frame))
        {
        fprintf(stderr, "MV read error: expected %c%d(%d), got %c%d(%d)\n",
                "IPBS"[curr_vop->prediction_type],
                curr_vop->frame, prev_vop->frame,
                "IPBS"[type], pic, ref);
        exit(1);
        }                        /* It should be 'vop->fcode_for' */
      if (rad > (8 << f_code))
        {
        fprintf(stderr, "ME radius from file (%d) > current ME radius (%d)\n",
                rad, (8 << f_code));
        exit(1);
        }
      break;
    case 2:
      (*curr_vop->mvlinenop)++;
      fprintf(curr_vop->mvfile, "Pic %d %d %d %d\n", curr_vop->frame,
              curr_vop->prediction_type, prev_vop->frame,
              curr_vop->sr_for);
      break;
    }
  MotionEstimatePicture(
             curr_vop,
             (SInt*)GetImageData(pr_rec) + edge * (vop_width + 2*edge + 1),
             (SInt*)GetImageData(alpha),
             prev_orig,                  /* Y padd with edge */
             prev_ipol,                  /* Y interpolated (from pi) */
             prev_rec_vop->hor_spat_ref+edge,
             prev_rec_vop->ver_spat_ref+edge,
             vop_width,vop_height,
             mosaic_edge,
             (Vop*)GetVopSprite(curr_vop),
             no_of_sprite_points, warpingaccuracy,
             offset_mos_x0, offset_mos_y0,
             enable_8x8_mv,
             edge,
             f_code,
             br_x,br_y,                  /* pos. of the bounding rectangle */
             br_width,br_height,         /* dims. of the bounding rectangle */
             curr_vop->ref_sel_code,    /*for Scalability (UPS)*/
             (Float*)GetImageData(mv16_w),
             (Float*)GetImageData(mv16_h),
             (Float*)GetImageData(mv8_w),
             (Float*)GetImageData(mv8_h),
             (Float*)GetImageData(mvF_w),
             (Float*)GetImageData(mvF_h),
             (SInt*) GetImageData(mode16),
             (SInt*) GetImageData(*alpha_sub)
             );

  /*
   * We don't generate here the predicted vop from MV or mv_x/mv_y
   */

  /* Convert output to MOMUSYS format
   *
   */

  if (GetMotionImages(mv16_w, mv16_h, mv8_w, mv8_h,
      mvF_w, mvF_h, mode16, mot_x, mot_y, mode) != 1)
    {
    fprintf(stderr,"GetMotionImages failed\n");
    return;
    }

  /* Deallocate dynamic memory
   */

  FreeImage(mv16_w); FreeImage(mv16_h);
  FreeImage(mv8_w);  FreeImage(mv8_h);
  FreeImage(mvF_w);  FreeImage(mvF_h);
  FreeImage(mode16);
  FreeImage(pi);
  if (GetVopSpriteUsage (curr_vop) == ONLINE_SPRITE)
    {
    FreeImage(mosaic_edge);
    FreeImage(alphamosaic_edge);
    }
}

/***********************************************************CommentBegin******
 *
 * --Load_ref_block--
 *
 * Author :
 *      HHI - Minhua Zhou
 *
 * Created :
 *      08.08.97
 *
 * Purpose :
 *      Loads the reference area for MC
 *
 * Arguments in :
 *      SInt *prev_Y_padded  reference luminance
 *      Int  prev_width      width of prev_Y_padded
 *      Int  prev_height     height of prev_Y_padded
 *      Int  Sr_x,Int Sr_y   search range
 *      Int  npix            width of ref_blk
 *      Int  col,Int row     spatial position of current block wrt prev_Y_padded
 *
 * Arguments in/out :
 *
 *
 * Arguments out :
 *      SInt *ref_blk        loaded reference area for the current MB
  *
 * Return values :
 *      none
 *
 * Side effects :
 *
 *
 * Description :
 *
 *
 * See also :
 *
 *
 * Modified :
 *
 *
 ***********************************************************CommentEnd********/

Void
Load_ref_block(
   SInt   *prev_Y_padded,
   Int    prev_width,
   Int    prev_height,
   Int    Sr_x,
   Int    Sr_y,
   Int    npix,
   Int    col,
   Int    row,
   SInt   *ref_blk
   )
{
  Int   k,l,index0,index,offset1,offset;
  Int   Sr_x_low,Sr_x_upper,Sr_y_low,Sr_y_upper;
  Int   col_ref,row_ref;

  if (Sr_x+16+col<=0)
    Sr_x_low = Sr_x_upper = Sr_x+16;
  else if (-Sr_x+col>=prev_width)
    Sr_x_low = Sr_x_upper = -Sr_x;
  else
    {
    if (-Sr_x+col<0)
      Sr_x_low = -col;
    else
      Sr_x_low =-Sr_x;
    if (Sr_x+col+16>prev_width)
      Sr_x_upper =prev_width-col;
    else
      Sr_x_upper=Sr_x+16;
    }

  if (Sr_y+16+row<=0)
    Sr_y_low = Sr_y_upper = Sr_y+16;
  else if (-Sr_y+row>=prev_height)
    Sr_y_low = Sr_y_upper = -Sr_y;
  else
    {
    if (-Sr_y+row<0)
      Sr_y_low = -row;
    else
      Sr_y_low =-Sr_y;
    if (Sr_y+row+16>prev_height)
      Sr_y_upper =prev_height-row;
    else
      Sr_y_upper=Sr_y+16;
    }

  offset = Sr_y*npix+Sr_x;

  for (k=-Sr_y;k<Sr_y_low;k++)
    for (l=-Sr_x;l<Sr_x_low;l++)
      {
      index0=k*npix+l+offset;
      ref_blk[index0] = prev_Y_padded[0];
      }

  for (k=-Sr_y;k<Sr_y_low;k++)
    for (l=Sr_x_low;l<Sr_x_upper;l++)
      {
      col_ref=l+col;
      index0=k*npix+l+offset;
      ref_blk[index0] = prev_Y_padded[col_ref];
      }

  index =prev_width-1;
  for (k=-Sr_y;k<Sr_y_low;k++)
    for (l=Sr_x_upper;l<Sr_x+16;l++)
      {
      index0=k*npix+l+offset;
      ref_blk[index0] = prev_Y_padded[index];
      }

  for (k=Sr_y_low;k<Sr_y_upper;k++)
    for (l=-Sr_x;l<Sr_x_low;l++)
      {
      row_ref=k+row;
      index0=k*npix+l+offset;
      index = row_ref*prev_width;
      ref_blk[index0] = prev_Y_padded[index];
      }

  for (k=Sr_y_low;k<Sr_y_upper;k++)
    for (l=Sr_x_low;l<Sr_x_upper;l++)
      {
      col_ref=l+col;
      row_ref=k+row;
      index0=k*npix+l+offset;
      index = row_ref*prev_width+col_ref;
      ref_blk[index0] = prev_Y_padded[index];
      }

  offset1 = prev_width-1;
  for (k=Sr_y_low;k<Sr_y_upper;k++)
    for (l=Sr_x_upper;l<Sr_x+16;l++)
      {
      row_ref=k+row;
      index0=k*npix+l+offset;
      index = row_ref*prev_width+offset1;
      ref_blk[index0] = prev_Y_padded[index];
      }

  offset1 = (prev_height-1)*prev_width;
  for (k=Sr_y_upper;k<Sr_y+16;k++)
    for (l=-Sr_x;l<Sr_x_low;l++)
      {
      index0=k*npix+l+offset;
      ref_blk[index0] = prev_Y_padded[offset1];
      }


  for (k=Sr_y_upper;k<Sr_y+16;k++)
    for (l=Sr_x_low;l<Sr_x_upper;l++)
      {
      col_ref=l+col;
      row_ref=k+row;
      index0=k*npix+l+offset;
      index = offset1+col_ref;
      ref_blk[index0] = prev_Y_padded[index];
      }

  index = offset1+prev_width-1;
  for (k=Sr_y_upper;k<Sr_y+16;k++)
    for (l=-Sr_x;l<Sr_x_low;l++)
      {
      index0=k*npix+l+offset;
      ref_blk[index0] = prev_Y_padded[index];
      }
}

/***********************************************************CommentBegin******
 *
 * --B_MotionEstimation--  MotionEstimation for B-VOPs
 *
 * Author :
 *      HHI - Minhua Zhou
 *
 * Created :
 *      05.03.97
 *
 * Purpose :
 *      Estimates motion vectors for B-VOPs
 *
 * Arguments in :
 *      Vop    *curr_vop,      current original Vop
 *      Vop    *curr_rec_vop   current origibal Vop with reconstructed shape
 *      Vop    *prev_rec_vop   previous reconstructed reference Vop  (forward)
 *      Vop    *next_rec_vop   next reconstructed reference Vop (backward)
 *      Int    f_code_for      forward f_code
 *      Int    f_code_back     backward f_code
 *      Int    TRB             forward temporal distance
 *      Int    TRD             temporal distance between two consecutive P-VOPs
 *      Image  *mot_x_P        horizontal MV coordinates of the next P-VOP
 *      Image  *mot_y_P        Vertical MV coordinates of the next P-VOP
 *      Image  *MB_decisions_P
 *
 * Arguments in/out :
 *
 *
 * Arguments out :
 *      Image   **mot_x         horizontal MV coordinates (Float) for B-VOP
 *      Image   **mot_y         vertical MV coordinates   (Float) for B-VOP
 *      Image   **mode          prediction modes for each MB          (SInt)
 *      Image   **alpha_sub     subsampled alpha plane (blocks units) (SInt)
 *
 * Return values :
 *      none
 *
 * Side effects :
 *
 *
 * Description :
 *      1)  memory is allocated for these images.
 *      2)
 *          Mode = 0  <----->  forward_mode;
 *          Mode = 1  <----->  backward;
 *          Mode = 2  <----->  interpolative;
 *          Mode = 3  <----->  direct mode
 *
 * See also :
 *
 *
 * Modified :
 *      14.05.97 Fernando Jaureguizar: B_MotionEstimation_after_padding() call
 *               modified to cope with enlarged f_code range (VM7.0).
 *      05.08.97 Minhua Zhou: modified to have MC carried out on the rec. VOPs
 *      12.11.97 Bob Eifrig: added explicit direct mode search radius
 *
 ***********************************************************CommentEnd********/

Void
B_MotionEstimation (
   Vop    *curr_vop,     /* <-- current Vop (for luminance)                  */
   Vop    *curr_rec_vop, /* <-- reconstructed current  Vop (for shape)       */
   Vop    *prev_rec_vop, /* <-- reference Vop (reconstructed) (1/2 pel)      */
   Vop    *next_rec_vop, /* <-- reference Vop (reconstructed) (1/2 pel)      */
   Int    f_code_for,    /* <-- MV search range 1/2 pel: 1=32,2=64,...,7=2048*/
   Int    f_code_back,   /* <---MV search range 1/2 pel: 1=32,2=64,...,7=2048*/
   Int    TRB,
   Int    TRD,
   Image  *mot_x_P,      /* <-- horizontal MV coordinates of the next P-VOP  */
   Image  *mot_y_P,      /* <-- Vertical MV coordinates of the next P-VOP    */
   Image  *MB_decisions_P,
   Image  **mot_x,       /* --> horizontal MV coordinates                    */
   Image  **mot_y,       /* --> vertical MV coordinates                      */
   Image  **mode,        /* --> modes for each MB                            */
   Image  **alpha_sub    /* --> subsampled alpha plane (blocks units)        */
   )
{
  Image   *prev_Y_padded,*prev_rec_Y_padded;
  Image   *next_Y_padded,*next_rec_Y_padded;
  Int     edge;

  *mode=AllocImage(curr_rec_vop->width/16,curr_rec_vop->height/16,SHORT_TYPE);

  /*  storage format of B-VOP vectors

      0 | 1          0: forward vector   1: backward vector
     ---|---
      2 | 3          2: not used         3: MVD
  */

  *mot_x=AllocImage(curr_rec_vop->width/8,curr_rec_vop->height/8,FLOAT_TYPE);
  *mot_y=AllocImage(curr_rec_vop->width/8,curr_rec_vop->height/8,FLOAT_TYPE);
  SetConstantImage (*mot_x,+0.0);
  SetConstantImage (*mot_y,+0.0);

  *alpha_sub=AllocImage(curr_rec_vop->width/8,curr_rec_vop->height/8,
             SHORT_TYPE);
  subsamp_alpha((SInt*)GetImageData(GetVopA(curr_rec_vop)),
                curr_rec_vop->width,curr_rec_vop->height,
                1 /* block-based */, (SInt*)GetImageData(*alpha_sub));

  edge=16;

  /*   Padding of the reference VOPs */

  prev_rec_Y_padded = GetVopY(prev_rec_vop);
  prev_Y_padded = GetVopY(prev_rec_vop);
  next_rec_Y_padded = GetVopY(next_rec_vop);
  next_Y_padded = GetVopY(next_rec_vop);

  B_MotionEstimation_after_padding(
               curr_rec_vop->hor_spat_ref,
               curr_rec_vop->ver_spat_ref,
               curr_rec_vop->width,
               curr_rec_vop->height,
               (SInt*)GetImageData(GetVopY(curr_vop)),
               (SInt*)GetImageData(GetVopA(curr_rec_vop)),
               (SInt*)GetImageData(*alpha_sub),
               prev_rec_vop->hor_spat_ref,
               prev_rec_vop->ver_spat_ref,
               prev_rec_vop->width,
               prev_rec_vop->height,
               (Int)pow(2.0,(Double)(f_code_for+3)),
               (SInt*)GetImageData(prev_Y_padded),
               (SInt*)GetImageData(prev_rec_Y_padded),
               curr_vop->scalability, /* for spatial scal (UPS) */
               GetVopRefSelCode(curr_vop),
               next_rec_vop->hor_spat_ref,
               next_rec_vop->ver_spat_ref,
               next_rec_vop->width,
               next_rec_vop->height,
               (Int)pow(2.0,(Double)(f_code_back+3)),
               (SInt*)GetImageData(next_Y_padded),
               (SInt*)GetImageData(next_rec_Y_padded),
               TRB,
               TRD,
               edge,
               (Float*)GetImageData(mot_x_P),
               (Float*)GetImageData(mot_y_P),
               (SInt *)GetImageData(MB_decisions_P),
               (Float*)GetImageData(*mot_x),
               (Float*)GetImageData(*mot_y),
               (SInt*) GetImageData(*mode),
               GetVopEnhanceType(curr_vop),
               curr_vop->sr_direct);
}

/***********************************************************CommentBegin******
 *
 * --B_MotionEstimation_after_padding-- MotionEstimation for B-VOPs after
 *      padding
 *
 * Author :
 *      HHI - Minhua Zhou
 *
 * Created :
 *      05.03.97
 *
 * Purpose :
 *      Estimates motion vectors for B-VOPs
 *
 * Arguments in :
 *      Int   curr_hor_spat_ref   horizont. spatial reference of the current VOP
 *      Int   curr_ver_spat_ref   vertical spatial reference of the current VOP
 *      Int   curr_vop_width      horizontal resolution of the current VOP
 *      Int   curr_vop_height     vertical resolution of the current VOP
 *      SInt  *curr_Y             luminance component of the current VOP
 *      SInt  *curr_alpha         alpha plane of the current VOP (reconstructed)
 *      SInt  *alpha_sub          sub-sampled alpha plane of the current VOP
 *                                (reconstructed)
 *      Int   prev_hor_spat_ref   horizontal spatial reference of the previous
 *                                padded VOP
 *      Int   prev_ver_spat_ref   vertical spatial reference of the previous
 *                                padded VOP
 *      Int   prev_vop_width      horizontal resolution of the previous padded
 *                                VOP
 *      Int   prev_vop_height     vertical resolution of the previous padded VOP
 *      Int   Sr_for              forward search range
 *      SInt  *prev_Y_padded      luminance component of the previous padded VOP
 *      SInt  *prev_rec_Y_padded  luminance component of the previous padded VOP
 *                                (reconstructed)
 *      Int   scal                scalability flag (for scal (UPS)
 *      Int   ref_sel_code        ref_sel_code:0 in BVop spatial scal
 *      Int   next_hor_spat_ref   horizontal spatial reference of the next
 *                                padded VOP
 *      Int   next_ver_spat_ref   vertical spatial reference of the next padded
 *                                VOP
 *      Int   next_vop_width      horizontal resolution of the next padded VOP
 *      Int   next_vop_height     vertical resolution of the next padded VOP
 *      Int   Sr_back             backward search range
 *      SInt  *next_Y_padded      luminance component of the next padded VOP
 *      SInt  *next_rec_Y_padded  luminance component of the next padded VOP
 *                                (reconstructed)
 *      Int   TRB                 forward temporal distance
 *      Int   TRD                 temporal distance between two consecutive
 *                                P-VOPs
 *      Int   edge                spatial offset for padding
 *      Float *mot_x_P            horizontal MV coordinates of the next P-VOP
 *      Float *mot_y_P            vertical MV coordinates of the next P-VOP
 *
 * Arguments in/out :
 *
 *
 * Arguments out :
 *      Image **mot_x             horizontal MV coordinates (Float) for B-VOP
 *      Image **mot_y             vertical MV coordinates   (Float) for B-VOP
 *      Image **mode              prediction modes for each MB         (SInt)
 *
 * Return values :
 *      none
 *
 * Side effects :
 *
 *
 * Description :
 *      Mode = 0  <----->  forward_mode;
 *      Mode = 1  <----->  backward;
 *      Mode = 2  <----->  interpolative;
 *      Mode = 3  <----->  direct mode
 *
 * See also :
 *
 *
 * Modified :
 *      12.11.97 Bob Eifrig: added explicit direct mode search radius
 *
 ************************************************************CommentEnd********/

Void
B_MotionEstimation_after_padding(
   Int   curr_hor_spat_ref,
   Int   curr_ver_spat_ref,
   Int   curr_vop_width,
   Int   curr_vop_height,
   SInt  *curr_Y,
   SInt  *curr_alpha,
   SInt  *alpha_sub,
   Int   prev_hor_spat_ref,
   Int   prev_ver_spat_ref,
   Int   prev_vop_width,
   Int   prev_vop_height,
   Int   Sr_for,
   SInt  *prev_Y_padded,
   SInt  *prev_rec_Y_padded,
   Int   scal, /* for scal (UPS)*/
   Int   ref_sel_code, /*end of scal (UPS)*/
   Int   next_hor_spat_ref,
   Int   next_ver_spat_ref,
   Int   next_vop_width,
   Int   next_vop_height,
   Int   Sr_back,
   SInt  *next_Y_padded,
   SInt  *next_rec_Y_padded,
   Int   TRB,
   Int   TRD,
   Int   edge,
   Float *mot_x_P,
   Float *mot_y_P,
   SInt  *mode_P,
   Float *mot_x,
   Float *mot_y,
   SInt  *mode,
   Int   enhancement_type,
   Int   sr_direct)
{
  Int   i,j,MB_num_x,MB_num_y;
  SInt  vx,vy,vx_hpel[4],vy_hpel[4];
  Int   k,l,err[4],row,col,npix_for,npix_back,index,index0;
  SInt  blk16[256];
  SInt  mask16[256],pred_for[256],pred_back[256];
  SInt  *ref_blk_for,*ref_blk_back;
  SInt  vx_P[4],vy_P[4],MVD_x,MVD_y;
  Int   Sr_for_x,Sr_for_y,Sr_back_x,Sr_back_y;
  Int   spat_offset_x_for,spat_offset_y_for,NB,sum_min;
  Int   spat_offset_x_back,spat_offset_y_back,col1,row1;
  Int   offset,co_located_MB_mode;
  SInt  sign[5];

  Sr_for_x=Sr_for_y=Sr_for;
  Sr_back_x=Sr_back_y=Sr_back;


  MB_num_x=curr_vop_width/16;
  MB_num_y=curr_vop_height/16;
  npix_for=18+2*Sr_for_x;
  ref_blk_for=(SInt *)malloc(npix_for*(18+2*Sr_for_y)*sizeof(ref_blk_for));
  npix_back=18+2*Sr_back_x;
  ref_blk_back=(SInt *)malloc(npix_back*(18+2*Sr_back_y)*sizeof(ref_blk_back));
  spat_offset_x_for=curr_hor_spat_ref-prev_hor_spat_ref;
  spat_offset_y_for=curr_ver_spat_ref-prev_ver_spat_ref;
  spat_offset_x_back=curr_hor_spat_ref-next_hor_spat_ref;
  spat_offset_y_back=curr_ver_spat_ref-next_ver_spat_ref;

  for (i=0;i<MB_num_y;i++)
    {
    row= i*16;
    for (j=0;j<MB_num_x;j++)
      {
      col=j*16; NB = 0;
      index = (i+i)*MB_num_x*2+j+j;
      mode[i*MB_num_x+j] = -1;

      sign[0] = alpha_sub[index];
      sign[1] = alpha_sub[index+1];
      sign[2] = alpha_sub[index+2*MB_num_x];
      sign[3] = alpha_sub[index+2*MB_num_x+1];
      if ((sign[0]==MBM_TRANSPARENT)&&(sign[1]==MBM_TRANSPARENT)&&
          (sign[2]==MBM_TRANSPARENT)&&(sign[3]==MBM_TRANSPARENT))
        sign[4] = MBM_TRANSPARENT;
      else if ((sign[0]==MBM_OPAQUE)&&(sign[1]==MBM_OPAQUE)&&
               (sign[2]==MBM_OPAQUE)&&(sign[3]==MBM_OPAQUE))
        sign[4] = MBM_OPAQUE;
      else
        sign[4] = MBM_BOUNDARY;

      if (sign[4]!=MBM_TRANSPARENT)
        {
        co_located_MB_mode = serach_for_vectors_of_co_located_MB(
                  mot_x_P,mot_y_P,mode_P,
                  col/*+spat_offset_x_back-edge*/,
                  row/*+spat_offset_y_back-edge*/,
                  next_vop_width-2*edge, next_vop_height-2*edge,
                  vx_P,vy_P);

        if ((co_located_MB_mode==MBM_SKIPPED)&&
            (!scal || (scal == 1 && (ref_sel_code == 3 &&
                                     enhancement_type == 0))))
          {
          /* in this case, the current B-MB is also skipped */
          index = i*MB_num_x+j;
          index0 =4*i*MB_num_x+j+j;
          mode[index] = 4;     /* forward mode, but skipped */
          mot_x[index0] = 0;
          mot_y[index0] = 0;
          }
        else
          {
          offset =row*curr_vop_width+col;
          for (k=0;k<16;k++)
            for (l=0;l<16;l++)
              {
              index = k*curr_vop_width+l+offset;
              index0=k*16+l;
              blk16[index0]=curr_Y[index];
              mask16[index0]=curr_alpha[index];
              }
          if (sign[4]==MBM_OPAQUE)
            NB=256;
          else
            for (k=0;k<256;k++)
              if (mask16[k])
                NB++;

          /* forward motion estimation */

          col1=col+spat_offset_x_for;
          row1=row+spat_offset_y_for;

          Load_ref_block(prev_Y_padded,prev_vop_width,
                         prev_vop_height,Sr_for_x+1,Sr_for_y+1,npix_for,
                         col1,row1,
                         ref_blk_for);

          MB_pel_motion_estimation(blk16,mask16,ref_blk_for,
                                   -Sr_for_x,Sr_for_x,
                                   -Sr_for_y,Sr_for_y,
                                   Sr_for_x,Sr_for_y,npix_for,NB,4,sign[4],
                                   &vx,&vy,&err[0]);
          MB_halfpel_motion_estimation(blk16,mask16,ref_blk_for,
                                       npix_for,Sr_for_x+1,Sr_for_y+1,
                                       Sr_for_x+Sr_for_x,Sr_for_y+Sr_for_y,
                                       4,sign[4],vx,vy,
                                       &vx_hpel[0],&vy_hpel[0],&err[0]);

          /* backward motion estimation */
          col1=col+spat_offset_x_back;
          row1=row+spat_offset_y_back;


          if(scal==1&&ref_sel_code==0)  /* B-VOP spatial scalability (UPS) */
            {
            Load_ref_block(next_Y_padded,next_vop_width,
                           next_vop_height,0,0,npix_back,
                           col1,row1,
                           ref_blk_back);
            err[1] =0;
            for(k=0;k<256;k++)
              if (mask16[k])
                err[1]+=abs(blk16[k]-ref_blk_back[k]);
            vx_hpel[1]=vy_hpel[1]=0;
            }
          else
            {
            Load_ref_block(next_Y_padded,next_vop_width,
                           next_vop_height,Sr_back_x+1,Sr_back_y+1,npix_back,
                           col1,row1,
                           ref_blk_back);
            MB_pel_motion_estimation(blk16,mask16,ref_blk_back,
                                     -Sr_back_x,Sr_back_x,
                                     -Sr_back_y,Sr_back_y,
                                     Sr_back_x,Sr_back_y,npix_back,NB,4,sign[4],
                                     &vx,&vy,&err[1]);

            MB_halfpel_motion_estimation(blk16,mask16,ref_blk_back,
                                         npix_back,Sr_back_x+1,Sr_back_y+1,
                                         Sr_back_x+Sr_back_x,
                                         Sr_back_y+Sr_back_y,
                                         4,sign[4],vx,vy,
                                         &vx_hpel[1],&vy_hpel[1],&err[1]);
            }

          /* interpolative motion estimation */

          halfpel_motion_compensation(prev_rec_Y_padded,vx_hpel[0],vy_hpel[0],
                                      col+spat_offset_x_for,
                                      row+spat_offset_y_for,
                                      prev_vop_width,prev_vop_height,
                                      16,16,2,0,16,pred_for);
          halfpel_motion_compensation(next_rec_Y_padded,vx_hpel[1],vy_hpel[1],
                                      col+spat_offset_x_back,
                                      row+spat_offset_y_back,
                                      next_vop_width,next_vop_height,
                                      16,16,2,0,16,pred_back);
          err[2] = 0;
          for(k=0;k<256;k++)
            if (mask16[k])
              {
              pred_for[k]=(pred_for[k]+pred_back[k]+1)>>1;
              err[2]+=abs(blk16[k]-pred_for[k]);
              }


          /* Motion Estimation for Direct Mode */
          /* reorgnization of mask16 and blk16 */

          if ( ( !(scal == 1) ||                  /* base B-Vop (UPS) */
               (scal == 1 && ref_sel_code==3) )   /* tps  B-Vop Case1 (TPS) */
               &&(co_located_MB_mode==MBM_INTER16))
            {
            reorgnization(mask16);reorgnization(blk16);

            err[3] = 200000;
            for (MVD_y = -sr_direct; MVD_y <= sr_direct; MVD_y++)
              for (MVD_x = -sr_direct; MVD_x <= sr_direct; MVD_x++)
                {
                Int sum=0;
                for (k=0;k<4;k++)
                  {
                  vx = (TRB*vx_P[k])/TRD+MVD_x;
                  vy = (TRB*vy_P[k])/TRD+MVD_y;
                  halfpel_motion_compensation (
                      prev_rec_Y_padded,vx,vy,
                      col+spat_offset_x_for+(k%2)*8,
                      row+spat_offset_y_for+(k/2)*8,
                      prev_vop_width,prev_vop_height,
                      8,8,2,0,16,&pred_for[k*64]);
                  vx = (MVD_x)? vx-vx_P[k]:((TRB-TRD)*vx_P[k])/TRD;
                  vy = (MVD_y)? vy-vy_P[k]:((TRB-TRD)*vy_P[k])/TRD;
                  halfpel_motion_compensation (
                      next_rec_Y_padded,vx,vy,
                      col+spat_offset_x_back+(k%2)*8,
                      row+spat_offset_y_back+(k/2)*8,
                      next_vop_width,next_vop_height,
                      8,8,2,0,16,&pred_back[k*64]);
                  }
                for (k=0;k<256;k++)
                  if (mask16[k])
                    {
                    pred_for[k] = (pred_for[k]+pred_back[k]+1)>>1;
                    sum+= abs(blk16[k]-pred_for[k]);
                    }
                if (sum<err[3])
                  {
                  err[3]=sum;
                  vx_hpel[3] = MVD_x;
                  vy_hpel[3] = MVD_y;
                  }
                else if (sum==err[3])
                  {
                  if (abs(MVD_x)+abs(MVD_y)<abs(vx_hpel[3])+abs(vy_hpel[3]))
                    {
                    vx_hpel[3] = MVD_x;
                    vy_hpel[3] = MVD_y;
                    }
                  }
                }
            }
          else
            err[3] = 200000;

          err[3] -= NB/2-1;
          sum_min = err[0];
          for (k=1;k<4;k++)
            if (sum_min>err[k])
              sum_min=err[k];

          index = i*MB_num_x+j;
          index0 =4*i*MB_num_x+j+j;

          if (err[3]==sum_min)
            {
            mode[index] = 3;   /* direct mode */
            mot_x[index0+2*MB_num_x+1] = (Float)vx_hpel[3]/2.0;
            mot_y[index0+2*MB_num_x+1] = (Float)vy_hpel[3]/2.0;
            }
          else if (err[2]==sum_min)
            {
            mode[index] = 2;   /* interpolate mode */
            mot_x[index0] = (Float)vx_hpel[0]/2.0;
            mot_y[index0] = (Float)vy_hpel[0]/2.0;
            mot_x[index0+1] = (Float)vx_hpel[1]/2.0;
            mot_y[index0+1] = (Float)vy_hpel[1]/2.0;
            }
          else if (err[1]==sum_min)
            {
            mode[index] = 1;  /* backward mode */
            mot_x[index0+1] = (Float)vx_hpel[1]/2.0;
            mot_y[index0+1] = (Float)vy_hpel[1]/2.0;
            }
          else
            {
            mode[index] = 0;     /* forward mode */
            mot_x[index0] = (Float)vx_hpel[0]/2.0;
            mot_y[index0] = (Float)vy_hpel[0]/2.0;
            }
          }
        }
      }
    }
  free(ref_blk_for);
  free(ref_blk_back);
}

/***********************************************************CommentBegin******
 * --- MB_pel_motion_estimation --- Motion Estimation (forward/backward) with
 *      pel accuracy for a 16x16 block
 *
 * Author :
 *      HHI - Minhua Zhou
 *
 * Created :
 *      05.03.97
 *
 * Purpose :
 *      Estimates motion vectors for a 16x16 block
 *
 * Arguments in :
 *      SInt *blk16     16x16 luminance block
 *      SInt *mask16    16x16 block of alpha plane
 *      SInt *ref_blk   reference block
 *      Int  Sr_x_low
 *      Int  Sr_x_upper
 *      Int  Sr_y_low
 *      Int  Sr_y_upper
 *      Int  Sr_x       horizontal search range
 *      Int  Sr_y       vertical serarch range
 *      Int  npix       horizontal resolution of the reefrence block
 *      Int  NB         Number of pixels belonging to VOP in the block
 *      Int  blk_nr
 *      SInt sign
 *
 * Arguments in/out :
 *
 *
 * Arguments out :
 *      SInt *vx_out    horizontal component of MV
 *      SInt *vy_out    vertical component of MV
 *      Int  *err
 *
 * Return values :
 *      none
 *
 * Side effects :
 *
 *
 * Description :
 *      The used method the full search block matching
 *
 * See also :
 *
 *
 * Modified :
 *
 *
 ***********************************************************CommentEnd********/

Void
MB_pel_motion_estimation(
   SInt *blk16,
   SInt *mask16,
   SInt *ref_blk,
   Int  Sr_x_low,
   Int  Sr_x_upper,
   Int  Sr_y_low,
   Int  Sr_y_upper,
   Int  Sr_x,
   Int  Sr_y,
   Int  npix,
   Int  NB,
   Int  blk_nr,
   SInt sign, /* 0 , 1,  2, 3,  4 <--> 16x16*/
   SInt *vx_out,
   SInt *vy_out,
   Int  *err
   )
{
  SInt         *i0,*j0,*k0,*k1;
  Int          index1;
  Int          sum,sum_min=0,offset;
  static Int   k1_start0[5]={0,8,128,136,0};
  static Int   k1_end0[5]  ={128,128,256,256,256};
  static Int   pos[5] ={0,0,8,8,0};
  SInt         *k1_start=k1_start0[blk_nr]+blk16;
  SInt         *k1_end=k1_end0[blk_nr]+blk16;
  SInt         *msk;
  Int          mask_offset;
  Int          i,j,vx,vy;
  Int          error;

  vx = vy =0;
  offset =(Sr_y+1)*npix+Sr_x+1;
  if (blk_nr==4)
    {
    for (i=0;i<16;i++)
      for (j=0;j<16;j++)
        {
        index1 = 16*i+j;
        if (mask16[index1])
          {
          error = blk16[index1]-ref_blk[i*npix+j+offset];
          if (error<0)
            sum_min-=error;else sum_min+=error;
          }
        }
    sum_min -= NB/2-1;
    }
  else
    sum_min = 2000000;

  mask_offset=mask16-blk16;
  i0 = (Sr_y+Sr_y_low+pos[blk_nr])*npix+Sr_x+Sr_x_low+(blk_nr%2)*8+ref_blk;
  if (blk_nr==4)
    {
    if (sign==MBM_OPAQUE)
      for (i=Sr_y_low;i<Sr_y_upper;i++)
        {
        i0+=npix;j0=i0;
        for (j=Sr_x_low;j<Sr_x_upper;j++)
          {
          j0++;
          sum=0;
          for (k0=j0,k1=k1_start;k1<k1_end;k0+=npix,k1+=16)
            {
            error = k0[0]-k1[0];
            if (error<0) sum-=error;else sum+=error;
            error = k0[1]-k1[1];
            if (error<0) sum-=error;else sum+=error;
            error = k0[2]-k1[2];
            if (error<0) sum-=error;else sum+=error;
            error = k0[3]-k1[3];
            if (error<0) sum-=error;else sum+=error;
            error = k0[4]-k1[4];
            if (error<0) sum-=error;else sum+=error;
            error = k0[5]-k1[5];
            if (error<0) sum-=error;else sum+=error;
            error = k0[6]-k1[6];
            if (error<0) sum-=error;else sum+=error;
            error = k0[7]-k1[7];
            if (error<0) sum-=error;else sum+=error;
            error = k0[8]-k1[8];
            if (error<0) sum-=error;else sum+=error;
            error = k0[9]-k1[9];
            if (error<0) sum-=error;else sum+=error;
            error = k0[10]-k1[10];
            if (error<0) sum-=error;else sum+=error;
            error = k0[11]-k1[11];
            if (error<0) sum-=error;else sum+=error;
            error = k0[12]-k1[12];
            if (error<0) sum-=error;else sum+=error;
            error = k0[13]-k1[13];
            if (error<0) sum-=error;else sum+=error;
            error = k0[14]-k1[14];
            if (error<0) sum-=error;else sum+=error;
            error = k0[15]-k1[15];
            if (error<0) sum-=error;else sum+=error;
            if (sum>sum_min) break;
            }

          if (sum<sum_min)
            {
            sum_min=sum;
            vx=j;vy=i;
            }
          else if (sum==sum_min)
            if (abs(vx)+abs(vy)>abs(i)+abs(j))
              {
              vx=j;vy=i;
              }
          }
        }
      else
        for (i=Sr_y_low;i<Sr_y_upper;i++)
          {
          i0+=npix;j0=i0;
          for (j=Sr_x_low;j<Sr_x_upper;j++)
            {
            j0++;
            sum=0;
            for (k0=j0,k1=k1_start;k1<k1_end;k0+=npix,k1+=16)
              {
              msk = k1+mask_offset;
              if(msk[0])
                {
                error = k0[0]-k1[0];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[1])
                {
                error = k0[1]-k1[1];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[2])
                {
                error = k0[2]-k1[2];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[3])
                {
                error = k0[3]-k1[3];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[4])
                {
                error = k0[4]-k1[4];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[5])
                {
                error = k0[5]-k1[5];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[6])
                {
                error = k0[6]-k1[6];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[7])
                {
                error = k0[7]-k1[7];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[8])
                {
                error = k0[8]-k1[8];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[9])
                {
                error = k0[9]-k1[9];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[10])
                {
                error = k0[10]-k1[10];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[11])
                {
                error = k0[11]-k1[11];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[12])
                {
                error = k0[12]-k1[12];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[13])
                {
                error = k0[13]-k1[13];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[14])
                {
                error = k0[14]-k1[14];
                if (error<0) sum-=error;else sum+=error;
                }
              if(msk[15])
                {
                error = k0[15]-k1[15];
                if (error<0) sum-=error;else sum+=error;
                }
              if (sum>sum_min) break;
              }

            if (sum<sum_min)
              {
              sum_min=sum;
              vx=j;vy=i;
              }
            else if (sum==sum_min)
              if (abs(vx)+abs(vy)>abs(i)+abs(j))
                {
                vx=j;vy=i;
                }
            }
          }
    }
  else
    {
    if (sign==MBM_OPAQUE)
      for (i=Sr_y_low;i<Sr_y_upper;i++)
        {
        i0+=npix;j0=i0;
        for (j=Sr_x_low;j<Sr_x_upper;j++)
          {
          j0++;
          sum=0;
          for (k0=j0,k1=k1_start;k1<k1_end;k0+=npix,k1+=16)
            {
            error = k0[0]-k1[0];
            if (error<0) sum-=error;else sum+=error;
            error = k0[1]-k1[1];
            if (error<0) sum-=error;else sum+=error;
            error = k0[2]-k1[2];
            if (error<0) sum-=error;else sum+=error;
            error = k0[3]-k1[3];
            if (error<0) sum-=error;else sum+=error;
            error = k0[4]-k1[4];
            if (error<0) sum-=error;else sum+=error;
            error = k0[5]-k1[5];
            if (error<0) sum-=error;else sum+=error;
            error = k0[6]-k1[6];
            if (error<0) sum-=error;else sum+=error;
            error = k0[7]-k1[7];
            if (error<0) sum-=error;else sum+=error;
            if (sum>sum_min) break;
            }

          if (sum<sum_min)
            {
            sum_min=sum;
            vx=j;vy=i;
            }
          else if (sum==sum_min)
            if (abs(vx)+abs(vy)>abs(i)+abs(j))
              {
              vx=j;vy=i;
              }
          }
        }
    else
      for (i=Sr_y_low;i<Sr_y_upper;i++)
        {
        i0+=npix;j0=i0;
        for (j=Sr_x_low;j<Sr_x_upper;j++)
          {
          j0++;
          sum=0;
          for (k0=j0,k1=k1_start;k1<k1_end;k0+=npix,k1+=16)
            {
            msk = k1+mask_offset;
            if(msk[0])
              {
              error = k0[0]-k1[0];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[1])
              {
              error = k0[1]-k1[1];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[2])
              {
              error = k0[2]-k1[2];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[3])
              {
              error = k0[3]-k1[3];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[4])
              {
              error = k0[4]-k1[4];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[5])
              {
              error = k0[5]-k1[5];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[6])
              {
              error = k0[6]-k1[6];
              if (error<0) sum-=error;else sum+=error;
              }
            if(msk[7])
              {
              error = k0[7]-k1[7];
              if (error<0) sum-=error;else sum+=error;
              }
            if (sum>sum_min) break;
            }

          if (sum<sum_min)
            {
            sum_min=sum;
            vx=j;vy=i;
            }
          else if (sum==sum_min)
            if (abs(vx)+abs(vy)>abs(i)+abs(j))
              {
              vx=j;vy=i;
              }
          }
        }
    }

  *vx_out = (SInt )vx;
  *vy_out = (SInt )vy;
  *err = sum_min;
}

/***********************************************************CommentBegin******
 * --- MB_halfpel_motion_estimation --- Motion Estimation (forward/backward)
 *      with half-pel accuracy for a block
 *
 * Author :
 *      HHI - Minhua Zhou
 *
 * Created :
 *      05.03.97
 *
 * Purpose :
 *      Estimates motion vectors for a block
 *
 * Arguments in :
 *      SInt *blk      a luminance block to be estimated
 *      SInt *mask     block of alpha plane
 *      SInt *ref_Y    luminance component of the padded reference VOP
 *                     (reconstructed)
 *      Int  npix      horizontal resolution of the padded reference VOP
 *      Int  offset_x  horizontal spatial offset between the block and the
 *                     reference VOP
 *      Int  offset_y  vertical spatial offset between the block and the
 *                     reference VOP
 *      Int  Sr_x      horizontal search range
 *      Int  Sr_y      vertical search range
 *      Int  bksize_x  block size of the luminance block (horizontal)
 *      Int  bksize_y  block size of the luminance block (vertical)
 *      Int  NB        number of pixels belonging to VOP in the block
 *      Int  blk_nr
 *      SInt sign
 *      SInt vx_pel    horizontal component of MV with pel accuarcy
 *      SInt vy_pel    vertical component of MV with pel accuarcy
 *
 * Arguments in/out :
 *
 *
 * Arguments out :
 *      SInt *vx_hpel_out  horizontal component of MV with half-pel accuarcy
 *      SInt *vy_hpel_out  vertical component of MV with half-pel accuarcy
 *      Int  *err          size of abosult prediction errors
 * Return values :
 *      none
 *
 * Side effects :
 *
 *
 * Description :
 *      The used method complies with the VM-definition
 *
 * See also :
 *
 *
 * Modified :
 *
 *
 ***********************************************************CommentEnd********/

Void
MB_halfpel_motion_estimation(
   SInt   *blk,
   SInt   *mask,
   SInt   *ref_Y,
   Int    npix,
   Int    offset_x,
   Int    offset_y,
   Int    Sr_x,
   Int    Sr_y,
   Int    blk_nr,
   SInt   sign,
   SInt   vx_pel,
   SInt   vy_pel,
   SInt   *vx_hpel_out,
   SInt   *vy_hpel_out,
   Int    *err
   )
{
  Int           i,k,l;
  static SInt   vyy[9]={-1,-1,-1, 0,0, 1,1,1,0};
  static SInt   vxx[9]={-1, 0, 1,-1,1,-1,0,1,0};
  static Int    k_start[5]={0,0, 8, 8, 0};
  static Int    k_end[5]=  {8,8,16,16,16};
  static Int    l_start[5]={0, 8,0, 8, 0};
  static Int    l_end[5]=  {8,16,8,16,16};

  SInt          vx[9],vy[9];
  Int           pel[9],int_pel[9],a18,a86;
  Int           index1,index2,sum[9],offset[9],error;
  Int           row0,col0;
  Int           blk_size=16;

  /*
    0 1 2
    3 8 4
    5 6 7
  */

  row0 = offset_y+vy_pel;
  col0 = offset_x+vx_pel;

  offset[0]=(row0-1)*npix+col0-1;
  offset[1]=offset[0]+1;
  offset[2]=offset[1]+1;
  offset[3]=offset[0]+npix;
  offset[8]=offset[3]+1;
  offset[4]=offset[8]+1;
  offset[5]=offset[3]+npix;
  offset[6]=offset[5]+1;
  offset[7]=offset[6]+1;

  for (i=0;i<8;i++)
    sum[i]=0; sum[8] = *err;
  if (sign==MBM_OPAQUE)
    for (k=k_start[blk_nr];k<k_end[blk_nr];k++)
      for (l=l_start[blk_nr];l<l_end[blk_nr];l++)
        {
        index1=k*blk_size+l;
        index2=k*npix+l;
        pel[0]=ref_Y[index2+offset[0]];
        pel[1]=ref_Y[index2+offset[1]];
        pel[2]=ref_Y[index2+offset[2]];
        pel[3]=ref_Y[index2+offset[3]];
        pel[4]=ref_Y[index2+offset[4]];
        pel[5]=ref_Y[index2+offset[5]];
        pel[6]=ref_Y[index2+offset[6]];
        pel[7]=ref_Y[index2+offset[7]];
        pel[8]=ref_Y[index2+offset[8]];

        a18=pel[1]+pel[8];
        a86=pel[8]+pel[6];
        int_pel[0]=(pel[0]+pel[3]+a18+2)>>2;
        int_pel[1]=(a18+1)>>1;
        int_pel[2]=(a18+pel[2]+pel[4]+2)>>2;
        int_pel[3]=(pel[3]+pel[8]+1)>>1;
        int_pel[4]=(pel[8]+pel[4]+1)>>1;
        int_pel[5]=(pel[3]+pel[5]+a86+2)>>2;
        int_pel[6]=(a86+1)>>1;
        int_pel[7]=(pel[4]+pel[7]+a86+2)>>2;
        int_pel[8]=pel[8];

        error = blk[index1]-int_pel[0];
        if (error<0) sum[0]-=error;else sum[0]+=error;
        error = blk[index1]-int_pel[1];
        if (error<0) sum[1]-=error;else sum[1]+=error;
        error = blk[index1]-int_pel[2];
        if (error<0) sum[2]-=error;else sum[2]+=error;
        error = blk[index1]-int_pel[3];
        if (error<0) sum[3]-=error;else sum[3]+=error;
        error = blk[index1]-int_pel[4];
        if (error<0) sum[4]-=error;else sum[4]+=error;
        error = blk[index1]-int_pel[5];
        if (error<0) sum[5]-=error;else sum[5]+=error;
        error = blk[index1]-int_pel[6];
        if (error<0) sum[6]-=error;else sum[6]+=error;
        error = blk[index1]-int_pel[7];
        if (error<0) sum[7]-=error;else sum[7]+=error;
        }
  else
    for (k=k_start[blk_nr];k<k_end[blk_nr];k++)
      for (l=l_start[blk_nr];l<l_end[blk_nr];l++)
        {
        index1=k*blk_size+l;
        if (mask[index1])
          {
          index2=k*npix+l;
          pel[0]=ref_Y[index2+offset[0]];
          pel[1]=ref_Y[index2+offset[1]];
          pel[2]=ref_Y[index2+offset[2]];
          pel[3]=ref_Y[index2+offset[3]];
          pel[4]=ref_Y[index2+offset[4]];
          pel[5]=ref_Y[index2+offset[5]];
          pel[6]=ref_Y[index2+offset[6]];
          pel[7]=ref_Y[index2+offset[7]];
          pel[8]=ref_Y[index2+offset[8]];

          a18=pel[1]+pel[8];
          a86=pel[8]+pel[6];
          int_pel[0]=(pel[0]+pel[3]+a18+2)>>2;
          int_pel[1]=(a18+1)>>1;
          int_pel[2]=(a18+pel[2]+pel[4]+2)>>2;
          int_pel[3]=(pel[3]+pel[8]+1)>>1;
          int_pel[4]=(pel[8]+pel[4]+1)>>1;
          int_pel[5]=(pel[3]+pel[5]+a86+2)>>2;
          int_pel[6]=(a86+1)>>1;
          int_pel[7]=(pel[4]+pel[7]+a86+2)>>2;
          int_pel[8]=pel[8];

          error = blk[index1]-int_pel[0];
          if (error<0) sum[0]-=error;else sum[0]+=error;
          error = blk[index1]-int_pel[1];
          if (error<0) sum[1]-=error;else sum[1]+=error;
          error = blk[index1]-int_pel[2];
          if (error<0) sum[2]-=error;else sum[2]+=error;
          error = blk[index1]-int_pel[3];
          if (error<0) sum[3]-=error;else sum[3]+=error;
          error = blk[index1]-int_pel[4];
          if (error<0) sum[4]-=error;else sum[4]+=error;
          error = blk[index1]-int_pel[5];
          if (error<0) sum[5]-=error;else sum[5]+=error;
          error = blk[index1]-int_pel[6];
          if (error<0) sum[6]-=error;else sum[6]+=error;
          error = blk[index1]-int_pel[7];
          if (error<0) sum[7]-=error;else sum[7]+=error;
          }
        }

  vx_pel += vx_pel;
  vy_pel += vy_pel;
  for (i=0;i<9;i++)
    {
    vx[i]=vx_pel+vxx[i];
    vy[i]=vy_pel+vyy[i];
    }

  for (i=0;i<8;i++)
    if ((vx[i]<-Sr_x)||(vx[i]>Sr_x-1)||(vy[i]<-Sr_y)||(vy[i]>Sr_y-1))
      sum[i]=sum[8]+1;
  k=0;
  for (i=1;i<9;i++)
    if (sum[i]<sum[k])
      k=i;
  *vx_hpel_out=vx[k];
  *vy_hpel_out=vy[k];
  *err = sum[k];
}
