09fe6626139a5c65469e4177e94e1078ecbc6e3c
[ffmpeg.git] / libavcodec / motion_est.c
1 /*
2  * Motion estimation 
3  * Copyright (c) 2000,2001 Gerard Lantau.
4  * 
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
19  */
20 #include <stdlib.h>
21 #include <stdio.h>
22 #include "avcodec.h"
23 #include "dsputil.h"
24 #include "mpegvideo.h"
25
26 static void halfpel_motion_search(MpegEncContext * s,
27                                   int *mx_ptr, int *my_ptr, int dmin,
28                                   int xmin, int ymin, int xmax, int ymax);
29
30 /* config it to test motion vector encoding (send random vectors) */
31 //#define CONFIG_TEST_MV_ENCODE
32
33 static int pix_sum(UINT8 * pix, int line_size)
34 {
35     int s, i, j;
36
37     s = 0;
38     for (i = 0; i < 16; i++) {
39         for (j = 0; j < 16; j += 8) {
40             s += pix[0];
41             s += pix[1];
42             s += pix[2];
43             s += pix[3];
44             s += pix[4];
45             s += pix[5];
46             s += pix[6];
47             s += pix[7];
48             pix += 8;
49         }
50         pix += line_size - 16;
51     }
52     return s;
53 }
54
55 static int pix_norm1(UINT8 * pix, int line_size)
56 {
57     int s, i, j;
58     UINT32 *sq = squareTbl + 256;
59
60     s = 0;
61     for (i = 0; i < 16; i++) {
62         for (j = 0; j < 16; j += 8) {
63             s += sq[pix[0]];
64             s += sq[pix[1]];
65             s += sq[pix[2]];
66             s += sq[pix[3]];
67             s += sq[pix[4]];
68             s += sq[pix[5]];
69             s += sq[pix[6]];
70             s += sq[pix[7]];
71             pix += 8;
72         }
73         pix += line_size - 16;
74     }
75     return s;
76 }
77
78 static int pix_norm(UINT8 * pix1, UINT8 * pix2, int line_size)
79 {
80     int s, i, j;
81     UINT32 *sq = squareTbl + 256;
82
83     s = 0;
84     for (i = 0; i < 16; i++) {
85         for (j = 0; j < 16; j += 8) {
86             s += sq[pix1[0] - pix2[0]];
87             s += sq[pix1[1] - pix2[1]];
88             s += sq[pix1[2] - pix2[2]];
89             s += sq[pix1[3] - pix2[3]];
90             s += sq[pix1[4] - pix2[4]];
91             s += sq[pix1[5] - pix2[5]];
92             s += sq[pix1[6] - pix2[6]];
93             s += sq[pix1[7] - pix2[7]];
94             pix1 += 8;
95             pix2 += 8;
96         }
97         pix1 += line_size - 16;
98         pix2 += line_size - 16;
99     }
100     return s;
101 }
102
103 static void no_motion_search(MpegEncContext * s,
104                              int *mx_ptr, int *my_ptr)
105 {
106     *mx_ptr = 16 * s->mb_x;
107     *my_ptr = 16 * s->mb_y;
108 }
109
110 static int full_motion_search(MpegEncContext * s,
111                               int *mx_ptr, int *my_ptr, int range,
112                               int xmin, int ymin, int xmax, int ymax)
113 {
114     int x1, y1, x2, y2, xx, yy, x, y;
115     int mx, my, dmin, d;
116     UINT8 *pix;
117
118     xx = 16 * s->mb_x;
119     yy = 16 * s->mb_y;
120     x1 = xx - range + 1;        /* we loose one pixel to avoid boundary pb with half pixel pred */
121     if (x1 < xmin)
122         x1 = xmin;
123     x2 = xx + range - 1;
124     if (x2 > xmax)
125         x2 = xmax;
126     y1 = yy - range + 1;
127     if (y1 < ymin)
128         y1 = ymin;
129     y2 = yy + range - 1;
130     if (y2 > ymax)
131         y2 = ymax;
132     pix = s->new_picture[0] + (yy * s->linesize) + xx;
133     dmin = 0x7fffffff;
134     mx = 0;
135     my = 0;
136     for (y = y1; y <= y2; y++) {
137         for (x = x1; x <= x2; x++) {
138             d = pix_abs16x16(pix, s->last_picture[0] + (y * s->linesize) + x,
139                              s->linesize, 16);
140             if (d < dmin ||
141                 (d == dmin &&
142                  (abs(x - xx) + abs(y - yy)) <
143                  (abs(mx - xx) + abs(my - yy)))) {
144                 dmin = d;
145                 mx = x;
146                 my = y;
147             }
148         }
149     }
150
151     *mx_ptr = mx;
152     *my_ptr = my;
153
154 #if 0
155     if (*mx_ptr < -(2 * range) || *mx_ptr >= (2 * range) ||
156         *my_ptr < -(2 * range) || *my_ptr >= (2 * range)) {
157         fprintf(stderr, "error %d %d\n", *mx_ptr, *my_ptr);
158     }
159 #endif
160     return dmin;
161 }
162
163
164 static int log_motion_search(MpegEncContext * s,
165                              int *mx_ptr, int *my_ptr, int range,
166                              int xmin, int ymin, int xmax, int ymax)
167 {
168     int x1, y1, x2, y2, xx, yy, x, y;
169     int mx, my, dmin, d;
170     UINT8 *pix;
171
172     xx = s->mb_x << 4;
173     yy = s->mb_y << 4;
174
175     /* Left limit */
176     x1 = xx - range;
177     if (x1 < xmin)
178         x1 = xmin;
179
180     /* Right limit */
181     x2 = xx + range;
182     if (x2 > xmax)
183         x2 = xmax;
184
185     /* Upper limit */
186     y1 = yy - range;
187     if (y1 < ymin)
188         y1 = ymin;
189
190     /* Lower limit */
191     y2 = yy + range;
192     if (y2 > ymax)
193         y2 = ymax;
194
195     pix = s->new_picture[0] + (yy * s->linesize) + xx;
196     dmin = 0x7fffffff;
197     mx = 0;
198     my = 0;
199
200     do {
201         for (y = y1; y <= y2; y += range) {
202             for (x = x1; x <= x2; x += range) {
203                 d = pix_abs16x16(pix, s->last_picture[0] + (y * s->linesize) + x, s->linesize, 16);
204                 if (d < dmin || (d == dmin && (abs(x - xx) + abs(y - yy)) < (abs(mx - xx) + abs(my - yy)))) {
205                     dmin = d;
206                     mx = x;
207                     my = y;
208                 }
209             }
210         }
211
212         range = range >> 1;
213
214         x1 = mx - range;
215         if (x1 < xmin)
216             x1 = xmin;
217
218         x2 = mx + range;
219         if (x2 > xmax)
220             x2 = xmax;
221
222         y1 = my - range;
223         if (y1 < ymin)
224             y1 = ymin;
225
226         y2 = my + range;
227         if (y2 > ymax)
228             y2 = ymax;
229
230     } while (range >= 1);
231
232 #ifdef DEBUG
233     fprintf(stderr, "log       - MX: %d\tMY: %d\n", mx, my);
234 #endif
235     *mx_ptr = mx;
236     *my_ptr = my;
237     return dmin;
238 }
239
240 static int phods_motion_search(MpegEncContext * s,
241                                int *mx_ptr, int *my_ptr, int range,
242                                int xmin, int ymin, int xmax, int ymax)
243 {
244     int x1, y1, x2, y2, xx, yy, x, y, lastx, d;
245     int mx, my, dminx, dminy;
246     UINT8 *pix;
247
248     xx = s->mb_x << 4;
249     yy = s->mb_y << 4;
250
251     /* Left limit */
252     x1 = xx - range;
253     if (x1 < xmin)
254         x1 = xmin;
255
256     /* Right limit */
257     x2 = xx + range;
258     if (x2 > xmax)
259         x2 = xmax;
260
261     /* Upper limit */
262     y1 = yy - range;
263     if (y1 < ymin)
264         y1 = ymin;
265
266     /* Lower limit */
267     y2 = yy + range;
268     if (y2 > ymax)
269         y2 = ymax;
270
271     pix = s->new_picture[0] + (yy * s->linesize) + xx;
272     mx = 0;
273     my = 0;
274
275     x = xx;
276     y = yy;
277     do {
278         dminx = 0x7fffffff;
279         dminy = 0x7fffffff;
280
281         lastx = x;
282         for (x = x1; x <= x2; x += range) {
283             d = pix_abs16x16(pix, s->last_picture[0] + (y * s->linesize) + x, s->linesize, 16);
284             if (d < dminx || (d == dminx && (abs(x - xx) + abs(y - yy)) < (abs(mx - xx) + abs(my - yy)))) {
285                 dminx = d;
286                 mx = x;
287             }
288         }
289
290         x = lastx;
291         for (y = y1; y <= y2; y += range) {
292             d = pix_abs16x16(pix, s->last_picture[0] + (y * s->linesize) + x, s->linesize, 16);
293             if (d < dminy || (d == dminy && (abs(x - xx) + abs(y - yy)) < (abs(mx - xx) + abs(my - yy)))) {
294                 dminy = d;
295                 my = y;
296             }
297         }
298
299         range = range >> 1;
300
301         x = mx;
302         y = my;
303         x1 = mx - range;
304         if (x1 < xmin)
305             x1 = xmin;
306
307         x2 = mx + range;
308         if (x2 > xmax)
309             x2 = xmax;
310
311         y1 = my - range;
312         if (y1 < ymin)
313             y1 = ymin;
314
315         y2 = my + range;
316         if (y2 > ymax)
317             y2 = ymax;
318
319     } while (range >= 1);
320
321 #ifdef DEBUG
322     fprintf(stderr, "phods     - MX: %d\tMY: %d\n", mx, my);
323 #endif
324
325     /* half pixel search */
326     *mx_ptr = mx;
327     *my_ptr = my;
328     return dminy;
329 }
330
331 /* The idea would be to make half pel ME after Inter/Intra decision to 
332    save time. */
333 static void halfpel_motion_search(MpegEncContext * s,
334                                   int *mx_ptr, int *my_ptr, int dmin,
335                                   int xmin, int ymin, int xmax, int ymax)
336 {
337     int mx, my, mx1, my1, d, xx, yy, dminh;
338     UINT8 *pix;
339
340     mx = *mx_ptr << 1;
341     my = *my_ptr << 1;
342
343     xx = 16 * s->mb_x;
344     yy = 16 * s->mb_y;
345
346     dminh = dmin;
347
348     /* Half pixel search */
349     mx1 = mx;
350     my1 = my;
351
352     pix = s->new_picture[0] + (yy * s->linesize) + xx;
353
354     if ((mx > (xmin << 1)) && mx < (xmax << 1) && 
355         (my > (ymin << 1)) && my < (ymax << 1)) {
356             int dx, dy, px, py;
357             UINT8 *ptr;
358         for (dy = -1; dy <= 1; dy++) {
359             for (dx = -1; dx <= 1; dx++) {
360                 if (dx != 0 || dy != 0) {
361                     px = mx1 + dx;
362                     py = my1 + dy;
363                     ptr = s->last_picture[0] + ((py >> 1) * s->linesize) + (px >> 1);
364                     switch (((py & 1) << 1) | (px & 1)) {
365                     default:
366                     case 0:
367                         d = pix_abs16x16(pix, ptr, s->linesize, 16);
368                         break;
369                     case 1:
370                         d = pix_abs16x16_x2(pix, ptr, s->linesize, 16);
371                         break;
372                     case 2:
373                         d = pix_abs16x16_y2(pix, ptr, s->linesize, 16);
374                         break;
375                     case 3:
376                         d = pix_abs16x16_xy2(pix, ptr, s->linesize, 16);
377                         break;
378                     }
379                     if (d < dminh) {
380                         dminh = d;
381                         mx = px;
382                         my = py;
383                     }
384                 }
385             }
386         }
387     }
388
389     *mx_ptr = mx - (xx << 1);
390     *my_ptr = my - (yy << 1);
391     //fprintf(stderr,"half  - MX: %d\tMY: %d\n",*mx_ptr ,*my_ptr);
392 }
393
394 #ifndef CONFIG_TEST_MV_ENCODE
395
396 int estimate_motion(MpegEncContext * s,
397                     int mb_x, int mb_y,
398                     int *mx_ptr, int *my_ptr)
399 {
400     UINT8 *pix, *ppix;
401     int sum, varc, vard, mx, my, range, dmin, xx, yy;
402     int xmin, ymin, xmax, ymax;
403
404     range = 8 * (1 << (s->f_code - 1));
405     /* XXX: temporary kludge to avoid overflow for msmpeg4 */
406     if (s->out_format == FMT_H263 && !s->h263_msmpeg4)
407         range = range * 2;
408
409     if (s->unrestricted_mv) {
410         xmin = -16;
411         ymin = -16;
412         xmax = s->width;
413         ymax = s->height;
414     } else {
415         xmin = 0;
416         ymin = 0;
417         xmax = s->width - 16;
418         ymax = s->height - 16;
419     }
420
421     switch(s->full_search) {
422     case ME_ZERO:
423     default:
424         no_motion_search(s, &mx, &my);
425         dmin = 0;
426         break;
427     case ME_FULL:
428         dmin = full_motion_search(s, &mx, &my, range, xmin, ymin, xmax, ymax);
429         break;
430     case ME_LOG:
431         dmin = log_motion_search(s, &mx, &my, range / 2, xmin, ymin, xmax, ymax);
432         break;
433     case ME_PHODS:
434         dmin = phods_motion_search(s, &mx, &my, range / 2, xmin, ymin, xmax, ymax);
435         break;
436     }
437 #ifdef CONFIG_MMX
438     if (mm_flags & MM_MMX)
439         emms();
440 #endif
441
442     /* intra / predictive decision */
443     xx = mb_x * 16;
444     yy = mb_y * 16;
445
446     pix = s->new_picture[0] + (yy * s->linesize) + xx;
447     /* At this point (mx,my) are full-pell and the absolute displacement */
448     ppix = s->last_picture[0] + (my * s->linesize) + mx;
449
450     sum = pix_sum(pix, s->linesize);
451     varc = pix_norm1(pix, s->linesize);
452     vard = pix_norm(pix, ppix, s->linesize);
453
454     vard = vard >> 8;
455     sum = sum >> 8;
456     varc = (varc >> 8) - (sum * sum);
457 #if 0
458     printf("varc=%d (sum=%d) vard=%d mx=%d my=%d\n",
459            varc, sum, vard, mx - xx, my - yy);
460 #endif
461     if (vard <= 64 || vard < varc) {
462         if (s->full_search != ME_ZERO) {
463             halfpel_motion_search(s, &mx, &my, dmin, xmin, ymin, xmax, ymax);
464         } else {
465             mx -= 16 * s->mb_x;
466             my -= 16 * s->mb_y;
467         }
468         *mx_ptr = mx;
469         *my_ptr = my;
470         return 0;
471     } else {
472         *mx_ptr = 0;
473         *my_ptr = 0;
474         return 1;
475     }
476 }
477
478 #else
479
480 /* test version which generates valid random vectors */
481 int estimate_motion(MpegEncContext * s,
482                     int mb_x, int mb_y,
483                     int *mx_ptr, int *my_ptr)
484 {
485     int xx, yy, x1, y1, x2, y2, range;
486
487     if ((random() % 10) >= 5) {
488         range = 8 * (1 << (s->f_code - 1));
489         if (s->out_format == FMT_H263 && !s->h263_msmpeg4)
490             range = range * 2;
491
492         xx = 16 * s->mb_x;
493         yy = 16 * s->mb_y;
494         x1 = xx - range;
495         if (x1 < 0)
496             x1 = 0;
497         x2 = xx + range - 1;
498         if (x2 > (s->width - 16))
499             x2 = s->width - 16;
500         y1 = yy - range;
501         if (y1 < 0)
502             y1 = 0;
503         y2 = yy + range - 1;
504         if (y2 > (s->height - 16))
505             y2 = s->height - 16;
506
507         *mx_ptr = (random() % (2 * (x2 - x1 + 1))) + 2 * (x1 - xx);
508         *my_ptr = (random() % (2 * (y2 - y1 + 1))) + 2 * (y1 - yy);
509         return 0;
510     } else {
511         *mx_ptr = 0;
512         *my_ptr = 0;
513         return 1;
514     }
515 }
516
517 #endif