featurealgrithm.cpp

影像融合与融合精度评价源码影像融合与融合精度评价源码

					//  g3     
					else 
					{ 
						g1 = pnMag[nPos-nWidth+1] ;
						g3 = pnMag[nPos+nWidth-1] ;
					}
				}

				// 下面利用g1-g4对梯度进行插值
				{
					dTmp1 = weight*g1 + (1-weight)*g2 ;
					dTmp2 = weight*g3 + (1-weight)*g4 ;
					
					// 当前象素的梯度是局部的最大值
					// 该点可能是个边界点
					if(dTmp>=dTmp1 && dTmp>=dTmp2)
					{
						pUnchRst[nPos] = 128 ;
					}
					else
					{
						// 不可能是边界点
						pUnchRst[nPos] = 0 ;
					}
				}
			} //else
		} // for

	}
}

void CFeatureAlgrithm::Hysteresis(int *pnMag, int nWidth, int nHeight, double dRatioLow,double dRatioHigh, unsigned char *pUnchEdge)
{
	/*\说明:
	   本函数实现类似“磁滞现象”的一个功能，也就是，先调用EstimateThreshold
		*   函数对经过non-maximum处理后的数据pUnchSpr估计一个高阈值，然后判断
		*   pUnchSpr中可能的边界象素(=128)的梯度是不是大于高阈值nThdHigh，如果比
		*   该阈值大，该点将作为一个边界的起点，调用TraceEdge函数，把对应该边界
		*   的所有象素找出来。最后，当整个搜索完毕时，如果还有象素没有被标志成
		*   边界点，那么就一定不是边界点。*/
	// 循环控制变量
	int y;
	int x;
	
	int nThdHigh ;
	int nThdLow  ;
	
	int nPos;
	
	// 估计TraceEdge需要的低阈值，以及Hysteresis函数使用的高阈值
	EstimateThreshold(pnMag, nWidth, nHeight, &nThdHigh, 
		&nThdLow, pUnchEdge,dRatioHigh, dRatioLow);
	
	// 这个循环用来寻找大于nThdHigh的点，这些点被用来当作边界点，然后用
	// TraceEdge函数来跟踪该点对应的边界
	for(y=0; y<nHeight; y++)
	{
		for(x=0; x<nWidth; x++)
		{
			nPos = y*nWidth + x ; 
			
			// 如果该象素是可能的边界点，并且梯度大于高阈值，该象素作为
			// 一个边界的起点
			if((pUnchEdge[nPos] == 128) && (pnMag[nPos] >= nThdHigh))
			{
				// 设置该点为边界点
				pUnchEdge[nPos] = 255;
				TraceEdge(y, x, nThdLow, pUnchEdge, pnMag, nWidth);
			}
		}
	}
	
	// 那些还没有被设置为边界点的象素已经不可能成为边界点
	for(y=0; y<nHeight; y++)
	{
		for(x=0; x<nWidth; x++)
		{
			nPos = y*nWidth + x ; 
			if(pUnchEdge[nPos] != 255)
			{
				// 设置为非边界点
				pUnchEdge[nPos] = 0 ;
			}
		}
	}
}

void CFeatureAlgrithm::EstimateThreshold(int *pnMag, int nWidth, int nHeight, int *pnThdHigh, int *pnThdLow,
										 unsigned char * pUnchEdge, double dRatioHigh,  double dRationLow)
{
/*	\输入参数:
		*   int *pnMag               - 梯度幅度图
		*	 int nWidth               - 图象数据宽度
		*	 int nHeight              - 图象数据高度
		*   int *pnThdHigh           - 高阈值
		*   int *pnThdLow            - 低阈值
		*	 double dRatioLow         - 低阈值和高阈值之间的比例
		*	 double dRatioHigh        - 高阈值占图象象素总数的比例
		*   unsigned char *pUnchEdge - 经过non-maximum处理后的数据
		\说明:
		*   经过non-maximum处理后的数据pUnchEdge，统计pnMag的直方图，确定阈值。
		*   本函数中只是统计pUnchEdge中可能为边界点的那些象素。然后利用直方图，
		*   根据dRatioHigh设置高阈值，存储到pnThdHigh。利用dRationLow和高阈值，
		*   设置低阈值，存储到*pnThdLow。dRatioHigh是一种比例：表明梯度小于
		*   *pnThdHigh的象素数目占象素总数目的比例。dRationLow表明*pnThdHigh
		*   和*pnThdLow的比例，这个比例在canny算法的原文里，作者给出了一个区间。
	*/
	// 循环控制变量
	int y;
	int x;
	int k;
	
	// 该数组的大小和梯度值的范围有关，如果采用本程序的算法，那么梯度的范围不会超过pow(2,10)
	int nHist[1024] ;
	
	// 可能的边界数目
	int nEdgeNb     ;
	
	// 最大梯度值
	int nMaxMag     ;
	
	int nHighCount  ;
	
	nMaxMag = 0     ; 
	
	// 初始化
	for(k=0; k<1024; k++) 
	{
		nHist[k] = 0; 
	}
	
	// 统计直方图，然后利用直方图计算阈值
	for(y=0; y<nHeight; y++)
	{
		for(x=0; x<nWidth; x++)
		{
			// 只是统计那些可能是边界点，并且还没有处理过的象素
			if(pUnchEdge[y*nWidth+x]==128)
			{
				nHist[ pnMag[y*nWidth+x] ]++;
			}
		}
	}
	
	nEdgeNb = nHist[0]  ;
	nMaxMag = 0         ;
	// 统计经过“非最大值抑止(non-maximum suppression)”后有多少象素
	for(k=1; k<1024; k++)
	{
		if(nHist[k] != 0)
		{
			// 最大梯度值
			nMaxMag = k;
		}
		
		// 梯度为0的点是不可能为边界点的
		// 经过non-maximum suppression后有多少象素
		nEdgeNb += nHist[k];
	}
	
	// 梯度比高阈值*pnThdHigh小的象素点总数目
	nHighCount = (int)(dRatioHigh * nEdgeNb +0.5);
	
	k = 1;
	nEdgeNb = nHist[1];
	
	// 计算高阈值
	while( (k<(nMaxMag-1)) && (nEdgeNb < nHighCount) )
	{
		k++;
		nEdgeNb += nHist[k];
	}
	
	// 设置高阈值
	*pnThdHigh = k ;
	
	// 设置低阈值
	*pnThdLow  = (int)((*pnThdHigh) * dRationLow+ 0.5);
}

