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