motion_est.c
上传用户:jxp0626
上传日期:2007-01-08
资源大小:102k
文件大小:14k
源码类别:

流媒体/Mpeg4/MP4

开发平台:

Unix_Linux

  1. /*
  2.  * Motion estimation 
  3.  * Copyright (c) 2000 Gerard Lantau.
  4.  *
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or
  8.  * (at your option) any later version.
  9.  *
  10.  * This program is distributed in the hope that it will be useful,
  11.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13.  * GNU General Public License for more details.
  14.  *
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19. #include <stdlib.h>
  20. #include <stdio.h>
  21. #include "avcodec.h"
  22. #include "dsputil.h"
  23. #include "mpegvideo.h"
  24. static int pix_sum(UINT8 *pix, int line_size)
  25. {
  26.     int s, i, j;
  27.     s = 0;
  28.     for(i=0;i<16;i++) {
  29.         for(j=0;j<16;j+=8) {
  30.             s += pix[0];
  31.             s += pix[1];
  32.             s += pix[2];
  33.             s += pix[3];
  34.             s += pix[4];
  35.             s += pix[5];
  36.             s += pix[6];
  37.             s += pix[7];
  38.             pix += 8;
  39.         }
  40.         pix += line_size - 16;
  41.     }
  42.     return s;
  43. }
  44. static int pix_norm1(UINT8 *pix, int line_size)
  45. {
  46.     int s, i, j;
  47.     UINT32 *sq = squareTbl + 256;
  48.     s = 0;
  49.     for(i=0;i<16;i++) {
  50.         for(j=0;j<16;j+=8) {
  51.             s += sq[pix[0]];
  52.             s += sq[pix[1]];
  53.             s += sq[pix[2]];
  54.             s += sq[pix[3]];
  55.             s += sq[pix[4]];
  56.             s += sq[pix[5]];
  57.             s += sq[pix[6]];
  58.             s += sq[pix[7]];
  59.             pix += 8;
  60.         }
  61.         pix += line_size - 16;
  62.     }
  63.     return s;
  64. }
  65. static int pix_norm(UINT8 *pix1, UINT8 *pix2, int line_size)
  66. {
  67.     int s, i, j;
  68.     UINT32 *sq = squareTbl + 256;
  69.     s = 0;
  70.     for(i=0;i<16;i++) {
  71.         for(j=0;j<16;j+=8) {
  72.             s += sq[pix1[0] - pix2[0]];
  73.             s += sq[pix1[1] - pix2[1]];
  74.             s += sq[pix1[2] - pix2[2]];
  75.             s += sq[pix1[3] - pix2[3]];
  76.             s += sq[pix1[4] - pix2[4]];
  77.             s += sq[pix1[5] - pix2[5]];
  78.             s += sq[pix1[6] - pix2[6]];
  79.             s += sq[pix1[7] - pix2[7]];
  80.             pix1 += 8;
  81.             pix2 += 8;
  82.         }
  83.         pix1 += line_size - 16;
  84.         pix2 += line_size - 16;
  85.     }
  86.     return s;
  87. }
  88. #define USE_MMX
  89. #ifdef USE_MMX
  90. /* original mmx code from mpeg2movie by Adam Williams */
  91. static void inline mmx_start_block(void)
  92. {
  93. asm(" 
  94. .align 8
  95. pxor %%mm7, %%mm7; 
  96. pxor %%mm6, %%mm6;
  97. " : : );
  98. }
  99. static void inline mmx_absdiff(unsigned char *a, unsigned char *b)
  100. {
  101. asm("
  102. .align 8
  103. movq (%%ebx), %%mm0;     // Get first half of row1
  104. movq (%%ecx), %%mm1;     // Get first half of row2
  105. movq %%mm0, %%mm2;     // Make a copy of row1 for absdiff operation
  106. movq 8(%%ebx), %%mm3;     // Get second half of row1
  107. psubusb %%mm1, %%mm0;     // Subtract the first halves one way
  108. psubusb %%mm2, %%mm1;     // Subtract the other way
  109. movq        8(%%ecx),   %%mm4;     // Get second half of row2
  110. por %%mm1,      %%mm0;     // Merge first half results
  111. movq %%mm3, %%mm5;     // Copy for absdiff operation
  112. movq %%mm0, %%mm1;     // Keep a copy
  113. psubusb %%mm4, %%mm3;     // Subtract second halves one way
  114. punpcklbw %%mm6, %%mm0;     // Unpack to higher precision for accumulation
  115. psubusb %%mm5, %%mm4;     // Subtract the other way
  116. psrlq $32, %%mm1;     // Shift registeres for accumulation
  117. por %%mm4, %%mm3;     // merge results of 2nd half
  118. punpcklbw %%mm6, %%mm1;     // unpack to higher precision for accumulation
  119. movq %%mm3, %%mm4;     // keep a copy
  120. punpcklbw %%mm6, %%mm3;     // unpack to higher precision for accumulation
  121. paddw %%mm0, %%mm7;     // accumulate difference
  122. psrlq $32, %%mm4;     // shift results for accumulation
  123. paddw %%mm1, %%mm7;     // accumulate difference
  124. punpcklbw %%mm6, %%mm4;     // unpack to higher precision for accumulation
  125. paddw %%mm3, %%mm7;     // accumulate difference
  126. paddw %%mm4, %%mm7;     // accumulate difference
  127. "
  128. : "b" (a), "c" (b) 
  129. : "st" );
  130. }
  131. static inline unsigned int mmx_accum_absdiff()
  132. {
  133. unsigned long long r = 0;
  134. asm("
  135. .align 8
  136. movq %%mm7, %%mm5;
  137. movq %%mm7, %%mm4;
  138. punpcklwd %%mm6, %%mm4;
  139. punpckhwd %%mm6, %%mm5;
  140. paddd %%mm5, %%mm4;
  141. movq %%mm4, %%mm5;
  142. punpckldq %%mm6, %%mm5;
  143. punpckhdq %%mm6, %%mm4;
  144. paddd %%mm5, %%mm4;
  145. movq %%mm4, (%%ebx);
  146. emms;
  147. "
  148. : :  "b" (&r));
  149. return r;
  150. }
  151. static unsigned int pix_abs16x16(UINT8 *pix1, UINT8 *pix2, int line_size)
  152. {
  153.     int i;
  154.     mmx_start_block();
  155.     for(i=0;i<16;i++) {
  156.         mmx_absdiff(pix1, pix2);
  157.         pix1 += line_size;
  158.         pix2 += line_size;
  159.     }
  160.     return mmx_accum_absdiff();
  161. }
  162. #else
  163. static int pix_abs16x16(UINT8 *pix1, UINT8 *pix2, int line_size)
  164. {
  165.     int s, i;
  166.     s = 0;
  167.     for(i=0;i<16;i++) {
  168.         s += abs(pix1[0] - pix2[0]);
  169.         s += abs(pix1[1] - pix2[1]);
  170.         s += abs(pix1[2] - pix2[2]);
  171.         s += abs(pix1[3] - pix2[3]);
  172.         s += abs(pix1[4] - pix2[4]);
  173.         s += abs(pix1[5] - pix2[5]);
  174.         s += abs(pix1[6] - pix2[6]);
  175.         s += abs(pix1[7] - pix2[7]);
  176.         s += abs(pix1[8] - pix2[8]);
  177.         s += abs(pix1[9] - pix2[9]);
  178.         s += abs(pix1[10] - pix2[10]);
  179.         s += abs(pix1[11] - pix2[11]);
  180.         s += abs(pix1[12] - pix2[12]);
  181.         s += abs(pix1[13] - pix2[13]);
  182.         s += abs(pix1[14] - pix2[14]);
  183.         s += abs(pix1[15] - pix2[15]);
  184.         pix1 += line_size;
  185.         pix2 += line_size;
  186.     }
  187.     return s;
  188. }
  189. #endif
  190. #define avg2(a,b) ((a+b+1)>>1)
  191. #define avg4(a,b,c,d) ((a+b+c+d+2)>>2)
  192. static int pix_abs16x16_x2(UINT8 *pix1, UINT8 *pix2, int line_size)
  193. {
  194.     int s, i;
  195.     s = 0;
  196.     for(i=0;i<16;i++) {
  197.         s += abs(pix1[0] - avg2(pix2[0], pix2[1]));
  198.         s += abs(pix1[1] - avg2(pix2[1], pix2[2]));
  199.         s += abs(pix1[2] - avg2(pix2[2], pix2[3]));
  200.         s += abs(pix1[3] - avg2(pix2[3], pix2[4]));
  201.         s += abs(pix1[4] - avg2(pix2[4], pix2[5]));
  202.         s += abs(pix1[5] - avg2(pix2[5], pix2[6]));
  203.         s += abs(pix1[6] - avg2(pix2[6], pix2[7]));
  204.         s += abs(pix1[7] - avg2(pix2[7], pix2[8]));
  205.         s += abs(pix1[8] - avg2(pix2[8], pix2[9]));
  206.         s += abs(pix1[9] - avg2(pix2[9], pix2[10]));
  207.         s += abs(pix1[10] - avg2(pix2[10], pix2[11]));
  208.         s += abs(pix1[11] - avg2(pix2[11], pix2[12]));
  209.         s += abs(pix1[12] - avg2(pix2[12], pix2[13]));
  210.         s += abs(pix1[13] - avg2(pix2[13], pix2[14]));
  211.         s += abs(pix1[14] - avg2(pix2[14], pix2[15]));
  212.         s += abs(pix1[15] - avg2(pix2[15], pix2[16]));
  213.         pix1 += line_size;
  214.         pix2 += line_size;
  215.     }
  216.     return s;
  217. }
  218. static int pix_abs16x16_y2(UINT8 *pix1, UINT8 *pix2, int line_size)
  219. {
  220.     int s, i;
  221.     UINT8 *pix3 = pix2 + line_size;
  222.     s = 0;
  223.     for(i=0;i<16;i++) {
  224.         s += abs(pix1[0] - avg2(pix2[0], pix3[0]));
  225.         s += abs(pix1[1] - avg2(pix2[1], pix3[1]));
  226.         s += abs(pix1[2] - avg2(pix2[2], pix3[2]));
  227.         s += abs(pix1[3] - avg2(pix2[3], pix3[3]));
  228.         s += abs(pix1[4] - avg2(pix2[4], pix3[4]));
  229.         s += abs(pix1[5] - avg2(pix2[5], pix3[5]));
  230.         s += abs(pix1[6] - avg2(pix2[6], pix3[6]));
  231.         s += abs(pix1[7] - avg2(pix2[7], pix3[7]));
  232.         s += abs(pix1[8] - avg2(pix2[8], pix3[8]));
  233.         s += abs(pix1[9] - avg2(pix2[9], pix3[9]));
  234.         s += abs(pix1[10] - avg2(pix2[10], pix3[10]));
  235.         s += abs(pix1[11] - avg2(pix2[11], pix3[11]));
  236.         s += abs(pix1[12] - avg2(pix2[12], pix3[12]));
  237.         s += abs(pix1[13] - avg2(pix2[13], pix3[13]));
  238.         s += abs(pix1[14] - avg2(pix2[14], pix3[14]));
  239.         s += abs(pix1[15] - avg2(pix2[15], pix3[15]));
  240.         pix1 += line_size;
  241.         pix2 += line_size;
  242.     }
  243.     return s;
  244. }
  245. static int pix_abs16x16_xy2(UINT8 *pix1, UINT8 *pix2, int line_size)
  246. {
  247.     int s, i;
  248.     UINT8 *pix3 = pix2 + line_size;
  249.     s = 0;
  250.     for(i=0;i<16;i++) {
  251.         s += abs(pix1[0] - avg4(pix2[0], pix2[1], pix3[0], pix3[1]));
  252.         s += abs(pix1[1] - avg4(pix2[1], pix2[2], pix3[1], pix3[2]));
  253.         s += abs(pix1[2] - avg4(pix2[2], pix2[3], pix3[2], pix3[3]));
  254.         s += abs(pix1[3] - avg4(pix2[3], pix2[4], pix3[3], pix3[4]));
  255.         s += abs(pix1[4] - avg4(pix2[4], pix2[5], pix3[4], pix3[5]));
  256.         s += abs(pix1[5] - avg4(pix2[5], pix2[6], pix3[5], pix3[6]));
  257.         s += abs(pix1[6] - avg4(pix2[6], pix2[7], pix3[6], pix3[7]));
  258.         s += abs(pix1[7] - avg4(pix2[7], pix2[8], pix3[7], pix3[8]));
  259.         s += abs(pix1[8] - avg4(pix2[8], pix2[9], pix3[8], pix3[9]));
  260.         s += abs(pix1[9] - avg4(pix2[9], pix2[10], pix3[9], pix3[10]));
  261.         s += abs(pix1[10] - avg4(pix2[10], pix2[11], pix3[10], pix3[11]));
  262.         s += abs(pix1[11] - avg4(pix2[11], pix2[12], pix3[11], pix3[12]));
  263.         s += abs(pix1[12] - avg4(pix2[12], pix2[13], pix3[12], pix3[13]));
  264.         s += abs(pix1[13] - avg4(pix2[13], pix2[14], pix3[13], pix3[14]));
  265.         s += abs(pix1[14] - avg4(pix2[14], pix2[15], pix3[14], pix3[15]));
  266.         s += abs(pix1[15] - avg4(pix2[15], pix2[16], pix3[15], pix3[16]));
  267.         pix1 += line_size;
  268.         pix2 += line_size;
  269.     }
  270.     return s;
  271. }
  272. static void no_motion_search(MpegEncContext *s, 
  273.                              int *mx_ptr, int *my_ptr)
  274. {
  275.     *mx_ptr = 0;
  276.     *my_ptr = 0;
  277. }
  278. static void full_motion_search(MpegEncContext *s, 
  279.                                int *mx_ptr, int *my_ptr)
  280. {
  281.     int range, x1, y1, x2, y2, xx, yy, x, y;
  282.     int mx, my, mx1, my1, dmin, d;
  283.     UINT8 *pix;
  284.     range = 8 * (1 << (s->f_code - 1));
  285.     if (s->out_format == FMT_H263)
  286.         range = range * 2;
  287.     /* do not allow vectors outside the image (will change one day for
  288.        h263) */
  289.     xx = 16 * s->mb_x;
  290.     yy = 16 * s->mb_y;
  291.     x1 = xx - range + 1; /* we loose one pixel to avoid boundary pb with half pixel pred */
  292.     if (x1 < 0)
  293.         x1 = 0;
  294.     x2 = xx + range - 1;
  295.     if (x2 > (s->width - 16))
  296.         x2 = s->width - 16;
  297.     y1 = yy - range + 1;
  298.     if (y1 < 0)
  299.         y1 = 0;
  300.     y2 = yy + range - 1;
  301.     if (y2 > (s->height - 16))
  302.         y2 = s->height - 16;
  303.     pix = s->new_picture[0] + (yy * s->width) + xx;
  304.     dmin = 0x7fffffff;
  305.     mx = 0;
  306.     my = 0;
  307.     for(y=y1;y<=y2;y++) {
  308.         for(x=x1;x<=x2;x++) {
  309.             d = pix_abs16x16(pix, s->last_picture[0] + (y * s->width) + x, 
  310.                              s->width);
  311.             if (d < dmin ||
  312.                 (d == dmin && 
  313.                  (abs(x - xx) + abs(y - yy)) < 
  314.                  (abs(mx - xx) + abs(my - yy)))) {
  315.                 dmin = d;
  316.                 mx = x;
  317.                 my = y;
  318.             }
  319.         }
  320.     }
  321.     /* half pixel search */
  322.     mx1 = mx = 2 * mx;
  323.     my1 = my = 2 * my;
  324.     if (mx > 0 && mx < 2 * (s->width - 16) &&
  325.         my > 0 && my < 2 * (s->height - 16)) {
  326.         int dx, dy, px, py;
  327.         UINT8 *ptr;
  328.         for(dy=-1;dy<=1;dy++) {
  329.             for(dx=-1;dx<=1;dx++) {
  330.                 if (dx != 0 || dy != 0) {
  331.                     px = mx1 + dx;
  332.                     py = my1 + dy;
  333.                     ptr = s->last_picture[0] + ((py >> 1) * s->width) + (px >> 1);
  334.                     switch(((py & 1) << 1) | (px & 1)) {
  335.                     default:
  336.                     case 0:
  337.                         d = pix_abs16x16(pix, ptr, s->width);
  338.                         break;
  339.                     case 1:
  340.                         d = pix_abs16x16_x2(pix, ptr, s->width);
  341.                         break;
  342.                     case 2:
  343.                         d = pix_abs16x16_y2(pix, ptr, s->width);
  344.                         break;
  345.                     case 3:
  346.                         d = pix_abs16x16_xy2(pix, ptr, s->width);
  347.                         break;
  348.                     }
  349.                     if (d < dmin) {
  350.                         dmin = d;
  351.                         mx = px;
  352.                         my = py;
  353.                     }
  354.                 }
  355.             }
  356.         }
  357.     }
  358.     *mx_ptr = mx - 2 * xx;
  359.     *my_ptr = my - 2 * yy;
  360. #if 0
  361.     if (*mx_ptr < -(2 * range) || *mx_ptr >= (2 * range) ||
  362.         *my_ptr < -(2 * range) || *my_ptr >= (2 * range)) {
  363.         fprintf(stderr, "error %d %dn", *mx_ptr, *my_ptr);
  364.     }
  365. #endif
  366. }
  367. #if 1
  368. int estimate_motion(MpegEncContext *s, 
  369.                     int mb_x, int mb_y,
  370.                     int *mx_ptr, int *my_ptr)
  371. {
  372.     UINT8 *pix, *ppix;
  373.     int sum, varc, vard, mx, my;
  374.     
  375.     if (s->full_search) {
  376.         full_motion_search(s, &mx, &my);
  377.     } else {
  378.         no_motion_search(s, &mx, &my);
  379.     }
  380.     /* intra / predictive decision */
  381.     pix = s->new_picture[0] + (mb_y * 16 * s->width) + mb_x * 16;
  382.     ppix = s->last_picture[0] + 
  383.         ((mb_y * 16 + (my >> 1)) * s->width) + (mb_x * 16 + (mx >> 1));
  384.     sum = pix_sum(pix, s->width);
  385.     varc = pix_norm1(pix, s->width);
  386.     vard = pix_norm(pix, ppix, s->width);
  387.     
  388.     vard = vard >> 8;
  389.     sum = sum >> 8;
  390.     varc = (varc >> 8)  - (sum * sum);
  391. #if 0
  392.     printf("varc=%d (sum=%d) vard=%d mx=%d my=%dn",
  393.            varc, sum, vard, mx, my);
  394. #endif
  395.     if (vard <= 64) {
  396.         *mx_ptr = mx;
  397.         *my_ptr = my;
  398. return 0;
  399.     } else if (vard < varc) {
  400.         *mx_ptr = mx;
  401.         *my_ptr = my;
  402. return 0;
  403.     } else {
  404.         *mx_ptr = 0;
  405.         *my_ptr = 0;
  406.         return 1;
  407.     }
  408. }
  409. #else
  410. /* test version which generates valid random vectors */
  411. int estimate_motion(MpegEncContext *s, 
  412.                     int mb_x, int mb_y,
  413.                     int *mx_ptr, int *my_ptr)
  414. {
  415.      int xx, yy, x1, y1, x2, y2, range;
  416.      if ((random() % 10) >= 5) {
  417.          range = 8 * (1 << (s->f_code - 1));
  418.          if (s->out_format == FMT_H263)
  419.              range = range * 2;
  420.          xx = 16 * s->mb_x;
  421.          yy = 16 * s->mb_y;
  422.          x1 = xx - range;
  423.          if (x1 < 0)
  424.              x1 = 0;
  425.          x2 = xx + range - 1;
  426.          if (x2 > (s->width - 16))
  427.              x2 = s->width - 16;
  428.          y1 = yy - range;
  429.          if (y1 < 0)
  430.              y1 = 0;
  431.          y2 = yy + range - 1;
  432.          if (y2 > (s->height - 16))
  433.              y2 = s->height - 16;
  434.          *mx_ptr = (random()%(2 * (x2 - x1 + 1))) + 2 * (x1 - xx);
  435.          *my_ptr = (random()%(2 * (y2 - y1 + 1))) + 2 * (y1 - yy);
  436.          return 0;
  437.     } else {
  438.         *mx_ptr = 0;
  439.         *my_ptr = 0;
  440.         return 1;
  441.     }
  442. }
  443. #endif