void CFeatureAlgrithm::TraceEdge(int y, int x, int nLowThd, unsigned char *pUnchEdge, int *pnMag, int nWidth)
{
	// 对8邻域象素进行查询
	int xNb[8] = {1, 1, 0,-1,-1,-1, 0, 1} ;
	int yNb[8] = {0, 1, 1, 1,0 ,-1,-1,-1} ;
	
	int yy ;
	int xx ;
	
	int k  ;
	
	for(k=0; k<8; k++)
	{
		yy = y + yNb[k] ;
		xx = x + xNb[k] ;
		// 如果该象素为可能的边界点，又没有处理过
		// 并且梯度大于阈值
		if(pUnchEdge[yy*nWidth+xx] == 128  && pnMag[yy*nWidth+xx]>=nLowThd)
		{
			// 把该点设置成为边界点
			pUnchEdge[yy*nWidth+xx] = 255 ;
			
			// 以该点为中心进行跟踪
			TraceEdge(yy, xx, nLowThd, pUnchEdge, pnMag, nWidth);
		}
	}
}

void CFeatureAlgrithm::SobelOperator(LPSTR pDib, double *pdGrad)
{
	// 图像的指针 LPSTR pDib;
	//图像数据的指针	m_image=(BYTE *)FindDIBBits(m_view);
	//  double * pdGrad	- 指向梯度数据的指针，含有图像的梯度信息
	// 遍历图象的纵坐标
	int y;
	
	// 遍历图象的横坐标
	int x;
	
	// 图象的长宽大小
	int nWidth			= ::DIBWidth(pDib)		;
	int nHeight			= ::DIBHeight(pDib)		;
	CSize sizeImage		= nWidth*nHeight        ;

	// 图像在计算机在存储中的实际大小
	int WIDE=((nWidth*8)+31)/32*4;
	CSize sizeImageSave	= WIDE*nHeight;
		//图像数据的指针
	BYTE* m_image=(BYTE *)FindDIBBits(pDib);
	

	
//	LPBYTE  lpImage = pDib->m_lpImage;
	
	// 初始化
	for(y=0; y<nHeight ; y++ )
		for(x=0 ; x<nWidth ; x++ )
		{
			*(pdGrad+y*nWidth+x)=0;
		}
		
		// 设置模板系数
		static int nWeight[2][3][3] ;
		nWeight[0][0][0] = -1 ;   
		nWeight[0][0][1] =  0 ;   
		nWeight[0][0][2] =  1 ;   
		nWeight[0][1][0] = -2 ;   
		nWeight[0][1][1] =  0 ;   
		nWeight[0][1][2] =  2 ;   
		nWeight[0][2][0] = -1 ;   
		nWeight[0][2][1] =  0 ;   
		nWeight[0][2][2] =  1 ;   
		
		nWeight[1][0][0] =  1 ;   
		nWeight[1][0][1] =  2 ;   
		nWeight[1][0][2] =  1 ;   
		nWeight[1][1][0] =  0 ;   
		nWeight[1][1][1] =  0 ;   
		nWeight[1][1][2] =  0 ;   
		nWeight[1][2][0] = -1 ;   
		nWeight[1][2][1] = -2 ;   
		nWeight[1][2][2] = -1 ;   
		
		
		
		//这个变量用来表示Laplacian算子象素值
		int nTmp[3][3];
		
		// 临时变量
		double dGrad   ;
		double dGradOne;
		double dGradTwo;
		
		// 模板循环控制变量
		int yy ;
		int xx ;
		
		
		// 下面开始利用Prewitt算子进行计算，为了保证计算所需要的
		// 的数据位于图像数据的内部，下面的两重循环的条件是
		// y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
		// 而不是x<nWidth
		for(y=1; y<nHeight-1 ; y++ )
			for(x=1 ; x<nWidth-1 ; x++ )
			{
				dGrad    = 0 ; 
				dGradOne = 0 ;
				dGradTwo = 0 ;
				// Laplacian算子需要的各点象素值
				
				// 模板第一行
				nTmp[0][0] = m_image[(y-1)*nWidth + x - 1 ] ; 
				nTmp[0][1] = m_image[(y-1)*nWidth + x     ] ; 
				nTmp[0][2] = m_image[(y-1)*nWidth + x + 1 ] ; 
				
				// 模板第二行
				nTmp[1][0] = m_image[y*nWidth + x - 1 ] ; 
				nTmp[1][1] = m_image[y*nWidth + x     ] ; 
				nTmp[1][2] = m_image[y*nWidth + x + 1 ] ; 
				
				// 模板第三行
				nTmp[2][0] = m_image[(y+1)*nWidth + x - 1 ] ; 
				nTmp[2][1] = m_image[(y+1)*nWidth + x     ] ; 
				nTmp[2][2] = m_image[(y+1)*nWidth + x + 1 ] ; 
				
				// 计算梯度
				for(yy=0; yy<3; yy++)
					for(xx=0; xx<3; xx++)
					{
						dGradOne += nTmp[yy][xx] * nWeight[0][yy][xx] ;
						dGradTwo += nTmp[yy][xx] * nWeight[1][yy][xx] ;
					}
					dGrad = dGradOne*dGradOne + dGradTwo*dGradTwo  ;
					dGrad = sqrt(dGrad) ;
					// 梯度值写入内存
					*(pdGrad+y*nWidth+x)=dGrad;
			}
}

