motion_est.c

用MPEG-4对YUV视频文件编码压缩成divx视频文件

/**************************************************************************
 *
 *  Modifications:
 *
 *	01.05.2002	updated MotionEstimationBVOP
 *	25.04.2002 partial prevMB conversion
 *  22.04.2002 remove some compile warning by chenm001 <[e-mail]>
 *  14.04.2002 added MotionEstimationBVOP()
 *  02.04.2002 add EPZS(^2) as ME algorithm, use PMV_USESQUARES to choose between 
 *             EPZS and EPZS^2
 *  08.02.2002 split up PMVfast into three routines: PMVFast, PMVFast_MainLoop
 *             PMVFast_Refine to support multiple searches with different start points
 *  07.01.2002 uv-block-based interpolation
 *  06.01.2002 INTER/INTRA-decision is now done before any SEARCH8 (speedup)
 *             changed INTER_BIAS to 150 (as suggested by suxen_drol)
 *             removed halfpel refinement step in PMVfastSearch8 + quality=5
 *             added new quality mode = 6 which performs halfpel refinement
 *             filesize difference between quality 5 and 6 is smaller than 1%
 *             (Isibaar)
 *  31.12.2001 PMVfastSearch16 and PMVfastSearch8 (gruel)
 *  30.12.2001 get_range/MotionSearchX simplified; blue/green bug fix
 *  22.12.2001 commented best_point==99 check
 *  19.12.2001 modified get_range (purple bug fix)
 *  15.12.2001 moved pmv displacement from mbprediction
 *  02.12.2001 motion estimation/compensation split (Isibaar)
 *  16.11.2001 rewrote/tweaked search algorithms; [e-mail]
 *  10.11.2001 support for sad16/sad8 functions
 *  28.08.2001 reactivated MODE_INTER4V for EXT_MODE
 *  24.08.2001 removed MODE_INTER4V_Q, disabled MODE_INTER4V for EXT_MODE
 *  22.08.2001 added MODE_INTER4V_Q			
 *  20.08.2001 added pragma to get rid of internal compiler error with VC6
 *             idea by Cyril. Thanks.
 *
 *  Michael Militzer <[e-mail]>
 *
 **************************************************************************/

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "../user_macro.h"
#include "../encoder.h"
#include "../utils/mbfunctions.h"
/*#include "../prediction/mbprediction.h"*/
#include "../global.h"
#include "../utils/timer.h"
#include "motion.h"
#include "sad.h"

/*******************************zjy********************************/
extern FILE *ftestlog;
/******************************************************************/

#define HalfPixel_Eight

/* very large value */
#define MV_MAX_ERROR	(4096 * 256)


/* sad16(0,0) bias; mpeg4 spec suggests nb/2+1 */
/* nb  = vop pixels * 2^(bpp-8) */
#define MV16_00_BIAS	(128+1)
#define MV8_00_BIAS	(0)

/* INTER bias for INTER/INTRA decision; mpeg4 spec suggests 2*nb */
#define MV16_INTER_BIAS	512

/* Parameters which control inter/inter4v decision */
#define IMV16X16			5

/* vector map (vlc delta size) smoother parameters */
#define NEIGH_TEND_16X16	2
#define NEIGH_TEND_8X8		2

/* fast ((A)/2)*2  */
#define EVEN(A)		(((A)<0?(A)+1:(A)) & ~1)

#define MVzero(A) ( ((A).x)==(0) && ((A).y)==(0) )
#define MVequal(A,B) ( ((A).x)==((B).x) && ((A).y)==((B).y) )

static int amvfastSearch16(
				   const uint8_t * const pRef,
                   const uint8_t * const pRefH,
				   const uint8_t * const pRefV,
                   const uint8_t * const pRefHV,
                   const IMAGE * const pCur,
				   const int x,
				   const int y,
				   const uint32_t MotionFlags,
				   const int32_t iQuant,
				   const uint32_t iFcode,
				   const MBParam * const pParam,
				   const MACROBLOCK * const pMBs,
				   const MACROBLOCK * const prevMBs,
				   VECTOR * const currMV,
				   VECTOR * const currPMV);
static int
amvfastSearch8(const uint8_t * const pRef,
			   const uint8_t * const pRefH,
			   const uint8_t * const pRefV,
			   const uint8_t * const pRefHV,
			   const IMAGE * const pCur,
			   const int x,
			   const int y,
			   const int start_x,
			   const int start_y,
			   const uint32_t MotionFlags,
			   const int32_t iQuant,
			   const uint32_t iFcode,
			   const MBParam * const pParam,
			   const MACROBLOCK * const pMBs,
			   const MACROBLOCK * const prevMBs,
			   VECTOR * const currMV,
			   VECTOR * const currPMV);


