match.cpp

图像处理软件,功能比较基础

// Match.cpp: implementation of the CMatch class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
//#include "Radar.h"
#include "Match.h"

#include "math.h"

#include "Comlib.h"
//#include "Mark.h"
//#include "Mask.h"

//#include "MyTimer.h"

#include "MatchPeak.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

#define  DSP_SDRAM_ADD0    0x03000000      // DSP外存起始地址

#define  F                     FALSE
#define  T                     TRUE
#define  SAVECORRELATION       T           // 是否保存相关面数据
#define  SAVEMASKFILE          F
#define  SAVEHISTOGRAM         F

#define  CORMAPFUSION          T           // DBS两部分作为整体计算相关面融合

#define  DBSMAPTOGETHER        (T && CORMAPFUSION)      // DBS两部分作为整体计算相关面

#define  MARKFIND              T           // 标记实际匹配位置，同MARKRIGHT
#define  MARKMIXED             T           // 标记融合相关面匹配位置
#define  MARKRIGHT             T           // 标记有特征边的匹配位置（反转）
#define  MARKOTHER             T           // 标记另外一边的匹配位置
#define  MARKLOCAL             T           // 标记理想匹配位置

#define  FULL_NORMALIZE_BOX    T           // 是否归一化互相关计算光学图匹配度量，可以变速
#define  FULL_NORMALIZE_REL    T           // 是否归一化互相关计算实孔径匹配度量
#define  FULL_NORMALIZE_DBS    T           // 是否归一化互相关计算DBS匹配度量

#define  SIG_ADD               T           // 相关面直接相加
#define  FUSIONMETHOD          ! SIG_ADD     // 采用模板融合得到匹配相关面

/* 以DSP板数据格式保存图象数据及参考图象预处理数据到文件 */
static CString DspFile[4]={
	"1.dat",    /* Ref.dat 参考图象*/
	"2.dat",    /* Rel.dat 实时图象*/
	"3.dat",    /* KRf.dat 参考图标准差*/
	"4.dat"};   /* Avg.dat 参考图均值*/

/* 图象数据储内存块首地址 */
static DWORD Addr0 = DSP_SDRAM_ADD0 ;

// 好于模板尺寸为5的情况
static int model3_09[3][3]={
	5,7,5,
	7,9,7,
	5,7,5
};

static int model3_13[3][3]={
	5,7,5,
	7,13,7,
	5,7,5
};

static int model3_15[3][3]={
	5,7,5,
	7,15,7,
	5,7,5
};

#define MODELSIZE  3

static int Total3_09 = 57;
static int Total3_13 = 61;
static int Total3_15 = 63;
static int nHSize = MODELSIZE/2;

static int model5[5][5]={
	1,2,3,2,1,
	2,5,7,5,2,
	3,7,9,7,3,
	2,5,7,5,2,
	1,2,3,2,1
};

static CString Memo[8] ={
	"C",
	"L",
	"R",
	"Tgt",
	"Add",
	"309",
	"313",
	"315"
};

inline BOOL Dist(NPOT a, NPOT b, float MaxSize)
{
	return ( fabs(a.x-b.x) <= MaxSize &&  fabs(a.y-b.y) <= MaxSize );
}

NPOT NPAvg(NPOT a, NPOT b, NPOT c)
{
	NPOT ret = a;
	ret.x += b.x;
	ret.x += c.x;
	ret.y += b.y;
	ret.y += c.y;

	ret.x = (ret.x+1.5)/3.0;
	ret.y = (ret.y+1.5)/3.0;

	return ret;
}

NPOT NPAvg(NPOT a, NPOT b)
{
	NPOT ret = a;
	ret.x += b.x;
	ret.y += b.y;

	ret.x = (ret.x+1)/2.0;
	ret.y = (ret.y+1)/2.0;

	ret.cor = 1;  //  only for flag

	return ret;
}

