motionES.cpp

//   MBSIZE : 簛E？榈拇笮�
//   P : 搜索步长参数

#define MBSIZE 16
#define P 9
#include "motionES.h"
#include <cmath>
#include "string.h"
//#include <stdafx.h>

const struct Point{int x;int y;};
const Point LDSP[9]={{0,-2},{-1,-1},{1,-1},{-2,0},{0,0},{2,0},{-1,1},{1,1},{0,2}};
const Point SDSP[5]={{0,-1},{-1,0},{0,0},{1,0},{0,1}};

//********************************************************************************************
// TSS算法——P越大，搜索窗口越大，计算量越大，MBSIZE愈大，搜索的运动矢量少，计算量小，效果阐靠
//
// 输葋E
//   imgI : 参考图蟻E
//   imgP : 要寻找运动矢量的图蟻E
//   imgComp : 补偿后图像的数据区的指諄E
//   imgWidth : 图像的窥胰
//   imgHeight: 图像的窥胰
// 输硜E
//   motionvector: 图像估计的运动矢量矩諄E
//   TSScomputations: TSS算法的计算量
//*********************************************************************************************
void motionEstTSS(unsigned char *imgP, unsigned char* imgI, int imgWidth, int imgHeight, int motionVect[2][MAX], double* TSScomputations)
{
//初始化参数	

     double costs[3][3]={{65537,65537,65537},{65537,65537,65537},{65537,65537,65537}};
     double computations = 0;
     int L = 0; 
     int stepMax = 0;
     int vectorsNum = 0;
	 // 计算所需要的畜步数
	 L= (int)floor(log10(P+1.0)/log10(2.0));  //p=7,s=4
	 stepMax = (int)pow(2.0,(L-1));
	 vectorsNum =(int)(imgWidth*imgHeight/pow((double)MBSIZE,2));
//动态分艁E
	 //初始化运动矢量为羴E
     int vectors[2][MAX] ={0};//int vectors[2][vectorsNum] ;
// 从图像左上角开始，说明像块的中心点在左上角处开始
	 //步长为MBSIZE，我们寻找的每个簛E？樗阉髌渖舷伦笥襭像素处的匹艁E
     int mbCount = 1;
     int i,j;
     for( i = 0; i< imgHeight;i+=MBSIZE)//防止溢硜E-图像必衼E荕BSIZE^2整数倍?周边的模块不去运藖E
	 {
        for(j = 0; j<imgWidth;j+=MBSIZE)
		{   
        //对于三步法，每一步我们评估周围的9部分找出猊配的
            double  pVector[3]={0,0,0};//----------[dx,dy,min]
            int  x = j;//图像的列
            int  y = i;//图像的行
			int  k,t;
			int  stepSize;
        // 为了避免对搜索中心的反复计算，我们存储其前一次运行得到的数值。
		//对于第一次祦E颐前汛酥荡娣诺窖分猓婧蟮牡丒颐前阎荡娣诺揭贫沟男轮行牡愦�

			//动态分配二维数讈E
			unsigned char** imgPtmp;
			unsigned char** imgItmp;
			imgPtmp = new unsigned char*[MBSIZE];
			imgItmp = new unsigned char*[MBSIZE];
			for (k=0;k<MBSIZE;k++)
			{		
				imgPtmp[k]=new unsigned char[MBSIZE];
			    memset(imgPtmp[k],0,MBSIZE*sizeof(unsigned char));//新建内存区莵E丒
			}

			for (k=0;k<MBSIZE;k++)
			{			
				imgItmp[k]=new unsigned char[MBSIZE];
                memset(imgItmp[k],0,MBSIZE*sizeof(unsigned char));//新建内存区莵E丒
			}

            for(int blkX=0;blkX<MBSIZE;blkX++)
			{
	    		for(int blkY=0;blkY<MBSIZE;blkY++)
				{
					int MAX1=0,MAX2=0;//防溢硜E
					//为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					MAX1 =((i+blkY)*imgWidth+(j+blkX)<=(imgWidth*imgHeight-1))?(i+blkY)*imgWidth+(j+blkX):(imgWidth*imgHeight-1); 
                    *(imgPtmp[blkX]+blkY) = *(imgP+MAX1);
		    	    *(imgItmp[blkX]+blkY) = *(imgI+MAX1);
				}
			}
			   //     costs(2,2) = costFuncMAD(imgP(i:i+MBSIZE-1,j:j+MBSIZE-1), ...
                             //       imgI(i:i+MBSIZE-1,j:j+MBSIZE-1),MBSIZE);
			//强制转换开销值为整数
			//循环前先计算参考中心处开蟻E
            costs[1][1] = costFuncMAD(imgPtmp,imgItmp,MBSIZE);
            computations = computations + 1;
            stepSize = stepMax;               

            while(stepSize >= 1)  
			{
                int m,n;
				int dx,dy;
	     		double min;
            // m 是行或者垂直搜索参数
            // n 是列或者水平搜索参数
                for( m = -stepSize;m <= stepSize ; m += stepSize)
				{
                   for( n = -stepSize;n <= stepSize ; n += stepSize)
				   {
				    	int refBlkVer,refBlkHor;
				    	int costRow,costCol;
                        refBlkVer = y + m;   // 参考块的行，垂直坐眮E
                        refBlkHor = x + n;   // 参考縼E牧校阶丒
               //         if (( refBlkVer < 1 )||( refBlkVer+MBSIZE-1 > imgHeight)|| (refBlkHor < 1 )|| (refBlkHor+MBSIZE-1)> imgWidth)
						if (( refBlkVer < 0 )||( refBlkVer+MBSIZE > imgHeight)|| (refBlkHor < 0 )|| (refBlkHor+MBSIZE)> imgWidth)
                           continue;
						//生成cost矩阵的行列坐眮E阉飨卤丒
                        costRow = m/stepSize + 1;//0~2  
                        costCol = n/stepSize + 1;//0~2
						//中心处跳硜E
                        if (costRow == 1 && costCol == 1)
                           continue;
				        for(int blkX2=0;blkX2<MBSIZE;blkX2++)
						{
	                     	for(int blkY2=0;blkY2<MBSIZE;blkY2++)
							{
							   int MAX1=0,MAX2=0;//防溢硜E
							//为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
							   MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				               MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                               *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                   *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
							}
						}
                        costs[costRow][costCol] = costFuncMAD(imgPtmp,imgItmp,MBSIZE);
     //                 costs(costRow, costCol ) = costFuncMAD(imgP(i:i+MBSIZE-1,j:j+MBSIZE-1), ...
     //                   imgI(refBlkVer:refBlkVer+MBSIZE-1, refBlkHor:refBlkHor+MBSIZE-1), MBSIZE);
                        computations = computations + 1;//运算量加一
				   }
				}
            
            // 到此我们找到臁开销并且存储该值
                minCost(costs,pVector);      // finds which macroblock in imgI gave us min Cost
                dx  = (int)(pVector[0]);//dx.dy都是0~2不是1~3
		     	dy  = (int)(pVector[1]);
			    min = pVector[2];
               // 将搜索窗口坐眮E丒碌阶】�
                x = x + (dx-1)*stepSize;
                y = y + (dy-1)*stepSize; 
            // Arohs thought: At this point we can check and see if the shifted co-ordinates are exactly the same 
			//as the root co-ordinates of the last step, then we check them against a preset threshold, 
			//and ifthe cost is less then that, than we can exit from teh loop right here. 
			//This way we can save more computations. However, as this is not implemented in the paper I am modeling, 
			//am not incorporating this test. May be later...as my own addition to the algorithm
            //  搜索步长减皝E
                stepSize = stepSize / 2;
			//  参考中心的开销值竵E挛壳白ヅ浯Φ�
                costs[1][1] = costs[dy][dx]; 
				
			}


			//判断是否写溢硜E
			if(mbCount>vectorsNum)
				mbCount = vectorsNum;
            vectors[0][mbCount-1] = y - i;    // 运动矢量的行坐眮E
            vectors[1][mbCount-1] = x - j;    // 运动矢量的列坐眮E
			//模块计数加一
            mbCount = mbCount + 1;
       //     costs[3][3]={{65537,65537,65537},{65537,65537,65537},{65537,65537,65537}};
			//将开销函数复原---有其他简单方法么？
			for(k=0;k<3;k++)  //memset(costs,65537,3*sizeof(int));//memset按字符操作导致65537截断成01,int 占2字节，产蓙E010101---16843009，结果不对！！！
			{
				for(t=0;t<3;t++)
					costs[k][t]=65537;
			}	 			//释放内磥E
     	    for (k=0;k<MBSIZE;k++)
            	delete []imgPtmp[k];
        	for (k=0;k<MBSIZE;k++)
            	delete []imgItmp[k];
	        delete []imgPtmp;
	        delete []imgItmp;
		}
	 }
 	 for(i=0;i<2;i++) 
		{
			for(j=0;j<MAX;j++)
			{
				if(j>=vectorsNum)
					motionVect[i][j]=0;
				else
				    motionVect[i][j]=vectors[i][j];
			}

		}
     *TSScomputations = computations/(mbCount - 1);    

}