static int32_t lambda_vec16[32] =	/* rounded values for lambda param for weight of motion bits as in modified H.26L */
{ 0, (int) (1.00235 + 0.5), (int) (1.15582 + 0.5), (int) (1.31976 + 0.5),
		(int) (1.49591 + 0.5), (int) (1.68601 + 0.5),
	(int) (1.89187 + 0.5), (int) (2.11542 + 0.5), (int) (2.35878 + 0.5),
		(int) (2.62429 + 0.5), (int) (2.91455 + 0.5),
	(int) (3.23253 + 0.5), (int) (3.58158 + 0.5), (int) (3.96555 + 0.5),
		(int) (4.38887 + 0.5), (int) (4.85673 + 0.5),
	(int) (5.37519 + 0.5), (int) (5.95144 + 0.5), (int) (6.59408 + 0.5),
		(int) (7.31349 + 0.5), (int) (8.12242 + 0.5),
	(int) (9.03669 + 0.5), (int) (10.0763 + 0.5), (int) (11.2669 + 0.5),
		(int) (12.6426 + 0.5), (int) (14.2493 + 0.5),
	(int) (16.1512 + 0.5), (int) (18.442 + 0.5), (int) (21.2656 + 0.5),
		(int) (24.8580 + 0.5), (int) (29.6436 + 0.5),
	(int) (36.4949 + 0.5)
};

static int32_t *lambda_vec8 = lambda_vec16;	/* same table for INTER and INTER4V for now */



/* mv.length table  */
static const uint32_t mvtab[33] = {
	1, 2, 3, 4, 6, 7, 7, 7,
	9, 9, 9, 10, 10, 10, 10, 10,
	10, 10, 10, 10, 10, 10, 10, 10,
	10, 11, 11, 11, 11, 11, 11, 12, 12
};


static __inline uint32_t
mv_bits(int32_t component,
		const uint32_t iFcode)
{
	if (component == 0)
		return 1;

	if (component < 0)
		component = -component;

	if (iFcode == 1) {
		if (component > 32)
			component = 32;

		return mvtab[component] + 1;
	}

	component += (1 << (iFcode - 1)) - 1;
	component >>= (iFcode - 1);

	if (component > 32)
		component = 32;

	return mvtab[component] + 1 + iFcode - 1;
}

/* replace by old  because of iQuant???? lxq 2002-7-12     */
static __inline uint32_t
calc_delta_16(const int32_t dx,
			  const int32_t dy,
			  const uint32_t iFcode,
			  const uint32_t iQuant)
{
	return NEIGH_TEND_16X16 * lambda_vec16[iQuant] * (mv_bits(dx, iFcode) +
													  mv_bits(dy, iFcode));
}

static __inline uint32_t
calc_delta_8(const int32_t dx,
			 const int32_t dy,
			 const uint32_t iFcode,
			 const uint32_t iQuant)
{
	return NEIGH_TEND_8X8 * lambda_vec8[iQuant] * (mv_bits(dx, iFcode) +
												   mv_bits(dy, iFcode));
}
/*
static __inline uint32_t calc_delta_16(const int32_t dx, const int32_t dy, const uint32_t iFcode,const uint32_t iQuant)
{
	return NEIGH_TEND_16X16 * (mv_bits(dx, iFcode) + mv_bits(dy, iFcode));
}

static __inline uint32_t calc_delta_8(const int32_t dx, const int32_t dy, const uint32_t iFcode,const uint32_t iQuant)

{
    return NEIGH_TEND_8X8 * (mv_bits(dx, iFcode) + mv_bits(dy, iFcode));
}
*/



#ifndef SEARCH16
#define SEARCH16    amvfastSearch16
//#define SEARCH16    fullSearch16
#endif

#ifndef SEARCH8
#define SEARCH8     amvfastSearch8
#endif
/*
bool
MotionEstimation(MBParam * const pParam,
				 FRAMEINFO * const current,
				 FRAMEINFO * const reference,
				 const IMAGE * const pRefH,
				 const IMAGE * const pRefV,
				 const IMAGE * const pRefHV,int * mad,
				 const uint32_t iLimit)
*/