// 匹配点融合规则
BOOL JudgMatchResult(SRESULT rlt[], NPOT pot[] , COORDINATE Local , float MaxDist )
{
	BOOL ret = 0;

	BOOL bLR = Dist(pot[1],pot[2], 2*MaxDist);
	BOOL bCL = Dist(pot[0],pot[1], 2*MaxDist);
	BOOL bCR = Dist(pot[0],pot[2], 2*MaxDist);

	if( bLR && bCL && bCR)            // C L R 分布集中
		pot[8] = NPAvg( pot[0] , pot[1] , pot[2] );
	else if( bLR )                    // L R集中，但和C远
		pot[8] = NPAvg( pot[1] , pot[2] );
	else if( !(bLR || bCL || bCR) )   // C L R 都相隔很远
		pot[8] = ( fabs(pot[1].cor) > fabs(pot[2].cor)) ? pot[1] : pot[2];
	else
		pot[8] = pot[0];

	int xx = Local.x;
	int yy = Local.y;
	float erry = pot[8].y - yy;
	float errx = pot[8].x - xx;
	
	rlt[8].uRow[0] += erry;
	rlt[8].uCol[0] += errx;
	rlt[8].dRow[0] += erry * erry;
	rlt[8].dCol[0] += errx * errx;
#if ERROR_RADIUS  
	// radius <= ErrSize
	rlt[8].nCorrects[0] += ( (erry * erry + errx * errx) <= MaxDist * MaxDist ) ? 1 : 0 ;
#else 
	rlt[8].nCorrects[0] += ( abs(errx) <= MaxDist  && abs(erry) <= MaxDist ) ? 1 : 0 ;
#endif

	return bLR;
}

// new here ! 2000/04/27
// 采用多子区匹配方法来融合得到匹配相关面
// 默认模板3,中心权重为13
void FusionCorMap(float **dRxy , float **dRxyL, float **dRxyR , int nr=65, int nc=65 , BOOL bVerL=0 , BOOL bVerR = 0, float fLWeight=0.5)
{
	//计算相关面:dRxy
	float fRWeight ;
	if(fLWeight < 0) fRWeight = - 1 - fLWeight;
	else fRWeight = 1 - fLWeight;

	int sgnL = (bVerL==1) ? -1 : 1;
	int sgnR = (bVerR==1) ? -1 : 1;

	for ( int i=0 ;i<nr;i++)
	for ( int j=0 ;j<nc;j++) 
	{
		dRxy[i][j] = 0;
		for (int m=0;m<MODELSIZE;m++)
		for (int n=0;n<MODELSIZE;n++)
		{
			int k = m+i-nHSize;
			int l = n+j-nHSize;
			if( k < 0 ) k *= -1;
			else if( k >= nr) k = 2*nr - k -1;
			if( l < 0 ) l *= -1;
			else if( l >= nc) l = 2*nc - l -1;

			dRxy[i][j] += sgnL*fLWeight * dRxyL[k][l]*model3_13[m][n];
			dRxy[i][j] += sgnR*fRWeight * dRxyR[k][l]*model3_13[m][n];
		}
		dRxy[i][j] /= Total3_13;
	}

	return ;
}
/////////////////////////////////////////////////////////////////////////////////////////
// 融合匹配相关面数据，包括方法：整体计算，子区相加，子区加权09，子区加权13，子区加权15
// 输入数据：
//           RxyL , RxyR , mr , mc
// 输出数据：
//           RxyTgt , RxyAdd , Rxy09 , Rxy13 , Rxy15
// 控制参数：
//           fLWeight  L子区权重
////////////////////////////////////////
void FusionCorMapALL(float **dRxyL, float **dRxyR ,  float **dRxyAdd , float **dRxy309, float **dRxy313 , float **dRxy315 ,int nr=65, int nc=65 , BOOL bVerL=0 , BOOL bVerR=0 , float fLWeight=0.5)
{
	//计算相关面:dRxy
	float fRWeight ;
	if(fLWeight < 0) fRWeight = - 1 - fLWeight;
	else fRWeight = 1 - fLWeight;

	int sgnL = (bVerL==0) ? 1 : -1;
	int sgnR = (bVerR==0) ? 1 : -1;

	for ( int i=0 ;i<nr;i++)
	for ( int j=0 ;j<nc;j++) 
	{
//		dRxyAdd[i][j] = (dRxyL[i][j] + dRxyR[i][j])/2;
		dRxyAdd[i][j] = sgnL*fLWeight*dRxyL[i][j] + sgnR*fRWeight * dRxyR[i][j] ;
	}

	for ( i=0 ;i<nr;i++)
	for ( int j=0 ;j<nc;j++) 
	{
		dRxy309[i][j] = 0;
		dRxy313[i][j] = 0;
		dRxy315[i][j] = 0;

		for (int m=0;m<MODELSIZE;m++)
		for (int n=0;n<MODELSIZE;n++)
		{
			int k = m+i-nHSize;
			int l = n+j-nHSize;
			if( k < 0 ) k *= -1;
			else if( k >= nr) k = 2*nr - k -1;
			if( l < 0 ) l *= -1;
			else if( l >= nc) l = 2*nc - l -1;

			dRxy309[i][j] += dRxyAdd[k][l]*model3_09[m][n];
			dRxy313[i][j] += dRxyAdd[k][l]*model3_13[m][n];
			dRxy315[i][j] += dRxyAdd[k][l]*model3_15[m][n];
/*
			dRxy309[i][j] += sgnL*fLWeight * dRxyL[k][l]*model3_09[m][n];
			dRxy309[i][j] += sgnR*fRWeight * dRxyR[k][l]*model3_09[m][n];

			dRxy313[i][j] += sgnL*fLWeight * dRxyL[k][l]*model3_13[m][n];
			dRxy313[i][j] += sgnR*fRWeight * dRxyR[k][l]*model3_13[m][n];

			dRxy315[i][j] += sgnL*fLWeight * dRxyL[k][l]*model3_15[m][n];
			dRxy315[i][j] += sgnR*fRWeight * dRxyR[k][l]*model3_15[m][n];
*/		}
		dRxy309[i][j] /= Total3_09;
		dRxy313[i][j] /= Total3_13;
		dRxy315[i][j] /= Total3_15;
	}
	return ;
}

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CMatch::CMatch()
{
}

