dlgreg.cpp

《Visual C++数字图像获取 处理及实践应用》一书的源代码

		pDbSMatrixT[count*m_nChsFeatureNum + 2] = 1;		
	}

	// 计算pDbSMatrix*pDbSMatrixT，并将结果放入pDbInvMatrix中
	CalMatProduct(pDbSMatrix,pDbSMatrixT,pDbInvMatrix,3,3,m_nChsFeatureNum);

	// 计算pDbInvMatrix的逆矩阵
	CalInvMatrix(pDbInvMatrix, 3);

	// 计算仿射变换系数
	CalMatProduct(pDbBMatrix, pDbSMatrixT, pDbTemp, 2, 3, m_nChsFeatureNum);
	CalMatProduct(pDbTemp, pDbInvMatrix, pDbAffPara, 2, 3, 3);

	// 释放内存
	delete[]pDbBMatrix;
	delete[]pDbSMatrix;
	delete[]pDbSMatrixT;
	delete[]pDbInvMatrix;
	delete[]pDbTemp;

}

/*************************************************************************
 *
 * \函数名称：
 *   CalMatProduct()
 *
 * \输入参数:
 *   double  *pDbSrc1	- 指向相乘矩阵1的内存
 *   double  *pDbSrc2	- 指向相乘矩阵2的内存
 *   double  *pDbDest   - 存放矩阵相乘运行结果的内存指针
 *   int     nX		- 矩阵的尺寸，具体参见函数说明
 *   int     nY		- 矩阵的尺寸，具体参见函数说明
 *   int     nZ		- 矩阵的尺寸，具体参见函数说明
 *
 * \返回值:
 *   无
 *
 * \说明:
 *   该函数计算两个矩阵的相乘，然后将相乘的结果存放在pDbDest中。其中pDbSrc1
 *的大小为nX*nZ，pDbSrc2的大小为nZ*nY，pDbDest的大小为nX*nY
 *
 *************************************************************************
 */
void CDlgReg::CalMatProduct(double* pDbSrc1, double *pDbSrc2, double *pDbDest, int y, int x, int z)
{
	// 循环变量
	//int i,j,m;
/*
	// 矩阵相乘
	for( i=0;i<nX;i++)
		for( j=0;j<nY;j++)
		{
			pDbDest[i*nY + j] = 0;
			for( m=0; m<nZ; m++)
				pDbDest[i*nY + j] += pDbSrc1[i*nZ + m]*pDbSrc2[m *nY + j];
		}
		*/
	for(int i=0;i<y;i++)
		for(int j=0;j<x;j++)
		{
			pDbDest[i*x + j] = 0;
			for(int m=0;m<z;m++)
				pDbDest[i*x + j] += pDbSrc1[i*z + m]*pDbSrc2[m*x + j];
		}
}

/*************************************************************************
 *
 * \函数名称：
 *   CalInvMatrix()
 *
 * \输入参数:
 *   double  *pDbSrc	- 指向矩阵的指针
 *   int     nLen	- 矩阵的尺寸 
 *
 * \返回值:
 *   无
 *
 * \说明:
 *   该函数计算矩阵pDbSrc的逆矩阵，其中pDbSrc的大小为nLen*nLen
 *
 *************************************************************************
 */