//*******************************************************************************************
//   计算补偿后图像的PSNR
// 输葋E
//   imgP: 原图
//   imgComp:补偿后图蟻E
//   imgHeight: 图像的高
//   imgWidth: 图像的縼E
//   n : 图像像素的赃峰峰值
// 输硜E
//   psnr:补偿图像的峰值信詠E�
//******************************************************************************************
double psnr(unsigned char* imgP,unsigned char* imgComp,int imgWidth,int imgHeight,int n)
{

    double err = 0;
    double mse,psnr;
	int i,j;
    for (i=0;i<imgHeight;i++)
	{
     	for(j=0;j<imgWidth;j++)
		{
     	  err = err + pow(double(*(imgP+i*imgWidth+j) - *(imgComp+i*imgWidth+j)),2);
		}
	}

    mse = err / (imgHeight*imgWidth);
    psnr = 10*log10(n*n/mse);
    return psnr;
	
}

//****************************************************************************
// 对两个块计算平均绝对误瞾E
// 输葋E
// currentBlk : 计算MAD值的像縼E
// refBlk : 计算MAD值的参考像縼E
// 输硜E
// cost : 两像块MAD值
//******************************************************************************
double costFuncMAD(unsigned char** currentBlk,unsigned char** refBlk,int mbSize)
{
	
	double err = 0;
	double cost;
	int i,j;
    for(i = 0;i<mbSize;i++)
	{
		for(j = 0;j<mbSize;j++)
		{
		  err = err + labs(*(currentBlk[i]+j) - *(refBlk[i]+j));
		}
	}
    cost = err / (mbSize*mbSize);
	return cost;
	
}