CMatch::~CMatch()
{
}

/*****************************************/
/* 计算匹配相关面数据（归一化相关系数）    */
/* 仅用于非实孔径图象为实时图象的匹配      */
/*****************************************/
float ** CMatch::CalCormap_Full_Spect(BYTE * * ref, BYTE * * rel, CSize s1, CSize s2)
{
	int r1 = s1.cy;
	int c1 = s1.cx;
	int r2 = s2.cy;
	int c2 = s2.cx;
	int mr = r1-r2 + 1;
	int mc = c1-c2 + 1;
	float SumX ,SumY , SumYY, SumXX , SumXY;
	
	float  ** ret = (float ** )CComlib::fspace_2d(mr,mc, sizeof(float));

	int  counts = r2*c2;
	SumX = SumXX = 0.0f;
	for (int m=0;m<r2;m++) 
	for (int n=0;n<c2;n++) 
	{
		SumX += rel[m][n];
		SumXX += rel[m][n]*rel[m][n];
	}

	for (int i=0;i<mr;i++)
	for (int j=0;j<mc;j++) 
	{
		SumY = SumYY = SumXY = 0.0f;
		for (int m=0;m<r2;m++) 
		for (int n=0;n<c2;n++) 
		{
			SumY  += ref[i+m][j+n];
			SumYY += ref[i+m][j+n] * ref[i+m][j+n];
			SumXY += ref[i+m][j+n] * rel[m][n];
		}
		
		float up = SumXY - SumX*SumY/counts;
		float down = (SumXX - SumX*SumX/counts) * (SumYY - SumY*SumY/counts);
		ret[i][j] = (float)  up/sqrt(down);
	}

	return ret;
}