BOOL CDlgReg::CalInvMatrix(double *pDbSrc, int nLen)
{
	int *is,*js,i,j,k;
	double d,p;
	is = new int[nLen];
	js = new int[nLen];
	for(k=0;k<nLen;k++)
	{
		d=0.0;
		for(i=k;i<nLen;i++)
			for(j=k;j<nLen;j++)
			{
				p=fabs(pDbSrc[i*nLen + j]);
				if(p>d)
				{
					d     = p; 
					is[k] = i;
					js[k] = j;
				}
			}
		if(d+1.0==1.0)
		{
			delete is;
			delete js;
			return FALSE;
		}
		if(is[k] != k)
			for(j=0;j<nLen;j++)
			{
				p = pDbSrc[k*nLen + j];
				pDbSrc[k*nLen + j] = pDbSrc[(is[k]*nLen) + j];
				pDbSrc[(is[k])*nLen + j] = p;
			}
		if(js[k] != k)
			for(i=0; i<nLen; i++)
			{
				p = pDbSrc[i*nLen + k];
				pDbSrc[i*nLen + k] = pDbSrc[i*nLen + (js[k])];
				pDbSrc[i*nLen + (js[k])] = p;
			}

		pDbSrc[k*nLen + k]=1.0/pDbSrc[k*nLen + k];
		for(j=0; j<nLen; j++)
			if(j != k)
			{
				pDbSrc[k*nLen + j]*=pDbSrc[k*nLen + k];
			}
		for(i=0; i<nLen; i++)
			if(i != k)
				for(j=0; j<nLen; j++)
					if(j!=k)
					{
						pDbSrc[i*nLen + j] -= pDbSrc[i*nLen + k]*pDbSrc[k*nLen + j];
					}
		for(i=0; i<nLen; i++)
			if(i != k)
			{
				pDbSrc[i*nLen + k] *= -pDbSrc[k*nLen + k];
			}
	}
	for(k=nLen-1; k>=0; k--)
	{
		if(js[k] != k)
			for(j=0; j<nLen; j++)
			{
				p = pDbSrc[k*nLen + j];
				pDbSrc[k*nLen + j] = pDbSrc[(js[k])*nLen + j];
				pDbSrc[(js[k])*nLen + j] = p;
			}
		if(is[k] != k)
			for(i=0; i<nLen; i++)
			{
				p = pDbSrc[i*nLen + k];
				pDbSrc[i*nLen + k] = pDbSrc[i*nLen +(is[k])];
				pDbSrc[i*nLen + (is[k])] = p;
			}
	}
	delete is;

	return TRUE;	
}

/*************************************************************************
 *
 * \函数名称：
 *   CalInvMatrix()
 *
 * \输入参数:
 *   double  *pDbAffPara	- 仿射变换系数矩阵
 *
 * \返回值:
 *   CRect			- 返回待配准图象经仿射变换后的区域
 *
 * \说明:
 *   该函数根据仿射变换系数，计算待配准图象仿射变换后的图象尺寸大小
 *
 *************************************************************************
 */
CRect CDlgReg::GetAftAffDim(double* pDbAffPara)
{
	// 基准图象的宽度和高度
	int nBaseImgWidth, nBaseImgHeight;
	nBaseImgWidth = m_pDibInit->m_lpBMIH->biWidth;
	nBaseImgHeight= m_pDibInit->m_lpBMIH->biHeight;

	// 待配准图象的宽度和高度
	int nSamplImgWidth, nSamplImgHeight;
	nSamplImgWidth = m_pDibSamp->m_lpBMIH->biWidth;
	nSamplImgHeight= m_pDibSamp->m_lpBMIH->biHeight;

	// 基准图象的原始区域
	CRect rectBase(0,0,nBaseImgWidth,nBaseImgHeight);

	// 临时变量
	CPoint pointTemp;
	double tx,ty;
	
	// 图象的端点
	pointTemp.x = 0; 
	pointTemp.y = 0;

	// 计算点pointTemp经过仿射变换后的坐标
	tx = pDbAffPara[0*3 +0]*pointTemp.x + 
		pDbAffPara[0*3 + 1]*pointTemp.y + pDbAffPara[0*3 + 2];
	ty = pDbAffPara[1*3 + 0]*pointTemp.x + 
		pDbAffPara[1*3 + 1]*pointTemp.y + pDbAffPara[1*3 + 2];

	// 判断pointTemp经过仿射变换后是否超出原来的大小
	if(tx<rectBase.left)
		rectBase.left = (int)tx;
	if(tx>rectBase.right)
		rectBase.right = (int)tx+1;
	if(ty<rectBase.top)
		rectBase.top = (int)ty;
	if(ty>rectBase.bottom)
		rectBase.bottom = (int)ty+1;

	// 计算端点(0, nSamplImgHeight)变换后的坐标
	pointTemp.x = 0; pointTemp.y = nSamplImgHeight;
	tx = pDbAffPara[0*3 + 0]*pointTemp.x + 
		pDbAffPara[0*3 + 1]*pointTemp.y + pDbAffPara[0*3 + 2];
	ty = pDbAffPara[1*3 +0]*pointTemp.x +
		pDbAffPara[1*3 + 1]*pointTemp.y + pDbAffPara[1*3 + 2];

	// 判断是否越界
	if(tx<rectBase.left)
		rectBase.left = (int)tx;
	if(tx>rectBase.right)
		rectBase.right = (int)tx+1;
	if(ty<rectBase.top)
		rectBase.top = (int)ty;
	if(ty>rectBase.bottom)
		rectBase.bottom = (int)ty+1;

	// 计算端点(nSamplImgWidth, nSamplImgHeight)变换后的坐标
	pointTemp.x = nSamplImgWidth; pointTemp.y = nSamplImgHeight;
	tx = pDbAffPara[0*3 + 0]*pointTemp.x +
		pDbAffPara[0*3 + 1]*pointTemp.y + pDbAffPara[0*3 + 2];
	ty = pDbAffPara[1*3 + 0]*pointTemp.x +
		pDbAffPara[1*3 + 1]*pointTemp.y + pDbAffPara[1*3 + 2];

	// 判断是否越界
	if(tx<rectBase.left)
		rectBase.left = (int)tx;
	if(tx>rectBase.right)
		rectBase.right = (int)tx+1;
	if(ty<rectBase.top)
		rectBase.top = (int)ty;
	if(ty>rectBase.bottom)
		rectBase.bottom = (int)ty+1;

	// 计算端点(nSamplImgWidth, 0)变换后的坐标
	pointTemp.x = nSamplImgWidth; pointTemp.y = 0;
	tx = pDbAffPara[0*3 + 0]*pointTemp.x +
		pDbAffPara[0*3 + 1]*pointTemp.y + pDbAffPara[0*3 + 2];
	ty = pDbAffPara[1*3 + 0]*pointTemp.x +
		pDbAffPara[1*3 + 1]*pointTemp.y + pDbAffPara[1*3 + 2];

	// 判断是否越界
	if(tx<rectBase.left)
		rectBase.left = (int)tx;
	if(tx>rectBase.right)
		rectBase.right = (int)tx+1;
	if(ty<rectBase.top)
		rectBase.top = (int)ty;
	if(ty>rectBase.bottom)
		rectBase.bottom = (int)ty+1;

	// 返回待配准图象变换后的区域大小
	return(rectBase);

}