//****************************************************************************
// 找到臁开销矩阵的臁值的下眮E
//
// 输葋E
//   costs : 簛E？榭卣丒
//         
// 输硜E
//   pVector: dx,dy,min输出臁块的列，行坐眮E妥≈�
//****************************************************************************
void minCost(double costs[3][3],double pVector[] )
{
//初始化
  double min = *(costs[1]+1);
  double dx=1,dy=1;
  int i,j;
  for(i =0; i< 3;i++) 
  {
	 for(j =0; j< 3;j++) 
	 {     
		 if (*(costs[i]+j) < min)
		 {
		    min = *(costs[i]+j);
            dx = j; dy = i;
		 }
	 }
  }
  pVector[0]=dx;
  pVector[1]=dy;
  pVector[2]=min;

}

//*****************************************************************************
// 根据运动估计矢量生成运动补偿图蟻E
//
// 输葋E
//   imgI : 参考图蟻E
//   motionVect : 估计运动矢量矩諄E
//   MBSIZE ： 簛E？榇笮�
///  imgWidth : 图像窥胰
//   imgHeight: 图像的高度
// 输硜E
//   imgComp :补偿图像数据指諄E
//*****************************************************************************
void motionComp(unsigned char* imgI,int mbSize,int imgWidth,int imgHeight,int motionVect[2][MAX],unsigned char*imgComp)
{

   int mbCount = 1;
   int i,j;
//指示是否输葋E枷袷荕BSIZE的整数倍
   int xExtra,yExtra;
   for(i = 0; i<imgHeight; i+= mbSize) //运动补偿图像的产生需要循环所有像素
   { 
	   for(j = 0; j<imgWidth; j+=mbSize) 
	   {
		// dy 行 垂直坐眮E
        // dx 列 水平坐眮E
		   xExtra = imgWidth-j;//残阐楷煮一次等于0则整数倍  图像的列
		   yExtra = imgHeight-i;

		   int dy,dx,refBlkVer,refBlkHor;
           dy = motionVect[0][mbCount-1];//读硜E858993460
           dx = motionVect[1][mbCount-1];
           refBlkVer = i + dy;
           refBlkHor = j + dx;
		   //像块运动补偿，中心点在像块左上角
		   for(int m=0;m<mbSize;m++)
		   {
		     for(int n=0;n<mbSize;n++)
			 {
				 int MAX1 = 0,MAX2 = 0;
				 //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
				 MAX1 = ((refBlkVer+m)*imgWidth+(refBlkHor+n)<=(imgWidth*imgHeight-1))?((refBlkVer+m)*imgWidth+(refBlkHor+n)):(imgWidth*imgHeight-1);
				 MAX2 = (((i+m)*imgWidth+(n+j))<=(imgWidth*imgHeight-1))?((i+m)*imgWidth+(n+j)):(imgWidth*imgHeight-1);
				 imgComp[MAX2]= imgI[MAX1];
			//	imageComp(i:i+MBSIZE-1,j:j+MBSIZE-1) = imgI(refBlkVer:refBlkVer+MBSIZE-1, refBlkHor:refBlkHor+MBSIZE-1);
			 }
		   }
           mbCount = mbCount + 1;
	   }
    }
    if(yExtra!=0)//非整数倍
		if(xExtra==0)
	     	for(i=imgHeight-1;i>=imgHeight-1-yExtra;i--)
		    	for(j=0;j<imgWidth;j++)
				{				
			    	 imgComp[i*imgWidth+j]= imgI[i*imgWidth+j];
				}
		else
		{
			for(i=imgHeight-1;i>=imgHeight-1-yExtra;i--)
		    	for(j=0;j<imgWidth;j++)
				{				
			    	 imgComp[i*imgWidth+j]= imgI[i*imgWidth+j];
				}
			for(i=0;i<imgHeight-yExtra;i++)
		    	for(j=imgWidth-1;j>=imgWidth-1-xExtra;j--)
				{				
			    	 imgComp[i*imgWidth+j]= imgI[i*imgWidth+j];
				}
		}
	else
	;


}