void 
MotionEstimation(MBParam * const pParam,
				 FRAMEINFO * const current,
				 FRAMEINFO * const reference,
				 const IMAGE * const pRefH,
				 const IMAGE * const pRefV,
				 const IMAGE * const pRefHV,
				 int * mad)
{
	const int32_t iWcount = pParam->mb_width;
	const int32_t iHcount = pParam->mb_height;
	MACROBLOCK *const pMBs = current->mbs;
	MACROBLOCK *const prevMBs = reference->mbs;
	const IMAGE *const pCurrent = &current->image;
	const IMAGE *const pRef = &reference->image;

	const VECTOR zeroMV = { 0, 0 };

	int32_t x, y;
	/*int32_t iIntra = 0;*/
	VECTOR pmv;

	/* fyh add code */
    MACROBLOCK *pMB;/* pointer to current macroblock struct. Add by fyh */
	int32_t deviation;
	int difference;
	int sad8;
	*mad=0;
	
	/* consider delete sadInit. Add by fyh */
/*	if (sadInit)
		(*sadInit) ();*/

	for (y = 0; y < iHcount; y++){
/*	//	if(y==1)
	//		printf("this is a test!");   */
		for (x = 0; x < iWcount; x++) {
			/* fyh comment code */
			/*MACROBLOCK *const pMB = &pMBs[x + y * iWcount];*/
            pMB = &pMBs[x + y * iWcount];

			pMB->sad16 =
				SEARCH16(pRef->y, pRefH->y, pRefV->y, pRefHV->y, pCurrent, x,
						 y, current->motion_flags, current->quant,
						 current->fcode, pParam, pMBs, prevMBs, &pMB->mv16,
						 &pMB->pmvs[0]);

			/*if (0 < (pMB->sad16 - MV16_INTER_BIAS)) {*/
			difference=pMB->sad16 - MV16_INTER_BIAS;
			if(difference) {

                /* fyh comment code */
				/*int32_t deviation;*/

				deviation =
					dev16(pCurrent->y + x * 16 + y * 16 * pParam->edged_width,
						  pParam->edged_width);

				/*if (deviation < (pMB->sad16 - MV16_INTER_BIAS)) {*/
                if (deviation < difference) {
					pMB->mode = MODE_INTRA;

					pMB->mv16 = pMB->mvs[0] = pMB->mvs[1] = pMB->mvs[2] =
						pMB->mvs[3] = zeroMV;
					pMB->sad16 = pMB->sad8[0] = pMB->sad8[1] = pMB->sad8[2] =
						pMB->sad8[3] = 0;
					*mad+=deviation; /* add by lxq   */

					/* iIntra del */
					/*iIntra++;
					if (iIntra >= iLimit)
						return 1;*/
					/* iIntra del */

					continue;
				}
			}

			pmv = pMB->pmvs[0];
			if (current->global_flags & XVID_INTER4V)
/*				if ((!(current->global_flags & XVID_LUMIMASKING) ||
					 pMB->dquant == NO_CHANGE)) */
					 {
     				/*int32_t sad8 = IMV16X16 * current->quant;*/
				    sad8 = IMV16X16 * current->quant;/* ???? fyh add */

					if (sad8 < pMB->sad16)

						sad8 += pMB->sad8[0] =
							SEARCH8(pRef->y, pRefH->y, pRefV->y, pRefHV->y,
									pCurrent, 2 * x, 2 * y, pMB->mv16.x,
									pMB->mv16.y, current->motion_flags,
									current->quant, current->fcode, pParam,
									pMBs, prevMBs, &pMB->mvs[0],
									&pMB->pmvs[0]);

					if (sad8 < pMB->sad16)
						sad8 += pMB->sad8[1] =
							SEARCH8(pRef->y, pRefH->y, pRefV->y, pRefHV->y,
									pCurrent, 2 * x + 1, 2 * y, pMB->mv16.x,
									pMB->mv16.y, current->motion_flags,
									current->quant, current->fcode, pParam,
									pMBs, prevMBs, &pMB->mvs[1],
									&pMB->pmvs[1]);

					if (sad8 < pMB->sad16)
						sad8 += pMB->sad8[2] =
							SEARCH8(pRef->y, pRefH->y, pRefV->y, pRefHV->y,
									pCurrent, 2 * x, 2 * y + 1, pMB->mv16.x,
									pMB->mv16.y, current->motion_flags,
									current->quant, current->fcode, pParam,
									pMBs, prevMBs, &pMB->mvs[2],
									&pMB->pmvs[2]);

					if (sad8 < pMB->sad16)
						sad8 += pMB->sad8[3] =
							SEARCH8(pRef->y, pRefH->y, pRefV->y, pRefHV->y,
									pCurrent, 2 * x + 1, 2 * y + 1,
									pMB->mv16.x, pMB->mv16.y,
									current->motion_flags, current->quant,
									current->fcode, pParam, pMBs, prevMBs,
									&pMB->mvs[3], &pMB->pmvs[3]);

					/* decide: MODE_INTER or MODE_INTER4V
					   mpeg4:   if (sad8 < pMB->sad16 - nb/2+1) use_inter4v
					 */

					if (sad8 < pMB->sad16) {
						pMB->mode = MODE_INTER4V;
						*mad+=pMB->sad8[0]+pMB->sad8[1]+pMB->sad8[2]+pMB->sad8[3];/*need to be test add by lxq*/
						
						/* 2002.2.26 fyh comment code */
						/*pMB->sad8[0] *= 4;
						pMB->sad8[1] *= 4;
						pMB->sad8[2] *= 4;
						pMB->sad8[3] *= 4;*/
						
						continue;
					}

				}

			pMB->mode = MODE_INTER;
			pMB->pmvs[0] = pmv;	/* pMB->pmvs[1] = pMB->pmvs[2] = pMB->pmvs[3]  are not needed for INTER */
			pMB->mvs[0] = pMB->mvs[1] = pMB->mvs[2] = pMB->mvs[3] = pMB->mv16;
			/*pMB->sad8[0] = pMB->sad8[1] = pMB->sad8[2] = pMB->sad8[3] =
				pMB->sad16;*/
			pMB->sad8[0] = pMB->sad8[1] = pMB->sad8[2] = pMB->sad8[3] =
				pMB->sad16>>2;
			*mad+=pMB->sad16; /*add by lxq*/
		}
		}
		*mad/=iHcount*iWcount;/*add by lxq  */
	/*return 0;*/
}