/*************************************************************************
 *
 * \函数名称：
 *   SetSampImgAftAff()
 *
 * \输入参数:
 *   double  *pDbAffPara	- 仿射变换系数矩阵
 *   CRect   rectNewImg		- 变换后图象的大小尺寸
 *
 * \返回值:
 *   LPBYTE			- 返回变换后的图象
 *
 * \说明:
 *   该函数根据仿射变换系数，计算待配准图象仿射变换后的图象。并返回此图象指针
 *此图象的大小为rectNewImg
 *
 *************************************************************************
 */
LPBYTE CDlgReg::SetSampImgAftAff(double* pDbAffPara, CRect rectNewImg)
{
	// pUnchSect是4*4大小的矩阵数组
	unsigned char *pUnchSect;
	pUnchSect = new unsigned char[4*4];

	// 新的图象宽度和高度
	int nNewImgWidth, nNewImgHeight;
	nNewImgWidth  = rectNewImg.right - rectNewImg.left;
	nNewImgHeight = rectNewImg.bottom- rectNewImg.top;

	// 待配准图象的宽度和高度
	int nSamplImgWidth, nSamplImgHeight;
	nSamplImgWidth = m_pDibSamp->m_lpBMIH->biWidth;
	nSamplImgHeight= m_pDibSamp->m_lpBMIH->biHeight;

	// 待配准图象的存储宽度
	int nSampSaveWidth;
	nSampSaveWidth = m_pDibSamp->GetDibSaveDim().cx;
	
	// pUnchAftAffSamp是一个大小为rectNewImg大小的图象，
	// 其中rectNewImg表示变换后的图象大小
	unsigned char *pUnchAftAffSamp;
	pUnchAftAffSamp = new unsigned char[nNewImgWidth * nNewImgHeight];
	
	double tx,ty;

	// 计算在变换后的图象的数据
	for(int i=0;i<rectNewImg.bottom-rectNewImg.top;i++)
		for(int j=0;j<rectNewImg.right-rectNewImg.left;j++)
		{
			tx = pDbAffPara[0*3 + 0]*(j+rectNewImg.left) +
				pDbAffPara[0*3 + 1]*(i+rectNewImg.top) + pDbAffPara[0*3 + 2];
			ty = pDbAffPara[1*3 + 0]*(j+rectNewImg.left) + 
				pDbAffPara[1*3 + 1]*(i+rectNewImg.top) + pDbAffPara[1*3 + 2];

			for(int m=(int)ty-1;m<=(int)ty+2;m++)
				for(int n=(int)tx-1;n<=(int)tx+2;n++)
				{
					if(m<0||m>=nSamplImgHeight||n<0||n>=nSamplImgWidth)
						pUnchSect[(m-(int)ty+1)*4 + (n-(int)tx+1)] = 0;
					else
						pUnchSect[(m-(int)ty+1)*4 + (n-(int)tx+1)] = 
							m_pDibSamp->m_lpImage[(nSamplImgHeight-m-1)*nSampSaveWidth + n];
				}

			// 确定变换的坐标
			ty = ty - (int)ty + 1;
			tx = tx - (int)tx + 1;

			// 确定变换后此坐标的数值
			pUnchAftAffSamp[i*nNewImgWidth + j] = CalSpline(pUnchSect,tx,ty);
		}
	
	
	// 是否内存
	delete[]pUnchSect;

	// 返回指针
	return (LPBYTE)pUnchAftAffSamp;
	
	
}