//***************************************************************************************************/
//计算一维矩阵的臁值和其所在的下眮E
//输葋E
//matrix:N个元素矩諄E
//N: 矩阵长度，为奇数
//minNum:矩阵中的臁值
//index:矩阵臁值的索引
//***************************************************************************************************
void minSearch(double* matrix,int N,double* minNum,int* index)
{
			//初始臁在中紒E
	int i,k;
	double value;
	k = (N-1)/2;
	value = *(matrix+k);
	for(i=0;i<N;i++)
	{
		if(*(matrix+i)<*(matrix+k))
		{
			k = i;
			value = *(matrix+i);
		}

	}
	*index = k;
	*minNum =value;
}

//********************************************************************************************
// DS运动估计算法---蛘匹配块设置在中心点为中心，P为皝E兜姆娇蚰冢孀臥增大计算量增大，效果变好，MBSIZE变小计算量增大，效果变好？
//
// 输葋E
//   imgI : 参考图蟻E
//   imgP : 要寻找运动矢量的图蟻E
//   imgComp : 补偿后图像的数据区的指諄E
//   imgWidth : 图像的窥胰
//   imgHeight: 图像的窥胰
// 输硜E
//   motionvector: 图像估计的运动矢量矩諄E
//   DScomputations: DS算法的计算量
//*********************************************************************************************