#define CHECK_MV16_CANDIDATE(X,Y) { \
  if ( ((X) <= max_dx) && ((X) >= min_dx) \
    && ((Y) <= max_dy) && ((Y) >= min_dy) ) \
  { \
    iSAD = sad16( cur, get_ref(pRef, pRefH, pRefV, pRefHV, x, y, 16, X, Y, iEdgedWidth),iEdgedWidth, iMinSAD); \
    iSAD += calc_delta_16((X) - pmv[0].x, (Y) - pmv[0].y, (uint8_t)iFcode, iQuant);\
    if (iSAD < iMinSAD) \
    {  iMinSAD=iSAD; currMV->x=(X); currMV->y=(Y); } } \
}

int32_t
Halfpel16_Refine(const uint8_t * const pRef,
				 const uint8_t * const pRefH,
				 const uint8_t * const pRefV,
				 const uint8_t * const pRefHV,
				 const uint8_t * const cur,
				 const int x,
				 const int y,
				 VECTOR * const currMV,
				 int32_t iMinSAD,
				 const VECTOR * const pmv,
				 const int32_t min_dx,
				 const int32_t max_dx,
				 const int32_t min_dy,
				 const int32_t max_dy,
				 const int32_t iFcode,
				 const int32_t iQuant,
				 const int32_t iEdgedWidth)
{
/* Do a half-pel refinement (or rather a "smallest possible amount" refinement) */

	int32_t iSAD;
	VECTOR backupMV = *currMV;

    CHECK_MV16_CANDIDATE(backupMV.x - 1, backupMV.y - 1);
    CHECK_MV16_CANDIDATE(backupMV.x, backupMV.y - 1);
    CHECK_MV16_CANDIDATE(backupMV.x + 1, backupMV.y - 1);
    CHECK_MV16_CANDIDATE(backupMV.x - 1, backupMV.y);
    CHECK_MV16_CANDIDATE(backupMV.x + 1, backupMV.y);
    CHECK_MV16_CANDIDATE(backupMV.x - 1, backupMV.y + 1);
    CHECK_MV16_CANDIDATE(backupMV.x, backupMV.y + 1);
    CHECK_MV16_CANDIDATE(backupMV.x + 1, backupMV.y + 1);

	return iMinSAD;
}

#define CHECK_MV8_CANDIDATE(X,Y) { \
  if ( ((X) <= max_dx) && ((X) >= min_dx) \
    && ((Y) <= max_dy) && ((Y) >= min_dy) ) \
  { \
    iSAD = sad8( cur, get_ref(pRef, pRefH, pRefV, pRefHV, x, y, 8, (X), (Y), iEdgedWidth),iEdgedWidth); \
    iSAD += calc_delta_8((X)-pmv[0].x, (Y)-pmv[0].y, (uint8_t)iFcode, iQuant);\
    if (iSAD < iMinSAD) \
    {  iMinSAD=iSAD; currMV->x=(X); currMV->y=(Y); } } \
}

