src/gdcmmpeg2/src/mpeg2enc/motion.c

   1 /* motion.c, motion estimation                                              */
   2
   3 /* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */
   4
   5 /*
   6  * Disclaimer of Warranty
   7  *
   8  * These software programs are available to the user without any license fee or
   9  * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
  10  * any and all warranties, whether express, implied, or statuary, including any
  11  * implied warranties or merchantability or of fitness for a particular
  12  * purpose.  In no event shall the copyright-holder be liable for any
  13  * incidental, punitive, or consequential damages of any kind whatsoever
  14  * arising from the use of these programs.
  15  *
  16  * This disclaimer of warranty extends to the user of these programs and user's
  17  * customers, employees, agents, transferees, successors, and assigns.
  18  *
  19  * The MPEG Software Simulation Group does not represent or warrant that the
  20  * programs furnished hereunder are free of infringement of any third-party
  21  * patents.
  22  *
  23  * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
  24  * are subject to royalty fees to patent holders.  Many of these patents are
  25  * general enough such that they are unavoidable regardless of implementation
  26  * design.
  27  *
  28  */
  29
  30 #include <stdio.h>
  31 #include "config.h"
  32 #include "global.h"
  33
  34 /* private prototypes */
  35
  36 static void frame_ME _ANSI_ARGS_((unsigned char *oldorg, unsigned char *neworg,
  37   unsigned char *oldref, unsigned char *newref, unsigned char *cur,
  38   int i, int j, int sxf, int syf, int sxb, int syb, struct mbinfo *mbi));
  39
  40 static void field_ME _ANSI_ARGS_((unsigned char *oldorg, unsigned char *neworg,
  41   unsigned char *oldref, unsigned char *newref, unsigned char *cur,
  42   unsigned char *curref, int i, int j, int sxf, int syf, int sxb, int syb,
  43   struct mbinfo *mbi, int secondfield, int ipflag));
  44
  45 static void frame_estimate _ANSI_ARGS_((unsigned char *org,
  46   unsigned char *ref, unsigned char *mb,
  47   int i, int j,
  48   int sx, int sy, int *iminp, int *jminp, int *imintp, int *jmintp,
  49   int *iminbp, int *jminbp, int *dframep, int *dfieldp,
  50   int *tselp, int *bselp, int imins[2][2], int jmins[2][2]));
  51
  52 static void field_estimate _ANSI_ARGS_((unsigned char *toporg,
  53   unsigned char *topref, unsigned char *botorg, unsigned char *botref,
  54   unsigned char *mb, int i, int j, int sx, int sy, int ipflag,
  55   int *iminp, int *jminp, int *imin8up, int *jmin8up, int *imin8lp,
  56   int *jmin8lp, int *dfieldp, int *d8p, int *selp, int *sel8up, int *sel8lp,
  57   int *iminsp, int *jminsp, int *dsp));
  58
  59 static void dpframe_estimate _ANSI_ARGS_((unsigned char *ref,
  60   unsigned char *mb, int i, int j, int iminf[2][2], int jminf[2][2],
  61   int *iminp, int *jminp, int *imindmvp, int *jmindmvp,
  62   int *dmcp, int *vmcp));
  63
  64 static void dpfield_estimate _ANSI_ARGS_((unsigned char *topref,
  65   unsigned char *botref, unsigned char *mb,
  66   int i, int j, int imins, int jmins, int *imindmvp, int *jmindmvp,
  67   int *dmcp, int *vmcp));
  68
  69 static int fullsearch _ANSI_ARGS_((unsigned char *org, unsigned char *ref,
  70   unsigned char *blk,
  71   int lx, int i0, int j0, int sx, int sy, int h, int xmax, int ymax,
  72   int *iminp, int *jminp));
  73
  74 static int dist1 _ANSI_ARGS_((unsigned char *blk1, unsigned char *blk2,
  75   int lx, int hx, int hy, int h, int distlim));
  76
  77 static int dist2 _ANSI_ARGS_((unsigned char *blk1, unsigned char *blk2,
  78   int lx, int hx, int hy, int h));
  79
  80 static int bdist1 _ANSI_ARGS_((unsigned char *pf, unsigned char *pb,
  81   unsigned char *p2, int lx, int hxf, int hyf, int hxb, int hyb, int h));
  82
  83 static int bdist2 _ANSI_ARGS_((unsigned char *pf, unsigned char *pb,
  84   unsigned char *p2, int lx, int hxf, int hyf, int hxb, int hyb, int h));
  85
  86 static int variance _ANSI_ARGS_((unsigned char *p, int lx));
  87
  88 /*
  89  * motion estimation for progressive and interlaced frame pictures
  90  *
  91  * oldorg: source frame for forward prediction (used for P and B frames)
  92  * neworg: source frame for backward prediction (B frames only)
  93  * oldref: reconstructed frame for forward prediction (P and B frames)
  94  * newref: reconstructed frame for backward prediction (B frames only)
  95  * cur:    current frame (the one for which the prediction is formed)
  96  * sxf,syf: forward search window (frame coordinates)
  97  * sxb,syb: backward search window (frame coordinates)
  98  * mbi:    pointer to macroblock info structure
  99  *
 100  * results:
 101  * mbi->
 102  *  mb_type: 0, MB_INTRA, MB_FORWARD, MB_BACKWARD, MB_FORWARD|MB_BACKWARD
 103  *  MV[][][]: motion vectors (frame format)
 104  *  mv_field_sel: top/bottom field (for field prediction)
 105  *  motion_type: MC_FRAME, MC_FIELD
 106  *
 107  * uses global vars: pict_type, frame_pred_dct
 108  */
 109 void motion_estimation(oldorg,neworg,oldref,newref,cur,curref,
 110   sxf,syf,sxb,syb,mbi,secondfield,ipflag)
 111 unsigned char *oldorg,*neworg,*oldref,*newref,*cur,*curref;
 112 int sxf,syf,sxb,syb;
 113 struct mbinfo *mbi;
 114 int secondfield,ipflag;
 115 {
 116   int i, j;
 117
 118   /* loop through all macroblocks of the picture */
 119   for (j=0; j<height2; j+=16)
 120   {
 121     for (i=0; i<width; i+=16)
 122     {
 123       if (pict_struct==FRAME_PICTURE)
 124         frame_ME(oldorg,neworg,oldref,newref,cur,i,j,sxf,syf,sxb,syb,mbi);
 125       else
 126         field_ME(oldorg,neworg,oldref,newref,cur,curref,i,j,sxf,syf,sxb,syb,
 127           mbi,secondfield,ipflag);
 128       mbi++;
 129     }
 130
 131     if (!quiet)
 132     {
 133       putc('.',stderr);
 134       fflush(stderr);
 135     }
 136   }
 137   if (!quiet)
 138     putc('\n',stderr);
 139 }
 140
 141 static void frame_ME(oldorg,neworg,oldref,newref,cur,i,j,sxf,syf,sxb,syb,mbi)
 142 unsigned char *oldorg,*neworg,*oldref,*newref,*cur;
 143 int i,j,sxf,syf,sxb,syb;
 144 struct mbinfo *mbi;
 145 {
 146   int imin,jmin,iminf,jminf,iminr,jminr;
 147   int imint,jmint,iminb,jminb;
 148   int imintf,jmintf,iminbf,jminbf;
 149   int imintr,jmintr,iminbr,jminbr;
 150   int var,v0;
 151   int dmc,dmcf,dmcr,dmci,vmc,vmcf,vmcr,vmci;
 152   int dmcfield,dmcfieldf,dmcfieldr,dmcfieldi;
 153   int tsel,bsel,tself,bself,tselr,bselr;
 154   unsigned char *mb;
 155   int imins[2][2],jmins[2][2];
 156   int imindp,jmindp,imindmv,jmindmv,dmc_dp,vmc_dp;
 157
 158   mb = cur + i + width*j;
 159
 160   var = variance(mb,width);
 161
 162   if (pict_type==I_TYPE)
 163     mbi->mb_type = MB_INTRA;
 164   else if (pict_type==P_TYPE)
 165   {
 166     if (frame_pred_dct)
 167     {
 168       dmc = fullsearch(oldorg,oldref,mb,
 169                        width,i,j,sxf,syf,16,width,height,&imin,&jmin);
 170       vmc = dist2(oldref+(imin>>1)+width*(jmin>>1),mb,
 171                   width,imin&1,jmin&1,16);
 172       mbi->motion_type = MC_FRAME;
 173     }
 174     else
 175     {
 176       frame_estimate(oldorg,oldref,mb,i,j,sxf,syf,
 177         &imin,&jmin,&imint,&jmint,&iminb,&jminb,
 178         &dmc,&dmcfield,&tsel,&bsel,imins,jmins);
 179
 180       if (M==1)
 181         dpframe_estimate(oldref,mb,i,j>>1,imins,jmins,
 182           &imindp,&jmindp,&imindmv,&jmindmv,&dmc_dp,&vmc_dp);
 183
 184       /* select between dual prime, frame and field prediction */
 185       if (M==1 && dmc_dp<dmc && dmc_dp<dmcfield)
 186       {
 187         mbi->motion_type = MC_DMV;
 188         dmc = dmc_dp;
 189         vmc = vmc_dp;
 190       }
 191       else if (dmc<=dmcfield)
 192       {
 193         mbi->motion_type = MC_FRAME;
 194         vmc = dist2(oldref+(imin>>1)+width*(jmin>>1),mb,
 195                     width,imin&1,jmin&1,16);
 196       }
 197       else
 198       {
 199         mbi->motion_type = MC_FIELD;
 200         dmc = dmcfield;
 201         vmc = dist2(oldref+(tsel?width:0)+(imint>>1)+(width<<1)*(jmint>>1),
 202                     mb,width<<1,imint&1,jmint&1,8);
 203         vmc+= dist2(oldref+(bsel?width:0)+(iminb>>1)+(width<<1)*(jminb>>1),
 204                     mb+width,width<<1,iminb&1,jminb&1,8);
 205       }
 206     }
 207
 208     /* select between intra or non-intra coding:
 209      *
 210      * selection is based on intra block variance (var) vs.
 211      * prediction error variance (vmc)
 212      *
 213      * blocks with small prediction error are always coded non-intra
 214      * even if variance is smaller (is this reasonable?)
 215      */
 216     if (vmc>var && vmc>=9*256)
 217       mbi->mb_type = MB_INTRA;
 218     else
 219     {
 220       /* select between MC / No-MC
 221        *
 222        * use No-MC if var(No-MC) <= 1.25*var(MC)
 223        * (i.e slightly biased towards No-MC)
 224        *
 225        * blocks with small prediction error are always coded as No-MC
 226        * (requires no motion vectors, allows skipping)
 227        */
 228       v0 = dist2(oldref+i+width*j,mb,width,0,0,16);
 229       if (4*v0>5*vmc && v0>=9*256)
 230       {
 231         /* use MC */
 232         var = vmc;
 233         mbi->mb_type = MB_FORWARD;
 234         if (mbi->motion_type==MC_FRAME)
 235         {
 236           mbi->MV[0][0][0] = imin - (i<<1);
 237           mbi->MV[0][0][1] = jmin - (j<<1);
 238         }
 239         else if (mbi->motion_type==MC_DMV)
 240         {
 241           /* these are FRAME vectors */
 242           /* same parity vector */
 243           mbi->MV[0][0][0] = imindp - (i<<1);
 244           mbi->MV[0][0][1] = (jmindp<<1) - (j<<1);
 245
 246           /* opposite parity vector */
 247           mbi->dmvector[0] = imindmv;
 248           mbi->dmvector[1] = jmindmv;
 249         }
 250         else
 251         {
 252           /* these are FRAME vectors */
 253           mbi->MV[0][0][0] = imint - (i<<1);
 254           mbi->MV[0][0][1] = (jmint<<1) - (j<<1);
 255           mbi->MV[1][0][0] = iminb - (i<<1);
 256           mbi->MV[1][0][1] = (jminb<<1) - (j<<1);
 257           mbi->mv_field_sel[0][0] = tsel;
 258           mbi->mv_field_sel[1][0] = bsel;
 259         }
 260       }
 261       else
 262       {
 263         /* No-MC */
 264         var = v0;
 265         mbi->mb_type = 0;
 266         mbi->motion_type = MC_FRAME;
 267         mbi->MV[0][0][0] = 0;
 268         mbi->MV[0][0][1] = 0;
 269       }
 270     }
 271   }
 272   else /* if (pict_type==B_TYPE) */
 273   {
 274     if (frame_pred_dct)
 275     {
 276       /* forward */
 277       dmcf = fullsearch(oldorg,oldref,mb,
 278                         width,i,j,sxf,syf,16,width,height,&iminf,&jminf);
 279       vmcf = dist2(oldref+(iminf>>1)+width*(jminf>>1),mb,
 280                    width,iminf&1,jminf&1,16);
 281
 282       /* backward */
 283       dmcr = fullsearch(neworg,newref,mb,
 284                         width,i,j,sxb,syb,16,width,height,&iminr,&jminr);
 285       vmcr = dist2(newref+(iminr>>1)+width*(jminr>>1),mb,
 286                    width,iminr&1,jminr&1,16);
 287
 288       /* interpolated (bidirectional) */
 289       vmci = bdist2(oldref+(iminf>>1)+width*(jminf>>1),
 290                     newref+(iminr>>1)+width*(jminr>>1),
 291                     mb,width,iminf&1,jminf&1,iminr&1,jminr&1,16);
 292
 293       /* decisions */
 294
 295       /* select between forward/backward/interpolated prediction:
 296        * use the one with smallest mean sqaured prediction error
 297        */
 298       if (vmcf<=vmcr && vmcf<=vmci)
 299       {
 300         vmc = vmcf;
 301         mbi->mb_type = MB_FORWARD;
 302       }
 303       else if (vmcr<=vmci)
 304       {
 305         vmc = vmcr;
 306         mbi->mb_type = MB_BACKWARD;
 307       }
 308       else
 309       {
 310         vmc = vmci;
 311         mbi->mb_type = MB_FORWARD|MB_BACKWARD;
 312       }
 313
 314       mbi->motion_type = MC_FRAME;
 315     }
 316     else
 317     {
 318       /* forward prediction */
 319       frame_estimate(oldorg,oldref,mb,i,j,sxf,syf,
 320         &iminf,&jminf,&imintf,&jmintf,&iminbf,&jminbf,
 321         &dmcf,&dmcfieldf,&tself,&bself,imins,jmins);
 322
 323       /* backward prediction */
 324       frame_estimate(neworg,newref,mb,i,j,sxb,syb,
 325         &iminr,&jminr,&imintr,&jmintr,&iminbr,&jminbr,
 326         &dmcr,&dmcfieldr,&tselr,&bselr,imins,jmins);
 327
 328       /* calculate interpolated distance */
 329       /* frame */
 330       dmci = bdist1(oldref+(iminf>>1)+width*(jminf>>1),
 331                     newref+(iminr>>1)+width*(jminr>>1),
 332                     mb,width,iminf&1,jminf&1,iminr&1,jminr&1,16);
 333
 334       /* top field */
 335       dmcfieldi = bdist1(
 336                     oldref+(imintf>>1)+(tself?width:0)+(width<<1)*(jmintf>>1),
 337                     newref+(imintr>>1)+(tselr?width:0)+(width<<1)*(jmintr>>1),
 338                     mb,width<<1,imintf&1,jmintf&1,imintr&1,jmintr&1,8);
 339
 340       /* bottom field */
 341       dmcfieldi+= bdist1(
 342                     oldref+(iminbf>>1)+(bself?width:0)+(width<<1)*(jminbf>>1),
 343                     newref+(iminbr>>1)+(bselr?width:0)+(width<<1)*(jminbr>>1),
 344                     mb+width,width<<1,iminbf&1,jminbf&1,iminbr&1,jminbr&1,8);
 345
 346       /* select prediction type of minimum distance from the
 347        * six candidates (field/frame * forward/backward/interpolated)
 348        */
 349       if (dmci<dmcfieldi && dmci<dmcf && dmci<dmcfieldf
 350           && dmci<dmcr && dmci<dmcfieldr)
 351       {
 352         /* frame, interpolated */
 353         mbi->mb_type = MB_FORWARD|MB_BACKWARD;
 354         mbi->motion_type = MC_FRAME;
 355         vmc = bdist2(oldref+(iminf>>1)+width*(jminf>>1),
 356                      newref+(iminr>>1)+width*(jminr>>1),
 357                      mb,width,iminf&1,jminf&1,iminr&1,jminr&1,16);
 358       }
 359       else if (dmcfieldi<dmcf && dmcfieldi<dmcfieldf
 360                && dmcfieldi<dmcr && dmcfieldi<dmcfieldr)
 361       {
 362         /* field, interpolated */
 363         mbi->mb_type = MB_FORWARD|MB_BACKWARD;
 364         mbi->motion_type = MC_FIELD;
 365         vmc = bdist2(oldref+(imintf>>1)+(tself?width:0)+(width<<1)*(jmintf>>1),
 366                      newref+(imintr>>1)+(tselr?width:0)+(width<<1)*(jmintr>>1),
 367                      mb,width<<1,imintf&1,jmintf&1,imintr&1,jmintr&1,8);
 368         vmc+= bdist2(oldref+(iminbf>>1)+(bself?width:0)+(width<<1)*(jminbf>>1),
 369                      newref+(iminbr>>1)+(bselr?width:0)+(width<<1)*(jminbr>>1),
 370                      mb+width,width<<1,iminbf&1,jminbf&1,iminbr&1,jminbr&1,8);
 371       }
 372       else if (dmcf<dmcfieldf && dmcf<dmcr && dmcf<dmcfieldr)
 373       {
 374         /* frame, forward */
 375         mbi->mb_type = MB_FORWARD;
 376         mbi->motion_type = MC_FRAME;
 377         vmc = dist2(oldref+(iminf>>1)+width*(jminf>>1),mb,
 378                     width,iminf&1,jminf&1,16);
 379       }
 380       else if (dmcfieldf<dmcr && dmcfieldf<dmcfieldr)
 381       {
 382         /* field, forward */
 383         mbi->mb_type = MB_FORWARD;
 384         mbi->motion_type = MC_FIELD;
 385         vmc = dist2(oldref+(tself?width:0)+(imintf>>1)+(width<<1)*(jmintf>>1),
 386                     mb,width<<1,imintf&1,jmintf&1,8);
 387         vmc+= dist2(oldref+(bself?width:0)+(iminbf>>1)+(width<<1)*(jminbf>>1),
 388                     mb+width,width<<1,iminbf&1,jminbf&1,8);
 389       }
 390       else if (dmcr<dmcfieldr)
 391       {
 392         /* frame, backward */
 393         mbi->mb_type = MB_BACKWARD;
 394         mbi->motion_type = MC_FRAME;
 395         vmc = dist2(newref+(iminr>>1)+width*(jminr>>1),mb,
 396                     width,iminr&1,jminr&1,16);
 397       }
 398       else
 399       {
 400         /* field, backward */
 401         mbi->mb_type = MB_BACKWARD;
 402         mbi->motion_type = MC_FIELD;
 403         vmc = dist2(newref+(tselr?width:0)+(imintr>>1)+(width<<1)*(jmintr>>1),
 404                     mb,width<<1,imintr&1,jmintr&1,8);
 405         vmc+= dist2(newref+(bselr?width:0)+(iminbr>>1)+(width<<1)*(jminbr>>1),
 406                     mb+width,width<<1,iminbr&1,jminbr&1,8);
 407       }
 408     }
 409
 410     /* select between intra or non-intra coding:
 411      *
 412      * selection is based on intra block variance (var) vs.
 413      * prediction error variance (vmc)
 414      *
 415      * blocks with small prediction error are always coded non-intra
 416      * even if variance is smaller (is this reasonable?)
 417      */
 418     if (vmc>var && vmc>=9*256)
 419       mbi->mb_type = MB_INTRA;
 420     else
 421     {
 422       var = vmc;
 423       if (mbi->motion_type==MC_FRAME)
 424       {
 425         /* forward */
 426         mbi->MV[0][0][0] = iminf - (i<<1);
 427         mbi->MV[0][0][1] = jminf - (j<<1);
 428         /* backward */
 429         mbi->MV[0][1][0] = iminr - (i<<1);
 430         mbi->MV[0][1][1] = jminr - (j<<1);
 431       }
 432       else
 433       {
 434         /* these are FRAME vectors */
 435         /* forward */
 436         mbi->MV[0][0][0] = imintf - (i<<1);
 437         mbi->MV[0][0][1] = (jmintf<<1) - (j<<1);
 438         mbi->MV[1][0][0] = iminbf - (i<<1);
 439         mbi->MV[1][0][1] = (jminbf<<1) - (j<<1);
 440         mbi->mv_field_sel[0][0] = tself;
 441         mbi->mv_field_sel[1][0] = bself;
 442         /* backward */
 443         mbi->MV[0][1][0] = imintr - (i<<1);
 444         mbi->MV[0][1][1] = (jmintr<<1) - (j<<1);
 445         mbi->MV[1][1][0] = iminbr - (i<<1);
 446         mbi->MV[1][1][1] = (jminbr<<1) - (j<<1);
 447         mbi->mv_field_sel[0][1] = tselr;
 448         mbi->mv_field_sel[1][1] = bselr;
 449       }
 450     }
 451   }
 452
 453   mbi->var = var;
 454 }
 455
 456 /*
 457  * motion estimation for field pictures
 458  *
 459  * oldorg: original frame for forward prediction (P and B frames)
 460  * neworg: original frame for backward prediction (B frames only)
 461  * oldref: reconstructed frame for forward prediction (P and B frames)
 462  * newref: reconstructed frame for backward prediction (B frames only)
 463  * cur:    current original frame (the one for which the prediction is formed)
 464  * curref: current reconstructed frame (to predict second field from first)
 465  * sxf,syf: forward search window (frame coordinates)
 466  * sxb,syb: backward search window (frame coordinates)
 467  * mbi:    pointer to macroblock info structure
 468  * secondfield: indicates second field of a frame (in P fields this means
 469  *              that reference field of opposite parity is in curref instead
 470  *              of oldref)
 471  * ipflag: indicates a P type field which is the second field of a frame
 472  *         in which the first field is I type (this restricts predictions
 473  *         to be based only on the opposite parity (=I) field)
 474  *
 475  * results:
 476  * mbi->
 477  *  mb_type: 0, MB_INTRA, MB_FORWARD, MB_BACKWARD, MB_FORWARD|MB_BACKWARD
 478  *  MV[][][]: motion vectors (field format)
 479  *  mv_field_sel: top/bottom field
 480  *  motion_type: MC_FIELD, MC_16X8
 481  *
 482  * uses global vars: pict_type, pict_struct
 483  */
 484 static void field_ME(oldorg,neworg,oldref,newref,cur,curref,i,j,
 485   sxf,syf,sxb,syb,mbi,secondfield,ipflag)
 486 unsigned char *oldorg,*neworg,*oldref,*newref,*cur,*curref;
 487 int i,j,sxf,syf,sxb,syb;
 488 struct mbinfo *mbi;
 489 int secondfield,ipflag;
 490 {
 491   int w2;
 492   unsigned char *mb, *toporg, *topref, *botorg, *botref;
 493   int var,vmc,v0,dmcfieldi,dmc8i;
 494   /* int dmc; */
 495   int imin,jmin,imin8u,jmin8u,imin8l,jmin8l,dmcfield,dmc8,sel,sel8u,sel8l;
 496   int iminf,jminf,imin8uf,jmin8uf,imin8lf,jmin8lf,dmcfieldf,dmc8f,self,sel8uf,sel8lf;
 497   int iminr,jminr,imin8ur,jmin8ur,imin8lr,jmin8lr,dmcfieldr,dmc8r,selr,sel8ur,sel8lr;
 498   int imins,jmins,ds,imindmv,jmindmv,vmc_dp,dmc_dp;
 499
 500   w2 = width<<1;
 501
 502   mb = cur + i + w2*j;
 503   if (pict_struct==BOTTOM_FIELD)
 504     mb += width;
 505
 506   var = variance(mb,w2);
 507
 508   if (pict_type==I_TYPE)
 509     mbi->mb_type = MB_INTRA;
 510   else if (pict_type==P_TYPE)
 511   {
 512     toporg = oldorg;
 513     topref = oldref;
 514     botorg = oldorg + width;
 515     botref = oldref + width;
 516
 517     if (secondfield)
 518     {
 519       /* opposite parity field is in same frame */
 520       if (pict_struct==TOP_FIELD)
 521       {
 522         /* current is top field */
 523         botorg = cur + width;
 524         botref = curref + width;
 525       }
 526       else
 527       {
 528         /* current is bottom field */
 529         toporg = cur;
 530         topref = curref;
 531       }
 532     }
 533
 534     field_estimate(toporg,topref,botorg,botref,mb,i,j,sxf,syf,ipflag,
 535                    &imin,&jmin,&imin8u,&jmin8u,&imin8l,&jmin8l,
 536                    &dmcfield,&dmc8,&sel,&sel8u,&sel8l,&imins,&jmins,&ds);
 537
 538     if (M==1 && !ipflag)  /* generic condition which permits Dual Prime */
 539       dpfield_estimate(topref,botref,mb,i,j,imins,jmins,&imindmv,&jmindmv,
 540         &dmc_dp,&vmc_dp);
 541
 542     /* select between dual prime, field and 16x8 prediction */
 543     if (M==1 && !ipflag && dmc_dp<dmc8 && dmc_dp<dmcfield)
 544     {
 545       /* Dual Prime prediction */
 546       mbi->motion_type = MC_DMV;
 547       /*dmc = dmc_dp;*/     /* L1 metric */
 548       vmc = vmc_dp;     /* we already calculated L2 error for Dual */
 549
 550     }
 551     else if (dmc8<dmcfield)
 552     {
 553       /* 16x8 prediction */
 554       mbi->motion_type = MC_16X8;
 555       /* upper half block */
 556       vmc = dist2((sel8u?botref:topref) + (imin8u>>1) + w2*(jmin8u>>1),
 557                   mb,w2,imin8u&1,jmin8u&1,8);
 558       /* lower half block */
 559       vmc+= dist2((sel8l?botref:topref) + (imin8l>>1) + w2*(jmin8l>>1),
 560                   mb+8*w2,w2,imin8l&1,jmin8l&1,8);
 561     }
 562     else
 563     {
 564       /* field prediction */
 565       mbi->motion_type = MC_FIELD;
 566       vmc = dist2((sel?botref:topref) + (imin>>1) + w2*(jmin>>1),
 567                   mb,w2,imin&1,jmin&1,16);
 568     }
 569
 570     /* select between intra and non-intra coding */
 571     if (vmc>var && vmc>=9*256)
 572       mbi->mb_type = MB_INTRA;
 573     else
 574     {
 575       /* zero MV field prediction from same parity ref. field
 576        * (not allowed if ipflag is set)
 577        */
 578       if (!ipflag)
 579         v0 = dist2(((pict_struct==BOTTOM_FIELD)?botref:topref) + i + w2*j,
 580                    mb,w2,0,0,16);
 581       if (ipflag || (4*v0>5*vmc && v0>=9*256))
 582       {
 583         var = vmc;
 584         mbi->mb_type = MB_FORWARD;
 585         if (mbi->motion_type==MC_FIELD)
 586         {
 587           mbi->MV[0][0][0] = imin - (i<<1);
 588           mbi->MV[0][0][1] = jmin - (j<<1);
 589           mbi->mv_field_sel[0][0] = sel;
 590         }
 591         else if (mbi->motion_type==MC_DMV)
 592         {
 593           /* same parity vector */
 594           mbi->MV[0][0][0] = imins - (i<<1);
 595           mbi->MV[0][0][1] = jmins - (j<<1);
 596
 597           /* opposite parity vector */
 598           mbi->dmvector[0] = imindmv;
 599           mbi->dmvector[1] = jmindmv;
 600         }
 601         else
 602         {
 603           mbi->MV[0][0][0] = imin8u - (i<<1);
 604           mbi->MV[0][0][1] = jmin8u - (j<<1);
 605           mbi->MV[1][0][0] = imin8l - (i<<1);
 606           mbi->MV[1][0][1] = jmin8l - ((j+8)<<1);
 607           mbi->mv_field_sel[0][0] = sel8u;
 608           mbi->mv_field_sel[1][0] = sel8l;
 609         }
 610       }
 611       else
 612       {
 613         /* No MC */
 614         var = v0;
 615         mbi->mb_type = 0;
 616         mbi->motion_type = MC_FIELD;
 617         mbi->MV[0][0][0] = 0;
 618         mbi->MV[0][0][1] = 0;
 619         mbi->mv_field_sel[0][0] = (pict_struct==BOTTOM_FIELD);
 620       }
 621     }
 622   }
 623   else /* if (pict_type==B_TYPE) */
 624   {
 625     /* forward prediction */
 626     field_estimate(oldorg,oldref,oldorg+width,oldref+width,mb,
 627                    i,j,sxf,syf,0,
 628                    &iminf,&jminf,&imin8uf,&jmin8uf,&imin8lf,&jmin8lf,
 629                    &dmcfieldf,&dmc8f,&self,&sel8uf,&sel8lf,&imins,&jmins,&ds);
 630
 631     /* backward prediction */
 632     field_estimate(neworg,newref,neworg+width,newref+width,mb,
 633                    i,j,sxb,syb,0,
 634                    &iminr,&jminr,&imin8ur,&jmin8ur,&imin8lr,&jmin8lr,
 635                    &dmcfieldr,&dmc8r,&selr,&sel8ur,&sel8lr,&imins,&jmins,&ds);
 636
 637     /* calculate distances for bidirectional prediction */
 638     /* field */
 639     dmcfieldi = bdist1(oldref + (self?width:0) + (iminf>>1) + w2*(jminf>>1),
 640                        newref + (selr?width:0) + (iminr>>1) + w2*(jminr>>1),
 641                        mb,w2,iminf&1,jminf&1,iminr&1,jminr&1,16);
 642
 643     /* 16x8 upper half block */
 644     dmc8i = bdist1(oldref + (sel8uf?width:0) + (imin8uf>>1) + w2*(jmin8uf>>1),
 645                    newref + (sel8ur?width:0) + (imin8ur>>1) + w2*(jmin8ur>>1),
 646                    mb,w2,imin8uf&1,jmin8uf&1,imin8ur&1,jmin8ur&1,8);
 647
 648     /* 16x8 lower half block */
 649     dmc8i+= bdist1(oldref + (sel8lf?width:0) + (imin8lf>>1) + w2*(jmin8lf>>1),
 650                    newref + (sel8lr?width:0) + (imin8lr>>1) + w2*(jmin8lr>>1),
 651                    mb+8*w2,w2,imin8lf&1,jmin8lf&1,imin8lr&1,jmin8lr&1,8);
 652
 653     /* select prediction type of minimum distance */
 654     if (dmcfieldi<dmc8i && dmcfieldi<dmcfieldf && dmcfieldi<dmc8f
 655         && dmcfieldi<dmcfieldr && dmcfieldi<dmc8r)
 656     {
 657       /* field, interpolated */
 658       mbi->mb_type = MB_FORWARD|MB_BACKWARD;
 659       mbi->motion_type = MC_FIELD;
 660       vmc = bdist2(oldref + (self?width:0) + (iminf>>1) + w2*(jminf>>1),
 661                    newref + (selr?width:0) + (iminr>>1) + w2*(jminr>>1),
 662                    mb,w2,iminf&1,jminf&1,iminr&1,jminr&1,16);
 663     }
 664     else if (dmc8i<dmcfieldf && dmc8i<dmc8f
 665              && dmc8i<dmcfieldr && dmc8i<dmc8r)
 666     {
 667       /* 16x8, interpolated */
 668       mbi->mb_type = MB_FORWARD|MB_BACKWARD;
 669       mbi->motion_type = MC_16X8;
 670
 671       /* upper half block */
 672       vmc = bdist2(oldref + (sel8uf?width:0) + (imin8uf>>1) + w2*(jmin8uf>>1),
 673                    newref + (sel8ur?width:0) + (imin8ur>>1) + w2*(jmin8ur>>1),
 674                    mb,w2,imin8uf&1,jmin8uf&1,imin8ur&1,jmin8ur&1,8);
 675
 676       /* lower half block */
 677       vmc+= bdist2(oldref + (sel8lf?width:0) + (imin8lf>>1) + w2*(jmin8lf>>1),
 678                    newref + (sel8lr?width:0) + (imin8lr>>1) + w2*(jmin8lr>>1),
 679                    mb+8*w2,w2,imin8lf&1,jmin8lf&1,imin8lr&1,jmin8lr&1,8);
 680     }
 681     else if (dmcfieldf<dmc8f && dmcfieldf<dmcfieldr && dmcfieldf<dmc8r)
 682     {
 683       /* field, forward */
 684       mbi->mb_type = MB_FORWARD;
 685       mbi->motion_type = MC_FIELD;
 686       vmc = dist2(oldref + (self?width:0) + (iminf>>1) + w2*(jminf>>1),
 687                   mb,w2,iminf&1,jminf&1,16);
 688     }
 689     else if (dmc8f<dmcfieldr && dmc8f<dmc8r)
 690     {
 691       /* 16x8, forward */
 692       mbi->mb_type = MB_FORWARD;
 693       mbi->motion_type = MC_16X8;
 694
 695       /* upper half block */
 696       vmc = dist2(oldref + (sel8uf?width:0) + (imin8uf>>1) + w2*(jmin8uf>>1),
 697                   mb,w2,imin8uf&1,jmin8uf&1,8);
 698
 699       /* lower half block */
 700       vmc+= dist2(oldref + (sel8lf?width:0) + (imin8lf>>1) + w2*(jmin8lf>>1),
 701                   mb+8*w2,w2,imin8lf&1,jmin8lf&1,8);
 702     }
 703     else if (dmcfieldr<dmc8r)
 704     {
 705       /* field, backward */
 706       mbi->mb_type = MB_BACKWARD;
 707       mbi->motion_type = MC_FIELD;
 708       vmc = dist2(newref + (selr?width:0) + (iminr>>1) + w2*(jminr>>1),
 709                   mb,w2,iminr&1,jminr&1,16);
 710     }
 711     else
 712     {
 713       /* 16x8, backward */
 714       mbi->mb_type = MB_BACKWARD;
 715       mbi->motion_type = MC_16X8;
 716
 717       /* upper half block */
 718       vmc = dist2(newref + (sel8ur?width:0) + (imin8ur>>1) + w2*(jmin8ur>>1),
 719                   mb,w2,imin8ur&1,jmin8ur&1,8);
 720
 721       /* lower half block */
 722       vmc+= dist2(newref + (sel8lr?width:0) + (imin8lr>>1) + w2*(jmin8lr>>1),
 723                   mb+8*w2,w2,imin8lr&1,jmin8lr&1,8);
 724     }
 725
 726     /* select between intra and non-intra coding */
 727     if (vmc>var && vmc>=9*256)
 728       mbi->mb_type = MB_INTRA;
 729     else
 730     {
 731       var = vmc;
 732       if (mbi->motion_type==MC_FIELD)
 733       {
 734         /* forward */
 735         mbi->MV[0][0][0] = iminf - (i<<1);
 736         mbi->MV[0][0][1] = jminf - (j<<1);
 737         mbi->mv_field_sel[0][0] = self;
 738         /* backward */
 739         mbi->MV[0][1][0] = iminr - (i<<1);
 740         mbi->MV[0][1][1] = jminr - (j<<1);
 741         mbi->mv_field_sel[0][1] = selr;
 742       }
 743       else /* MC_16X8 */
 744       {
 745         /* forward */
 746         mbi->MV[0][0][0] = imin8uf - (i<<1);
 747         mbi->MV[0][0][1] = jmin8uf - (j<<1);
 748         mbi->mv_field_sel[0][0] = sel8uf;
 749         mbi->MV[1][0][0] = imin8lf - (i<<1);
 750         mbi->MV[1][0][1] = jmin8lf - ((j+8)<<1);
 751         mbi->mv_field_sel[1][0] = sel8lf;
 752         /* backward */
 753         mbi->MV[0][1][0] = imin8ur - (i<<1);
 754         mbi->MV[0][1][1] = jmin8ur - (j<<1);
 755         mbi->mv_field_sel[0][1] = sel8ur;
 756         mbi->MV[1][1][0] = imin8lr - (i<<1);
 757         mbi->MV[1][1][1] = jmin8lr - ((j+8)<<1);
 758         mbi->mv_field_sel[1][1] = sel8lr;
 759       }
 760     }
 761   }
 762
 763   mbi->var = var;
 764 }
 765
 766 /*
 767  * frame picture motion estimation
 768  *
 769  * org: top left pel of source reference frame
 770  * ref: top left pel of reconstructed reference frame
 771  * mb:  macroblock to be matched
 772  * i,j: location of mb relative to ref (=center of search window)
 773  * sx,sy: half widths of search window
 774  * iminp,jminp,dframep: location and value of best frame prediction
 775  * imintp,jmintp,tselp: location of best field pred. for top field of mb
 776  * iminbp,jminbp,bselp: location of best field pred. for bottom field of mb
 777  * dfieldp: value of field prediction
 778  */
 779 static void frame_estimate(org,ref,mb,i,j,sx,sy,
 780   iminp,jminp,imintp,jmintp,iminbp,jminbp,dframep,dfieldp,tselp,bselp,
 781   imins,jmins)
 782 unsigned char *org,*ref,*mb;
 783 int i,j,sx,sy;
 784 int *iminp,*jminp;
 785 int *imintp,*jmintp,*iminbp,*jminbp;
 786 int *dframep,*dfieldp;
 787 int *tselp,*bselp;
 788 int imins[2][2],jmins[2][2];
 789 {
 790   int dt,db,dmint,dminb;
 791   int imint,iminb,jmint,jminb;
 792
 793   /* frame prediction */
 794   *dframep = fullsearch(org,ref,mb,width,i,j,sx,sy,16,width,height,
 795                         iminp,jminp);
 796
 797   /* predict top field from top field */
 798   dt = fullsearch(org,ref,mb,width<<1,i,j>>1,sx,sy>>1,8,width,height>>1,
 799                   &imint,&jmint);
 800
 801   /* predict top field from bottom field */
 802   db = fullsearch(org+width,ref+width,mb,width<<1,i,j>>1,sx,sy>>1,8,width,height>>1,
 803                   &iminb,&jminb);
 804
 805   imins[0][0] = imint;
 806   jmins[0][0] = jmint;
 807   imins[1][0] = iminb;
 808   jmins[1][0] = jminb;
 809
 810   /* select prediction for top field */
 811   if (dt<=db)
 812   {
 813     dmint=dt; *imintp=imint; *jmintp=jmint; *tselp=0;
 814   }
 815   else
 816   {
 817     dmint=db; *imintp=iminb; *jmintp=jminb; *tselp=1;
 818   }
 819
 820   /* predict bottom field from top field */
 821   dt = fullsearch(org,ref,mb+width,width<<1,i,j>>1,sx,sy>>1,8,width,height>>1,
 822                   &imint,&jmint);
 823
 824   /* predict bottom field from bottom field */
 825   db = fullsearch(org+width,ref+width,mb+width,width<<1,i,j>>1,sx,sy>>1,8,width,height>>1,
 826                   &iminb,&jminb);
 827
 828   imins[0][1] = imint;
 829   jmins[0][1] = jmint;
 830   imins[1][1] = iminb;
 831   jmins[1][1] = jminb;
 832
 833   /* select prediction for bottom field */
 834   if (db<=dt)
 835   {
 836     dminb=db; *iminbp=iminb; *jminbp=jminb; *bselp=1;
 837   }
 838   else
 839   {
 840     dminb=dt; *iminbp=imint; *jminbp=jmint; *bselp=0;
 841   }
 842
 843   *dfieldp=dmint+dminb;
 844 }
 845
 846 /*
 847  * field picture motion estimation subroutine
 848  *
 849  * toporg: address of original top reference field
 850  * topref: address of reconstructed top reference field
 851  * botorg: address of original bottom reference field
 852  * botref: address of reconstructed bottom reference field
 853  * mb:  macroblock to be matched
 854  * i,j: location of mb (=center of search window)
 855  * sx,sy: half width/height of search window
 856  *
 857  * iminp,jminp,selp,dfieldp: location and distance of best field prediction
 858  * imin8up,jmin8up,sel8up: location of best 16x8 pred. for upper half of mb
 859  * imin8lp,jmin8lp,sel8lp: location of best 16x8 pred. for lower half of mb
 860  * d8p: distance of best 16x8 prediction
 861  * iminsp,jminsp,dsp: location and distance of best same parity field
 862  *                    prediction (needed for dual prime, only valid if
 863  *                    ipflag==0)
 864  */
 865 static void field_estimate(toporg,topref,botorg,botref,mb,i,j,sx,sy,ipflag,
 866   iminp,jminp,imin8up,jmin8up,imin8lp,jmin8lp,dfieldp,d8p,selp,sel8up,sel8lp,
 867   iminsp,jminsp,dsp)
 868 unsigned char *toporg, *topref, *botorg, *botref, *mb;
 869 int i,j,sx,sy;
 870 int ipflag;
 871 int *iminp, *jminp;
 872 int *imin8up, *jmin8up, *imin8lp, *jmin8lp;
 873 int *dfieldp,*d8p;
 874 int *selp, *sel8up, *sel8lp;
 875 int *iminsp, *jminsp, *dsp;
 876 {
 877   int dt, db, imint, jmint, iminb, jminb, notop, nobot;
 878
 879   /* if ipflag is set, predict from field of opposite parity only */
 880   notop = ipflag && (pict_struct==TOP_FIELD);
 881   nobot = ipflag && (pict_struct==BOTTOM_FIELD);
 882
 883   /* field prediction */
 884
 885   /* predict current field from top field */
 886   if (notop)
 887     dt = 65536; /* infinity */
 888   else
 889     dt = fullsearch(toporg,topref,mb,width<<1,
 890                     i,j,sx,sy>>1,16,width,height>>1,
 891                     &imint,&jmint);
 892
 893   /* predict current field from bottom field */
 894   if (nobot)
 895     db = 65536; /* infinity */
 896   else
 897     db = fullsearch(botorg,botref,mb,width<<1,
 898                     i,j,sx,sy>>1,16,width,height>>1,
 899                     &iminb,&jminb);
 900
 901   /* same parity prediction (only valid if ipflag==0) */
 902   if (pict_struct==TOP_FIELD)
 903   {
 904     *iminsp = imint; *jminsp = jmint; *dsp = dt;
 905   }
 906   else
 907   {
 908     *iminsp = iminb; *jminsp = jminb; *dsp = db;
 909   }
 910
 911   /* select field prediction */
 912   if (dt<=db)
 913   {
 914     *dfieldp = dt; *iminp = imint; *jminp = jmint; *selp = 0;
 915   }
 916   else
 917   {
 918     *dfieldp = db; *iminp = iminb; *jminp = jminb; *selp = 1;
 919   }
 920
 921
 922   /* 16x8 motion compensation */
 923
 924   /* predict upper half field from top field */
 925   if (notop)
 926     dt = 65536;
 927   else
 928     dt = fullsearch(toporg,topref,mb,width<<1,
 929                     i,j,sx,sy>>1,8,width,height>>1,
 930                     &imint,&jmint);
 931
 932   /* predict upper half field from bottom field */
 933   if (nobot)
 934     db = 65536;
 935   else
 936     db = fullsearch(botorg,botref,mb,width<<1,
 937                     i,j,sx,sy>>1,8,width,height>>1,
 938                     &iminb,&jminb);
 939
 940   /* select prediction for upper half field */
 941   if (dt<=db)
 942   {
 943     *d8p = dt; *imin8up = imint; *jmin8up = jmint; *sel8up = 0;
 944   }
 945   else
 946   {
 947     *d8p = db; *imin8up = iminb; *jmin8up = jminb; *sel8up = 1;
 948   }
 949
 950   /* predict lower half field from top field */
 951   if (notop)
 952     dt = 65536;
 953   else
 954     dt = fullsearch(toporg,topref,mb+(width<<4),width<<1,
 955                     i,j+8,sx,sy>>1,8,width,height>>1,
 956                     &imint,&jmint);
 957
 958   /* predict lower half field from bottom field */
 959   if (nobot)
 960     db = 65536;
 961   else
 962     db = fullsearch(botorg,botref,mb+(width<<4),width<<1,
 963                     i,j+8,sx,sy>>1,8,width,height>>1,
 964                     &iminb,&jminb);
 965
 966   /* select prediction for lower half field */
 967   if (dt<=db)
 968   {
 969     *d8p += dt; *imin8lp = imint; *jmin8lp = jmint; *sel8lp = 0;
 970   }
 971   else
 972   {
 973     *d8p += db; *imin8lp = iminb; *jmin8lp = jminb; *sel8lp = 1;
 974   }
 975 }
 976
 977 static void dpframe_estimate(ref,mb,i,j,iminf,jminf,
 978   iminp,jminp,imindmvp, jmindmvp, dmcp, vmcp)
 979 unsigned char *ref, *mb;
 980 int i,j;
 981 int iminf[2][2], jminf[2][2];
 982 int *iminp, *jminp;
 983 int *imindmvp, *jmindmvp;
 984 int *dmcp,*vmcp;
 985 {
 986   int pref,ppred,delta_x,delta_y;
 987   int is,js,it,jt,ib,jb,it0,jt0,ib0,jb0;
 988   int imins,jmins,imint,jmint,iminb,jminb,imindmv,jmindmv;
 989   int vmc,local_dist;
 990
 991   /* Calculate Dual Prime distortions for 9 delta candidates
 992    * for each of the four minimum field vectors
 993    * Note: only for P pictures!
 994    */
 995
 996   /* initialize minimum dual prime distortion to large value */
 997   vmc = 1 << 30;
 998
 999   for (pref=0; pref<2; pref++)
1000   {
1001     for (ppred=0; ppred<2; ppred++)
1002     {
1003       /* convert Cartesian absolute to relative motion vector
1004        * values (wrt current macroblock address (i,j)
1005        */
1006       is = iminf[pref][ppred] - (i<<1);
1007       js = jminf[pref][ppred] - (j<<1);
1008
1009       if (pref!=ppred)
1010       {
1011         /* vertical field shift adjustment */
1012         if (ppred==0)
1013           js++;
1014         else
1015           js--;
1016
1017         /* mvxs and mvys scaling*/
1018         is<<=1;
1019         js<<=1;
1020         if (topfirst == ppred)
1021         {
1022           /* second field: scale by 1/3 */
1023           is = (is>=0) ? (is+1)/3 : -((-is+1)/3);
1024           js = (js>=0) ? (js+1)/3 : -((-js+1)/3);
1025         }
1026         else
1027           continue;
1028       }
1029
1030       /* vector for prediction from field of opposite 'parity' */
1031       if (topfirst)
1032       {
1033         /* vector for prediction of top field from bottom field */
1034         it0 = ((is+(is>0))>>1);
1035         jt0 = ((js+(js>0))>>1) - 1;
1036
1037         /* vector for prediction of bottom field from top field */
1038         ib0 = ((3*is+(is>0))>>1);
1039         jb0 = ((3*js+(js>0))>>1) + 1;
1040       }
1041       else
1042       {
1043         /* vector for prediction of top field from bottom field */
1044         it0 = ((3*is+(is>0))>>1);
1045         jt0 = ((3*js+(js>0))>>1) - 1;
1046
1047         /* vector for prediction of bottom field from top field */
1048         ib0 = ((is+(is>0))>>1);
1049         jb0 = ((js+(js>0))>>1) + 1;
1050       }
1051
1052       /* convert back to absolute half-pel field picture coordinates */
1053       is += i<<1;
1054       js += j<<1;
1055       it0 += i<<1;
1056       jt0 += j<<1;
1057       ib0 += i<<1;
1058       jb0 += j<<1;
1059
1060       if (is >= 0 && is <= (width-16)<<1 &&
1061           js >= 0 && js <= (height-16))
1062       {
1063         for (delta_y=-1; delta_y<=1; delta_y++)
1064         {
1065           for (delta_x=-1; delta_x<=1; delta_x++)
1066           {
1067             /* opposite field coordinates */
1068             it = it0 + delta_x;
1069             jt = jt0 + delta_y;
1070             ib = ib0 + delta_x;
1071             jb = jb0 + delta_y;
1072
1073             if (it >= 0 && it <= (width-16)<<1 &&
1074                 jt >= 0 && jt <= (height-16) &&
1075                 ib >= 0 && ib <= (width-16)<<1 &&
1076                 jb >= 0 && jb <= (height-16))
1077             {
1078               /* compute prediction error */
1079               local_dist = bdist2(
1080                 ref + (is>>1) + (width<<1)*(js>>1),
1081                 ref + width + (it>>1) + (width<<1)*(jt>>1),
1082                 mb,             /* current mb location */
1083                 width<<1,       /* adjacent line distance */
1084                 is&1, js&1, it&1, jt&1, /* half-pel flags */
1085                 8);             /* block height */
1086               local_dist += bdist2(
1087                 ref + width + (is>>1) + (width<<1)*(js>>1),
1088                 ref + (ib>>1) + (width<<1)*(jb>>1),
1089                 mb + width,     /* current mb location */
1090                 width<<1,       /* adjacent line distance */
1091                 is&1, js&1, ib&1, jb&1, /* half-pel flags */
1092                 8);             /* block height */
1093
1094               /* update delta with least distortion vector */
1095               if (local_dist < vmc)
1096               {
1097                 imins = is;
1098                 jmins = js;
1099                 imint = it;
1100                 jmint = jt;
1101                 iminb = ib;
1102                 jminb = jb;
1103                 imindmv = delta_x;
1104                 jmindmv = delta_y;
1105                 vmc = local_dist;
1106               }
1107             }
1108           }  /* end delta x loop */
1109         } /* end delta y loop */
1110       }
1111     }
1112   }
1113
1114   /* Compute L1 error for decision purposes */
1115   local_dist = bdist1(
1116     ref + (imins>>1) + (width<<1)*(jmins>>1),
1117     ref + width + (imint>>1) + (width<<1)*(jmint>>1),
1118     mb,
1119     width<<1,
1120     imins&1, jmins&1, imint&1, jmint&1,
1121     8);
1122   local_dist += bdist1(
1123     ref + width + (imins>>1) + (width<<1)*(jmins>>1),
1124     ref + (iminb>>1) + (width<<1)*(jminb>>1),
1125     mb + width,
1126     width<<1,
1127     imins&1, jmins&1, iminb&1, jminb&1,
1128     8);
1129
1130   *dmcp = local_dist;
1131   *iminp = imins;
1132   *jminp = jmins;
1133   *imindmvp = imindmv;
1134   *jmindmvp = jmindmv;
1135   *vmcp = vmc;
1136 }
1137
1138 static void dpfield_estimate(topref,botref,mb,i,j,imins,jmins,
1139   imindmvp, jmindmvp, dmcp, vmcp)
1140 unsigned char *topref, *botref, *mb;
1141 int i,j;
1142 int imins, jmins;
1143 int *imindmvp, *jmindmvp;
1144 int *dmcp,*vmcp;
1145 {
1146   unsigned char *sameref, *oppref;
1147   int io0,jo0,io,jo,delta_x,delta_y,mvxs,mvys,mvxo0,mvyo0;
1148   int imino,jmino,imindmv,jmindmv,vmc_dp,local_dist;
1149
1150   /* Calculate Dual Prime distortions for 9 delta candidates */
1151   /* Note: only for P pictures! */
1152
1153   /* Assign opposite and same reference pointer */
1154   if (pict_struct==TOP_FIELD)
1155   {
1156     sameref = topref;
1157     oppref = botref;
1158   }
1159   else
1160   {
1161     sameref = botref;
1162     oppref = topref;
1163   }
1164
1165   /* convert Cartesian absolute to relative motion vector
1166    * values (wrt current macroblock address (i,j)
1167    */
1168   mvxs = imins - (i<<1);
1169   mvys = jmins - (j<<1);
1170
1171   /* vector for prediction from field of opposite 'parity' */
1172   mvxo0 = (mvxs+(mvxs>0)) >> 1;  /* mvxs // 2 */
1173   mvyo0 = (mvys+(mvys>0)) >> 1;  /* mvys // 2 */
1174
1175   /* vertical field shift correction */
1176   if (pict_struct==TOP_FIELD)
1177     mvyo0--;
1178   else
1179     mvyo0++;
1180
1181   /* convert back to absolute coordinates */
1182   io0 = mvxo0 + (i<<1);
1183   jo0 = mvyo0 + (j<<1);
1184
1185   /* initialize minimum dual prime distortion to large value */
1186   vmc_dp = 1 << 30;
1187
1188   for (delta_y = -1; delta_y <= 1; delta_y++)
1189   {
1190     for (delta_x = -1; delta_x <=1; delta_x++)
1191     {
1192       /* opposite field coordinates */
1193       io = io0 + delta_x;
1194       jo = jo0 + delta_y;
1195
1196       if (io >= 0 && io <= (width-16)<<1 &&
1197           jo >= 0 && jo <= (height2-16)<<1)
1198       {
1199         /* compute prediction error */
1200         local_dist = bdist2(
1201           sameref + (imins>>1) + width2*(jmins>>1),
1202           oppref  + (io>>1)    + width2*(jo>>1),
1203           mb,             /* current mb location */
1204           width2,         /* adjacent line distance */
1205           imins&1, jmins&1, io&1, jo&1, /* half-pel flags */
1206           16);            /* block height */
1207
1208         /* update delta with least distortion vector */
1209         if (local_dist < vmc_dp)
1210         {
1211           imino = io;
1212           jmino = jo;
1213           imindmv = delta_x;
1214           jmindmv = delta_y;
1215           vmc_dp = local_dist;
1216         }
1217       }
1218     }  /* end delta x loop */
1219   } /* end delta y loop */
1220
1221   /* Compute L1 error for decision purposes */
1222   *dmcp = bdist1(
1223     sameref + (imins>>1) + width2*(jmins>>1),
1224     oppref  + (imino>>1) + width2*(jmino>>1),
1225     mb,             /* current mb location */
1226     width2,         /* adjacent line distance */
1227     imins&1, jmins&1, imino&1, jmino&1, /* half-pel flags */
1228     16);            /* block height */
1229
1230   *imindmvp = imindmv;
1231   *jmindmvp = jmindmv;
1232   *vmcp = vmc_dp;
1233 }
1234
1235 /*
1236  * full search block matching
1237  *
1238  * blk: top left pel of (16*h) block
1239  * h: height of block
1240  * lx: distance (in bytes) of vertically adjacent pels in ref,blk
1241  * org: top left pel of source reference picture
1242  * ref: top left pel of reconstructed reference picture
1243  * i0,j0: center of search window
1244  * sx,sy: half widths of search window
1245  * xmax,ymax: right/bottom limits of search area
1246  * iminp,jminp: pointers to where the result is stored
1247  *              result is given as half pel offset from ref(0,0)
1248  *              i.e. NOT relative to (i0,j0)
1249  */
1250 static int fullsearch(org,ref,blk,lx,i0,j0,sx,sy,h,xmax,ymax,iminp,jminp)
1251 unsigned char *org,*ref,*blk;
1252 int lx,i0,j0,sx,sy,h,xmax,ymax;
1253 int *iminp,*jminp;
1254 {
1255   int i,j,imin,jmin,ilow,ihigh,jlow,jhigh;
1256   int d,dmin;
1257   int k,l,sxy;
1258
1259   ilow = i0 - sx;
1260   ihigh = i0 + sx;
1261
1262   if (ilow<0)
1263     ilow = 0;
1264
1265   if (ihigh>xmax-16)
1266     ihigh = xmax-16;
1267
1268   jlow = j0 - sy;
1269   jhigh = j0 + sy;
1270
1271   if (jlow<0)
1272     jlow = 0;
1273
1274   if (jhigh>ymax-h)
1275     jhigh = ymax-h;
1276
1277   /* full pel search, spiraling outwards */
1278
1279   imin = i0;
1280   jmin = j0;
1281   dmin = dist1(org+imin+lx*jmin,blk,lx,0,0,h,65536);
1282
1283   sxy = (sx>sy) ? sx : sy;
1284
1285   for (l=1; l<=sxy; l++)
1286   {
1287     i = i0 - l;
1288     j = j0 - l;
1289     for (k=0; k<8*l; k++)
1290     {
1291       if (i>=ilow && i<=ihigh && j>=jlow && j<=jhigh)
1292       {
1293         d = dist1(org+i+lx*j,blk,lx,0,0,h,dmin);
1294
1295         if (d<dmin)
1296         {
1297           dmin = d;
1298           imin = i;
1299           jmin = j;
1300         }
1301       }
1302
1303       if      (k<2*l) i++;
1304       else if (k<4*l) j++;
1305       else if (k<6*l) i--;
1306       else            j--;
1307     }
1308   }
1309
1310   /* half pel */
1311   dmin = 65536;
1312   imin <<= 1;
1313   jmin <<= 1;
1314   ilow = imin - (imin>0);
1315   ihigh = imin + (imin<((xmax-16)<<1));
1316   jlow = jmin - (jmin>0);
1317   jhigh = jmin + (jmin<((ymax-h)<<1));
1318
1319   for (j=jlow; j<=jhigh; j++)
1320     for (i=ilow; i<=ihigh; i++)
1321     {
1322       d = dist1(ref+(i>>1)+lx*(j>>1),blk,lx,i&1,j&1,h,dmin);
1323
1324       if (d<dmin)
1325       {
1326         dmin = d;
1327         imin = i;
1328         jmin = j;
1329       }
1330     }
1331
1332   *iminp = imin;
1333   *jminp = jmin;
1334
1335   return dmin;
1336 }
1337
1338 /*
1339  * total absolute difference between two (16*h) blocks
1340  * including optional half pel interpolation of blk1 (hx,hy)
1341  * blk1,blk2: addresses of top left pels of both blocks
1342  * lx:        distance (in bytes) of vertically adjacent pels
1343  * hx,hy:     flags for horizontal and/or vertical interpolation
1344  * h:         height of block (usually 8 or 16)
1345  * distlim:   bail out if sum exceeds this value
1346  */
1347 static int dist1(blk1,blk2,lx,hx,hy,h,distlim)
1348 unsigned char *blk1,*blk2;
1349 int lx,hx,hy,h;
1350 int distlim;
1351 {
1352   unsigned char *p1,*p1a,*p2;
1353   int i,j;
1354   int s,v;
1355
1356   s = 0;
1357   p1 = blk1;
1358   p2 = blk2;
1359
1360   if (!hx && !hy)
1361     for (j=0; j<h; j++)
1362     {
1363       if ((v = p1[0]  - p2[0])<0)  v = -v; s+= v;
1364       if ((v = p1[1]  - p2[1])<0)  v = -v; s+= v;
1365       if ((v = p1[2]  - p2[2])<0)  v = -v; s+= v;
1366       if ((v = p1[3]  - p2[3])<0)  v = -v; s+= v;
1367       if ((v = p1[4]  - p2[4])<0)  v = -v; s+= v;
1368       if ((v = p1[5]  - p2[5])<0)  v = -v; s+= v;
1369       if ((v = p1[6]  - p2[6])<0)  v = -v; s+= v;
1370       if ((v = p1[7]  - p2[7])<0)  v = -v; s+= v;
1371       if ((v = p1[8]  - p2[8])<0)  v = -v; s+= v;
1372       if ((v = p1[9]  - p2[9])<0)  v = -v; s+= v;
1373       if ((v = p1[10] - p2[10])<0) v = -v; s+= v;
1374       if ((v = p1[11] - p2[11])<0) v = -v; s+= v;
1375       if ((v = p1[12] - p2[12])<0) v = -v; s+= v;
1376       if ((v = p1[13] - p2[13])<0) v = -v; s+= v;
1377       if ((v = p1[14] - p2[14])<0) v = -v; s+= v;
1378       if ((v = p1[15] - p2[15])<0) v = -v; s+= v;
1379
1380       if (s >= distlim)
1381         break;
1382
1383       p1+= lx;
1384       p2+= lx;
1385     }
1386   else if (hx && !hy)
1387     for (j=0; j<h; j++)
1388     {
1389       for (i=0; i<16; i++)
1390       {
1391         v = ((unsigned int)(p1[i]+p1[i+1]+1)>>1) - p2[i];
1392         if (v>=0)
1393           s+= v;
1394         else
1395           s-= v;
1396       }
1397       p1+= lx;
1398       p2+= lx;
1399     }
1400   else if (!hx && hy)
1401   {
1402     p1a = p1 + lx;
1403     for (j=0; j<h; j++)
1404     {
1405       for (i=0; i<16; i++)
1406       {
1407         v = ((unsigned int)(p1[i]+p1a[i]+1)>>1) - p2[i];
1408         if (v>=0)
1409           s+= v;
1410         else
1411           s-= v;
1412       }
1413       p1 = p1a;
1414       p1a+= lx;
1415       p2+= lx;
1416     }
1417   }
1418   else /* if (hx && hy) */
1419   {
1420     p1a = p1 + lx;
1421     for (j=0; j<h; j++)
1422     {
1423       for (i=0; i<16; i++)
1424       {
1425         v = ((unsigned int)(p1[i]+p1[i+1]+p1a[i]+p1a[i+1]+2)>>2) - p2[i];
1426         if (v>=0)
1427           s+= v;
1428         else
1429           s-= v;
1430       }
1431       p1 = p1a;
1432       p1a+= lx;
1433       p2+= lx;
1434     }
1435   }
1436
1437   return s;
1438 }
1439
1440 /*
1441  * total squared difference between two (16*h) blocks
1442  * including optional half pel interpolation of blk1 (hx,hy)
1443  * blk1,blk2: addresses of top left pels of both blocks
1444  * lx:        distance (in bytes) of vertically adjacent pels
1445  * hx,hy:     flags for horizontal and/or vertical interpolation
1446  * h:         height of block (usually 8 or 16)
1447  */
1448 static int dist2(blk1,blk2,lx,hx,hy,h)
1449 unsigned char *blk1,*blk2;
1450 int lx,hx,hy,h;
1451 {
1452   unsigned char *p1,*p1a,*p2;
1453   int i,j;
1454   int s,v;
1455
1456   s = 0;
1457   p1 = blk1;
1458   p2 = blk2;
1459   if (!hx && !hy)
1460     for (j=0; j<h; j++)
1461     {
1462       for (i=0; i<16; i++)
1463       {
1464         v = p1[i] - p2[i];
1465         s+= v*v;
1466       }
1467       p1+= lx;
1468       p2+= lx;
1469     }
1470   else if (hx && !hy)
1471     for (j=0; j<h; j++)
1472     {
1473       for (i=0; i<16; i++)
1474       {
1475         v = ((unsigned int)(p1[i]+p1[i+1]+1)>>1) - p2[i];
1476         s+= v*v;
1477       }
1478       p1+= lx;
1479       p2+= lx;
1480     }
1481   else if (!hx && hy)
1482   {
1483     p1a = p1 + lx;
1484     for (j=0; j<h; j++)
1485     {
1486       for (i=0; i<16; i++)
1487       {
1488         v = ((unsigned int)(p1[i]+p1a[i]+1)>>1) - p2[i];
1489         s+= v*v;
1490       }
1491       p1 = p1a;
1492       p1a+= lx;
1493       p2+= lx;
1494     }
1495   }
1496   else /* if (hx && hy) */
1497   {
1498     p1a = p1 + lx;
1499     for (j=0; j<h; j++)
1500     {
1501       for (i=0; i<16; i++)
1502       {
1503         v = ((unsigned int)(p1[i]+p1[i+1]+p1a[i]+p1a[i+1]+2)>>2) - p2[i];
1504         s+= v*v;
1505       }
1506       p1 = p1a;
1507       p1a+= lx;
1508       p2+= lx;
1509     }
1510   }
1511
1512   return s;
1513 }
1514
1515 /*
1516  * absolute difference error between a (16*h) block and a bidirectional
1517  * prediction
1518  *
1519  * p2: address of top left pel of block
1520  * pf,hxf,hyf: address and half pel flags of forward ref. block
1521  * pb,hxb,hyb: address and half pel flags of backward ref. block
1522  * h: height of block
1523  * lx: distance (in bytes) of vertically adjacent pels in p2,pf,pb
1524  */
1525 static int bdist1(pf,pb,p2,lx,hxf,hyf,hxb,hyb,h)
1526 unsigned char *pf,*pb,*p2;
1527 int lx,hxf,hyf,hxb,hyb,h;
1528 {
1529   unsigned char *pfa,*pfb,*pfc,*pba,*pbb,*pbc;
1530   int i,j;
1531   int s,v;
1532
1533   pfa = pf + hxf;
1534   pfb = pf + lx*hyf;
1535   pfc = pfb + hxf;
1536
1537   pba = pb + hxb;
1538   pbb = pb + lx*hyb;
1539   pbc = pbb + hxb;
1540
1541   s = 0;
1542
1543   for (j=0; j<h; j++)
1544   {
1545     for (i=0; i<16; i++)
1546     {
1547       v = ((((unsigned int)(*pf++ + *pfa++ + *pfb++ + *pfc++ + 2)>>2) +
1548             ((unsigned int)(*pb++ + *pba++ + *pbb++ + *pbc++ + 2)>>2) + 1)>>1)
1549            - *p2++;
1550       if (v>=0)
1551         s+= v;
1552       else
1553         s-= v;
1554     }
1555     p2+= lx-16;
1556     pf+= lx-16;
1557     pfa+= lx-16;
1558     pfb+= lx-16;
1559     pfc+= lx-16;
1560     pb+= lx-16;
1561     pba+= lx-16;
1562     pbb+= lx-16;
1563     pbc+= lx-16;
1564   }
1565
1566   return s;
1567 }
1568
1569 /*
1570  * squared error between a (16*h) block and a bidirectional
1571  * prediction
1572  *
1573  * p2: address of top left pel of block
1574  * pf,hxf,hyf: address and half pel flags of forward ref. block
1575  * pb,hxb,hyb: address and half pel flags of backward ref. block
1576  * h: height of block
1577  * lx: distance (in bytes) of vertically adjacent pels in p2,pf,pb
1578  */
1579 static int bdist2(pf,pb,p2,lx,hxf,hyf,hxb,hyb,h)
1580 unsigned char *pf,*pb,*p2;
1581 int lx,hxf,hyf,hxb,hyb,h;
1582 {
1583   unsigned char *pfa,*pfb,*pfc,*pba,*pbb,*pbc;
1584   int i,j;
1585   int s,v;
1586
1587   pfa = pf + hxf;
1588   pfb = pf + lx*hyf;
1589   pfc = pfb + hxf;
1590
1591   pba = pb + hxb;
1592   pbb = pb + lx*hyb;
1593   pbc = pbb + hxb;
1594
1595   s = 0;
1596
1597   for (j=0; j<h; j++)
1598   {
1599     for (i=0; i<16; i++)
1600     {
1601       v = ((((unsigned int)(*pf++ + *pfa++ + *pfb++ + *pfc++ + 2)>>2) +
1602             ((unsigned int)(*pb++ + *pba++ + *pbb++ + *pbc++ + 2)>>2) + 1)>>1)
1603           - *p2++;
1604       s+=v*v;
1605     }
1606     p2+= lx-16;
1607     pf+= lx-16;
1608     pfa+= lx-16;
1609     pfb+= lx-16;
1610     pfc+= lx-16;
1611     pb+= lx-16;
1612     pba+= lx-16;
1613     pbb+= lx-16;
1614     pbc+= lx-16;
1615   }
1616
1617   return s;
1618 }
1619
1620 /*
1621  * variance of a (16*16) block, multiplied by 256
1622  * p:  address of top left pel of block
1623  * lx: distance (in bytes) of vertically adjacent pels
1624  */
1625 static int variance(p,lx)
1626 unsigned char *p;
1627 int lx;
1628 {
1629   int i,j;
1630   unsigned int v,s,s2;
1631
1632   s = s2 = 0;
1633
1634   for (j=0; j<16; j++)
1635   {
1636     for (i=0; i<16; i++)
1637     {
1638       v = *p++;
1639       s+= v;
1640       s2+= v*v;
1641     }
1642     p+= lx-16;
1643   }
1644   return s2 - (s*s)/256;
1645 }