void motionEstDS(unsigned char *imgP, unsigned char* imgI, int imgWidth, int imgHeight, int motionVect[2][MAX], double* DScomputations)
{

//初始化参数	
     int vectorsNum = 0;
     double costs[9];
     double computations = 0;
     int L = 0; 
// 计算所需要的畜步数
	 L= (int)floor(log10(P+1.0)/log10(2.0));  //p=7
	 vectorsNum =(int)(imgWidth*imgHeight/pow((double)MBSIZE,2));
//动态分艁E
	 //初始化运动矢量为羴E
     int vectors[2][MAX] ={0};//int vectors[2][vectorsNum] ;
// 从图像左上角开始，说明像块的中心点在左上角处开始
	 //步长为MBSIZE，我们寻找的每个簛E？樗阉髌渖舷伦笥襭像素处的匹艁E
     int mbCount = 1;
     int i=0,j=0;
	 //初始化开销矩諄E
	 for(j=0;j<9;j++)
	 {
		 costs[j]=65537;
	 }
     int SDSPFlag;//大菱形中心眮E�
	 int cornerFlag;//大菱形角点眮E�

	 // 从图像左上角开始，说明像块的中心点在左上角处开始
	 //步长为MBSIZE，我们寻找的每个簛E？樗阉髌渖舷伦笥襭像素处的匹艁E
     for( i = 0; i< imgHeight;i+=MBSIZE)//防止溢硜E-图像必衼E荕BSIZE^2整数倍?周边的模块不去运藖E
	 {
//		if(i == 128)
//		{	
//			i = i;
//		}
		for(j = 0; j<imgWidth;j+=MBSIZE)
		{   
//			if(j == 42)
//			{
//				j = j;
//			}
			double costsSDSP[5]={65537,65537,65537,65537,65537};
			double cost[1]={65535};//costs 中臁值
            int    point[1]={4};// costs 中臁值的坐眮E
			double costSDSP[1]={65535};
			int    pointSDSP[1]={4};
			int    xLast,yLast;//存储上一次的中心坐眮E
        //对于三步法，每一步我们评估周围的9部分找出猊配的
            int    x = j;//图像的列
            int    y = i;//图像的行
			int    k;
        // 为了避免对搜索中心的反复计算，我们存储其前一次运行得到的数值。
		//对于第一次祦E颐前汛酥荡娣诺窖分猓婧蟮牡丒颐前阎荡娣诺揭贫沟男轮行牡愦�

			//动态分配二维数讈E
			unsigned char** imgPtmp;
			unsigned char** imgItmp;
			imgPtmp = new unsigned char*[MBSIZE];
			imgItmp = new unsigned char*[MBSIZE];
			for (k=0;k<MBSIZE;k++)
			{		
				imgPtmp[k]=new unsigned char[MBSIZE];
			    memset(imgPtmp[k],0,MBSIZE*sizeof(unsigned char));//新建内存区莵E丒
			}

			for (k=0;k<MBSIZE;k++)
			{			
				imgItmp[k]=new unsigned char[MBSIZE];
                memset(imgItmp[k],0,MBSIZE*sizeof(unsigned char));//新建内存区莵E丒
			}

            for(int blkX=0;blkX<MBSIZE;blkX++)
			{
	    		for(int blkY=0;blkY<MBSIZE;blkY++)
				{
					int MAX1=0;
					//为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					MAX1 =((i+blkY)*imgWidth+(j+blkX)<=(imgWidth*imgHeight-1))?(i+blkY)*imgWidth+(j+blkX):(imgWidth*imgHeight-1); 
                    *(imgPtmp[blkX]+blkY) = *(imgP+MAX1);
		    	    *(imgItmp[blkX]+blkY) = *(imgI+MAX1);
				}
			}
			   //     costs(5) = costFuncMAD(imgP(i:i+MBSIZE-1,j:j+MBSIZE-1), ...
                             //       imgI(i:i+MBSIZE-1,j:j+MBSIZE-1),MBSIZE);
			//循环前先计算参考中心处开蟻E
            costs[4] = costFuncMAD(imgPtmp,imgItmp,MBSIZE);
            computations = computations + 1;

    //**************************************LDSP初次搜藖E***************************************//
		    for(k=0;k<9;k++)
			{
				int refBlkVer,refBlkHor;
				refBlkVer = y + LDSP[k].y;//参考块的行、垂直坐眮E
				refBlkHor = x + LDSP[k].x;//参考块的列、水平坐眮E
				if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
				continue;
				if (k == 4)
                continue;
			    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
				{
	               for(int blkY2=0;blkY2<MBSIZE;blkY2++)
				   {
					  int MAX1=0,MAX2=0;
					//为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					  MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?((i+blkY2)*imgWidth+(j+blkX2)):(imgWidth*imgHeight-1); 
				      MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                      *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			          *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
					}
				}
                costs[k] = costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                computations = computations + 1;
			}
//******************************一次全LDSP搜索结蕘E***************************************//

	//第一次调用计算初始搜索@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
			minSearch(costs,9,cost,point);
            // 菱形中褋EDSPFlag 眮E局�1---历捋上没有位移的话应该蛠E睾希琒DSPFlag始终为1,point[0]不是4就是2
            if (point[0] == 4)
                SDSPFlag = 1;
            else
			{
			    SDSPFlag = 0;
                if ( abs(LDSP[point[0]].x) == abs(LDSP[point[0]].y) )
                   cornerFlag = 0;
                else
				   cornerFlag = 1; //角点眮E�
                xLast = x;
                yLast = y;
			    //搜索中心移到新的中心
                x = x + LDSP[point[0]].x;
                y = y + LDSP[point[0]].y;
                for(int index=0;index<9;index++)
				{
				    costs[index]=65537;
				}
					 //中心值变成臁值
                costs[4] = cost[0];
			}
		//**********************************LDSP循环搜藖E****************************************//
            while (SDSPFlag == 0)//LDSP搜索至SDSPFlag =1为止,分成角点和非角点两纴E
			{
				if (cornerFlag == 1)
					//排除非角点和中心点，计算角点处的开蟻E
                for(int k=0;k<9;k++)
				{      
					int refBlkVer,refBlkHor;
					refBlkVer = y + LDSP[k].y;   
                    refBlkHor = x + LDSP[k].x;   
                    if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                        continue;
                    if (k == 4)
                        continue;      
                    if ( refBlkHor >= xLast - 1  && refBlkHor <= xLast + 1 && refBlkVer >= yLast - 1  && refBlkVer <= yLast + 1 )//非角点组成的正方形区域内
                        continue;// 移除正方形区域内的三个祦E
                    else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P || refBlkVer > i+P)//搜索窗口之蛠E
                        continue;
                    else
					{   
			            for(int blkX2=0;blkX2<MBSIZE;blkX2++)
						{
	                        for(int blkY2=0;blkY2<MBSIZE;blkY2++)
							{
							   int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					           MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				               MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                               *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                   *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
							}
						}
                        costs[k] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);//一般情况计藖E祦E
		            	 //   costs[k] = costFuncMAD(imgP(i:i+MBSIZE-1,j:j+MBSIZE-1), imgI(refBlkVer:refBlkVer+MBSIZE-1, refBlkHor:refBlkHor+MBSIZE-1), MBSIZE);
                        computations = computations + 1;
					}
				}     
                else//非角点，4种情况讨论
				{      
					switch (point[0])// 1  2  6  7
                    {
				    	case '1':      // just 3 block 1 2 4
						{
							int refBlkVer,refBlkHor;
						    refBlkVer = y + LDSP[0].y;  
                            refBlkHor = x + LDSP[0].x;   
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P || refBlkVer > i+P)
                            ;
                            else 
							{                      
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
										int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[0] = (int) costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
                            }
                            refBlkVer = y + LDSP[1].y;   
                            refBlkHor = x + LDSP[1].x;   
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;  //图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ;  //搜索窗口之外，跳硜E
                            else
							{       
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
										int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[1] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
                            }
                            refBlkVer = y + LDSP[3].y;   
                            refBlkHor = x + LDSP[3].x;  
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ; //图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ; //搜索窗口之外，跳硜E
                            else
							{   
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                               for(int blkY2=0;blkY2<MBSIZE;blkY2++)
								   {
									    int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
								   }
								}
                                costs[3] = costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
                                break;
							}
						}
                        case '2': // three blocks 1  3 6
						{
							int refBlkVer,refBlkHor;
					    	refBlkVer = y + LDSP[0].y;  
                            refBlkHor = x + LDSP[0].x;   
                            if ( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight|| refBlkHor < 0 || refBlkHor+MBSIZE> imgWidth)
                            ;
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P || refBlkVer > i+P)
                            ;
                            else 
							{    
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                               for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
									    int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
								   }
								}
                                costs[0] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}

                            refBlkVer = y + LDSP[2].y;   
                            refBlkHor = x + LDSP[2].x;   
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;//图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ; //搜索窗口之外，跳硜E
                            else
							{        
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
									   int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                   MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                       MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                       *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                           *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[2] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}
                        
                            refBlkVer = y + LDSP[5].y;   
                            refBlkHor = x + LDSP[5].x;  
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;//图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ;//搜索窗口之外，跳硜E
                            else
							{        
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
										int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[5] = costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}
				        	break;
						}
                        case '6': // three blocks 4 7 9
						{
							int refBlkVer,refBlkHor;
					        refBlkVer = y + LDSP[3].y;  
                            refBlkHor = x + LDSP[3].x;   
                            if ( refBlkVer < 0 || refBlkVer+MBSIZE> imgHeight|| refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P || refBlkVer > i+P)
                            ;
                            else 
							{         
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                               for(int blkY2=0;blkY2<MBSIZE;blkY2++)
								   {
									   int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                   MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                       MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                       *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                           *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
								   }
								}
                                costs[3] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}

                            refBlkVer = y + LDSP[6].y;   
                            refBlkHor = x + LDSP[6].x;   
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;//图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ;//搜索窗口之外，跳硜E
                            else
							{      
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
										int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                       *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[6] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}
                        
                            refBlkVer = y + LDSP[8].y;   
                            refBlkHor = x + LDSP[8].x;  
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;//图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ;//搜索窗口之外，跳硜E
                            else
							{        
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
										int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[8] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}
                            break;
						}
                        case '7': // three blocks 6 8 9
						{
							int refBlkVer,refBlkHor;
                            refBlkVer = y + LDSP[5].y;  
                            refBlkHor = x + LDSP[5].x;   
                            if ( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight|| refBlkHor < 0 || refBlkHor+MBSIZE> imgWidth)
                            ;
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P || refBlkVer > i+P)
                            ;
                            else 
							{            
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                               for(int blkY2=0;blkY2<MBSIZE;blkY2++)
								   {
									    int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
								   }
								}
                                costs[5] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}

                            refBlkVer = y + LDSP[7].y;   
                            refBlkHor = x + LDSP[7].x;   
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ; //图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ;//搜索窗口之外，跳硜E
                            else
							{              
			                    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
								{
	                                for(int blkY2=0;blkY2<MBSIZE;blkY2++)
									{
										int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                    MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                        MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                        *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                            *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
									}
								}
                                costs[7] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                                computations = computations + 1;
							}
                         
                            refBlkVer = y + LDSP[8].y;   
                            refBlkHor = x + LDSP[8].x;  
                            if( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight || refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                            ;//图像边缘之外，跳硜E
                            else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P  || refBlkVer > i+P)
                            ;//搜索窗口之外，跳硜E
                            else
							{             
			                   for(int blkX2=0;blkX2<MBSIZE;blkX2++)
							   {
	                              for(int blkY2=0;blkY2<MBSIZE;blkY2++)
								  {
									   int MAX1=0,MAX2=0;
					       //为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					                   MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				                       MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                                       *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			                           *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
								  }
							   }
                               costs[8] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                               computations = computations + 1;
							}
				        	break;
						}
				    	default:
						{
				        	break;
						}
					}//switch

				}//else
  
	//第二次调用计算LDSP运算蛠E葽@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//	
		     	minSearch(costs,9,cost,point);
                if (point[0] == 4)
                   SDSPFlag = 1;
                else
				{
			       SDSPFlag = 0;
                   if ( abs(LDSP[point[0]].x) == abs(LDSP[point[0]].y) )
                      cornerFlag = 0;
                   else
					   cornerFlag = 1; //角点眮E�
                   xLast = x;
                   yLast = y;
                   x = x + LDSP[point[0]].x;
                   y = y + LDSP[point[0]].y;
                   for(int index=0;index<9;index++)
				   {
					  costs[index]=65537;
				   }
                   costs[4] = cost[0];

				}
			}//while 循环结蕘E
