常见边缘检测算子极其实现(基于Python和OpenCV)
Posted Hello World
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了常见边缘检测算子极其实现(基于Python和OpenCV)相关的知识，希望对你有一定的参考价值。
1.边缘检测类型和基本原理
在图像处理中，图像边缘常包括三种模型
1).台阶模型：相邻两个像素的灰度值快速变化；如：在离散灰度图像中，灰度值为...0,0,255,255,...就可视为台阶模型；
2).斜坡模型：图像中从亮到暗（或暗到亮）呈现一个类似于斜坡，如在离散灰度图像灰度：...,0,50,100,150,200,255,255...;
3).屋顶模型：...0,0,80,160,255,160,80,0,0,0...
如下图所示
基本原理
边缘检测简单来说，就是对图像进行按模板进行卷积；边缘在图像上一般就表示为明暗交接的地方，从上面的图中看出，若我们将图像看为一个连续可导的函数，那么，对边缘检测，就是求突变的位置，在数学上，我们对这样连续可导图像是通过求导，因此，边缘也是通过同样的原理。但由于图像是离散矩阵形式，求导采用微分的形式。边缘检测中常见的数学方法如下图所示
2.常见的边缘检测算子
2.1Roberts算子
当对对角线方向的边缘感兴趣时，需要一个二维模板，Roberts交叉梯度算子是最早尝试使用这样的具有对角优势的二维算子，具体的函数实现如下
#传入灰度图像 #返回处理后图像 def robert(self,img): #获取边界模板 roberX=roberY=np.zeros((2,2)) roberX[0][0]=-1;roberX[1][1]=1 roberY[0][1]=-1;roberY[1][0]=1 #计算图像大小 m,n=np.shape(img) gx=0 gy=0 M=img[:] for i in range(m): for j in range(n): gx=np.sum(img[i:i+2,j:j+2]*roberX) gy=np.sum(img[i:i+2,j:j+2]*roberY) M[i][j]=abs(gx)+abs(gy) return M
2.2Prewitt算子
#传入灰度图像 #返回处理后图像 def prewitt(self,img): prewitX=[[-1,-1,-1], [0,0,0], [1,1,1]] prewitY=[[-1,0,1], [-1,0,1], [-1,0,1]] m,n=np.shape(img) M=img[:] gx=0 gy=0 for i in range(m-3): for j in range(n-3): gx=np.sum(img[i:i+3,j:j+3]*prewitX) gy=np.sum(img[i:i+3,j:j+3]*prewitY) M[i][j]=abs(gx)+abs(gy) return M
2.3Sobel算子
可以看出两个算子计算相差并不大，
#传入灰度图像 #返回处理后图像 def sobel(self,img): #计算边界模板 sobelX=np.array([[-1,-2,-1], [0,0,0], [1,2,1]]) sobelY=np.array([[-1,0,1], [-2,0,2], [-1,0,1]]) m,n=np.shape(img) gx=0 gy=0 #初始化与原图像相同数组 M=np.zeros_like(img) for i in range(1,m-1): for j in range(1,n-1): gx=np.sum(img[i-1:i+2,j-1:j+2]*sobelX) gy=np.sum(img[i-1:i+2,j-1:j+2]*sobelY) M[i][j]=abs(gx)+abs(gy) return M
3.LOG算子
LOG边缘检测步骤：
1.先用n*n的高斯进行平滑处理；
2.对平滑后的图像再使用拉普拉斯算子进行边缘检测；
3.再找出图像的零交叉点，进行边界连接。
传统的边缘检测算法，实质就是一些模板（又称为算子）在图像上滑动，LOG算子也不例外，其模板计算公式如下图所示：
对于参数δ是要输入的参数，因此，对于模板大小n>=6δ原则，以上函数图像如下，因长的像草帽，又叫墨西哥草帽算子。
图像零交叉点的查找是该方法的核心部分之一，从上述公式可以看出，这是对原始图像进行二阶导数处理过程。若我们考虑上面提到的边缘检测的斜坡模型，求导数的图像如下所示
最左边的为原始图像，中间部分为一阶导数图像，经过一阶导后，斜坡部分为一个常数，最右边为二级求导后图像，二阶处理后斜坡部分变为0，但其进入斜坡的两个端点处，却不为零，而且方向相反（一正一负）。变化部分如图像中箭头所示。因此。我们可以基于该原理，进行边界提取。即通过3*3模板进行边界查找，判断经过上面模板处理后图像为0像素的上下、左右，两个对角4个方向是否相反，若相反，则标记为边缘。对于台阶模型边缘只需其上下、左右4个方向上符号相反即可。
#MarrHildreth边界检测 class marrHildreth(): def edgesMarrHildreth(self,img,sigma): #动态确定滤波器大小，大于6sigma的奇数 size=int(2*(np.ceil(3*sigma))+1) x,y=np.meshgrid(np.arange(-size/2-1,size/2+1),np.arange(-size/2-1,size/2+1)) # normal=1/(2.0*np.pi*sigma**2) #计算LoG滤波器 kernel=((x**2+y**2-(2.0*sigma**2))/sigma**4)*np.exp(-(x**2+y**2)/(2.0*sigma**2)) kern_size=kernel.shape[0] m,n=np.shape(img) log=np.zeros_like(img,dtype=float) t=int(kern_size/2)#定义起始点 for i in range(t,m-t): for j in range(t,n-t): window=img[i-t:i+(t+1),j-t:j+(t+1)]*kernel log[i,j]=np.sum(window) log=log.astype(np.int64,copy=False) zero_crossing=np.zeros_like(log) #寻找零交叉点 for i in range(1,log.shape[0]-1): for j in range(1,log.shape[1]-1): #斜坡型，使用3*3模板，零交叉点处，至少相邻俩个方向二阶导数符号相反， #上/下；左/右，两个对角方向 if log[i][j]==0: if(log[i][j-1] < 0 and log[i][j+1] > 0) \\ or (log[i][j-1] > 0 and log[i][j+1] < 0) \\ or (log[i-1][j] < 0 and log[i+1][j] > 0) \\ or (log[i-1][j] > 0 and log[i+1][j] < 0) \\ or ((log[i - 1][i - 1] > 0) and (log[i + 1][j + 1] < 0)) \\ or ((log[i - 1][i - 1] < 0) and (log[i + 1][j + 1] > 0)): zero_crossing[i][j]=255 #阶梯型，只要两个相邻（4领域）中符号相反 if log[i][j]<0: if(log[i][j-1]>0)or(log[i][j+1]>0)or(log[i-1][j]>0)or(log[i+1][j]>0): zero_crossing[i][j]=255 return zero_crossing
可以看出LOG结果并不理想，这应该是程序设计有问题！
4.Canny算子，最经典的一种算子之一
基本步骤
1.用二维高斯平滑图像；由于二维高斯过程计算比较复杂，通常可以将其近似看为沿x，y两个方向上计算；
2.计算图像梯度；
3.非最大抑制；
4.双阈值进行边界提取
高斯平滑、计算梯度和角度的公式如上所示。这里还有一个重要的性质是梯度的方向和边缘的方向相互垂直。
#将RGB图像转化为灰度图像 def gray(self,img_path): img=plt.imread(img_path) img_rgb=cv2.cvtColor(img,cv2.COLOR_BGR2RGB) img_gray=np.dot(img_rgb[...,:3],[0.299,0.587,0.114]) return img_gray #高斯滤波平滑图像 def smooth(self,img_gray,sigmal1=1.4,sigmal2=1.4): sigmal1=sigmal1 sigmal2=sigmal2 gau_sum=0 gaussian=np.zeros([5,5]) for i in range(5): for j in range(5): gaussian[i,j]=math.exp((-1/(2*sigmal2*sigmal1))*(np.square(i-3)+ np.square(j-3)))/(2*math.pi*sigmal1*sigmal2) gau_sum=gau_sum+gaussian[i,j] gaussian=gaussian/gau_sum W,H=img_gray.shape new_gray=np.zeros([W-5,H-5]) for i in range(W-5): for j in range(H-5): new_gray[i,j]=np.sum(img_gray[i:i+5,j:j+5]*gaussian) return new_gray #计算dx,dy,M,thana def gradients(self,new_gray): W,H=new_gray.shape dx=np.zeros([W-1,H-1]) dy=np.ones([W-1,H-1]) M=np.zeros([W-1,H-1]) theta=np.zeros([W-1,H-1]) for i in range(W-1): for j in range(H-1): dx[i,j]=new_gray[i+1,j]-new_gray[i,j] dy[i,j]=new_gray[i,j+1]-new_gray[i,j] M[i,j]=np.sqrt(np.square(dx[i,j])+np.square(dy[i,j])) theta[i,j]=math.atan(dx[i,j]/(dy[i,j]+0.0000001)) return dx,dy,M,theta
非最大值抑制：
核心：梯度是有方向矢量，通常在进行非最大抑制时使用3*3模板进行查找该方向线上的两个向量，要求给定像元梯度要大于这两个方向线上的梯度值。
下面介绍两种非最大值抑制的方法
1)基于灰度图像插值的方法
从上图中可以看到，如前面介绍，梯度方向与边缘相互垂直，即，图中短线代表边缘，长的带箭头的表示梯度方向，在图中所示，上面部分在b中，但离a较近；下面部分在g中，但其离h较近，这时，我们不能简单的将上面部分看为b，下半部视为g，因该通过插值的方法通过a、b计算上面部分值，同样通过下半部分g、h计算下面值。
将以上过程简化为上图所示，图中虚线代表梯度方向，中间两个图均有一个共同的性质，即沿y轴方向导数大于沿x轴方向导数。同理，最下面两个图沿x轴方向导数大于沿y方向导数。首先我们要明确我们的目标是插值计算出图中g,h值。以中间的两个图为例，前面介绍了中间两个图沿y轴方向导数大于x轴；那么我们就可以明确在该条件下首先可以获取到上下两个值。之后还要进一步判断，中间两个图中左图：两个方向（x,y）上的梯度相乘大于0（判断方法：梯度有方向，首先假设梯度方向，再看x、y方向）；中间两个图中右图，两个方向相乘小于零。这样就可以判断出要进行插值的另外两个点。
#非最大抑制方案2 def NMS(self,M,dx,dy): #M:梯度 #dx:x方向导数 #dy:y方向导数 d=np.copy(M) W,H=M.shape NMS=np.zeros((W,H)) NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0 for i in range(1,W-1): for j in range(1,H-1): #找到可能边缘点 if M[i,j]!=0: #x方向导数 gradx=dx[i,j] #y方向导数 grady=dy[i,j] #当前梯度 gradTemp=M[i,j] #如果 x 方向梯度值比较大，说明导数方向趋向于 x 分量 if np.abs(gradx)>np.abs(grady): #权重 weight=np.abs(grady)/np.abs(gradx) grad2=M[i-1,j] grad4=M[i+1,j] if gradx*grady>0: # 如果 x, y 方向导数符号一致 # 像素点位置关系 # g1 # g2 c g4 # g3 grad1=M[i-1,j-1] grad3=M[i+1,j+1] else: # 如果 x，y 方向导数符号相反 # 像素点位置关系 # g3 # g2 c g4 # g1 grad1=M[i-1,j+1] grad3=M[i+1,j-1] else: weight=np.abs(gradx)/np.abs(grady) grad2=M[i,j-1] grad4=M[i,j+1] # 如果 x, y 方向导数符号相反 # 像素点位置关系 # g2 g1 # c # g3 g4 if gradx*grady<0: grad1=M[i+1,j-1] grad3=M[i-1,j+1] # 如果 x，y 方向导数符号一致 # 像素点位置关系 # g1 g2 # c # g4 g3 else: grad1=M[i-1,j-1] grad3=M[i+1,j+1] #分别计算两个方向的灰度差值 gradTemp1=weight*grad2+(1-weight)*grad1 gradTemp2=weight*grad4+(1-weight)*grad3 #若该梯度均大于差值的灰度值，则保留，否则赋值为零，抑制 if gradTemp>=gradTemp1 and gradTemp>=gradTemp2: NMS[i,j]=gradTemp return NMS
2)基于角度非最大值抑制方法
前面提到图像边界与梯度方向相互垂直。
基本流程：
1).令d1、d2、d3、d4分别表示水平、垂直、-45°、+45°角度范围。
2).同样采用3*3模板进行非极大值抑制，但现在使用的是角度。
3).首先根据角度判断是在什么方向。
4).中间像素要大于相邻俩个方向上的值。
#非最大抑制方案1 def NMS1(self,M,dx,dy,theta): #定义常数，将弧度转化为度 gree=180.0/math.pi W,H=np.shape(M) #获取角度值 d = np.copy(theta) NMS=np.zeros((W,H)) NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0 for i in range(1,W-1): for j in range(1,H-1): if M[i,j]!=0: gradTemp=d[i][j]*gree #水平边缘 if((-157.5<=gradTemp<=157.5) or (-22.5<=gradTemp<=22.5)): if(M[i][j]>M[i-1][j] and M[i][j]>M[i+1][j]): NMS[i][j]=M[i][j] #45°边缘 if((112.5<=gradTemp<=157.5) or (-67.5<=gradTemp<=-22.5)): if(M[i][j]>M[i+1][j-1] and(M[i][j]>M[i-1][j+1])): NMS[i][j]=M[i][j] #垂直边缘 if((-112.5<=gradTemp<=-67.5) or (67.5<=gradTemp<=122.5)): if(M[i][j]>M[i-1][j]and(M[i][j]>M[i][j+1])): NMS[i][j]=M[i][j] #-45°边缘 if((-157.5<=gradTemp<=112.5) or (22.5<=gradTemp<=67.5)): if(M[i][j]>M[i-1][j-1]and(M[i][j]>M[i+1][j+1])): NMS[i][j]=M[i][j] return NMS
双阈值边界提取
1).Canny算子采用的是高低俩个阈值来进行边界连接，通常高低阈值取值比为2:1或3:1，即TH:TL=2:1(或3:1);
2).高于TH高阈值的点认为是边界点，而低于低阈值点认为不可能为边界点，
3).高于低阈值，低于高阈值间的点，要进行判断（这里采用3*3模板）在该点的8个邻域中是否有高于高阈值(TH)点，若有，则将该点视为边界点
#双阈值处理 def double_threshold(self,NMS): W,H=NMS.shape DT=np.zeros([W,H]) #高低阈值比为：1：3 TL=0.1*np.max(NMS) TH=0.3*np.max(NMS) for i in range(1,W-1): for j in range(1,H-1): #将低于最小阈值都设为0 if(NMS[i,j]<TL): DT[i,j]=0 #将大于最高阈值都设为1 elif(NMS[i,j]>=TH): DT[i,j]=255 #用8连通方法得到高低阈值间的边缘像素 elif((NMS[i-1,:]>=TH).any()) or ((NMS[i+1,:]>=TH).any()) or ((NMS[i,[j-1,j+1]]>=TH).any()): DT[i,j]=255 # elif((NMS[i-1,[j-1,j+1]>TH].any())or((NMS[i+1,[j-1,j+1]]>TH).any())or((NMS[i,[j-1,j+1]]>TH).any())): # DT[i,j]=255 return DT
Canny算子总代码
#Canny边缘检测 class canny(): #将RGB图像转化为灰度图像 def gray(self,img_path): img=plt.imread(img_path) img_rgb=cv2.cvtColor(img,cv2.COLOR_BGR2RGB) img_gray=np.dot(img_rgb[...,:3],[0.299,0.587,0.114]) return img_gray #高斯滤波平滑图像 def smooth(self,img_gray,sigmal1=1.4,sigmal2=1.4): sigmal1=sigmal1 sigmal2=sigmal2 gau_sum=0 gaussian=np.zeros([5,5]) for i in range(5): for j in range(5): gaussian[i,j]=math.exp((-1/(2*sigmal2*sigmal1))*(np.square(i-3)+ np.square(j-3)))/(2*math.pi*sigmal1*sigmal2) gau_sum=gau_sum+gaussian[i,j] gaussian=gaussian/gau_sum W,H=img_gray.shape new_gray=np.zeros([W-5,H-5]) for i in range(W-5): for j in range(H-5): new_gray[i,j]=np.sum(img_gray[i:i+5,j:j+5]*gaussian) return new_gray #计算dx,dy,M,thana def gradients(self,new_gray): W,H=new_gray.shape dx=np.zeros([W-1,H-1]) dy=np.ones([W-1,H-1]) M=np.zeros([W-1,H-1]) theta=np.zeros([W-1,H-1]) for i in range(W-1): for j in range(H-1): dx[i,j]=new_gray[i+1,j]-new_gray[i,j] dy[i,j]=new_gray[i,j+1]-new_gray[i,j] M[i,j]=np.sqrt(np.square(dx[i,j])+np.square(dy[i,j])) theta[i,j]=math.atan(dx[i,j]/(dy[i,j]+0.0000001)) return dx,dy,M,theta #非最大抑制方案1 def NMS1(self,M,dx,dy,theta): #定义常数，将弧度转化为度 gree=180.0/math.pi W,H=np.shape(M) #获取角度值 d = np.copy(theta) NMS=np.zeros((W,H)) NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0 for i in range(1,W-1): for j in range(1,H-1): if M[i,j]!=0: gradTemp=d[i][j]*gree #水平边缘 if((-157.5<=gradTemp<=157.5) or (-22.5<=gradTemp<=22.5)): if(M[i][j]>M[i-1][j] and M[i][j]>M[i+1][j]): NMS[i][j]=M[i][j] #45°边缘 if((112.5<=gradTemp<=157.5) or (-67.5<=gradTemp<=-22.5)): if(M[i][j]>M[i+1][j-1] and(M[i][j]>M[i-1][j+1])): NMS[i][j]=M[i][j] #垂直边缘 if((-112.5<=gradTemp<=-67.5) or (67.5<=gradTemp<=122.5)): if(M[i][j]>M[i-1][j]and(M[i][j]>M[i][j+1])): NMS[i][j]=M[i][j] #-45°边缘 if((-157.5<=gradTemp<=112.5) or (22.5<=gradTemp<=67.5)): if(M[i][j]>M[i-1][j-1]and(M[i][j]>M[i+1][j+1])): NMS[i][j]=M[i][j] return NMS #非最大抑制方案2 def NMS(self,M,dx,dy): #M:梯度 #dx:x方向导数 #dy:y方向导数 d=np.copy(M) W,H=M.shape NMS=np.zeros((W,H)) NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0 for i in range(1,W-1): for j in range(1,H-1): #找到可能边缘点 if M[i,j]!=0: #x方向导数 gradx=dx[i,j] #y方向导数 grady=dy[i,j] #当前梯度 gradTemp=M[i,j] #如果 x 方向梯度值比较大，说明导数方向趋向于 x 分量 if np.abs(gradx)>np.abs(grady): #权重 weight=np.abs(grady)/np.abs(gradx) grad2=M[i-1,j] grad4=M[i+1,j] if gradx*grady>0: # 如果 x, y 方向导数符号一致 # 像素点位置关系 # g1 # g2 c g4 # g3 grad1=M[i-1,j-1] grad3=M[i+1,j+1] else: # 如果 x，y 方向导数符号相反 # 像素点位置关系 # g3 # g2 c g4 # g1 grad1=M[i-1,j+1] grad3=M[i+1,j-1] else: weight=np.abs(gradx)/np.abs(grady) grad2=M[i,j-1] grad4=M[i,j+1] # 如果 x, y 方向导数符号相反 # 像素点位置关系 # g2 g1 # c # g3 g4 if gradx*grady<0: grad1=M[i+1,j-1] grad3=M[i-1,j+1] # 如果 x，y 方向导数符号一致 # 像素点位置关系 # g1 g2 # c # g4 g3 else: grad1=M[i-1,j-1] grad3=M[i+1,j+1] #分别计算两个方向的灰度差值 gradTemp1=weight*grad2+(1-weight)*grad1 gradTemp2=weight*grad4+(1-weight)*grad3 #若该梯度均大于差值的灰度值，则保留，否则赋值为零，抑制 if gradTemp>=gradTemp1 and gradTemp>=gradTemp2: NMS[i,j]=gradTemp return NMS #双阈值处理 def double_threshold(self,NMS): W,H=NMS.shape DT=np.zeros([W,H]) #高低阈值比为：1：3 TL=0.1*np.max(NMS) TH=0.3*np.max(NMS) for i in range(1,W-1): for j in range(1,H-1): #将低于最小阈值都设为0 if(NMS[i,j]<TL): DT[i,j]=0 #将大于最高阈值都设为1 elif(NMS[i,j]>=TH): DT[i,j]=255 #用8连通方法得到高低阈值间的边缘像素 elif((NMS[i-1,:]>=TH).any()) or ((NMS[i+1,:]>=TH).any()) or ((NMS[i,[j-1,j+1]]>=TH).any()): DT[i,j]=255 # elif((NMS[i-1,[j-1,j+1]>TH].any())or((NMS[i+1,[j-1,j+1]]>TH).any())or((NMS[i,[j-1,j+1]]>TH).any())): # DT[i,j]=255 return DT
以上是关于常见边缘检测算子极其实现(基于Python和OpenCV)的主要内容，如果未能解决你的问题，请参考以下文章
常见边缘检测对比（Roberts算子Prewitt算子Sobel算子Laplacian算子Canny算子）
常见边缘检测对比（Roberts算子Prewitt算子Sobel算子Laplacian算子Canny算子）
常见边缘检测对比（Roberts算子Prewitt算子Sobel算子Laplacian算子Canny算子）