RD subpel motion estimation (--subme 7)
[x262.git] / encoder / me.c
1 /*****************************************************************************
2  * me.c: h264 encoder library (Motion Estimation)
3  *****************************************************************************
4  * Copyright (C) 2003 Laurent Aimar
5  * $Id: me.c,v 1.1 2004/06/03 19:27:08 fenrir Exp $
6  *
7  * Authors: Laurent Aimar <fenrir@via.ecp.fr>
8  *          Loren Merritt <lorenm@u.washington.edu>
9  *
10  * This program is free software; you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation; either version 2 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111, USA.
23  *****************************************************************************/
24
25 #include <stdio.h>
26 #include <string.h>
27
28 #include "common/common.h"
29 #include "me.h"
30
31 /* presets selected from good points on the speed-vs-quality curve of several test videos
32  * subpel_iters[i_subpel_refine] = { refine_hpel, refine_qpel, me_hpel, me_qpel }
33  * where me_* are the number of EPZS iterations run on all candidate block types,
34  * and refine_* are run only on the winner. */
35 static const int subpel_iterations[][4] = 
36    {{1,0,0,0},
37     {1,1,0,0},
38     {0,1,1,0},
39     {0,2,1,0},
40     {0,2,1,1},
41     {0,2,1,2},
42     {0,0,2,2},
43     {0,0,2,2}};
44
45 static void refine_subpel( x264_t *h, x264_me_t *m, int hpel_iters, int qpel_iters, int *p_halfpel_thresh, int b_refine_qpel );
46
47 #define COST_MV_INT( mx, my, bd, d ) \
48 { \
49     int cost = h->pixf.sad[i_pixel]( m->p_fenc[0], FENC_STRIDE, \
50                    &p_fref[(my)*m->i_stride[0]+(mx)], m->i_stride[0] ) \
51              + p_cost_mvx[ (mx)<<2 ]  \
52              + p_cost_mvy[ (my)<<2 ]; \
53     if( cost < bcost ) \
54     {                  \
55         bcost = cost;  \
56         bmx = mx;      \
57         bmy = my;      \
58         if( bd ) \
59             dir = d; \
60     } \
61 }
62 #define COST_MV( mx, my )         COST_MV_INT( mx, my, 0, 0 )
63 #define COST_MV_DIR( mx, my, d )  COST_MV_INT( mx, my, 1, d )
64
65 #define COST_MV_PDE( mx, my ) \
66 { \
67     int cost = h->pixf.sad_pde[i_pixel]( m->p_fenc[0], FENC_STRIDE, \
68                    &p_fref[(my)*m->i_stride[0]+(mx)], m->i_stride[0], \
69                    bcost - p_cost_mvx[ (mx)<<2 ] - p_cost_mvy[ (my)<<2 ] ); \
70     if( cost < bcost - p_cost_mvx[ (mx)<<2 ] - p_cost_mvy[ (my)<<2 ] ) \
71     {                  \
72         bcost = cost + p_cost_mvx[ (mx)<<2 ] + p_cost_mvy[ (my)<<2 ];  \
73         bmx = mx;      \
74         bmy = my;      \
75     } \
76 }
77
78 #define DIA1_ITER( mx, my )\
79     {\
80         omx = mx; omy = my;\
81         COST_MV( omx  , omy-1 );/*  1  */\
82         COST_MV( omx  , omy+1 );/* 101 */\
83         COST_MV( omx-1, omy   );/*  1  */\
84         COST_MV( omx+1, omy   );\
85     }
86
87 #define DIA2 \
88     {\
89         COST_MV( omx  , omy-2 );\
90         COST_MV( omx-1, omy-1 );/*   1   */\
91         COST_MV( omx+1, omy-1 );/*  1 1  */\
92         COST_MV( omx-2, omy   );/* 1 0 1 */\
93         COST_MV( omx+2, omy   );/*  1 1  */\
94         COST_MV( omx-1, omy+1 );/*   1   */\
95         COST_MV( omx+1, omy+1 );\
96         COST_MV( omx  , omy+2 );\
97     }\
98
99 #define OCT2 \
100     {\
101         COST_MV( omx-1, omy-2 );\
102         COST_MV( omx+1, omy-2 );/*  1 1  */\
103         COST_MV( omx-2, omy-1 );/* 1   1 */\
104         COST_MV( omx+2, omy-1 );/*   0   */\
105         COST_MV( omx-2, omy+1 );/* 1   1 */\
106         COST_MV( omx+2, omy+1 );/*  1 1  */\
107         COST_MV( omx-1, omy+2 );\
108         COST_MV( omx+1, omy+2 );\
109     }
110
111 #define CROSS( start, x_max, y_max ) \
112     { \
113         for( i = start; i < x_max; i+=2 ) \
114         { \
115             if( omx + i <= mv_x_max ) \
116                 COST_MV( omx + i, omy ); \
117             if( omx - i >= mv_x_min ) \
118                 COST_MV( omx - i, omy ); \
119         } \
120         for( i = start; i < y_max; i+=2 ) \
121         { \
122             if( omy + i <= mv_y_max ) \
123                 COST_MV( omx, omy + i ); \
124             if( omy - i >= mv_y_min ) \
125                 COST_MV( omx, omy - i ); \
126         } \
127     }
128
129
130 void x264_me_search_ref( x264_t *h, x264_me_t *m, int (*mvc)[2], int i_mvc, int *p_halfpel_thresh )
131 {
132     const int i_pixel = m->i_pixel;
133     int i_me_range = h->param.analyse.i_me_range;
134     int bmx, bmy, bcost;
135     int omx, omy, pmx, pmy;
136     uint8_t *p_fref = m->p_fref[0];
137     int i, j;
138     int dir;
139
140     int mv_x_min = h->mb.mv_min_fpel[0];
141     int mv_y_min = h->mb.mv_min_fpel[1];
142     int mv_x_max = h->mb.mv_max_fpel[0];
143     int mv_y_max = h->mb.mv_max_fpel[1];
144
145     const int16_t *p_cost_mvx = m->p_cost_mv - m->mvp[0];
146     const int16_t *p_cost_mvy = m->p_cost_mv - m->mvp[1];
147
148     if( h->mb.i_me_method == X264_ME_UMH )
149     {
150         /* clamp mvp to inside frame+padding, so that we don't have to check it each iteration */
151         p_cost_mvx = m->p_cost_mv - x264_clip3( m->mvp[0], h->mb.mv_min_spel[0], h->mb.mv_max_spel[0] );
152         p_cost_mvy = m->p_cost_mv - x264_clip3( m->mvp[1], h->mb.mv_min_spel[1], h->mb.mv_max_spel[1] );
153     }
154
155     bmx = pmx = x264_clip3( ( m->mvp[0] + 2 ) >> 2, mv_x_min, mv_x_max );
156     bmy = pmy = x264_clip3( ( m->mvp[1] + 2 ) >> 2, mv_y_min, mv_y_max );
157     bcost = COST_MAX;
158     COST_MV( pmx, pmy );
159     /* I don't know why this helps */
160     bcost -= p_cost_mvx[ bmx<<2 ] + p_cost_mvy[ bmy<<2 ];
161
162     /* try extra predictors if provided */
163     for( i = 0; i < i_mvc; i++ )
164     {
165         const int mx = x264_clip3( ( mvc[i][0] + 2 ) >> 2, mv_x_min, mv_x_max );
166         const int my = x264_clip3( ( mvc[i][1] + 2 ) >> 2, mv_y_min, mv_y_max );
167         if( mx != bmx || my != bmy )
168             COST_MV( mx, my );
169     }
170     
171     COST_MV( 0, 0 );
172
173     mv_x_max += 8;
174     mv_y_max += 8;
175     mv_x_min -= 8;
176     mv_y_min -= 8;
177
178     switch( h->mb.i_me_method )
179     {
180     case X264_ME_DIA:
181         /* diamond search, radius 1 */
182         for( i = 0; i < i_me_range; i++ )
183         {
184             DIA1_ITER( bmx, bmy );
185             if( bmx == omx && bmy == omy )
186                 break;
187         }
188         break;
189
190     case X264_ME_HEX:
191 me_hex2:
192         /* hexagon search, radius 2 */
193 #if 0
194         for( i = 0; i < i_me_range/2; i++ )
195         {
196             omx = bmx; omy = bmy;
197             COST_MV( omx-2, omy   );
198             COST_MV( omx-1, omy+2 );
199             COST_MV( omx+1, omy+2 );
200             COST_MV( omx+2, omy   );
201             COST_MV( omx+1, omy-2 );
202             COST_MV( omx-1, omy-2 );
203             if( bmx == omx && bmy == omy )
204                 break;
205         }
206 #else
207         /* equivalent to the above, but eliminates duplicate candidates */
208         dir = -1;
209         omx = bmx; omy = bmy;
210         COST_MV_DIR( omx-2, omy,   0 );
211         COST_MV_DIR( omx-1, omy+2, 1 );
212         COST_MV_DIR( omx+1, omy+2, 2 );
213         COST_MV_DIR( omx+2, omy,   3 );
214         COST_MV_DIR( omx+1, omy-2, 4 );
215         COST_MV_DIR( omx-1, omy-2, 5 );
216         if( dir != -1 )
217         {
218             for( i = 1; i < i_me_range/2; i++ )
219             {
220                 static const int hex2[8][2] = {{-1,-2}, {-2,0}, {-1,2}, {1,2}, {2,0}, {1,-2}, {-1,-2}, {-2,0}};
221                 static const int mod6[8] = {5,0,1,2,3,4,5,0};
222                 const int odir = mod6[dir+1];
223                 omx = bmx; omy = bmy;
224                 COST_MV_DIR( omx + hex2[odir+0][0], omy + hex2[odir+0][1], odir-1 );
225                 COST_MV_DIR( omx + hex2[odir+1][0], omy + hex2[odir+1][1], odir   );
226                 COST_MV_DIR( omx + hex2[odir+2][0], omy + hex2[odir+2][1], odir+1 );
227                 if( bmx == omx && bmy == omy )
228                     break;
229             }
230         }
231 #endif
232         /* square refine */
233         DIA1_ITER( bmx, bmy );
234         COST_MV( omx-1, omy-1 );
235         COST_MV( omx-1, omy+1 );
236         COST_MV( omx+1, omy-1 );
237         COST_MV( omx+1, omy+1 );
238         break;
239
240     case X264_ME_UMH:
241         {
242             /* Uneven-cross Multi-Hexagon-grid Search
243              * as in JM, except with different early termination */
244
245             static const int x264_pixel_size_shift[7] = { 0, 1, 1, 2, 3, 3, 4 };
246
247             int ucost1, ucost2;
248             int cross_start = 1;
249
250             /* refine predictors */
251             ucost1 = bcost;
252             DIA1_ITER( pmx, pmy );
253             if( pmx || pmy )
254                 DIA1_ITER( 0, 0 );
255
256             if(i_pixel == PIXEL_4x4)
257                 goto me_hex2;
258
259             ucost2 = bcost;
260             if( (bmx || bmy) && (bmx!=pmx || bmy!=pmy) )
261                 DIA1_ITER( bmx, bmy );
262             if( bcost == ucost2 )
263                 cross_start = 3;
264             omx = bmx; omy = bmy;
265
266             /* early termination */
267 #define SAD_THRESH(v) ( bcost < ( v >> x264_pixel_size_shift[i_pixel] ) )
268             if( bcost == ucost2 && SAD_THRESH(2000) )
269             {
270                 DIA2;
271                 if( bcost == ucost1 && SAD_THRESH(500) )
272                     break;
273                 if( bcost == ucost2 )
274                 {
275                     int range = (i_me_range>>1) | 1;
276                     CROSS( 3, range, range );
277                     OCT2;
278                     if( bcost == ucost2 )
279                         break;
280                     cross_start = range + 2;
281                 }
282             }
283
284             /* adaptive search range */
285             if( i_mvc ) 
286             {
287                 /* range multipliers based on casual inspection of some statistics of
288                  * average distance between current predictor and final mv found by ESA.
289                  * these have not been tuned much by actual encoding. */
290                 static const int range_mul[4][4] =
291                 {
292                     { 3, 3, 4, 4 },
293                     { 3, 4, 4, 4 },
294                     { 4, 4, 4, 5 },
295                     { 4, 4, 5, 6 },
296                 };
297                 int mvd;
298                 int sad_ctx, mvd_ctx;
299
300                 if( i_mvc == 1 )
301                 {
302                     if( i_pixel == PIXEL_16x16 )
303                         /* mvc is probably the same as mvp, so the difference isn't meaningful.
304                          * but prediction usually isn't too bad, so just use medium range */
305                         mvd = 25;
306                     else
307                         mvd = abs( m->mvp[0] - mvc[0][0] )
308                             + abs( m->mvp[1] - mvc[0][1] );
309                 }
310                 else
311                 {
312                     /* calculate the degree of agreement between predictors. */
313                     /* in 16x16, mvc includes all the neighbors used to make mvp,
314                      * so don't count mvp separately. */
315                     int i_denom = i_mvc - 1;
316                     mvd = 0;
317                     if( i_pixel != PIXEL_16x16 )
318                     {
319                         mvd = abs( m->mvp[0] - mvc[0][0] )
320                             + abs( m->mvp[1] - mvc[0][1] );
321                         i_denom++;
322                     }
323                     for( i = 0; i < i_mvc-1; i++ )
324                         mvd += abs( mvc[i][0] - mvc[i+1][0] )
325                              + abs( mvc[i][1] - mvc[i+1][1] );
326                     mvd /= i_denom; //FIXME idiv
327                 }
328
329                 sad_ctx = SAD_THRESH(1000) ? 0
330                         : SAD_THRESH(2000) ? 1
331                         : SAD_THRESH(4000) ? 2 : 3;
332                 mvd_ctx = mvd < 10 ? 0
333                         : mvd < 20 ? 1
334                         : mvd < 40 ? 2 : 3;
335
336                 i_me_range = i_me_range * range_mul[mvd_ctx][sad_ctx] / 4;
337             }
338
339             /* FIXME if the above DIA2/OCT2/CROSS found a new mv, it has not updated omx/omy.
340              * we are still centered on the same place as the DIA2. is this desirable? */
341             CROSS( cross_start, i_me_range, i_me_range/2 );
342
343             /* 5x5 ESA */
344             omx = bmx; omy = bmy;
345             for( i = (bcost == ucost2) ? 4 : 0; i < 24; i++ )
346             {
347                 static const int square2[24][2] = {
348                     { 1, 0}, { 0, 1}, {-1, 0}, { 0,-1},
349                     { 1, 1}, {-1, 1}, {-1,-1}, { 1,-1},
350                     { 2,-1}, { 2, 0}, { 2, 1}, { 2, 2},
351                     { 1, 2}, { 0, 2}, {-1, 2}, {-2, 2},
352                     {-2, 1}, {-2, 0}, {-2,-1}, {-2,-2},
353                     {-1,-2}, { 0,-2}, { 1,-2}, { 2,-2}
354                 };
355                 COST_MV( omx + square2[i][0], omy + square2[i][1] );
356             }
357
358             /* hexagon grid */
359             omx = bmx; omy = bmy;
360             for( i = 1; i <= i_me_range/4; i++ )
361             {
362                 static const int hex4[16][2] = {
363                     {-4, 2}, {-4, 1}, {-4, 0}, {-4,-1}, {-4,-2},
364                     { 4,-2}, { 4,-1}, { 4, 0}, { 4, 1}, { 4, 2},
365                     { 2, 3}, { 0, 4}, {-2, 3},
366                     {-2,-3}, { 0,-4}, { 2,-3},
367                 };
368                 const int bounds_check = 4*i > X264_MIN4( mv_x_max-omx, mv_y_max-omy, omx-mv_x_min, omy-mv_y_min );
369
370                 if( h->pixf.sad_pde[i_pixel] )
371                 {
372                     for( j = 0; j < 16; j++ )
373                     {
374                         int mx = omx + hex4[j][0]*i;
375                         int my = omy + hex4[j][1]*i;
376                         if( !bounds_check || ( mx >= mv_x_min && mx <= mv_x_max
377                                             && my >= mv_y_min && my <= mv_y_max ) )
378                             COST_MV_PDE( mx, my );
379                     }
380                 }
381                 else
382                 {
383                     for( j = 0; j < 16; j++ )
384                     {
385                         int mx = omx + hex4[j][0]*i;
386                         int my = omy + hex4[j][1]*i;
387                         if( !bounds_check || ( mx >= mv_x_min && mx <= mv_x_max
388                                             && my >= mv_y_min && my <= mv_y_max ) )
389                             COST_MV( mx, my );
390                     }
391                 }
392             }
393             goto me_hex2;
394         }
395
396     case X264_ME_ESA:
397         {
398             const int min_x = X264_MAX( bmx - i_me_range, mv_x_min);
399             const int min_y = X264_MAX( bmy - i_me_range, mv_y_min);
400             const int max_x = X264_MIN( bmx + i_me_range, mv_x_max);
401             const int max_y = X264_MIN( bmy + i_me_range, mv_y_max);
402             int mx, my;
403 #if 0
404             /* plain old exhaustive search */
405             for( my = min_y; my <= max_y; my++ )
406                 for( mx = min_x; mx <= max_x; mx++ )
407                     COST_MV( mx, my );
408 #else
409             /* successive elimination by comparing DC before a full SAD,
410              * because sum(abs(diff)) >= abs(diff(sum)). */
411             const int stride = m->i_stride[0];
412             const int dw = x264_pixel_size[i_pixel].w;
413             const int dh = x264_pixel_size[i_pixel].h * stride;
414             static uint8_t zero[16*16] = {0,};
415             const int enc_dc = h->pixf.sad[i_pixel]( m->p_fenc[0], FENC_STRIDE, zero, 16 );
416             const uint16_t *integral_base = &m->integral[ -1 - 1*stride ];
417
418             if( h->pixf.sad_pde[i_pixel] )
419             {
420                 for( my = min_y; my <= max_y; my++ )
421                     for( mx = min_x; mx <= max_x; mx++ )
422                     {
423                         const uint16_t *integral = &integral_base[ mx + my * stride ];
424                         const uint16_t ref_dc = integral[  0 ] + integral[ dh + dw ]
425                                               - integral[ dw ] - integral[ dh ];
426                         const int bsad = bcost - p_cost_mvx[ (mx)<<2 ] - p_cost_mvy[ (my)<<2 ];
427                         if( abs( ref_dc - enc_dc ) < bsad )
428                             COST_MV_PDE( mx, my );
429                     }
430             }
431             else
432             {
433                 for( my = min_y; my <= max_y; my++ )
434                     for( mx = min_x; mx <= max_x; mx++ )
435                     {
436                         const uint16_t *integral = &integral_base[ mx + my * stride ];
437                         const uint16_t ref_dc = integral[  0 ] + integral[ dh + dw ]
438                                               - integral[ dw ] - integral[ dh ];
439                         const int bsad = bcost - p_cost_mvx[ (mx)<<2 ] - p_cost_mvy[ (my)<<2 ];
440                         if( abs( ref_dc - enc_dc ) < bsad )
441                             COST_MV( mx, my );
442                     }
443             }
444 #endif
445         }
446         break;
447     }
448
449     /* -> qpel mv */
450     m->mv[0] = bmx << 2;
451     m->mv[1] = bmy << 2;
452
453     /* compute the real cost */
454     m->cost_mv = p_cost_mvx[ m->mv[0] ] + p_cost_mvy[ m->mv[1] ];
455     m->cost = bcost;
456     if( bmx == pmx && bmy == pmy )
457         m->cost += m->cost_mv;
458     
459     /* subpel refine */
460     if( h->mb.i_subpel_refine >= 2 )
461     {
462         int hpel = subpel_iterations[h->mb.i_subpel_refine][2];
463         int qpel = subpel_iterations[h->mb.i_subpel_refine][3];
464         refine_subpel( h, m, hpel, qpel, p_halfpel_thresh, 0 );
465     }
466 }
467 #undef COST_MV
468
469 void x264_me_refine_qpel( x264_t *h, x264_me_t *m )
470 {
471     int hpel = subpel_iterations[h->mb.i_subpel_refine][0];
472     int qpel = subpel_iterations[h->mb.i_subpel_refine][1];
473
474     if( m->i_pixel <= PIXEL_8x8 && h->sh.i_type == SLICE_TYPE_P )
475         m->cost -= m->i_ref_cost;
476         
477     refine_subpel( h, m, hpel, qpel, NULL, 1 );
478 }
479
480 #define COST_MV_SAD( mx, my, dir ) \
481 if( b_refine_qpel || (dir^1) != odir ) \
482 { \
483     int stride = 16; \
484     uint8_t *src = h->mc.get_ref( m->p_fref, m->i_stride[0], pix, &stride, mx, my, bw, bh ); \
485     int cost = h->pixf.sad[i_pixel]( m->p_fenc[0], FENC_STRIDE, src, stride ) \
486              + p_cost_mvx[ mx ] + p_cost_mvy[ my ]; \
487     if( cost < bcost ) \
488     {                  \
489         bcost = cost;  \
490         bmx = mx;      \
491         bmy = my;      \
492         bdir = dir;    \
493     } \
494 }
495
496 #define COST_MV_SATD( mx, my, dir ) \
497 if( b_refine_qpel || (dir^1) != odir ) \
498 { \
499     int stride = 16; \
500     uint8_t *src = h->mc.get_ref( m->p_fref, m->i_stride[0], pix, &stride, mx, my, bw, bh ); \
501     int cost = h->pixf.mbcmp[i_pixel]( m->p_fenc[0], FENC_STRIDE, src, stride ) \
502              + p_cost_mvx[ mx ] + p_cost_mvy[ my ]; \
503     if( b_chroma_me && cost < bcost ) \
504     { \
505         h->mc.mc_chroma( m->p_fref[4], m->i_stride[1], pix, 8, mx, my, bw/2, bh/2 ); \
506         cost += h->pixf.mbcmp[i_pixel+3]( m->p_fenc[1], FENC_STRIDE, pix, 8 ); \
507         if( cost < bcost ) \
508         { \
509             h->mc.mc_chroma( m->p_fref[5], m->i_stride[1], pix, 8, mx, my, bw/2, bh/2 ); \
510             cost += h->pixf.mbcmp[i_pixel+3]( m->p_fenc[2], FENC_STRIDE, pix, 8 ); \
511         } \
512     } \
513     if( cost < bcost ) \
514     {                  \
515         bcost = cost;  \
516         bmx = mx;      \
517         bmy = my;      \
518         bdir = dir;    \
519     } \
520 }
521
522 static void refine_subpel( x264_t *h, x264_me_t *m, int hpel_iters, int qpel_iters, int *p_halfpel_thresh, int b_refine_qpel )
523 {
524     const int bw = x264_pixel_size[m->i_pixel].w;
525     const int bh = x264_pixel_size[m->i_pixel].h;
526     const int16_t *p_cost_mvx = m->p_cost_mv - m->mvp[0];
527     const int16_t *p_cost_mvy = m->p_cost_mv - m->mvp[1];
528     const int i_pixel = m->i_pixel;
529     const int b_chroma_me = h->mb.b_chroma_me && i_pixel <= PIXEL_8x8;
530
531     DECLARE_ALIGNED( uint8_t, pix[16*16], 16 );
532     int omx, omy;
533     int i;
534
535     int bmx = m->mv[0];
536     int bmy = m->mv[1];
537     int bcost = m->cost;
538     int odir = -1, bdir;
539
540
541     /* try the subpel component of the predicted mv */
542     if( hpel_iters )
543     {
544         int mx = x264_clip3( m->mvp[0], h->mb.mv_min_spel[0], h->mb.mv_max_spel[0] );
545         int my = x264_clip3( m->mvp[1], h->mb.mv_min_spel[1], h->mb.mv_max_spel[1] );
546         if( mx != bmx || my != bmy )
547             COST_MV_SAD( mx, my, -1 );
548     }
549     
550     /* hpel search */
551     bdir = -1;
552     for( i = hpel_iters; i > 0; i-- )
553     {
554         odir = bdir;
555         omx = bmx;
556         omy = bmy;
557         COST_MV_SAD( omx, omy - 2, 0 );
558         COST_MV_SAD( omx, omy + 2, 1 );
559         COST_MV_SAD( omx - 2, omy, 2 );
560         COST_MV_SAD( omx + 2, omy, 3 );
561         if( bmx == omx && bmy == omy )
562             break;
563     }
564     
565     if( !b_refine_qpel )
566     {
567         bcost = COST_MAX;
568         COST_MV_SATD( bmx, bmy, -1 );
569     }
570     
571     /* early termination when examining multiple reference frames */
572     if( p_halfpel_thresh )
573     {
574         if( (bcost*7)>>3 > *p_halfpel_thresh )
575         {
576             m->cost = bcost;
577             m->mv[0] = bmx;
578             m->mv[1] = bmy;
579             // don't need cost_mv
580             return;
581         }
582         else if( bcost < *p_halfpel_thresh )
583             *p_halfpel_thresh = bcost;
584     }
585
586     /* qpel search */
587     bdir = -1;
588     for( i = qpel_iters; i > 0; i-- )
589     {
590         odir = bdir;
591         omx = bmx;
592         omy = bmy;
593         COST_MV_SATD( omx, omy - 1, 0 );
594         COST_MV_SATD( omx, omy + 1, 1 );
595         COST_MV_SATD( omx - 1, omy, 2 );
596         COST_MV_SATD( omx + 1, omy, 3 );
597         if( bmx == omx && bmy == omy )
598             break;
599     }
600
601     m->cost = bcost;
602     m->mv[0] = bmx;
603     m->mv[1] = bmy;
604     m->cost_mv = p_cost_mvx[ bmx ] + p_cost_mvy[ bmy ];
605 }
606
607 #define BIME_CACHE( dx, dy ) \
608 { \
609     int i = 4 + 3*dx + dy; \
610     h->mc.mc_luma( m0->p_fref, m0->i_stride[0], pix0[i], bw, om0x+dx, om0y+dy, bw, bh ); \
611     h->mc.mc_luma( m1->p_fref, m1->i_stride[0], pix1[i], bw, om1x+dx, om1y+dy, bw, bh ); \
612 }
613
614 #define BIME_CACHE2(a,b) \
615     BIME_CACHE(a,b) \
616     BIME_CACHE(-(a),-(b))
617
618 #define COST_BIMV_SATD( m0x, m0y, m1x, m1y ) \
619 if( pass == 0 || !visited[(m0x)&7][(m0y)&7][(m1x)&7][(m1y)&7] ) \
620 { \
621     int cost; \
622     int i0 = 4 + 3*(m0x-om0x) + (m0y-om0y); \
623     int i1 = 4 + 3*(m1x-om1x) + (m1y-om1y); \
624     visited[(m0x)&7][(m0y)&7][(m1x)&7][(m1y)&7] = 1; \
625     memcpy( pix, pix0[i0], bs ); \
626     if( i_weight == 32 ) \
627         h->mc.avg[i_pixel]( pix, bw, pix1[i1], bw ); \
628     else \
629         h->mc.avg_weight[i_pixel]( pix, bw, pix1[i1], bw, i_weight ); \
630     cost = h->pixf.mbcmp[i_pixel]( m0->p_fenc[0], FENC_STRIDE, pix, bw ) \
631          + p_cost_m0x[ m0x ] + p_cost_m0y[ m0y ] \
632          + p_cost_m1x[ m1x ] + p_cost_m1y[ m1y ]; \
633     if( cost < bcost ) \
634     {                  \
635         bcost = cost;  \
636         bm0x = m0x;    \
637         bm0y = m0y;    \
638         bm1x = m1x;    \
639         bm1y = m1y;    \
640     } \
641 }
642
643 #define CHECK_BIDIR(a,b,c,d) \
644     COST_BIMV_SATD(om0x+a, om0y+b, om1x+c, om1y+d)
645
646 #define CHECK_BIDIR2(a,b,c,d) \
647     CHECK_BIDIR(a,b,c,d) \
648     CHECK_BIDIR(-(a),-(b),-(c),-(d))
649
650 #define CHECK_BIDIR8(a,b,c,d) \
651     CHECK_BIDIR2(a,b,c,d) \
652     CHECK_BIDIR2(b,c,d,a) \
653     CHECK_BIDIR2(c,d,a,b) \
654     CHECK_BIDIR2(d,a,b,c)
655
656 int x264_me_refine_bidir( x264_t *h, x264_me_t *m0, x264_me_t *m1, int i_weight )
657 {
658     const int i_pixel = m0->i_pixel;
659     const int bw = x264_pixel_size[i_pixel].w;
660     const int bh = x264_pixel_size[i_pixel].h;
661     const int bs = bw*bh;
662     const int16_t *p_cost_m0x = m0->p_cost_mv - x264_clip3( m0->mvp[0], h->mb.mv_min_spel[0], h->mb.mv_max_spel[0] );
663     const int16_t *p_cost_m0y = m0->p_cost_mv - x264_clip3( m0->mvp[1], h->mb.mv_min_spel[0], h->mb.mv_max_spel[0] );
664     const int16_t *p_cost_m1x = m1->p_cost_mv - x264_clip3( m1->mvp[0], h->mb.mv_min_spel[0], h->mb.mv_max_spel[0] );
665     const int16_t *p_cost_m1y = m1->p_cost_mv - x264_clip3( m1->mvp[1], h->mb.mv_min_spel[0], h->mb.mv_max_spel[0] );
666     DECLARE_ALIGNED( uint8_t, pix0[9][16*16], 16 );
667     DECLARE_ALIGNED( uint8_t, pix1[9][16*16], 16 );
668     DECLARE_ALIGNED( uint8_t, pix[16*16], 16 );
669     int bm0x = m0->mv[0], om0x = bm0x;
670     int bm0y = m0->mv[1], om0y = bm0y;
671     int bm1x = m1->mv[0], om1x = bm1x;
672     int bm1y = m1->mv[1], om1y = bm1y;
673     int bcost = COST_MAX;
674     int pass = 0;
675     uint8_t visited[8][8][8][8];
676     memset( visited, 0, sizeof(visited) );
677
678     BIME_CACHE( 0, 0 );
679     CHECK_BIDIR( 0, 0, 0, 0 );
680
681     for( pass = 0; pass < 8; pass++ )
682     {
683         /* check all mv pairs that differ in at most 2 components from the current mvs. */
684         /* doesn't do chroma ME. this probably doesn't matter, as the gains
685          * from bidir ME are the same with and without chroma ME. */
686
687         BIME_CACHE2( 1, 0 );
688         BIME_CACHE2( 0, 1 );
689         BIME_CACHE2( 1, 1 );
690         BIME_CACHE2( 1,-1 );
691
692         CHECK_BIDIR8( 0, 0, 0, 1 );
693         CHECK_BIDIR8( 0, 0, 1, 1 );
694         CHECK_BIDIR2( 0, 1, 0, 1 );
695         CHECK_BIDIR2( 1, 0, 1, 0 );
696         CHECK_BIDIR8( 0, 0,-1, 1 );
697         CHECK_BIDIR2( 0,-1, 0, 1 );
698         CHECK_BIDIR2(-1, 0, 1, 0 );
699
700         if( om0x == bm0x && om0y == bm0y && om1x == bm1x && om1y == bm1y )
701             break;
702
703         om0x = bm0x;
704         om0y = bm0y;
705         om1x = bm1x;
706         om1y = bm1y;
707         BIME_CACHE( 0, 0 );
708     }
709
710     m0->mv[0] = bm0x;
711     m0->mv[1] = bm0y;
712     m1->mv[0] = bm1x;
713     m1->mv[1] = bm1y;
714     return bcost;
715 }
716
717 #define COST_MV_RD( mx, my, dir ) \
718 { \
719     if( (dir^1) != odir && (dir<0 || !p_visited[(mx)+(my)*16]) ) \
720     { \
721         int cost; \
722         cache_mv[0] = cache_mv2[0] = mx; \
723         cache_mv[1] = cache_mv2[1] = my; \
724         cost = x264_rd_cost_part( h, i_lambda2, i8, m->i_pixel ); \
725         if( cost < bcost ) \
726         {                  \
727             bcost = cost;  \
728             bmx = mx;      \
729             bmy = my;      \
730         } \
731         if(dir>=0) p_visited[(mx)+(my)*16] = 1; \
732     } \
733 }
734
735 void x264_me_refine_qpel_rd( x264_t *h, x264_me_t *m, int i_lambda2, int i8 )
736 {
737     // don't have to fill the whole mv cache rectangle
738     static const int pixel_mv_offs[] = { 0, 4, 4*8, 0 };
739     int16_t *cache_mv = h->mb.cache.mv[0][x264_scan8[i8*4]];
740     int16_t *cache_mv2 = cache_mv + pixel_mv_offs[m->i_pixel];
741     const int bw = x264_pixel_size[m->i_pixel].w>>2;
742     const int bh = x264_pixel_size[m->i_pixel].h>>2;
743
744     int bcost = m->i_pixel == PIXEL_16x16 ? m->cost : COST_MAX;
745     int bmx = m->mv[0]; 
746     int bmy = m->mv[1];
747     int omx, omy, i;
748     int odir = -1, bdir;
749
750     int visited[16*13] = {0}; // only need 13x13, but 16 is more convenient
751     int *p_visited = &visited[6+6*16];
752
753     if( m->i_pixel != PIXEL_16x16 )
754     {
755         COST_MV_RD( bmx, bmy, -1 );
756         x264_mb_predict_mv( h, 0, i8*4, bw, m->mvp );
757     }
758
759     /* check the predicted mv */
760     if( bmx != m->mvp[0] || bmy != m->mvp[1] )
761         COST_MV_RD( m->mvp[0], m->mvp[1], -1 );
762
763     /* mark mv and mvp as visited */
764     p_visited[0] = 1;
765     p_visited -= bmx + bmy*16;
766     {
767         int mx = bmx ^ m->mv[0] ^ m->mvp[0];
768         int my = bmy ^ m->mv[1] ^ m->mvp[1];
769         if( abs(mx-bmx) < 7 && abs(my-bmy) < 7 )
770             p_visited[mx + my*16] = 1;
771     }
772
773     /* hpel */  
774     bdir = -1;
775     for( i = 0; i < 2; i++ )
776     {
777          omx = bmx;
778          omy = bmy;
779          odir = bdir;
780          COST_MV_RD( omx, omy - 2, 0 );
781          COST_MV_RD( omx, omy + 2, 1 );
782          COST_MV_RD( omx - 2, omy, 2 );
783          COST_MV_RD( omx + 2, omy, 3 );
784          if( bmx == omx && bmy == omy )
785             break;
786     }
787     
788     /* qpel */
789     bdir = -1;
790     for( i = 0; i < 2; i++ )
791     {
792          omx = bmx;
793          omy = bmy;
794          odir = bdir;
795          COST_MV_RD( omx, omy - 1, 0 );
796          COST_MV_RD( omx, omy + 1, 1 );
797          COST_MV_RD( omx - 1, omy, 2 );
798          COST_MV_RD( omx + 1, omy, 3 );
799          if( bmx == omx && bmy == omy )
800             break;
801     }
802
803     m->cost = bcost;
804     m->mv[0] = bmx;
805     m->mv[1] = bmy;
806
807     x264_macroblock_cache_mv ( h, 2*(i8&1), i8&2, bw, bh, 0, bmx, bmy );
808     x264_macroblock_cache_mvd( h, 2*(i8&1), i8&2, bw, bh, 0, bmx - m->mvp[0], bmy - m->mvp[1] );
809 }
810