//************************************* 进葋ESDSP 计藖E**********************************************//
			costsSDSP[2]=cost[0];
            for(k=0;k<5;k++)
			{        
				int refBlkVer,refBlkHor;
				refBlkVer = y + SDSP[k].y;   
                refBlkHor = x + SDSP[k].x;   
                if ( refBlkVer < 0 || refBlkVer+MBSIZE > imgHeight|| refBlkHor < 0 || refBlkHor+MBSIZE > imgWidth)
                continue; 
                else if (refBlkHor < j-P || refBlkHor > j+P || refBlkVer < i-P || refBlkVer > i+P)
                continue;  
			    if (k == 2)
                continue;
			    for(int blkX2=0;blkX2<MBSIZE;blkX2++)
				{
	               for(int blkY2=0;blkY2<MBSIZE;blkY2++)
				   {
					  int MAX1=0,MAX2=0;
					//为了避免溢硜E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!过了则补上边缘
					  MAX1 =((i+blkY2)*imgWidth+(j+blkX2)<=(imgWidth*imgHeight-1))?(i+blkY2)*imgWidth+(j+blkX2):(imgWidth*imgHeight-1); 
				      MAX2 =(((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2))<=(imgWidth*imgHeight-1))?((refBlkVer+blkY2)*imgWidth+(refBlkHor+blkX2)):(imgWidth*imgHeight-1); 
                      *(imgPtmp[blkX2]+blkY2) = *(imgP+MAX1);
			          *(imgItmp[blkX2]+blkY2) = *(imgI+MAX2);
					}
				}
                costsSDSP[k] =  costFuncMAD(imgPtmp,imgItmp,MBSIZE);
                computations = computations + 1;
                   
            }
		//第三次调用用于SDSP@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//		
            minSearch(costsSDSP,5,cost,point);
            x = x + SDSP[point[0]].x;
            y = y + SDSP[point[0]].y;  
         	//判断是否写溢硜E
			if(mbCount>vectorsNum)
				mbCount = vectorsNum;
            vectors[0][mbCount-1] = y - i;    // 运动矢量的行坐眮E
            vectors[1][mbCount-1] = x - j;    // 运动矢量的列坐眮E
			//模块计数加一
            mbCount = mbCount + 1;
       //     costs[3][3]={{65537,65537,65537},{65537,65537,65537},{65537,65537,65537}};
			//将开销函数复原---有其他简单方法么？
			for(k=0;k<5;k++)  
			{
				costsSDSP[k]=65537;
			}	 		
            //还原costs值
		    for(k=0;k<9;k++)
			{
		        costs[k]=65537;
			}
			//释放内磥E
     	    for (k=0;k<MBSIZE;k++)
            	delete []imgPtmp[k];
        	for (k=0;k<MBSIZE;k++)
            	delete []imgItmp[k];
	        delete []imgPtmp;
	        delete []imgItmp;
		}
	 }


	for(i=0;i<2;i++) 
		{
			for(j=0;j<MAX;j++)
			{
				if(j>=vectorsNum)
					motionVect[i][j]=0;
				else
				    motionVect[i][j]=vectors[i][j];
			}
		}
    *DScomputations = computations/(mbCount - 1);    

}