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