home *** CD-ROM | disk | FTP | other *** search
/ MPEG Toolkit / MPEG Toolkit.iso / win / mpg2w11b / mpeg2enc / motion.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-01  |  47.2 KB  |  1,648 lines
  1. /* motion.c, motion estimation                                              */
  2.  
  3. /* Copyright (C) 1994, 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(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(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. #ifdef _MPEG_WINDOWS_
  124.       myYield();
  125. #endif
  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,sxf,syf,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.       if ((v = p1[0]  - p2[0])<0)  v = -v; s+= v;
  1366.       if ((v = p1[1]  - p2[1])<0)  v = -v; s+= v;
  1367.       if ((v = p1[2]  - p2[2])<0)  v = -v; s+= v;
  1368.       if ((v = p1[3]  - p2[3])<0)  v = -v; s+= v;
  1369.       if ((v = p1[4]  - p2[4])<0)  v = -v; s+= v;
  1370.       if ((v = p1[5]  - p2[5])<0)  v = -v; s+= v;
  1371.       if ((v = p1[6]  - p2[6])<0)  v = -v; s+= v;
  1372.       if ((v = p1[7]  - p2[7])<0)  v = -v; s+= v;
  1373.       if ((v = p1[8]  - p2[8])<0)  v = -v; s+= v;
  1374.       if ((v = p1[9]  - p2[9])<0)  v = -v; s+= v;
  1375.       if ((v = p1[10] - p2[10])<0) v = -v; s+= v;
  1376.       if ((v = p1[11] - p2[11])<0) v = -v; s+= v;
  1377.       if ((v = p1[12] - p2[12])<0) v = -v; s+= v;
  1378.       if ((v = p1[13] - p2[13])<0) v = -v; s+= v;
  1379.       if ((v = p1[14] - p2[14])<0) v = -v; s+= v;
  1380.       if ((v = p1[15] - p2[15])<0) v = -v; s+= v;
  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.       for (i=0; i<16; i++)
  1392.       {
  1393.         v = ((unsigned int)(p1[i]+p1[i+1]+1)>>1) - p2[i];
  1394.         if (v>=0)
  1395.           s+= v;
  1396.         else
  1397.           s-= v;
  1398.       }
  1399.       p1+= lx;
  1400.       p2+= lx;
  1401.     }
  1402.   else if (!hx && hy)
  1403.   {
  1404.     p1a = p1 + lx;
  1405.     for (j=0; j<h; j++)
  1406.     {
  1407.       for (i=0; i<16; i++)
  1408.       {
  1409.         v = ((unsigned int)(p1[i]+p1a[i]+1)>>1) - p2[i];
  1410.         if (v>=0)
  1411.           s+= v;
  1412.         else
  1413.           s-= v;
  1414.       }
  1415.       p1 = p1a;
  1416.       p1a+= lx;
  1417.       p2+= lx;
  1418.     }
  1419.   }
  1420.   else /* if (hx && hy) */
  1421.   {
  1422.     p1a = p1 + lx;
  1423.     for (j=0; j<h; j++)
  1424.     {
  1425.       for (i=0; i<16; i++)
  1426.       {
  1427.         v = ((unsigned int)(p1[i]+p1[i+1]+p1a[i]+p1a[i+1]+2)>>2) - p2[i];
  1428.         if (v>=0)
  1429.           s+= v;
  1430.         else
  1431.           s-= v;
  1432.       }
  1433.       p1 = p1a;
  1434.       p1a+= lx;
  1435.       p2+= lx;
  1436.     }
  1437.   }
  1438.  
  1439.   return s;
  1440. }
  1441.  
  1442. /*
  1443.  * total squared difference between two (16*h) blocks
  1444.  * including optional half pel interpolation of blk1 (hx,hy)
  1445.  * blk1,blk2: addresses of top left pels of both blocks
  1446.  * lx:        distance (in bytes) of vertically adjacent pels
  1447.  * hx,hy:     flags for horizontal and/or vertical interpolation
  1448.  * h:         height of block (usually 8 or 16)
  1449.  */
  1450. static int dist2(blk1,blk2,lx,hx,hy,h)
  1451. unsigned char *blk1,*blk2;
  1452. int lx,hx,hy,h;
  1453. {
  1454.   unsigned char *p1,*p1a,*p2;
  1455.   int i,j;
  1456.   int s,v;
  1457.  
  1458.   s = 0;
  1459.   p1 = blk1;
  1460.   p2 = blk2;
  1461.   if (!hx && !hy)
  1462.     for (j=0; j<h; j++)
  1463.     {
  1464.       for (i=0; i<16; i++)
  1465.       {
  1466.         v = p1[i] - p2[i];
  1467.         s+= v*v;
  1468.       }
  1469.       p1+= lx;
  1470.       p2+= lx;
  1471.     }
  1472.   else if (hx && !hy)
  1473.     for (j=0; j<h; j++)
  1474.     {
  1475.       for (i=0; i<16; i++)
  1476.       {
  1477.         v = ((unsigned int)(p1[i]+p1[i+1]+1)>>1) - p2[i];
  1478.         s+= v*v;
  1479.       }
  1480.       p1+= lx;
  1481.       p2+= lx;
  1482.     }
  1483.   else if (!hx && hy)
  1484.   {
  1485.     p1a = p1 + lx;
  1486.     for (j=0; j<h; j++)
  1487.     {
  1488.       for (i=0; i<16; i++)
  1489.       {
  1490.         v = ((unsigned int)(p1[i]+p1a[i]+1)>>1) - p2[i];
  1491.         s+= v*v;
  1492.       }
  1493.       p1 = p1a;
  1494.       p1a+= lx;
  1495.       p2+= lx;
  1496.     }
  1497.   }
  1498.   else /* if (hx && hy) */
  1499.   {
  1500.     p1a = p1 + lx;
  1501.     for (j=0; j<h; j++)
  1502.     {
  1503.       for (i=0; i<16; i++)
  1504.       {
  1505.         v = ((unsigned int)(p1[i]+p1[i+1]+p1a[i]+p1a[i+1]+2)>>2) - p2[i];
  1506.         s+= v*v;
  1507.       }
  1508.       p1 = p1a;
  1509.       p1a+= lx;
  1510.       p2+= lx;
  1511.     }
  1512.   }
  1513.  
  1514.   return s;
  1515. }
  1516.  
  1517. /*
  1518.  * absolute difference error between a (16*h) block and a bidirectional
  1519.  * prediction
  1520.  *
  1521.  * p2: address of top left pel of block
  1522.  * pf,hxf,hyf: address and half pel flags of forward ref. block
  1523.  * pb,hxb,hyb: address and half pel flags of backward ref. block
  1524.  * h: height of block
  1525.  * lx: distance (in bytes) of vertically adjacent pels in p2,pf,pb
  1526.  */
  1527. static int bdist1(pf,pb,p2,lx,hxf,hyf,hxb,hyb,h)
  1528. unsigned char *pf,*pb,*p2;
  1529. int lx,hxf,hyf,hxb,hyb,h;
  1530. {
  1531.   unsigned char *pfa,*pfb,*pfc,*pba,*pbb,*pbc;
  1532.   int i,j;
  1533.   int s,v;
  1534.  
  1535.   pfa = pf + hxf;
  1536.   pfb = pf + lx*hyf;
  1537.   pfc = pfb + hxf;
  1538.  
  1539.   pba = pb + hxb;
  1540.   pbb = pb + lx*hyb;
  1541.   pbc = pbb + hxb;
  1542.  
  1543.   s = 0;
  1544.  
  1545.   for (j=0; j<h; j++)
  1546.   {
  1547.     for (i=0; i<16; i++)
  1548.     {
  1549.       v = ((((unsigned int)(*pf++ + *pfa++ + *pfb++ + *pfc++ + 2)>>2) +
  1550.             ((unsigned int)(*pb++ + *pba++ + *pbb++ + *pbc++ + 2)>>2) + 1)>>1)
  1551.            - *p2++;
  1552.       if (v>=0)
  1553.         s+= v;
  1554.       else
  1555.         s-= v;
  1556.     }
  1557.     p2+= lx-16;
  1558.     pf+= lx-16;
  1559.     pfa+= lx-16;
  1560.     pfb+= lx-16;
  1561.     pfc+= lx-16;
  1562.     pb+= lx-16;
  1563.     pba+= lx-16;
  1564.     pbb+= lx-16;
  1565.     pbc+= lx-16;
  1566.   }
  1567.  
  1568.   return s;
  1569. }
  1570.  
  1571. /*
  1572.  * squared error between a (16*h) block and a bidirectional
  1573.  * prediction
  1574.  *
  1575.  * p2: address of top left pel of block
  1576.  * pf,hxf,hyf: address and half pel flags of forward ref. block
  1577.  * pb,hxb,hyb: address and half pel flags of backward ref. block
  1578.  * h: height of block
  1579.  * lx: distance (in bytes) of vertically adjacent pels in p2,pf,pb
  1580.  */
  1581. static int bdist2(pf,pb,p2,lx,hxf,hyf,hxb,hyb,h)
  1582. unsigned char *pf,*pb,*p2;
  1583. int lx,hxf,hyf,hxb,hyb,h;
  1584. {
  1585.   unsigned char *pfa,*pfb,*pfc,*pba,*pbb,*pbc;
  1586.   int i,j;
  1587.   int s,v;
  1588.  
  1589.   pfa = pf + hxf;
  1590.   pfb = pf + lx*hyf;
  1591.   pfc = pfb + hxf;
  1592.  
  1593.   pba = pb + hxb;
  1594.   pbb = pb + lx*hyb;
  1595.   pbc = pbb + hxb;
  1596.  
  1597.   s = 0;
  1598.  
  1599.   for (j=0; j<h; j++)
  1600.   {
  1601.     for (i=0; i<16; i++)
  1602.     {
  1603.       v = ((((unsigned int)(*pf++ + *pfa++ + *pfb++ + *pfc++ + 2)>>2) +
  1604.             ((unsigned int)(*pb++ + *pba++ + *pbb++ + *pbc++ + 2)>>2) + 1)>>1)
  1605.           - *p2++;
  1606.       s+=v*v;
  1607.     }
  1608.     p2+= lx-16;
  1609.     pf+= lx-16;
  1610.     pfa+= lx-16;
  1611.     pfb+= lx-16;
  1612.     pfc+= lx-16;
  1613.     pb+= lx-16;
  1614.     pba+= lx-16;
  1615.     pbb+= lx-16;
  1616.     pbc+= lx-16;
  1617.   }
  1618.  
  1619.   return s;
  1620. }
  1621.  
  1622. /*
  1623.  * variance of a (16*16) block, multiplied by 256
  1624.  * p:  address of top left pel of block
  1625.  * lx: distance (in bytes) of vertically adjacent pels
  1626.  */
  1627. static int variance(p,lx)
  1628. unsigned char *p;
  1629. int lx;
  1630. {
  1631.   int i,j;
  1632.   unsigned int v,s,s2;
  1633.  
  1634.   s = s2 = 0;
  1635.  
  1636.   for (j=0; j<16; j++)
  1637.   {
  1638.     for (i=0; i<16; i++)
  1639.     {
  1640.       v = *p++;
  1641.       s+= v;
  1642.       s2+= v*v;
  1643.     }
  1644.     p+= lx-16;
  1645.   }
  1646.   return s2 - (s*s)/256;
  1647. }
  1648.