图像处理_Canny算法

Canny算子求边缘点具体算法步骤如下：

1. 用高斯滤波器平滑图像．

2. 用一阶偏导有限差分计算梯度幅值和方向.

3. 对梯度幅值进行非极大值抑制.

4. 用双阈值算法检测和连接边缘．

4.1.双阈值的选取

4.2.滞后边界跟踪

第一步：高斯滤波 

int TempltExcuteAsh(BYTE** imageBuf0, int w, int h, int* templt, int tw, int x, int y)
{
	int i, j;
	int m = 0;
	int px, py;
	//依次对邻域中每个像素进行运算
	for (i = 0; i < tw; i++){
		for (j = 0; j < tw; j++){
			//计算对应模板上位置的像素在源图像中的位置, 此处tw是int，则tw/2也是int
			py = y - (int)tw / 2 + i;
			px = x - (int)tw / 2 + j;
			//加权求和
			m += imageBuf0[py][px] * templt[i * tw + j];
		}
	}
	return m;
}

static void SmoothGaussAsh(BYTE* image0, BYTE* image1, UINT w, UINT h)
{
	BYTE** imageBuf0 = CreateImage(image0, w, h);
	BYTE** imageBuf1 = CreateImage(image1, w, h);
	int templt[9] = { 1, 2, 1, 2, 4, 2, 1, 2, 1 };
	int x, y;
	int a;
	int scale;//衰减因子
	scale = 16;
	for (y = 1; y < h - 1; y++){
		for (x = 1; x < w - 1; x++){
			a = TempltExcuteAsh(imageBuf0, w, h, templt, 3, x, y);
			a /= scale;
			//过限处理
			a = a> 255 ? 255 : a;
			a = a < 0 ? 0 : a;
			imageBuf1[y][x] = a;
		}
	}
	free(imageBuf0);
	free(imageBuf1);
}

        模板中每一个点的高斯系数可以由上面的公式计算，这样得到的是不是最终的模板呢？答案不是，需要归一化，也即是每一个点的系数要除以所有系数之和，这样才是最终的二维高斯模板。

      这个里面有个小知识点，要想计算上面的系数，需要知道高斯函数的标准差σ （sigma），还需要知道选3*3还是5*5的模板，也就是模板要多大，实际应用的时候，这两者是有关系的，根据数理统计的知识，高斯分布的特点就是数值分布在（μ—3σ,μ+3σ)中的概率为0.9974，也就是模板的大小其实就是6σ这么大就OK了，但是6σ可能不是奇数，因为我们一定要保证有核心。所以模板窗口的大小一般采用1+2*ceil(3*nSigma) ceil是向上取整函数，例如ceil(0.6)=1。

该处引用https://blog.csdn.net/fengye2two/article/details/79190759

归一化用模板数值的和，取倒数。即 对templt[9]内元素求和，最终得16.

模板窗口取3*3。

第二步：计算梯度值和方向 

/***********
4->A0 A1
   A2 A3
8->N0 N1 N2
   N7 C  N3
   N6 N5 N4
************/
void SideSobel(BYTE* image0, BYTE* image1, UINT w, UINT h, bool type)
{
	int x, y ,a;
	int N0, N1, N2, N3, N4, N5, N6, N7,N;
	int n1, n2;
	double scale = 0.5;
	for (int y = 1; y < h - 1; y++){
		for (int x = 1; x < w - 1; x++){
			N0 = image0[(y - 1) * w + x - 1];
			N1 = image0[(y - 1) * w + x];
			N2 = image0[(y - 1) * w + x + 1];

			N3 = image0[y * w + x - 1];
			N7 = image0[y * w + x + 1];

			N4 = image0[(y + 1) * w + x + 1];
			N5 = image0[(y + 1) * w];
			N6 = image0[(y + 1) * w + x - 1];
			if (type == true)
				N = (abs((N0 + N1 * 2 + N2) - (N4 + N5 * 2 + N6)) +
					 abs((N0 + N7 * 2 + N6) - (N2 + N3 * 2 + N4))) * scale;
			else
			{
				n1 = abs((N0 + N1 * 2 + N2) - (N4 + N5 * 2 + N6));
				n2 = abs((N0 + N7 * 2 + N6) - (N2 + N3 * 2 + N4));
				if (n1 > n2)
					N = n1;
				else
					N = n2;
			}
			a = N / 2 * scale;
			a = a > 255 ? 255 : a;
			image1[y * w + x] = a;
		}
	}
}

第四步：非极大值抑制

void NonMaxSuppress(int *pMag, int* pGradX, int* pGradY, UINT w, UINT h, LPBYTE pNSRst)
{
	LONG x, y;
	int nPos;
	//梯度方向
	int gx, gy;
	//临时变量
	int g1, g2, g3, g4;//4个方向 0，45,90,135
	double weight;
	double dTemp, dTemp1, dTemp2;
	for (x = 0; x < w; x++){
		pNSRst[x] = 0;
		pNSRst[(h - 1)*w + x] = 0;
	}
		
	for (y = 0; y < h; y++){
		pNSRst[y * w] = 0;
		pNSRst[y * w + w - 1] = 0;
	}
	for (y = 1; y < h - 1; y++){
		for (x = 1; x < w - 1; x++){
			nPos = y * w + x;
			if (pMag[nPos] == 0)
				pNSRst[nPos] = 0;
			else{
				//梯度图索引
				dTemp = pMag[nPos];
				//x,y方向导数
				gx = pGradX[nPos];
				gy = pGradY[nPos];
				//如果方向导数y分量比x分量大，说明导数方向趋向于y分量
				if (abs(gy) > abs(gx)){
					weight = (double)(abs(gx)) / abs(gy);
					g2 = pMag[nPos - w];
					g4 = pMag[nPos + w];
					if (gx * gy > 0){   //gx 和 gy不同符号
						g1 = pMag[nPos - w - 1];
						g3 = pMag[nPos + w + 1];
					}
					else{
						g1 = pMag[nPos - w + 1];
						g3 = pMag[nPos + w - 1];
					}
				}
				else{
					weight = (double)(abs(gy)) / abs(gx);
					g2 = pMag[nPos + 1];
					g4 = pMag[nPos - 1];
					if (gx * gy > 0){
						g1 = pMag[nPos + w + 1];
						g3 = pMag[nPos - w - 1];
					}
					else{
						g1 = pMag[nPos - w + 1];
						g3 = pMag[nPos + w - 1];
					}
				}
				dTemp1 = weight * g1 + (1 - weight) * g2;
				dTemp2 = weight * g3 + (1 - weight) * g4;
				//当前像素的梯度是局部最大值，该点可能是边界点
				if (dTemp >= dTemp1 && dTemp >= dTemp2){
					pNSRst[nPos] = 128;
				}
				else{
					pNSRst[nPos] = 0;
				}
			}
		}
	}
}