void CFeatureAlgrithm::RobertsOperator(LPSTR pDib, double *pdGrad)
{
	//LPSTR  pDib		- 指向原始图象信息
	//	*   double * pdGrad	- 指向梯度数据的指针，含有图像的梯度信息
	// 遍历图象的纵坐标
	int y;
	
	// 遍历图象的横坐标
	int x;
	
	// 图象的长宽大小
	int nWidth			= ::DIBWidth(pDib)		;
	int nHeight			= ::DIBHeight(pDib)		;
	CSize sizeImage		= nWidth*nHeight        ;
	
	// 图像在计算机在存储中的实际大小
	int WIDE=((nWidth*8)+31)/32*4;
	CSize sizeImageSave	= WIDE*nHeight;

	// 图像数据的指针
    BYTE* m_image=(BYTE *)FindDIBBits(pDib);
	
	// 初始化
	for(y=0; y<nHeight ; y++ )
		for(x=0 ; x<nWidth ; x++ )
		{
			*(pdGrad+y*nWidth+x)=0;
		}
		
		
		// 下面开始利用Roberts算子进行计算，为了保证计算所需要的
		// 的数据位于图像数据的内部，下面的两重循环的条件是
		// y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
		// 而不是x<nWidth
		
		//这两个变量用来表示Roberts算子第一个模板的两个象素值
		int nUpLeft;
		int nDownRight;
		
		// 这两个变量用来表示Roberts算子第二个模板的两个象素值
		int nUpRight;
		int nDownLeft;
		
		// 这两个变量用来表示Roberts算子计算的结果
		int nValueOne;
		int nValueTwo;
		
		// 临时变量
		double dGrad;
		
		for(y=0; y<nHeight-1 ; y++ )
			for(x=0 ; x<nWidth-1 ; x++ )
			{
				// Roberts算子第一个模板需要的象素值
				nUpLeft		=*(m_image+y*WIDE+x) ; 
				nDownRight	=*(	m_image+(y+1)*WIDE+x+1	);
				nDownRight	*=-1;
				
				//Roberts算子的第一个模板计算结果
				nValueOne	=nUpLeft+nDownRight	;
				
				// Roberts算子第二个模板需要的象素值
				nUpRight	=*(	m_image+y*WIDE+x+1	)		;
				nDownLeft	=*(	m_image+(y+1)*WIDE+x	);
				nDownLeft	*=-1;
				
				// Roberts算子的第二个模板计算结果
				nValueTwo	=nUpRight+nDownLeft;
				
				// 计算两个偏导数的平方和
				dGrad=nValueOne*nValueOne + nValueTwo*nValueTwo;
				
				// 开方
				dGrad=pow(dGrad,0.5);
				
				// 范数采用欧式距离
				*(pdGrad+y*nWidth+x)=dGrad;
			}
}