/*************************************************************************
 *
 * \函数名称：
 *   CalSpline()
 *
 * \输入参数:
 *   unsigned char	*pUnchCorr	- 插值的点
 *   double		dX		- X坐标
 *   double		dY		- Y坐标
 *
 * \返回值:
 *   unsigned char			- 插值后的值
 *
 * \说明:
 *   该函数根据邻近位置的数值进行插值。
 *此图象的大小为rectNewImg
 *
 *************************************************************************
 */
unsigned char CDlgReg::CalSpline(unsigned char *pUnchCorr, double x, double y)
{
	double ret=0, Cx, Cy;
	double Temp;

	for(int i=0;i<4;i++)
		for(int j=0;j<4;j++)
		{
			Temp = pUnchCorr[i*4 + j];
			if(fabs(y-i)<1)
				Cy = 1-2*fabs(y-i)*fabs(y-i)+fabs(y-i)*fabs(y-i)*fabs(y-i);
			if(fabs(y-i)>=1)
				Cy = 4-8*fabs(y-i)+5*fabs(y-i)*fabs(y-i)-fabs(y-i)*fabs(y-i)*fabs(y-i);
			if(fabs(x-j)<1)
				Cx = 1-2*fabs(x-j)*fabs(x-j)+fabs(x-j)*fabs(x-j)*fabs(x-j);
			if(fabs(x-j)>=1)
				Cx = 4-8*fabs(x-j)+5*fabs(x-j)*fabs(x-j)-fabs(x-j)*fabs(x-j)*fabs(x-j);
			ret += Temp*Cy*Cx;
		}
	if(ret<0)
		ret=0;
	if(ret>255)
		ret=255;

	return (unsigned char)ret;

}

/*************************************************************************
 *
 * \函数名称：
 *   SetBaseImgAftAff()
 *
 * \输入参数:
 *   double  *pDbAffPara	- 仿射变换系数矩阵
 *
 * \返回值:
 *   无
 *
 * \说明:
 *   该函数根据仿射变换系数，计算基准图象仿射变换后的图象，并返回存放此
 *数据的指针
 *
 *************************************************************************
 */
LPBYTE CDlgReg::SetBaseImgAftAff(CRect rectNewImg)
{
	// 新图象的大小
	int nNewImgWidth, nNewImgHeight;
	nNewImgWidth  = rectNewImg.right  - rectNewImg.left;
	nNewImgHeight = rectNewImg.bottom - rectNewImg.top;	

	// 变换后图象
	unsigned char *pUnchAftAffBase;
	pUnchAftAffBase = new unsigned char[nNewImgWidth*nNewImgHeight];

	// 基准图象的高度和宽度
	int nBaseWidth, nBaseHeight;
	nBaseWidth  = m_pDibInit->m_lpBMIH->biWidth;
	nBaseHeight = m_pDibInit->m_lpBMIH->biHeight;

	// 基准图象的存储宽度
	int nBaseSaveWidth;
	nBaseSaveWidth = m_pDibInit->GetDibSaveDim().cx;
	
	for(int i=0;i<rectNewImg.bottom-rectNewImg.top;i++)
		for(int j=0;j<rectNewImg.right-rectNewImg.left;j++)
		{
			if(i<-rectNewImg.top||i>=-rectNewImg.top+nBaseHeight||j<-rectNewImg.left||j>=-rectNewImg.left+nBaseWidth)
				pUnchAftAffBase[i*nNewImgWidth + j] = 0;
			else
				pUnchAftAffBase[i*nNewImgWidth + j] = m_pDibInit->m_lpImage[(nBaseHeight - (i+rectNewImg.top) - 1)*nBaseSaveWidth + (j+rectNewImg.left)];
		}

	// 返回
	return (LPBYTE)pUnchAftAffBase;

}