/*******************************************************/
/* 计算匹配相关面数据并保存图象相关数据到DSP板格式文件 */
/* 仅用于实孔径图象为实时图象的匹配                    */
/*******************************************************/
float ** CMatch::CalCormap_SaveToDsp(BYTE * * ref, BYTE * * rel, BYTE ** Mask,CSize s1,CSize s2)
{
	int r1 = s1.cy;
	int c1 = s1.cx;
	int r2 = s2.cy;
	int c2 = s2.cx;
	int mr = r1-r2 + 1;
	int mc = c1-c2 + 1;
	long SumX ,SumY , SumYY,  SumXY;
	
	float  ** ret = (float ** )CComlib::fspace_2d(mr,mc, sizeof(float));
	long  * KRef = new long[mr * mc];
	long  * AvgY = new long[mr * mc];

	SumX = 0;
	int counts = 0;
	for (int i=0;i<r2;i++) 
	for (int j=0;j<c2;j++) 
	{
		if ( Mask[i][j] == 255 )
		{
			SumX += rel[i][j];
			counts ++;
		}
	}

	for ( i=0;i<mr;i++)
	for (int j=0;j<mc;j++) 
	{
		SumY = SumYY = SumXY = 0;
		for (int m=0;m<r2;m++) 
		for (int n=0;n<c2;n++) 
		{
			SumXY += ref[i+m][j+n] * rel[m][n];

			if( Mask[m][n] == 255 ) 
			{
				SumY  += ref[i+m][j+n];
				SumYY += ref[i+m][j+n] * ref[i+m][j+n];
			}
		}

		KRef[i*mc+j] =(int) (sqrt(SumYY-SumY*SumY/counts) +0.5);
		AvgY[i*mc+j] = (int) (SumY/counts+0.5);

		long  tmp = (long) (SumXY- SumX * AvgY[i*mc+j]) / KRef[i*mc+j];

		ret[i][j] = (float) tmp;
	}

	// 图象数据储内存块首地址
	DWORD Addr = Addr0 ;

	// 存储参考图象数据到文件
	BYTE * pref = CComlib::Trans2To1(ref , CSize(c1,r1));
	BYTE * prel = CComlib::Trans2To1(rel , CSize(c2,r2));
	Addr += sizeof(DWORD) * mc*mr;
	CComlib::SaveAsDspDataFile(pref, r1*c1, Addr, DspFile[0] , 'b');
	
	// 存储实时图象数据到文件
	Addr += sizeof(BYTE) * r1*c1;
	CComlib::SaveAsDspDataFile(prel, r2*c2, Addr, DspFile[1] , 'b');
	
	// 存储参考图象方差数据到文件
	Addr += sizeof(BYTE) * r2*c2;
	CComlib::SaveAsDspDataFile(KRef, mr*mc, Addr, DspFile[2] , 'l');

	// 存储参考图象均值数据到文件
	Addr += sizeof(DWORD) * mr*mc;
	CComlib::SaveAsDspDataFile(KRef, mr*mc, Addr, DspFile[3] , 'l');

	delete KRef;
	delete AvgY;
	delete pref;
	delete prel;

	return ret;
}

// 将DBS实时图象L和R作为整体计算相关面
float ** CMatch::CalTwoPartsCorMap(BYTE **ref1,BYTE ** rel1,  BYTE ** Mask1, 
	BYTE **ref2,BYTE ** rel2,  BYTE ** Mask2, 
	CSize s1,CSize s2)
{
	int r1 = s1.cy;
	int c1 = s1.cx;
	int r2 = s2.cy;
	int c2 = s2.cx;
	int mr = r1-r2 + 1;
	int mc = c1-c2 + 1;
	float SumX ,SumY , SumYY, SumXX , SumXY;
	
	float  ** ret = (float ** )CComlib::fspace_2d(mr,mc, sizeof(float));

	SumX = SumXX = 0;
	int counts = 0;

	for (int i=0;i<r2;i++) 
	for (int j=0;j<c2;j++) 
	{