int32_t
Halfpel8_Refine(const uint8_t * const pRef,
				const uint8_t * const pRefH,
				const uint8_t * const pRefV,
				const uint8_t * const pRefHV,
				const uint8_t * const cur,
				const int x,
				const int y,
				VECTOR * const currMV,
				int32_t iMinSAD,
				const VECTOR * const pmv,
				const int32_t min_dx,
				const int32_t max_dx,
				const int32_t min_dy,
				const int32_t max_dy,
				const int32_t iFcode,
				const int32_t iQuant,
				const int32_t iEdgedWidth)
{
/* Do a half-pel refinement (or rather a "smallest possible amount" refinement) */

	int32_t iSAD;
	VECTOR backupMV = *currMV;

	CHECK_MV8_CANDIDATE(backupMV.x - 1, backupMV.y - 1);
	CHECK_MV8_CANDIDATE(backupMV.x, backupMV.y - 1);
	CHECK_MV8_CANDIDATE(backupMV.x + 1, backupMV.y - 1);
	CHECK_MV8_CANDIDATE(backupMV.x - 1, backupMV.y);
	CHECK_MV8_CANDIDATE(backupMV.x + 1, backupMV.y);
	CHECK_MV8_CANDIDATE(backupMV.x - 1, backupMV.y + 1);
	CHECK_MV8_CANDIDATE(backupMV.x, backupMV.y + 1);
	CHECK_MV8_CANDIDATE(backupMV.x + 1, backupMV.y + 1);

	return iMinSAD;
}
/* This is somehow a copy of get_pmv, but returning all MVs instead of only Median MV */

static __inline int/* 1:means left\top\righttop block motion vector eauql */ 
get_pmvdata_amv(const MACROBLOCK * const pMBs,/* pointer to curent frame marcblok struct */
			const uint32_t x,/* macroblock coordinate */ 
			const uint32_t y,
			const uint32_t x_dim,/* image width by 16 pixel */
			const uint32_t block,/* block index in a macroblock */
			VECTOR * const pmv/* store predict motion vectors */
            )
{
	/*
	 * pmv are filled with: 
	 *  [0]: Median (or whatever is correct in a special case)
	 *  [1]: left neighbour
	 *  [2]: top neighbour
	 *  [3]: topright neighbour
	 */

	int xin1, xin2, xin3;
	int yin1, yin2, yin3;
	int vec1, vec2, vec3;

	uint32_t index = x + y * x_dim;
	const VECTOR zeroMV = { 0, 0 };

	/* first row of blocks (special case)*/
	if (y == 0 && (block == 0 || block == 1)) {
		if ((x == 0) && (block == 0)) 	/* first column, first block*/
		{
			pmv[0] = pmv[1] = pmv[2] = pmv[3] = zeroMV;
			return 0;
		}
		if (block == 1)			/* second block; has only a left neighbour*/
		{
			pmv[0] = pmv[1] = pMBs[index].mvs[0];
			pmv[2] = pmv[3] = zeroMV;
			return 0;
		} else {				/* block==0, but x!=0, so again, there is a left neighbour */

			pmv[0] = pmv[1] = pMBs[index - 1].mvs[1];
			pmv[2] = pmv[3] = zeroMV;
			return 0;
		}
	}

	/*
	 * MODE_INTER, vm18 page 48
	 * MODE_INTER4V vm18 page 51
	 *
	 *  (x,y-1)      (x+1,y-1)
	 *  [   |   ]    [   |   ]
	 *  [ 2 | 3 ]    [ 2 |   ]
	 *
	 *  (x-1,y)      (x,y)        (x+1,y)
	 *  [   | 1 ]    [ 0 | 1 ]    [ 0 |   ]
	 *  [   | 3 ]    [ 2 | 3 ]    [   |   ]
	 */

	switch (block) {
	case 0:
		xin1 = x - 1;
		yin1 = y;
		vec1 = 1;				/* left */
		xin2 = x;
		yin2 = y - 1;
		vec2 = 2;				/* top */
		xin3 = x + 1;
		yin3 = y - 1;
		vec3 = 2;				/* top right */
		break;
	case 1:
		xin1 = x;
		yin1 = y;
		vec1 = 0;
		xin2 = x;
		yin2 = y - 1;
		vec2 = 3;
		xin3 = x + 1;
		yin3 = y - 1;
		vec3 = 2;
		break;
	case 2: