git.sesse.net Git - ffmpeg/blob - libavcodec/motion_est.c

   1 /*
   2  * Motion estimation
   3  * Copyright (c) 2000,2001 Fabrice Bellard.
   4  * Copyright (c) 2002 Michael Niedermayer
   5  *
   6  *
   7  * This library is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2 of the License, or (at your option) any later version.
  11  *
  12  * This library is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with this library; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  20  *
  21  * new Motion Estimation (X1/EPZS) by Michael Niedermayer <michaelni@gmx.at>
  22  */
  23 #include <stdlib.h>
  24 #include <stdio.h>
  25 #include "avcodec.h"
  26 #include "dsputil.h"
  27 #include "mpegvideo.h"
  28
  29 //#undef NDEBUG
  30 //#include <assert.h>
  31
  32 #define SQ(a) ((a)*(a))
  33
  34 #define P_LEFT P[1]
  35 #define P_TOP P[2]
  36 #define P_TOPRIGHT P[3]
  37 #define P_MEDIAN P[4]
  38 #define P_MV1 P[9]
  39
  40 static inline int sad_hpel_motion_search(MpegEncContext * s,
  41                                   int *mx_ptr, int *my_ptr, int dmin,
  42                                   int xmin, int ymin, int xmax, int ymax,
  43                                   int pred_x, int pred_y, Picture *picture,
  44                                   int n, int size, uint16_t * const mv_penalty);
  45
  46 static inline int update_map_generation(MpegEncContext * s)
  47 {
  48     s->me.map_generation+= 1<<(ME_MAP_MV_BITS*2);
  49     if(s->me.map_generation==0){
  50         s->me.map_generation= 1<<(ME_MAP_MV_BITS*2);
  51         memset(s->me.map, 0, sizeof(uint32_t)*ME_MAP_SIZE);
  52     }
  53     return s->me.map_generation;
  54 }
  55
  56 /* shape adaptive search stuff */
  57 typedef struct Minima{
  58     int height;
  59     int x, y;
  60     int checked;
  61 }Minima;
  62
  63 static int minima_cmp(const void *a, const void *b){
  64     Minima *da = (Minima *) a;
  65     Minima *db = (Minima *) b;
  66
  67     return da->height - db->height;
  68 }
  69
  70 /* SIMPLE */
  71 #define RENAME(a) simple_ ## a
  72
  73 #define CMP(d, x, y, size)\
  74 d = cmp(s, src_y, (ref_y) + (x) + (y)*(stride), stride);
  75
  76 #define CMP_HPEL(d, dx, dy, x, y, size)\
  77 {\
  78     const int dxy= (dx) + 2*(dy);\
  79     hpel_put[0][dxy](s->me.scratchpad, (ref_y) + (x) + (y)*(stride), stride, (16>>size));\
  80     d = cmp_sub(s, s->me.scratchpad, src_y, stride);\
  81 }
  82
  83 #define CMP_QPEL(d, dx, dy, x, y, size)\
  84 {\
  85     const int dxy= (dx) + 4*(dy);\
  86     qpel_put[0][dxy](s->me.scratchpad, (ref_y) + (x) + (y)*(stride), stride);\
  87     d = cmp_sub(s, s->me.scratchpad, src_y, stride);\
  88 }
  89
  90 #include "motion_est_template.c"
  91 #undef RENAME
  92 #undef CMP
  93 #undef CMP_HPEL
  94 #undef CMP_QPEL
  95 #undef INIT
  96
  97 /* SIMPLE CHROMA */
  98 #define RENAME(a) simple_chroma_ ## a
  99
 100 #define CMP(d, x, y, size)\
 101 d = cmp(s, src_y, (ref_y) + (x) + (y)*(stride), stride);\
 102 if(chroma_cmp){\
 103     int dxy= ((x)&1) + 2*((y)&1);\
 104     int c= ((x)>>1) + ((y)>>1)*uvstride;\
 105 \
 106     chroma_hpel_put[0][dxy](s->me.scratchpad, ref_u + c, uvstride, 8);\
 107     d += chroma_cmp(s, s->me.scratchpad, src_u, uvstride);\
 108     chroma_hpel_put[0][dxy](s->me.scratchpad, ref_v + c, uvstride, 8);\
 109     d += chroma_cmp(s, s->me.scratchpad, src_v, uvstride);\
 110 }
 111
 112 #define CMP_HPEL(d, dx, dy, x, y, size)\
 113 {\
 114     const int dxy= (dx) + 2*(dy);\
 115     hpel_put[0][dxy](s->me.scratchpad, (ref_y) + (x) + (y)*(stride), stride, (16>>size));\
 116     d = cmp_sub(s, s->me.scratchpad, src_y, stride);\
 117     if(chroma_cmp_sub){\
 118         int cxy= (dxy) | ((x)&1) | (2*((y)&1));\
 119         int c= ((x)>>1) + ((y)>>1)*uvstride;\
 120         chroma_hpel_put[0][cxy](s->me.scratchpad, ref_u + c, uvstride, 8);\
 121         d += chroma_cmp_sub(s, s->me.scratchpad, src_u, uvstride);\
 122         chroma_hpel_put[0][cxy](s->me.scratchpad, ref_v + c, uvstride, 8);\
 123         d += chroma_cmp_sub(s, s->me.scratchpad, src_v, uvstride);\
 124     }\
 125 }
 126
 127 #define CMP_QPEL(d, dx, dy, x, y, size)\
 128 {\
 129     const int dxy= (dx) + 4*(dy);\
 130     qpel_put[0][dxy](s->me.scratchpad, (ref_y) + (x) + (y)*(stride), stride);\
 131     d = cmp_sub(s, s->me.scratchpad, src_y, stride);\
 132     if(chroma_cmp_sub){\
 133         int cxy, c;\
 134         int cx= (4*(x) + (dx))/2;\
 135         int cy= (4*(y) + (dy))/2;\
 136         cx= (cx>>1)|(cx&1);\
 137         cy= (cy>>1)|(cy&1);\
 138         cxy= (cx&1) + 2*(cy&1);\
 139         c= ((cx)>>1) + ((cy)>>1)*uvstride;\
 140         chroma_hpel_put[0][cxy](s->me.scratchpad, ref_u + c, uvstride, 8);\
 141         d += chroma_cmp_sub(s, s->me.scratchpad, src_u, uvstride);\
 142         chroma_hpel_put[0][cxy](s->me.scratchpad, ref_v + c, uvstride, 8);\
 143         d += chroma_cmp_sub(s, s->me.scratchpad, src_v, uvstride);\
 144     }\
 145 }
 146
 147 #include "motion_est_template.c"
 148 #undef RENAME
 149 #undef CMP
 150 #undef CMP_HPEL
 151 #undef CMP_QPEL
 152 #undef INIT
 153
 154 /* SIMPLE DIRECT HPEL */
 155 #define RENAME(a) simple_direct_hpel_ ## a
 156 //FIXME precalc divisions stuff
 157
 158 #define CMP_DIRECT(d, dx, dy, x, y, size, cmp_func)\
 159 if((x) >= xmin && 2*(x) + (dx) <= 2*xmax && (y) >= ymin && 2*(y) + (dy) <= 2*ymax){\
 160     const int hx= 2*(x) + (dx);\
 161     const int hy= 2*(y) + (dy);\
 162     if(s->mv_type==MV_TYPE_8X8){\
 163         int i;\
 164         for(i=0; i<4; i++){\
 165             int fx = s->me.direct_basis_mv[i][0] + hx;\
 166             int fy = s->me.direct_basis_mv[i][1] + hy;\
 167             int bx = hx ? fx - s->me.co_located_mv[i][0] : s->me.co_located_mv[i][0]*(time_pb - time_pp)/time_pp + (i &1)*16;\
 168             int by = hy ? fy - s->me.co_located_mv[i][1] : s->me.co_located_mv[i][1]*(time_pb - time_pp)/time_pp + (i>>1)*16;\
 169             int fxy= (fx&1) + 2*(fy&1);\
 170             int bxy= (bx&1) + 2*(by&1);\
 171 \
 172             uint8_t *dst= s->me.scratchpad + 8*(i&1) + 8*stride*(i>>1);\
 173             hpel_put[1][fxy](dst, (ref_y ) + (fx>>1) + (fy>>1)*(stride), stride, 8);\
 174             hpel_avg[1][bxy](dst, (ref2_y) + (bx>>1) + (by>>1)*(stride), stride, 8);\
 175         }\
 176     }else{\
 177         int fx = s->me.direct_basis_mv[0][0] + hx;\
 178         int fy = s->me.direct_basis_mv[0][1] + hy;\
 179         int bx = hx ? fx - s->me.co_located_mv[0][0] : s->me.co_located_mv[0][0]*(time_pb - time_pp)/time_pp;\
 180         int by = hy ? fy - s->me.co_located_mv[0][1] : s->me.co_located_mv[0][1]*(time_pb - time_pp)/time_pp;\
 181         int fxy= (fx&1) + 2*(fy&1);\
 182         int bxy= (bx&1) + 2*(by&1);\
 183 \
 184         hpel_put[0][fxy](s->me.scratchpad, (ref_y ) + (fx>>1) + (fy>>1)*(stride), stride, 16);\
 185         hpel_avg[0][bxy](s->me.scratchpad, (ref2_y) + (bx>>1) + (by>>1)*(stride), stride, 16);\
 186     }\
 187     d = cmp_func(s, s->me.scratchpad, src_y, stride);\
 188 }else\
 189     d= 256*256*256*32;
 190
 191
 192 #define CMP_HPEL(d, dx, dy, x, y, size)\
 193     CMP_DIRECT(d, dx, dy, x, y, size, cmp_sub)
 194
 195 #define CMP(d, x, y, size)\
 196     CMP_DIRECT(d, 0, 0, x, y, size, cmp)
 197
 198 #include "motion_est_template.c"
 199 #undef RENAME
 200 #undef CMP
 201 #undef CMP_HPEL
 202 #undef CMP_QPEL
 203 #undef INIT
 204 #undef CMP_DIRECT
 205
 206 /* SIMPLE DIRECT QPEL */
 207 #define RENAME(a) simple_direct_qpel_ ## a
 208
 209 #define CMP_DIRECT(d, dx, dy, x, y, size, cmp_func)\
 210 if((x) >= xmin && 4*(x) + (dx) <= 4*xmax && (y) >= ymin && 4*(y) + (dy) <= 4*ymax){\
 211     const int qx= 4*(x) + (dx);\
 212     const int qy= 4*(y) + (dy);\
 213     if(s->mv_type==MV_TYPE_8X8){\
 214         int i;\
 215         for(i=0; i<4; i++){\
 216             int fx = s->me.direct_basis_mv[i][0] + qx;\
 217             int fy = s->me.direct_basis_mv[i][1] + qy;\
 218             int bx = qx ? fx - s->me.co_located_mv[i][0] : s->me.co_located_mv[i][0]*(time_pb - time_pp)/time_pp + (i &1)*16;\
 219             int by = qy ? fy - s->me.co_located_mv[i][1] : s->me.co_located_mv[i][1]*(time_pb - time_pp)/time_pp + (i>>1)*16;\
 220             int fxy= (fx&3) + 4*(fy&3);\
 221             int bxy= (bx&3) + 4*(by&3);\
 222 \
 223             uint8_t *dst= s->me.scratchpad + 8*(i&1) + 8*stride*(i>>1);\
 224             qpel_put[1][fxy](dst, (ref_y ) + (fx>>2) + (fy>>2)*(stride), stride);\
 225             qpel_avg[1][bxy](dst, (ref2_y) + (bx>>2) + (by>>2)*(stride), stride);\
 226         }\
 227     }else{\
 228         int fx = s->me.direct_basis_mv[0][0] + qx;\
 229         int fy = s->me.direct_basis_mv[0][1] + qy;\
 230         int bx = qx ? fx - s->me.co_located_mv[0][0] : s->me.co_located_mv[0][0]*(time_pb - time_pp)/time_pp;\
 231         int by = qy ? fy - s->me.co_located_mv[0][1] : s->me.co_located_mv[0][1]*(time_pb - time_pp)/time_pp;\
 232         int fxy= (fx&3) + 4*(fy&3);\
 233         int bxy= (bx&3) + 4*(by&3);\
 234 \
 235         qpel_put[0][fxy](s->me.scratchpad, (ref_y ) + (fx>>2) + (fy>>2)*(stride), stride);\
 236         qpel_avg[0][bxy](s->me.scratchpad, (ref2_y) + (bx>>2) + (by>>2)*(stride), stride);\
 237     }\
 238     d = cmp_func(s, s->me.scratchpad, src_y, stride);\
 239 }else\
 240     d= 256*256*256*32;
 241
 242
 243 #define CMP_QPEL(d, dx, dy, x, y, size)\
 244     CMP_DIRECT(d, dx, dy, x, y, size, cmp_sub)
 245
 246 #define CMP(d, x, y, size)\
 247     CMP_DIRECT(d, 0, 0, x, y, size, cmp)
 248
 249 #include "motion_est_template.c"
 250 #undef RENAME
 251 #undef CMP
 252 #undef CMP_HPEL
 253 #undef CMP_QPEL
 254 #undef INIT
 255 #undef CMP__DIRECT
 256
 257
 258 static int zero_cmp(void *s, uint8_t *a, uint8_t *b, int stride){
 259     return 0;
 260 }
 261
 262 static void set_cmp(MpegEncContext *s, me_cmp_func *cmp, int type){
 263     DSPContext* c= &s->dsp;
 264     int i;
 265
 266     memset(cmp, 0, sizeof(void*)*11);
 267
 268     switch(type&0xFF){
 269     case FF_CMP_SAD:
 270         cmp[0]= c->sad[0];
 271         cmp[1]= c->sad[1];
 272         break;
 273     case FF_CMP_SATD:
 274         cmp[0]= c->hadamard8_diff[0];
 275         cmp[1]= c->hadamard8_diff[1];
 276         break;
 277     case FF_CMP_SSE:
 278         cmp[0]= c->sse[0];
 279         cmp[1]= c->sse[1];
 280         break;
 281     case FF_CMP_DCT:
 282         cmp[0]= c->dct_sad[0];
 283         cmp[1]= c->dct_sad[1];
 284         break;
 285     case FF_CMP_PSNR:
 286         cmp[0]= c->quant_psnr[0];
 287         cmp[1]= c->quant_psnr[1];
 288         break;
 289     case FF_CMP_ZERO:
 290         for(i=0; i<7; i++){
 291             cmp[i]= zero_cmp;
 292         }
 293         break;
 294     default:
 295         fprintf(stderr,"internal error in cmp function selection\n");
 296     }
 297 };
 298
 299 static inline int get_penalty_factor(MpegEncContext *s, int type){
 300
 301     switch(type){
 302     default:
 303     case FF_CMP_SAD:
 304         return s->qscale;
 305     case FF_CMP_SSE:
 306 //        return s->qscale*8;
 307     case FF_CMP_DCT:
 308     case FF_CMP_SATD:
 309         return s->qscale*8;
 310     }
 311 }
 312
 313 void ff_init_me(MpegEncContext *s){
 314     set_cmp(s, s->dsp.me_pre_cmp, s->avctx->me_pre_cmp);
 315     set_cmp(s, s->dsp.me_cmp, s->avctx->me_cmp);
 316     set_cmp(s, s->dsp.me_sub_cmp, s->avctx->me_sub_cmp);
 317     set_cmp(s, s->dsp.mb_cmp, s->avctx->mb_cmp);
 318
 319     if(s->flags&CODEC_FLAG_QPEL){
 320         if(s->avctx->me_sub_cmp&FF_CMP_CHROMA)
 321             s->me.sub_motion_search= simple_chroma_qpel_motion_search;
 322         else
 323             s->me.sub_motion_search= simple_qpel_motion_search;
 324     }else{
 325         if(s->avctx->me_sub_cmp&FF_CMP_CHROMA)
 326             s->me.sub_motion_search= simple_chroma_hpel_motion_search;
 327         else if(s->avctx->me_sub_cmp == FF_CMP_SAD && s->avctx->me_cmp == FF_CMP_SAD)
 328             s->me.sub_motion_search= sad_hpel_motion_search;
 329         else
 330             s->me.sub_motion_search= simple_hpel_motion_search;
 331     }
 332
 333     if(s->avctx->me_cmp&FF_CMP_CHROMA){
 334         s->me.motion_search[0]= simple_chroma_epzs_motion_search;
 335         s->me.motion_search[1]= simple_chroma_epzs_motion_search4;
 336     }else{
 337         s->me.motion_search[0]= simple_epzs_motion_search;
 338         s->me.motion_search[1]= simple_epzs_motion_search4;
 339     }
 340
 341     if(s->avctx->me_pre_cmp&FF_CMP_CHROMA){
 342         s->me.pre_motion_search= simple_chroma_epzs_motion_search;
 343     }else{
 344         s->me.pre_motion_search= simple_epzs_motion_search;
 345     }
 346 }
 347
 348 static int pix_dev(UINT8 * pix, int line_size, int mean)
 349 {
 350     int s, i, j;
 351
 352     s = 0;
 353     for (i = 0; i < 16; i++) {
 354         for (j = 0; j < 16; j += 8) {
 355             s += ABS(pix[0]-mean);
 356             s += ABS(pix[1]-mean);
 357             s += ABS(pix[2]-mean);
 358             s += ABS(pix[3]-mean);
 359             s += ABS(pix[4]-mean);
 360             s += ABS(pix[5]-mean);
 361             s += ABS(pix[6]-mean);
 362             s += ABS(pix[7]-mean);
 363             pix += 8;
 364         }
 365         pix += line_size - 16;
 366     }
 367     return s;
 368 }
 369
 370 static inline void no_motion_search(MpegEncContext * s,
 371                                     int *mx_ptr, int *my_ptr)
 372 {
 373     *mx_ptr = 16 * s->mb_x;
 374     *my_ptr = 16 * s->mb_y;
 375 }
 376
 377 static int full_motion_search(MpegEncContext * s,
 378                               int *mx_ptr, int *my_ptr, int range,
 379                               int xmin, int ymin, int xmax, int ymax, uint8_t *ref_picture)
 380 {
 381     int x1, y1, x2, y2, xx, yy, x, y;
 382     int mx, my, dmin, d;
 383     UINT8 *pix;
 384
 385     xx = 16 * s->mb_x;
 386     yy = 16 * s->mb_y;
 387     x1 = xx - range + 1;        /* we loose one pixel to avoid boundary pb with half pixel pred */
 388     if (x1 < xmin)
 389         x1 = xmin;
 390     x2 = xx + range - 1;
 391     if (x2 > xmax)
 392         x2 = xmax;
 393     y1 = yy - range + 1;
 394     if (y1 < ymin)
 395         y1 = ymin;
 396     y2 = yy + range - 1;
 397     if (y2 > ymax)
 398         y2 = ymax;
 399     pix = s->new_picture.data[0] + (yy * s->linesize) + xx;
 400     dmin = 0x7fffffff;
 401     mx = 0;
 402     my = 0;
 403     for (y = y1; y <= y2; y++) {
 404         for (x = x1; x <= x2; x++) {
 405             d = s->dsp.pix_abs16x16(pix, ref_picture + (y * s->linesize) + x,
 406                              s->linesize);
 407             if (d < dmin ||
 408                 (d == dmin &&
 409                  (abs(x - xx) + abs(y - yy)) <
 410                  (abs(mx - xx) + abs(my - yy)))) {
 411                 dmin = d;
 412                 mx = x;
 413                 my = y;
 414             }
 415         }
 416     }
 417
 418     *mx_ptr = mx;
 419     *my_ptr = my;
 420
 421 #if 0
 422     if (*mx_ptr < -(2 * range) || *mx_ptr >= (2 * range) ||
 423         *my_ptr < -(2 * range) || *my_ptr >= (2 * range)) {
 424         fprintf(stderr, "error %d %d\n", *mx_ptr, *my_ptr);
 425     }
 426 #endif
 427     return dmin;
 428 }
 429
 430
 431 static int log_motion_search(MpegEncContext * s,
 432                              int *mx_ptr, int *my_ptr, int range,
 433                              int xmin, int ymin, int xmax, int ymax, uint8_t *ref_picture)
 434 {
 435     int x1, y1, x2, y2, xx, yy, x, y;
 436     int mx, my, dmin, d;
 437     UINT8 *pix;
 438
 439     xx = s->mb_x << 4;
 440     yy = s->mb_y << 4;
 441
 442     /* Left limit */
 443     x1 = xx - range;
 444     if (x1 < xmin)
 445         x1 = xmin;
 446
 447     /* Right limit */
 448     x2 = xx + range;
 449     if (x2 > xmax)
 450         x2 = xmax;
 451
 452     /* Upper limit */
 453     y1 = yy - range;
 454     if (y1 < ymin)
 455         y1 = ymin;
 456
 457     /* Lower limit */
 458     y2 = yy + range;
 459     if (y2 > ymax)
 460         y2 = ymax;
 461
 462     pix = s->new_picture.data[0] + (yy * s->linesize) + xx;
 463     dmin = 0x7fffffff;
 464     mx = 0;
 465     my = 0;
 466
 467     do {
 468         for (y = y1; y <= y2; y += range) {
 469             for (x = x1; x <= x2; x += range) {
 470                 d = s->dsp.pix_abs16x16(pix, ref_picture + (y * s->linesize) + x, s->linesize);
 471                 if (d < dmin || (d == dmin && (abs(x - xx) + abs(y - yy)) < (abs(mx - xx) + abs(my - yy)))) {
 472                     dmin = d;
 473                     mx = x;
 474                     my = y;
 475                 }
 476             }
 477         }
 478
 479         range = range >> 1;
 480
 481         x1 = mx - range;
 482         if (x1 < xmin)
 483             x1 = xmin;
 484
 485         x2 = mx + range;
 486         if (x2 > xmax)
 487             x2 = xmax;
 488
 489         y1 = my - range;
 490         if (y1 < ymin)
 491             y1 = ymin;
 492
 493         y2 = my + range;
 494         if (y2 > ymax)
 495             y2 = ymax;
 496
 497     } while (range >= 1);
 498
 499 #ifdef DEBUG
 500     fprintf(stderr, "log       - MX: %d\tMY: %d\n", mx, my);
 501 #endif
 502     *mx_ptr = mx;
 503     *my_ptr = my;
 504     return dmin;
 505 }
 506
 507 static int phods_motion_search(MpegEncContext * s,
 508                                int *mx_ptr, int *my_ptr, int range,
 509                                int xmin, int ymin, int xmax, int ymax, uint8_t *ref_picture)
 510 {
 511     int x1, y1, x2, y2, xx, yy, x, y, lastx, d;
 512     int mx, my, dminx, dminy;
 513     UINT8 *pix;
 514
 515     xx = s->mb_x << 4;
 516     yy = s->mb_y << 4;
 517
 518     /* Left limit */
 519     x1 = xx - range;
 520     if (x1 < xmin)
 521         x1 = xmin;
 522
 523     /* Right limit */
 524     x2 = xx + range;
 525     if (x2 > xmax)
 526         x2 = xmax;
 527
 528     /* Upper limit */
 529     y1 = yy - range;
 530     if (y1 < ymin)
 531         y1 = ymin;
 532
 533     /* Lower limit */
 534     y2 = yy + range;
 535     if (y2 > ymax)
 536         y2 = ymax;
 537
 538     pix = s->new_picture.data[0] + (yy * s->linesize) + xx;
 539     mx = 0;
 540     my = 0;
 541
 542     x = xx;
 543     y = yy;
 544     do {
 545         dminx = 0x7fffffff;
 546         dminy = 0x7fffffff;
 547
 548         lastx = x;
 549         for (x = x1; x <= x2; x += range) {
 550             d = s->dsp.pix_abs16x16(pix, ref_picture + (y * s->linesize) + x, s->linesize);
 551             if (d < dminx || (d == dminx && (abs(x - xx) + abs(y - yy)) < (abs(mx - xx) + abs(my - yy)))) {
 552                 dminx = d;
 553                 mx = x;
 554             }
 555         }
 556
 557         x = lastx;
 558         for (y = y1; y <= y2; y += range) {
 559             d = s->dsp.pix_abs16x16(pix, ref_picture + (y * s->linesize) + x, s->linesize);
 560             if (d < dminy || (d == dminy && (abs(x - xx) + abs(y - yy)) < (abs(mx - xx) + abs(my - yy)))) {
 561                 dminy = d;
 562                 my = y;
 563             }
 564         }
 565
 566         range = range >> 1;
 567
 568         x = mx;
 569         y = my;
 570         x1 = mx - range;
 571         if (x1 < xmin)
 572             x1 = xmin;
 573
 574         x2 = mx + range;
 575         if (x2 > xmax)
 576             x2 = xmax;
 577
 578         y1 = my - range;
 579         if (y1 < ymin)
 580             y1 = ymin;
 581
 582         y2 = my + range;
 583         if (y2 > ymax)
 584             y2 = ymax;
 585
 586     } while (range >= 1);
 587
 588 #ifdef DEBUG
 589     fprintf(stderr, "phods     - MX: %d\tMY: %d\n", mx, my);
 590 #endif
 591
 592     /* half pixel search */
 593     *mx_ptr = mx;
 594     *my_ptr = my;
 595     return dminy;
 596 }
 597
 598
 599 #define Z_THRESHOLD 256
 600
 601 #define CHECK_SAD_HALF_MV(suffix, x, y) \
 602 {\
 603     d= pix_abs_ ## suffix(pix, ptr+((x)>>1), s->linesize);\
 604     d += (mv_penalty[pen_x + x] + mv_penalty[pen_y + y])*penalty_factor;\
 605     COPY3_IF_LT(dminh, d, dx, x, dy, y)\
 606 }
 607
 608 static inline int sad_hpel_motion_search(MpegEncContext * s,
 609                                   int *mx_ptr, int *my_ptr, int dmin,
 610                                   int xmin, int ymin, int xmax, int ymax,
 611                                   int pred_x, int pred_y, Picture *picture,
 612                                   int n, int size, uint16_t * const mv_penalty)
 613 {
 614     uint8_t *ref_picture= picture->data[0];
 615     uint32_t *score_map= s->me.score_map;
 616     const int penalty_factor= s->me.sub_penalty_factor;
 617     int mx, my, xx, yy, dminh;
 618     UINT8 *pix, *ptr;
 619     op_pixels_abs_func pix_abs_x2;
 620     op_pixels_abs_func pix_abs_y2;
 621     op_pixels_abs_func pix_abs_xy2;
 622
 623     if(size==0){
 624         pix_abs_x2 = s->dsp.pix_abs16x16_x2;
 625         pix_abs_y2 = s->dsp.pix_abs16x16_y2;
 626         pix_abs_xy2= s->dsp.pix_abs16x16_xy2;
 627     }else{
 628         pix_abs_x2 = s->dsp.pix_abs8x8_x2;
 629         pix_abs_y2 = s->dsp.pix_abs8x8_y2;
 630         pix_abs_xy2= s->dsp.pix_abs8x8_xy2;
 631     }
 632
 633     if(s->me.skip){
 634 //    printf("S");
 635         *mx_ptr = 0;
 636         *my_ptr = 0;
 637         return dmin;
 638     }
 639 //    printf("N");
 640
 641     xx = 16 * s->mb_x + 8*(n&1);
 642     yy = 16 * s->mb_y + 8*(n>>1);
 643     pix =  s->new_picture.data[0] + (yy * s->linesize) + xx;
 644
 645     mx = *mx_ptr;
 646     my = *my_ptr;
 647     ptr = ref_picture + ((yy + my) * s->linesize) + (xx + mx);
 648
 649     dminh = dmin;
 650
 651     if (mx > xmin && mx < xmax &&
 652         my > ymin && my < ymax) {
 653         int dx=0, dy=0;
 654         int d, pen_x, pen_y;
 655         const int index= (my<<ME_MAP_SHIFT) + mx;
 656         const int t= score_map[(index-(1<<ME_MAP_SHIFT))&(ME_MAP_SIZE-1)];
 657         const int l= score_map[(index- 1               )&(ME_MAP_SIZE-1)];
 658         const int r= score_map[(index+ 1               )&(ME_MAP_SIZE-1)];
 659         const int b= score_map[(index+(1<<ME_MAP_SHIFT))&(ME_MAP_SIZE-1)];
 660         mx<<=1;
 661         my<<=1;
 662
 663
 664         pen_x= pred_x + mx;
 665         pen_y= pred_y + my;
 666
 667         ptr-= s->linesize;
 668         if(t<=b){
 669             CHECK_SAD_HALF_MV(y2 , 0, -1)
 670             if(l<=r){
 671                 CHECK_SAD_HALF_MV(xy2, -1, -1)
 672                 if(t+r<=b+l){
 673                     CHECK_SAD_HALF_MV(xy2, +1, -1)
 674                     ptr+= s->linesize;
 675                 }else{
 676                     ptr+= s->linesize;
 677                     CHECK_SAD_HALF_MV(xy2, -1, +1)
 678                 }
 679                 CHECK_SAD_HALF_MV(x2 , -1,  0)
 680             }else{
 681                 CHECK_SAD_HALF_MV(xy2, +1, -1)
 682                 if(t+l<=b+r){
 683                     CHECK_SAD_HALF_MV(xy2, -1, -1)
 684                     ptr+= s->linesize;
 685                 }else{
 686                     ptr+= s->linesize;
 687                     CHECK_SAD_HALF_MV(xy2, +1, +1)
 688                 }
 689                 CHECK_SAD_HALF_MV(x2 , +1,  0)
 690             }
 691         }else{
 692             if(l<=r){
 693                 if(t+l<=b+r){
 694                     CHECK_SAD_HALF_MV(xy2, -1, -1)
 695                     ptr+= s->linesize;
 696                 }else{
 697                     ptr+= s->linesize;
 698                     CHECK_SAD_HALF_MV(xy2, +1, +1)
 699                 }
 700                 CHECK_SAD_HALF_MV(x2 , -1,  0)
 701                 CHECK_SAD_HALF_MV(xy2, -1, +1)
 702             }else{
 703                 if(t+r<=b+l){
 704                     CHECK_SAD_HALF_MV(xy2, +1, -1)
 705                     ptr+= s->linesize;
 706                 }else{
 707                     ptr+= s->linesize;
 708                     CHECK_SAD_HALF_MV(xy2, -1, +1)
 709                 }
 710                 CHECK_SAD_HALF_MV(x2 , +1,  0)
 711                 CHECK_SAD_HALF_MV(xy2, +1, +1)
 712             }
 713             CHECK_SAD_HALF_MV(y2 ,  0, +1)
 714         }
 715         mx+=dx;
 716         my+=dy;
 717
 718     }else{
 719         mx<<=1;
 720         my<<=1;
 721     }
 722
 723     *mx_ptr = mx;
 724     *my_ptr = my;
 725     return dminh;
 726 }
 727
 728 static inline void set_p_mv_tables(MpegEncContext * s, int mx, int my, int mv4)
 729 {
 730     const int xy= s->mb_x + 1 + (s->mb_y + 1)*(s->mb_width + 2);
 731
 732     s->p_mv_table[xy][0] = mx;
 733     s->p_mv_table[xy][1] = my;
 734
 735     /* has allready been set to the 4 MV if 4MV is done */
 736     if(mv4){
 737         int mot_xy= s->block_index[0];
 738
 739         s->motion_val[mot_xy  ][0]= mx;
 740         s->motion_val[mot_xy  ][1]= my;
 741         s->motion_val[mot_xy+1][0]= mx;
 742         s->motion_val[mot_xy+1][1]= my;
 743
 744         mot_xy += s->block_wrap[0];
 745         s->motion_val[mot_xy  ][0]= mx;
 746         s->motion_val[mot_xy  ][1]= my;
 747         s->motion_val[mot_xy+1][0]= mx;
 748         s->motion_val[mot_xy+1][1]= my;
 749     }
 750 }
 751
 752 static inline void get_limits(MpegEncContext *s, int *range, int *xmin, int *ymin, int *xmax, int *ymax, int f_code)
 753 {
 754     *range = 8 * (1 << (f_code - 1));
 755     /* XXX: temporary kludge to avoid overflow for msmpeg4 */
 756     if (s->out_format == FMT_H263 && !s->h263_msmpeg4)
 757         *range *= 2;
 758
 759     if (s->unrestricted_mv) {
 760         *xmin = -16;
 761         *ymin = -16;
 762         if (s->h263_plus)
 763             *range *= 2;
 764         if(s->avctx->codec->id!=CODEC_ID_MPEG4){
 765             *xmax = s->mb_width*16;
 766             *ymax = s->mb_height*16;
 767         }else {
 768             *xmax = s->width;
 769             *ymax = s->height;
 770         }
 771     } else {
 772         *xmin = 0;
 773         *ymin = 0;
 774         *xmax = s->mb_width*16 - 16;
 775         *ymax = s->mb_height*16 - 16;
 776     }
 777 }
 778
 779 static inline int mv4_search(MpegEncContext *s, int xmin, int ymin, int xmax, int ymax, int mx, int my, int shift)
 780 {
 781     int block;
 782     int P[10][2];
 783     uint8_t *ref_picture= s->last_picture.data[0];
 784     int dmin_sum=0;
 785     uint16_t * const mv_penalty= s->me.mv_penalty[s->f_code] + MAX_MV;
 786
 787     for(block=0; block<4; block++){
 788         int mx4, my4;
 789         int pred_x4, pred_y4;
 790         int dmin4;
 791         static const int off[4]= {2, 1, 1, -1};
 792         const int mot_stride = s->block_wrap[0];
 793         const int mot_xy = s->block_index[block];
 794 //        const int block_x= (block&1);
 795 //        const int block_y= (block>>1);
 796 #if 1 // this saves us a bit of cliping work and shouldnt affect compression in a negative way
 797         const int rel_xmin4= xmin;
 798         const int rel_xmax4= xmax;
 799         const int rel_ymin4= ymin;
 800         const int rel_ymax4= ymax;
 801 #else
 802         const int rel_xmin4= xmin - block_x*8;
 803         const int rel_xmax4= xmax - block_x*8 + 8;
 804         const int rel_ymin4= ymin - block_y*8;
 805         const int rel_ymax4= ymax - block_y*8 + 8;
 806 #endif
 807         P_LEFT[0] = s->motion_val[mot_xy - 1][0];
 808         P_LEFT[1] = s->motion_val[mot_xy - 1][1];
 809
 810         if(P_LEFT[0]       > (rel_xmax4<<shift)) P_LEFT[0]       = (rel_xmax4<<shift);
 811
 812         /* special case for first line */
 813         if (s->mb_y == 0 && block<2) {
 814             pred_x4= P_LEFT[0];
 815             pred_y4= P_LEFT[1];
 816         } else {
 817             P_TOP[0]      = s->motion_val[mot_xy - mot_stride             ][0];
 818             P_TOP[1]      = s->motion_val[mot_xy - mot_stride             ][1];
 819             P_TOPRIGHT[0] = s->motion_val[mot_xy - mot_stride + off[block]][0];
 820             P_TOPRIGHT[1] = s->motion_val[mot_xy - mot_stride + off[block]][1];
 821             if(P_TOP[1]      > (rel_ymax4<<shift)) P_TOP[1]     = (rel_ymax4<<shift);
 822             if(P_TOPRIGHT[0] < (rel_xmin4<<shift)) P_TOPRIGHT[0]= (rel_xmin4<<shift);
 823             if(P_TOPRIGHT[0] > (rel_xmax4<<shift)) P_TOPRIGHT[0]= (rel_xmax4<<shift);
 824             if(P_TOPRIGHT[1] > (rel_ymax4<<shift)) P_TOPRIGHT[1]= (rel_ymax4<<shift);
 825
 826             P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
 827             P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
 828
 829             if(s->out_format == FMT_H263){
 830                 pred_x4 = P_MEDIAN[0];
 831                 pred_y4 = P_MEDIAN[1];
 832             }else { /* mpeg1 at least */
 833                 pred_x4= P_LEFT[0];
 834                 pred_y4= P_LEFT[1];
 835             }
 836         }
 837         P_MV1[0]= mx;
 838         P_MV1[1]= my;
 839
 840         dmin4 = s->me.motion_search[1](s, block, &mx4, &my4, P, pred_x4, pred_y4, rel_xmin4, rel_ymin4, rel_xmax4, rel_ymax4,
 841                                        &s->last_picture, s->p_mv_table, (1<<16)>>shift, mv_penalty);
 842
 843         dmin4= s->me.sub_motion_search(s, &mx4, &my4, dmin4, rel_xmin4, rel_ymin4, rel_xmax4, rel_ymax4,
 844                                           pred_x4, pred_y4, &s->last_picture, block, 1, mv_penalty);
 845
 846         s->motion_val[ s->block_index[block] ][0]= mx4;
 847         s->motion_val[ s->block_index[block] ][1]= my4;
 848         dmin_sum+= dmin4;
 849     }
 850     return dmin_sum;
 851 }
 852
 853 void ff_estimate_p_frame_motion(MpegEncContext * s,
 854                                 int mb_x, int mb_y)
 855 {
 856     UINT8 *pix, *ppix;
 857     int sum, varc, vard, mx, my, range, dmin, xx, yy;
 858     int xmin, ymin, xmax, ymax;
 859     int rel_xmin, rel_ymin, rel_xmax, rel_ymax;
 860     int pred_x=0, pred_y=0;
 861     int P[10][2];
 862     const int shift= 1+s->quarter_sample;
 863     int mb_type=0;
 864     uint8_t *ref_picture= s->last_picture.data[0];
 865     Picture * const pic= &s->current_picture;
 866     uint16_t * const mv_penalty= s->me.mv_penalty[s->f_code] + MAX_MV;
 867
 868     assert(s->quarter_sample==0 || s->quarter_sample==1);
 869
 870     s->me.penalty_factor    = get_penalty_factor(s, s->avctx->me_cmp);
 871     s->me.sub_penalty_factor= get_penalty_factor(s, s->avctx->me_sub_cmp);
 872
 873     get_limits(s, &range, &xmin, &ymin, &xmax, &ymax, s->f_code);
 874     rel_xmin= xmin - mb_x*16;
 875     rel_xmax= xmax - mb_x*16;
 876     rel_ymin= ymin - mb_y*16;
 877     rel_ymax= ymax - mb_y*16;
 878     s->me.skip=0;
 879
 880     switch(s->me_method) {
 881     case ME_ZERO:
 882     default:
 883         no_motion_search(s, &mx, &my);
 884         mx-= mb_x*16;
 885         my-= mb_y*16;
 886         dmin = 0;
 887         break;
 888     case ME_FULL:
 889         dmin = full_motion_search(s, &mx, &my, range, xmin, ymin, xmax, ymax, ref_picture);
 890         mx-= mb_x*16;
 891         my-= mb_y*16;
 892         break;
 893     case ME_LOG:
 894         dmin = log_motion_search(s, &mx, &my, range / 2, xmin, ymin, xmax, ymax, ref_picture);
 895         mx-= mb_x*16;
 896         my-= mb_y*16;
 897         break;
 898     case ME_PHODS:
 899         dmin = phods_motion_search(s, &mx, &my, range / 2, xmin, ymin, xmax, ymax, ref_picture);
 900         mx-= mb_x*16;
 901         my-= mb_y*16;
 902         break;
 903     case ME_X1:
 904     case ME_EPZS:
 905        {
 906             const int mot_stride = s->block_wrap[0];
 907             const int mot_xy = s->block_index[0];
 908
 909             P_LEFT[0]       = s->motion_val[mot_xy - 1][0];
 910             P_LEFT[1]       = s->motion_val[mot_xy - 1][1];
 911
 912             if(P_LEFT[0]       > (rel_xmax<<shift)) P_LEFT[0]       = (rel_xmax<<shift);
 913
 914             if(mb_y) {
 915                 P_TOP[0]      = s->motion_val[mot_xy - mot_stride    ][0];
 916                 P_TOP[1]      = s->motion_val[mot_xy - mot_stride    ][1];
 917                 P_TOPRIGHT[0] = s->motion_val[mot_xy - mot_stride + 2][0];
 918                 P_TOPRIGHT[1] = s->motion_val[mot_xy - mot_stride + 2][1];
 919                 if(P_TOP[1]      > (rel_ymax<<shift)) P_TOP[1]     = (rel_ymax<<shift);
 920                 if(P_TOPRIGHT[0] < (rel_xmin<<shift)) P_TOPRIGHT[0]= (rel_xmin<<shift);
 921                 if(P_TOPRIGHT[1] > (rel_ymax<<shift)) P_TOPRIGHT[1]= (rel_ymax<<shift);
 922
 923                 P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
 924                 P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
 925
 926                 if(s->out_format == FMT_H263){
 927                     pred_x = P_MEDIAN[0];
 928                     pred_y = P_MEDIAN[1];
 929                 }else { /* mpeg1 at least */
 930                     pred_x= P_LEFT[0];
 931                     pred_y= P_LEFT[1];
 932                 }
 933             }else{
 934                 pred_x= P_LEFT[0];
 935                 pred_y= P_LEFT[1];
 936             }
 937
 938         }
 939         dmin = s->me.motion_search[0](s, 0, &mx, &my, P, pred_x, pred_y, rel_xmin, rel_ymin, rel_xmax, rel_ymax,
 940                                       &s->last_picture, s->p_mv_table, (1<<16)>>shift, mv_penalty);
 941
 942         break;
 943     }
 944
 945     /* intra / predictive decision */
 946     xx = mb_x * 16;
 947     yy = mb_y * 16;
 948
 949     pix = s->new_picture.data[0] + (yy * s->linesize) + xx;
 950     /* At this point (mx,my) are full-pell and the relative displacement */
 951     ppix = ref_picture + ((yy+my) * s->linesize) + (xx+mx);
 952
 953     sum = s->dsp.pix_sum(pix, s->linesize);
 954
 955     varc = (s->dsp.pix_norm1(pix, s->linesize) - (((unsigned)(sum*sum))>>8) + 500 + 128)>>8;
 956     vard = (s->dsp.sse[0](NULL, pix, ppix, s->linesize)+128)>>8;
 957
 958 //printf("%d %d %d %X %X %X\n", s->mb_width, mb_x, mb_y,(int)s, (int)s->mb_var, (int)s->mc_mb_var); fflush(stdout);
 959     pic->mb_var   [s->mb_width * mb_y + mb_x] = varc;
 960     pic->mc_mb_var[s->mb_width * mb_y + mb_x] = vard;
 961     pic->mb_mean  [s->mb_width * mb_y + mb_x] = (sum+128)>>8;
 962     pic->mb_var_sum    += varc;
 963     pic->mc_mb_var_sum += vard;
 964 //printf("E%d %d %d %X %X %X\n", s->mb_width, mb_x, mb_y,(int)s, (int)s->mb_var, (int)s->mc_mb_var); fflush(stdout);
 965
 966 #if 0
 967     printf("varc=%4d avg_var=%4d (sum=%4d) vard=%4d mx=%2d my=%2d\n",
 968            varc, s->avg_mb_var, sum, vard, mx - xx, my - yy);
 969 #endif
 970     if(s->flags&CODEC_FLAG_HQ){
 971         if (vard <= 64 || vard < varc)
 972             s->scene_change_score+= ff_sqrt(vard) - ff_sqrt(varc);
 973         else
 974             s->scene_change_score+= s->qscale;
 975
 976         if (vard*2 + 200 > varc)
 977             mb_type|= MB_TYPE_INTRA;
 978         if (varc*2 + 200 > vard){
 979             mb_type|= MB_TYPE_INTER;
 980             s->me.sub_motion_search(s, &mx, &my, dmin, rel_xmin, rel_ymin, rel_xmax, rel_ymax,
 981                                    pred_x, pred_y, &s->last_picture, 0, 0, mv_penalty);
 982         }else{
 983             mx <<=shift;
 984             my <<=shift;
 985         }
 986         if((s->flags&CODEC_FLAG_4MV)
 987            && !s->me.skip && varc>50 && vard>10){
 988             mv4_search(s, rel_xmin, rel_ymin, rel_xmax, rel_ymax, mx, my, shift);
 989             mb_type|=MB_TYPE_INTER4V;
 990
 991             set_p_mv_tables(s, mx, my, 0);
 992         }else
 993             set_p_mv_tables(s, mx, my, 1);
 994     }else{
 995         if (vard <= 64 || vard < varc) {
 996 //        if (sadP <= 32 || sadP < sadI + 500) {
 997             s->scene_change_score+= ff_sqrt(vard) - ff_sqrt(varc);
 998             mb_type|= MB_TYPE_INTER;
 999             if (s->me_method != ME_ZERO) {
1000                 dmin= s->me.sub_motion_search(s, &mx, &my, dmin, rel_xmin, rel_ymin, rel_xmax, rel_ymax,
1001                                             pred_x, pred_y, &s->last_picture, 0, 0, mv_penalty);
1002                 if((s->flags&CODEC_FLAG_4MV)
1003                    && !s->me.skip && varc>50 && vard>10){
1004                     int dmin4= mv4_search(s, rel_xmin, rel_ymin, rel_xmax, rel_ymax, mx, my, shift);
1005                     if(dmin4 + 128 <dmin)
1006                         mb_type= MB_TYPE_INTER4V;
1007                 }
1008                 set_p_mv_tables(s, mx, my, mb_type!=MB_TYPE_INTER4V);
1009
1010             } else {
1011                 mx <<=shift;
1012                 my <<=shift;
1013             }
1014 #if 0
1015             if (vard < 10) {
1016                 skip++;
1017                 fprintf(stderr,"\nEarly skip: %d vard: %2d varc: %5d dmin: %d",
1018                                 skip, vard, varc, dmin);
1019             }
1020 #endif
1021         }else{
1022             s->scene_change_score+= 20;
1023             mb_type|= MB_TYPE_INTRA;
1024             mx = 0;
1025             my = 0;
1026         }
1027     }
1028
1029     s->mb_type[mb_y*s->mb_width + mb_x]= mb_type;
1030 }
1031
1032 int ff_pre_estimate_p_frame_motion(MpegEncContext * s,
1033                                     int mb_x, int mb_y)
1034 {
1035     int mx, my, range, dmin;
1036     int xmin, ymin, xmax, ymax;
1037     int rel_xmin, rel_ymin, rel_xmax, rel_ymax;
1038     int pred_x=0, pred_y=0;
1039     int P[10][2];
1040     const int shift= 1+s->quarter_sample;
1041     uint16_t * const mv_penalty= s->me.mv_penalty[s->f_code] + MAX_MV;
1042     const int mv_stride= s->mb_width + 2;
1043     const int xy= mb_x + 1 + (mb_y + 1)*mv_stride;
1044
1045     assert(s->quarter_sample==0 || s->quarter_sample==1);
1046
1047     s->me.pre_penalty_factor    = get_penalty_factor(s, s->avctx->me_pre_cmp);
1048
1049     get_limits(s, &range, &xmin, &ymin, &xmax, &ymax, s->f_code);
1050     rel_xmin= xmin - mb_x*16;
1051     rel_xmax= xmax - mb_x*16;
1052     rel_ymin= ymin - mb_y*16;
1053     rel_ymax= ymax - mb_y*16;
1054     s->me.skip=0;
1055
1056     P_LEFT[0]       = s->p_mv_table[xy + 1][0];
1057     P_LEFT[1]       = s->p_mv_table[xy + 1][1];
1058
1059     if(P_LEFT[0]       < (rel_xmin<<shift)) P_LEFT[0]       = (rel_xmin<<shift);
1060
1061     /* special case for first line */
1062     if (mb_y == s->mb_height-1) {
1063         pred_x= P_LEFT[0];
1064         pred_y= P_LEFT[1];
1065         P_TOP[0]= P_TOPRIGHT[0]= P_MEDIAN[0]=
1066         P_TOP[1]= P_TOPRIGHT[1]= P_MEDIAN[1]= 0; //FIXME
1067     } else {
1068         P_TOP[0]      = s->p_mv_table[xy + mv_stride    ][0];
1069         P_TOP[1]      = s->p_mv_table[xy + mv_stride    ][1];
1070         P_TOPRIGHT[0] = s->p_mv_table[xy + mv_stride - 1][0];
1071         P_TOPRIGHT[1] = s->p_mv_table[xy + mv_stride - 1][1];
1072         if(P_TOP[1]      < (rel_ymin<<shift)) P_TOP[1]     = (rel_ymin<<shift);
1073         if(P_TOPRIGHT[0] > (rel_xmax<<shift)) P_TOPRIGHT[0]= (rel_xmax<<shift);
1074         if(P_TOPRIGHT[1] < (rel_ymin<<shift)) P_TOPRIGHT[1]= (rel_ymin<<shift);
1075
1076         P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
1077         P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
1078
1079         pred_x = P_MEDIAN[0];
1080         pred_y = P_MEDIAN[1];
1081     }
1082     dmin = s->me.pre_motion_search(s, 0, &mx, &my, P, pred_x, pred_y, rel_xmin, rel_ymin, rel_xmax, rel_ymax,
1083                                    &s->last_picture, s->p_mv_table, (1<<16)>>shift, mv_penalty);
1084
1085     s->p_mv_table[xy][0] = mx<<shift;
1086     s->p_mv_table[xy][1] = my<<shift;
1087
1088     return dmin;
1089 }
1090
1091 int ff_estimate_motion_b(MpegEncContext * s,
1092                        int mb_x, int mb_y, int16_t (*mv_table)[2], Picture *picture, int f_code)
1093 {
1094     int mx, my, range, dmin;
1095     int xmin, ymin, xmax, ymax;
1096     int rel_xmin, rel_ymin, rel_xmax, rel_ymax;
1097     int pred_x=0, pred_y=0;
1098     int P[10][2];
1099     const int shift= 1+s->quarter_sample;
1100     const int mot_stride = s->mb_width + 2;
1101     const int mot_xy = (mb_y + 1)*mot_stride + mb_x + 1;
1102     uint8_t * const ref_picture= picture->data[0];
1103     uint16_t * const mv_penalty= s->me.mv_penalty[f_code] + MAX_MV;
1104     int mv_scale;
1105
1106     s->me.penalty_factor    = get_penalty_factor(s, s->avctx->me_cmp);
1107     s->me.sub_penalty_factor= get_penalty_factor(s, s->avctx->me_sub_cmp);
1108
1109     get_limits(s, &range, &xmin, &ymin, &xmax, &ymax, f_code);
1110     rel_xmin= xmin - mb_x*16;
1111     rel_xmax= xmax - mb_x*16;
1112     rel_ymin= ymin - mb_y*16;
1113     rel_ymax= ymax - mb_y*16;
1114
1115     switch(s->me_method) {
1116     case ME_ZERO:
1117     default:
1118         no_motion_search(s, &mx, &my);
1119         dmin = 0;
1120         mx-= mb_x*16;
1121         my-= mb_y*16;
1122         break;
1123     case ME_FULL:
1124         dmin = full_motion_search(s, &mx, &my, range, xmin, ymin, xmax, ymax, ref_picture);
1125         mx-= mb_x*16;
1126         my-= mb_y*16;
1127         break;
1128     case ME_LOG:
1129         dmin = log_motion_search(s, &mx, &my, range / 2, xmin, ymin, xmax, ymax, ref_picture);
1130         mx-= mb_x*16;
1131         my-= mb_y*16;
1132         break;
1133     case ME_PHODS:
1134         dmin = phods_motion_search(s, &mx, &my, range / 2, xmin, ymin, xmax, ymax, ref_picture);
1135         mx-= mb_x*16;
1136         my-= mb_y*16;
1137         break;
1138     case ME_X1:
1139     case ME_EPZS:
1140        {
1141             P_LEFT[0]        = mv_table[mot_xy - 1][0];
1142             P_LEFT[1]        = mv_table[mot_xy - 1][1];
1143
1144             if(P_LEFT[0]       > (rel_xmax<<shift)) P_LEFT[0]       = (rel_xmax<<shift);
1145
1146             /* special case for first line */
1147             if (mb_y) {
1148                 P_TOP[0] = mv_table[mot_xy - mot_stride             ][0];
1149                 P_TOP[1] = mv_table[mot_xy - mot_stride             ][1];
1150                 P_TOPRIGHT[0] = mv_table[mot_xy - mot_stride + 1         ][0];
1151                 P_TOPRIGHT[1] = mv_table[mot_xy - mot_stride + 1         ][1];
1152                 if(P_TOP[1] > (rel_ymax<<shift)) P_TOP[1]= (rel_ymax<<shift);
1153                 if(P_TOPRIGHT[0] < (rel_xmin<<shift)) P_TOPRIGHT[0]= (rel_xmin<<shift);
1154                 if(P_TOPRIGHT[1] > (rel_ymax<<shift)) P_TOPRIGHT[1]= (rel_ymax<<shift);
1155
1156                 P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
1157                 P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
1158             }
1159             pred_x= P_LEFT[0];
1160             pred_y= P_LEFT[1];
1161         }
1162
1163         if(mv_table == s->b_forw_mv_table){
1164             mv_scale= (s->pb_time<<16) / (s->pp_time<<shift);
1165         }else{
1166             mv_scale= ((s->pb_time - s->pp_time)<<16) / (s->pp_time<<shift);
1167         }
1168
1169         dmin = s->me.motion_search[0](s, 0, &mx, &my, P, pred_x, pred_y, rel_xmin, rel_ymin, rel_xmax, rel_ymax,
1170                                       picture, s->p_mv_table, mv_scale, mv_penalty);
1171
1172         break;
1173     }
1174
1175     dmin= s->me.sub_motion_search(s, &mx, &my, dmin, rel_xmin, rel_ymin, rel_xmax, rel_ymax,
1176                                    pred_x, pred_y, picture, 0, 0, mv_penalty);
1177 //printf("%d %d %d %d//", s->mb_x, s->mb_y, mx, my);
1178 //    s->mb_type[mb_y*s->mb_width + mb_x]= mb_type;
1179     mv_table[mot_xy][0]= mx;
1180     mv_table[mot_xy][1]= my;
1181
1182     return dmin;
1183 }
1184
1185 static inline int check_bidir_mv(MpegEncContext * s,
1186                    int mb_x, int mb_y,
1187                    int motion_fx, int motion_fy,
1188                    int motion_bx, int motion_by,
1189                    int pred_fx, int pred_fy,
1190                    int pred_bx, int pred_by)
1191 {
1192     //FIXME optimize?
1193     //FIXME move into template?
1194     //FIXME better f_code prediction (max mv & distance)
1195     UINT16 *mv_penalty= s->me.mv_penalty[s->f_code] + MAX_MV; // f_code of the prev frame
1196     uint8_t *dest_y = s->me.scratchpad;
1197     uint8_t *ptr;
1198     int dxy;
1199     int src_x, src_y;
1200     int fbmin;
1201
1202     if(s->quarter_sample){
1203         dxy = ((motion_fy & 3) << 2) | (motion_fx & 3);
1204         src_x = mb_x * 16 + (motion_fx >> 2);
1205         src_y = mb_y * 16 + (motion_fy >> 2);
1206         assert(src_x >=-16 && src_x<=s->width);
1207         assert(src_y >=-16 && src_y<=s->height);
1208
1209         ptr = s->last_picture.data[0] + (src_y * s->linesize) + src_x;
1210         s->dsp.put_qpel_pixels_tab[0][dxy](dest_y    , ptr    , s->linesize);
1211
1212         dxy = ((motion_by & 3) << 2) | (motion_bx & 3);
1213         src_x = mb_x * 16 + (motion_bx >> 2);
1214         src_y = mb_y * 16 + (motion_by >> 2);
1215         assert(src_x >=-16 && src_x<=s->width);
1216         assert(src_y >=-16 && src_y<=s->height);
1217
1218         ptr = s->next_picture.data[0] + (src_y * s->linesize) + src_x;
1219         s->dsp.avg_qpel_pixels_tab[0][dxy](dest_y    , ptr    , s->linesize);
1220     }else{
1221         dxy = ((motion_fy & 1) << 1) | (motion_fx & 1);
1222         src_x = mb_x * 16 + (motion_fx >> 1);
1223         src_y = mb_y * 16 + (motion_fy >> 1);
1224         assert(src_x >=-16 && src_x<=s->width);
1225         assert(src_y >=-16 && src_y<=s->height);
1226
1227         ptr = s->last_picture.data[0] + (src_y * s->linesize) + src_x;
1228         s->dsp.put_pixels_tab[0][dxy](dest_y    , ptr    , s->linesize, 16);
1229
1230         dxy = ((motion_by & 1) << 1) | (motion_bx & 1);
1231         src_x = mb_x * 16 + (motion_bx >> 1);
1232         src_y = mb_y * 16 + (motion_by >> 1);
1233         assert(src_x >=-16 && src_x<=s->width);
1234         assert(src_y >=-16 && src_y<=s->height);
1235
1236         ptr = s->next_picture.data[0] + (src_y * s->linesize) + src_x;
1237         s->dsp.avg_pixels_tab[0][dxy](dest_y    , ptr    , s->linesize, 16);
1238     }
1239
1240     fbmin = (mv_penalty[motion_fx-pred_fx] + mv_penalty[motion_fy-pred_fy])*s->me.sub_penalty_factor
1241            +(mv_penalty[motion_bx-pred_bx] + mv_penalty[motion_by-pred_by])*s->me.sub_penalty_factor;
1242            + s->dsp.me_sub_cmp[0](s, s->new_picture.data[0] + mb_x*16 + mb_y*16*s->linesize, dest_y, s->linesize);
1243
1244     return fbmin;
1245 }
1246
1247 /* refine the bidir vectors in hq mode and return the score in both lq & hq mode*/
1248 static inline int bidir_refine(MpegEncContext * s,
1249                                   int mb_x, int mb_y)
1250 {
1251     const int mot_stride = s->mb_width + 2;
1252     const int xy = (mb_y + 1)*mot_stride + mb_x + 1;
1253     int fbmin;
1254     int pred_fx= s->b_bidir_forw_mv_table[xy-1][0];
1255     int pred_fy= s->b_bidir_forw_mv_table[xy-1][1];
1256     int pred_bx= s->b_bidir_back_mv_table[xy-1][0];
1257     int pred_by= s->b_bidir_back_mv_table[xy-1][1];
1258     int motion_fx= s->b_bidir_forw_mv_table[xy][0]= s->b_forw_mv_table[xy][0];
1259     int motion_fy= s->b_bidir_forw_mv_table[xy][1]= s->b_forw_mv_table[xy][1];
1260     int motion_bx= s->b_bidir_back_mv_table[xy][0]= s->b_back_mv_table[xy][0];
1261     int motion_by= s->b_bidir_back_mv_table[xy][1]= s->b_back_mv_table[xy][1];
1262
1263     //FIXME do refinement and add flag
1264
1265     fbmin= check_bidir_mv(s, mb_x, mb_y,
1266                           motion_fx, motion_fy,
1267                           motion_bx, motion_by,
1268                           pred_fx, pred_fy,
1269                           pred_bx, pred_by);
1270
1271    return fbmin;
1272 }
1273
1274 static inline int direct_search(MpegEncContext * s,
1275                                 int mb_x, int mb_y)
1276 {
1277     int P[10][2];
1278     const int mot_stride = s->mb_width + 2;
1279     const int mot_xy = (mb_y + 1)*mot_stride + mb_x + 1;
1280     const int shift= 1+s->quarter_sample;
1281     int dmin, i;
1282     const int time_pp= s->pp_time;
1283     const int time_pb= s->pb_time;
1284     int mx, my, xmin, xmax, ymin, ymax;
1285     int16_t (*mv_table)[2]= s->b_direct_mv_table;
1286     uint16_t * const mv_penalty= s->me.mv_penalty[1] + MAX_MV;
1287
1288     ymin= xmin=(-32)>>shift;
1289     ymax= xmax=   31>>shift;
1290
1291     if(s->co_located_type_table[mb_x + mb_y*s->mb_width]==CO_LOCATED_TYPE_4MV){
1292         s->mv_type= MV_TYPE_8X8;
1293     }else{
1294         s->mv_type= MV_TYPE_16X16;
1295     }
1296
1297     for(i=0; i<4; i++){
1298         int index= s->block_index[i];
1299         int min, max;
1300
1301         s->me.co_located_mv[i][0]= s->motion_val[index][0];
1302         s->me.co_located_mv[i][1]= s->motion_val[index][1];
1303         s->me.direct_basis_mv[i][0]= s->me.co_located_mv[i][0]*time_pb/time_pp + ((i& 1)<<(shift+3));
1304         s->me.direct_basis_mv[i][1]= s->me.co_located_mv[i][1]*time_pb/time_pp + ((i>>1)<<(shift+3));
1305 //        s->me.direct_basis_mv[1][i][0]= s->me.co_located_mv[i][0]*(time_pb - time_pp)/time_pp + ((i &1)<<(shift+3);
1306 //        s->me.direct_basis_mv[1][i][1]= s->me.co_located_mv[i][1]*(time_pb - time_pp)/time_pp + ((i>>1)<<(shift+3);
1307
1308         max= FFMAX(s->me.direct_basis_mv[i][0], s->me.direct_basis_mv[i][0] - s->me.co_located_mv[i][0])>>shift;
1309         min= FFMIN(s->me.direct_basis_mv[i][0], s->me.direct_basis_mv[i][0] - s->me.co_located_mv[i][0])>>shift;
1310         max+= (2*mb_x + (i& 1))*8 - 1; // +-1 is for the simpler rounding
1311         min+= (2*mb_x + (i& 1))*8 + 1;
1312         xmax= FFMIN(xmax, s->width - max);
1313         xmin= FFMAX(xmin, - 16     - min);
1314
1315         max= FFMAX(s->me.direct_basis_mv[i][1], s->me.direct_basis_mv[i][1] - s->me.co_located_mv[i][1])>>shift;
1316         min= FFMIN(s->me.direct_basis_mv[i][1], s->me.direct_basis_mv[i][1] - s->me.co_located_mv[i][1])>>shift;
1317         max+= (2*mb_y + (i>>1))*8 - 1; // +-1 is for the simpler rounding
1318         min+= (2*mb_y + (i>>1))*8 + 1;
1319         ymax= FFMIN(ymax, s->height - max);
1320         ymin= FFMAX(ymin, - 16      - min);
1321
1322         if(s->mv_type == MV_TYPE_16X16) break;
1323     }
1324
1325     assert(xmax <= 15 && ymax <= 15 && xmin >= -16 && ymin >= -16);
1326
1327     if(xmax < 0 || xmin >0 || ymax < 0 || ymin > 0){
1328         s->b_direct_mv_table[mot_xy][0]= 0;
1329         s->b_direct_mv_table[mot_xy][1]= 0;
1330
1331         return 256*256*256*64;
1332     }
1333
1334     P_LEFT[0]        = clip(mv_table[mot_xy - 1][0], xmin<<shift, xmax<<shift);
1335     P_LEFT[1]        = clip(mv_table[mot_xy - 1][1], ymin<<shift, ymax<<shift);
1336
1337     /* special case for first line */
1338     if (mb_y) {
1339         P_TOP[0]      = clip(mv_table[mot_xy - mot_stride             ][0], xmin<<shift, xmax<<shift);
1340         P_TOP[1]      = clip(mv_table[mot_xy - mot_stride             ][1], ymin<<shift, ymax<<shift);
1341         P_TOPRIGHT[0] = clip(mv_table[mot_xy - mot_stride + 1         ][0], xmin<<shift, xmax<<shift);
1342         P_TOPRIGHT[1] = clip(mv_table[mot_xy - mot_stride + 1         ][1], ymin<<shift, ymax<<shift);
1343
1344         P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
1345         P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
1346     }
1347
1348     if(s->flags&CODEC_FLAG_QPEL){
1349         dmin = simple_direct_qpel_epzs_motion_search(s, 0, &mx, &my, P, 0, 0, xmin, ymin, xmax, ymax,
1350                                                      &s->last_picture, mv_table, 1<<14, mv_penalty);
1351         dmin = simple_direct_qpel_qpel_motion_search(s, &mx, &my, dmin, xmin, ymin, xmax, ymax,
1352                                                 0, 0, &s->last_picture, 0, 0, mv_penalty);
1353     }else{
1354         dmin = simple_direct_hpel_epzs_motion_search(s, 0, &mx, &my, P, 0, 0, xmin, ymin, xmax, ymax,
1355                                                      &s->last_picture, mv_table, 1<<15, mv_penalty);
1356         dmin = simple_direct_hpel_hpel_motion_search(s, &mx, &my, dmin, xmin, ymin, xmax, ymax,
1357                                                 0, 0, &s->last_picture, 0, 0, mv_penalty);
1358     }
1359
1360     s->b_direct_mv_table[mot_xy][0]= mx;
1361     s->b_direct_mv_table[mot_xy][1]= my;
1362     return dmin;
1363 }
1364
1365 void ff_estimate_b_frame_motion(MpegEncContext * s,
1366                              int mb_x, int mb_y)
1367 {
1368     const int penalty_factor= s->me.penalty_factor;
1369     int fmin, bmin, dmin, fbmin;
1370     int type=0;
1371
1372     dmin= direct_search(s, mb_x, mb_y);
1373
1374     fmin= ff_estimate_motion_b(s, mb_x, mb_y, s->b_forw_mv_table, &s->last_picture, s->f_code);
1375     bmin= ff_estimate_motion_b(s, mb_x, mb_y, s->b_back_mv_table, &s->next_picture, s->b_code) - penalty_factor;
1376 //printf(" %d %d ", s->b_forw_mv_table[xy][0], s->b_forw_mv_table[xy][1]);
1377
1378     fbmin= bidir_refine(s, mb_x, mb_y);
1379
1380     {
1381         int score= dmin;
1382         type=MB_TYPE_DIRECT;
1383
1384         if(fmin<score){
1385             score=fmin;
1386             type= MB_TYPE_FORWARD;
1387         }
1388         if(bmin<score){
1389             score=bmin;
1390             type= MB_TYPE_BACKWARD;
1391         }
1392         if(fbmin<score){
1393             score=fbmin;
1394             type= MB_TYPE_BIDIR;
1395         }
1396         score= ((unsigned)(score*score + 128*256))>>16;
1397         s->current_picture.mc_mb_var_sum += score;
1398         s->current_picture.mc_mb_var[mb_y*s->mb_width + mb_x] = score; //FIXME use SSD
1399     }
1400
1401     if(s->flags&CODEC_FLAG_HQ){
1402         type= MB_TYPE_FORWARD | MB_TYPE_BACKWARD | MB_TYPE_BIDIR | MB_TYPE_DIRECT; //FIXME something smarter
1403         if(dmin>256*256*16) type&= ~MB_TYPE_DIRECT; //dont try direct mode if its invalid for this MB
1404     }
1405
1406     s->mb_type[mb_y*s->mb_width + mb_x]= type;
1407 }
1408
1409 /* find best f_code for ME which do unlimited searches */
1410 int ff_get_best_fcode(MpegEncContext * s, int16_t (*mv_table)[2], int type)
1411 {
1412     if(s->me_method>=ME_EPZS){
1413         int score[8];
1414         int i, y;
1415         UINT8 * fcode_tab= s->fcode_tab;
1416         int best_fcode=-1;
1417         int best_score=-10000000;
1418
1419         for(i=0; i<8; i++) score[i]= s->mb_num*(8-i);
1420
1421         for(y=0; y<s->mb_height; y++){
1422             int x;
1423             int xy= (y+1)* (s->mb_width+2) + 1;
1424             i= y*s->mb_width;
1425             for(x=0; x<s->mb_width; x++){
1426                 if(s->mb_type[i] & type){
1427                     int fcode= FFMAX(fcode_tab[mv_table[xy][0] + MAX_MV],
1428                                      fcode_tab[mv_table[xy][1] + MAX_MV]);
1429                     int j;
1430
1431                     for(j=0; j<fcode && j<8; j++){
1432                         if(s->pict_type==B_TYPE || s->current_picture.mc_mb_var[i] < s->current_picture.mb_var[i])
1433                             score[j]-= 170;
1434                     }
1435                 }
1436                 i++;
1437                 xy++;
1438             }
1439         }
1440
1441         for(i=1; i<8; i++){
1442             if(score[i] > best_score){
1443                 best_score= score[i];
1444                 best_fcode= i;
1445             }
1446 //            printf("%d %d\n", i, score[i]);
1447         }
1448
1449 //    printf("fcode: %d type: %d\n", i, s->pict_type);
1450         return best_fcode;
1451 /*        for(i=0; i<=MAX_FCODE; i++){
1452             printf("%d ", mv_num[i]);
1453         }
1454         printf("\n");*/
1455     }else{
1456         return 1;
1457     }
1458 }
1459
1460 void ff_fix_long_p_mvs(MpegEncContext * s)
1461 {
1462     const int f_code= s->f_code;
1463     int y;
1464     UINT8 * fcode_tab= s->fcode_tab;
1465 //int clip=0;
1466 //int noclip=0;
1467     /* clip / convert to intra 16x16 type MVs */
1468     for(y=0; y<s->mb_height; y++){
1469         int x;
1470         int xy= (y+1)* (s->mb_width+2)+1;
1471         int i= y*s->mb_width;
1472         for(x=0; x<s->mb_width; x++){
1473             if(s->mb_type[i]&MB_TYPE_INTER){
1474                 if(   fcode_tab[s->p_mv_table[xy][0] + MAX_MV] > f_code
1475                    || fcode_tab[s->p_mv_table[xy][0] + MAX_MV] == 0
1476                    || fcode_tab[s->p_mv_table[xy][1] + MAX_MV] > f_code
1477                    || fcode_tab[s->p_mv_table[xy][1] + MAX_MV] == 0 ){
1478                     s->mb_type[i] &= ~MB_TYPE_INTER;
1479                     s->mb_type[i] |= MB_TYPE_INTRA;
1480                     s->p_mv_table[xy][0] = 0;
1481                     s->p_mv_table[xy][1] = 0;
1482 //clip++;
1483                 }
1484 //else
1485 //  noclip++;
1486             }
1487             xy++;
1488             i++;
1489         }
1490     }
1491 //printf("%d no:%d %d//\n", clip, noclip, f_code);
1492     if(s->flags&CODEC_FLAG_4MV){
1493         const int wrap= 2+ s->mb_width*2;
1494
1495         /* clip / convert to intra 8x8 type MVs */
1496         for(y=0; y<s->mb_height; y++){
1497             int xy= (y*2 + 1)*wrap + 1;
1498             int i= y*s->mb_width;
1499             int x;
1500
1501             for(x=0; x<s->mb_width; x++){
1502                 if(s->mb_type[i]&MB_TYPE_INTER4V){
1503                     int block;
1504                     for(block=0; block<4; block++){
1505                         int off= (block& 1) + (block>>1)*wrap;
1506                         int mx= s->motion_val[ xy + off ][0];
1507                         int my= s->motion_val[ xy + off ][1];
1508
1509                         if(   fcode_tab[mx + MAX_MV] > f_code
1510                            || fcode_tab[mx + MAX_MV] == 0
1511                            || fcode_tab[my + MAX_MV] > f_code
1512                            || fcode_tab[my + MAX_MV] == 0 ){
1513                             s->mb_type[i] &= ~MB_TYPE_INTER4V;
1514                             s->mb_type[i] |= MB_TYPE_INTRA;
1515                         }
1516                     }
1517                 }
1518                 xy+=2;
1519                 i++;
1520             }
1521         }
1522     }
1523 }
1524
1525 void ff_fix_long_b_mvs(MpegEncContext * s, int16_t (*mv_table)[2], int f_code, int type)
1526 {
1527     int y;
1528     UINT8 * fcode_tab= s->fcode_tab;
1529
1530     /* clip / convert to intra 16x16 type MVs */
1531     for(y=0; y<s->mb_height; y++){
1532         int x;
1533         int xy= (y+1)* (s->mb_width+2)+1;
1534         int i= y*s->mb_width;
1535         for(x=0; x<s->mb_width; x++){
1536             if(   fcode_tab[mv_table[xy][0] + MAX_MV] > f_code
1537                || fcode_tab[mv_table[xy][0] + MAX_MV] == 0){
1538                 if(mv_table[xy][0]>0) mv_table[xy][0]=  (16<<f_code)-1;
1539                 else                  mv_table[xy][0]= -(16<<f_code);
1540             }
1541             if(   fcode_tab[mv_table[xy][1] + MAX_MV] > f_code
1542                || fcode_tab[mv_table[xy][1] + MAX_MV] == 0){
1543                 if(mv_table[xy][1]>0) mv_table[xy][1]=  (16<<f_code)-1;
1544                 else                  mv_table[xy][1]= -(16<<f_code);
1545             }
1546             xy++;
1547             i++;
1548         }
1549     }
1550 }