cinelerra-5.1/mpeg2enc/motion.c.table

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