if(abs(gy) > abs(gx)){ //计算插值比例
weight = fabs(gx)/fabs(gy); g2 = pMag[nPos-sz.cx]; g4 = pMag[nPos+sz.cx];
//如果x,y两个方向导数的符号相同，C 为当前像素，与g1-g4 的位置关系为：
// g1 g2 // C // g4 g3
if(gx*gy>0){
g1 = pMag[nPos-sz.cx-1]; g3 = pMag[nPos+sz.cx+1]; }
//如果x,y两个方向的方向导数方向相反，C是当前像素，与g1-g4的关系为：
// g2 g1 // C // g3 g4
else {
g1 = pMag[nPos-sz.cx+1]; g3 = pMag[nPos+sz.cx-1]; } }
//如果方向导数x分量比y分量大，说明导数的方向趋向于x分量
else{ //插值比例
weight = fabs(gy)/fabs(gx); g2 = pMag[nPos+1]; g4 = pMag[nPos-1];
//如果x,y两个方向的方向导数符号相同，当前像素C与 g1-g4的关系为 // g3 // g4 C g2 // g1
if(gx * gy > 0){
g1 = pMag[nPos+sz.cx+1]; g3 = pMag[nPos-sz.cx-1]; }
// 如果x,y两个方向导数的方向相反，C与g1-g4的关系为 // g1 // g4 C g2 // g3
else{
g1 = pMag[nPos-sz.cx+1]; g3 = pMag[nPos+sz.cx-1]; } }
{//利用 g1-g4 对梯度进行插值
dTmp1 = weight*g1 + (1-weight)*g2; dTmp2 = weight*g3 + (1-weight)*g4;
//当前像素的梯度是局部的最大值，该点可能是边界点
if(dTmp>=dTmp1 && dTmp>=dTmp2) {
pNSRst[nPos] = 128; } else {
pNSRst[nPos] = 0; //不可能是边界点 } } } } } }
// 统计pMag的直方图，判定阈值
void EstimateThreshold(int *pMag, SIZE sz, int *pThrHigh,
int *pThrLow, LPBYTE pGray, double dRatHigh, double dRatLow) {
LONG y,x,k;
//该数组的大小和梯度值的范围有关，如果采用本程序的算法 //那么梯度的范围不会超过pow(2,10)
int nHist[256]; //可能边界数 int nEdgeNum; int nMaxMag; int nHighCount; nMaxMag = 0; for(k=0;k<256;k++){
nHist[k] = 0; }
//最大梯度数
//初始化
//统计直方图,利用直方图计算阈值
for(y=0;y for(x=0;x if(pGray[y*sz.cx+x]==128){ nHist[pMag[y*sz.cx+x]]++; } } } nEdgeNum = nHist[0]; nMaxMag = 0; //统计经过“非最大值抑制”后有多少像素 for(k=1;k<256;k++){ if(nHist[k] != 0){ nMaxMag = k; } //梯度为0的点是不可能为边界点的经过non-maximum suppression后有多少像素 nEdgeNum += nHist[k]; } //梯度比高阈值*pThrHigh 小的像素点总书目 nHighCount = (int)(dRatHigh * nEdgeNum + 0.5); k=1; nEdgeNum = nHist[1]; //计算高阈值 while((k<(nMaxMag-1)) && (nEdgeNum < nHighCount)){ k++; nEdgeNum += nHist[k]; } *pThrHigh = k; //低阈值 *pThrLow = (int)((*pThrHigh) * dRatLow + 0.5); } //利用函数寻找边界起点 void Hysteresis(int *pMag, SIZE sz, double dRatLow, double dRatHigh, LPBYTE pResult) { LONG y,x; int nThrHigh,nThrLow; int nPos; //估计TraceEdge 函数需要的低阈值，以及Hysteresis函数使用的高阈值 EstimateThreshold(pMag,sz,&nThrHigh,&nThrLow,pResult, dRatHigh,dRatLow); //寻找大于dThrHigh的点，这些点用来当作边界点，然后用TraceEdge函数跟踪该点对应的边界 for(y=0;y for(x=0;x nPos = y*sz.cx + x; //如果该像素是可能的边界点，并且梯度大于高阈值，该像素作为一个边界的起点 if((pResult[nPos]==128) && (pMag[nPos] >= nThrHigh)){ pResult[nPos] = 255; //设置该点为边界点 TraceEdge(y,x,nThrLow,pResult,pMag,sz); } } } for(y=0;y //其他点已经不可能为边界点 for(x=0;x nPos = y*sz.cx + x; if(pResult[nPos] != 255){ pResult[nPos] = 0; } } } } /*根据Hysteresis 执行的结果，从一个像素点开始搜索，搜索以该像素点为边界起点的一条边界的一条边界的所有边界点，函数采用了递归算法。从（x,y)坐标出发，进行边界点的跟踪，跟踪只考虑pResult中没有处理并且可能是边界点的像素（=128），像素值为0表明该点不可能是边界点，像素值为255表明该点已经是边界点*/ void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz) { int xNum[8] = {1,1,0,-1,-1,-1,0,1}; int yNum[8] = {0,1,1,1,0,-1,-1,-1}; LONG yy,xx,k; for(k=0;k<8;k++){ yy = y+yNum[k]; xx = x+xNum[k]; if(pResult[yy*sz.cx+xx]==128 && pMag[yy*sz.cx+xx]>=nThrLow){ pResult[yy*sz.cx+xx]=255;//该点设为边界点，以该点为中心再进行跟踪 TraceEdge(yy,xx,nThrLow,pResult,pMag,sz); } } //对8邻域像素进行查询