改进的canny算法

import pic_process as pre
img = plt.imread('many.jpg')


image=img.copy()
sigma1 = sigma2 = 1
sum = 0

gaussian = np.zeros([5, 5])
for i in range(5):
    for j in range(5):
        gaussian[i, j] = math.exp(-1 / 2 * (np.square(i - 3) / np.square(sigma1)  # 生成二维高斯分布矩阵
                                            + (np.square(j - 3) / np.square(sigma2)))) / (2 * math.pi * sigma1 * sigma2)
        sum = sum + gaussian[i, j]

gaussian = gaussian / sum


# print(gaussian)

def rgb2gray(rgb):
    return np.dot(rgb[..., :3], [0.299, 0.587, 0.114])


# step1.高斯滤波
#gray = rgb2gray(img)
gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
#gray = cv2.equalizeHist(gray)
W, H = gray.shape
new_gray = np.zeros([W - 5, H - 5])
for i in range(W - 5):
    for j in range(H - 5):
        gray[i, j] = np.sum(gray[i:i + 5, j:j + 5] * gaussian)  # 与高斯矩阵卷积实现滤波
W, H = 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(gray[i:i + 5, j:j + 5] * gaussian)  # 与高斯矩阵卷积实现滤波

# plt.imshow(new_gray, cmap="gray")


# step2.增强 通过求梯度幅值
W1, H1 = new_gray.shape
dx = np.zeros([W1 - 1, H1 - 1])
dy = np.zeros([W1 - 1, H1 - 1])
dx2 = np.zeros([W1 - 1, H1 - 1])
dy2 = np.zeros([W1 - 1, H1 - 1])
d = np.zeros([W1 - 1, H1 - 1])
for i in range(W1 - 1):
    for j in range(H1 - 1):
        # dx[i, j] = (2*new_gray[i, j + 1] - 2*new_gray[i, j-1] + gray[i -1, j + 1] - gray[i -1, j-1]
        #             +new_gray[i+1, j + 1] - new_gray[i+1, j-1])
        # dy[i, j] = (2*new_gray[i+1, j ] - 2*new_gray[i-1, j] + gray[i +1, j +1] - gray[i -1, j+1]
        #             +new_gray[i+1, j -1] - new_gray[i-1, j-1])
        # dx[i, j] = new_gray[i, j + 1] - new_gray[i, j]
        # dy[i, j] = new_gray[i + 1, j] - new_gray[i, j]
        dx[i, j] = new_gray[i, j + 1] - new_gray[i, j]+new_gray[i+1, j + 1] - new_gray[i+1, j]
        dy[i, j] = new_gray[i + 1, j] - new_gray[i, j]+new_gray[i+1, j + 1] - new_gray[i, j+1]
        #
        # dx2[i, j] = new_gray[i+1, j + 1] - new_gray[i, j]
        # dy2[i, j] = new_gray[i -1, j+1] - new_gray[i, j]

        # dx[i, j] = new_gray[i, j + 1] -2* new_gray[i, j]+new_gray[i, j - 1]
        # dy[i, j] = new_gray[i + 1, j] - 2*new_gray[i, j]+new_gray[i - 1, j]
        d[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j]))  # 图像梯度幅值作为图像强度值
        #d[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j])+np.square(dx2[i, j]) + np.square(dy2[i, j]))  # 图像梯度幅值作为图像强度值

# plt.imshow(d, cmap="gray")


# setp3.非极大值抑制 NMS
W2, H2 = d.shape
NMS = np.copy(d)
NMS[0, :] = NMS[W2 - 1, :] = NMS[:, 0] = NMS[:, H2 - 1] = 0
p=0.5
for i in range(1, W2 - 1):
    for j in range(1, H2 - 1):

        if d[i, j] == 0:
            NMS[i, j] = 0
        else:
            gradX = dx[i, j]
            gradY = dy[i, j]
            gradTemp = d[i, j]

            # 如果Y方向幅度值较大
            if np.abs(gradY) > np.abs(gradX):
                weight = np.abs(gradX) / np.abs(gradY)
                grad2 = d[i,j]+p*(d[i-1,j]-d[i,j])
                grad4 = d[i, j] + p * (d[i + 1, j] - d[i, j])
                # 如果x,y方向梯度符号相同
                if gradX * gradY > 0:

                    grad1 = d[i,j-1]+p*(d[i-1,j-1]-d[i,j-1])
                    grad3 = d[i, j + 1] + p * (d[i + 1, j + 1] - d[i, j + 1])

                # 如果x,y方向梯度符号相反
                else:

                    grad1 = d[i , j + 1]+p*(d[i-1,j+1]-d[i , j + 1])
                    grad3 = d[i , j - 1]+p*(d[i+1,j-1]-d[i , j - 1])


            # 如果X方向幅度值较大
            else:
                weight = np.abs(gradY) / np.abs(gradX)
                grad2 = d[i , j ]+p*(d[i,j-1]-d[i , j ])
                grad4 = d[i , j ]+p*(d[i,j+1]-d[i , j ])
                # 如果x,y方向梯度符号相同
                if gradX * gradY > 0:
                    grad1 = d[i+1,j]+p*(d[i - 1, j + 1]-d[i+1,j])
                    grad3 = d[i-1,j]+p*(d[i - 1, j + 1]-d[i-1,j])
                # 如果x,y方向梯度符号相反
                else:
                    grad1 = d[i-1,j]+p*(d[i - 1, j - 1]-d[i-1,j])
                    grad3 = d[i+1,j]+p*(d[i + 1, j + 1]-d[i+1,j])

            gradTemp1 = weight *p* grad1 + (1 - weight*p) * grad2
            gradTemp2 = weight *p* grad3 + (1 - weight*p) * grad4
            if gradTemp >= gradTemp1 and gradTemp >= gradTemp2:
                NMS[i, j] = gradTemp
            else:
                NMS[i, j] = 0

# plt.imshow(NMS, cmap = "gray")
#hist = cv2.calcHist([NMS],[0],None,[256],[0,256])
#util.image_read(NMS)

arr=NMS.flatten()
arrs=[]
for i in arr:
    if i!=0 and i!=255:
        arrs.append(i)
#arrs=arr

print("hk",NMS.shape)


n, bins, patches = plt.hist(arrs, bins=256, normed=1, facecolor='green', alpha=0.75)
a,b=pre.OTSU(arrs)
print(a,"  ",b)
plt.show()
#print(arr)

# step4. 双阈值算法检测、连接边缘
W3, H3 = NMS.shape
DT = np.zeros([W3, H3])
# 定义高低阈值
TL = 0.15 * np.max(NMS)
TH = 0.3* np.max(NMS)
print("THTH",TH)
TL=(2/3)*b
TH=b

for i in range(1, W3 - 1):
    for j in range(1, H3 - 1):
        if (NMS[i, j] < TL):
            DT[i, j] = 0
        elif (NMS[i, j] > TH):
            DT[i, j] = 255
        elif ((NMS[i - 1, j - 1:j + 1] < TH).any() or (NMS[i + 1, j - 1:j + 1]).any()
              or (NMS[i, [j - 1, j + 1]] < TH).any()):
            DT[i, j] = 255


#plt.imshow(DT, cmap="gray")
util.image_read("DT",DT)
cv2.imwrite("images/xkgslotus_jh_many.jpg",DT)

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
import my_util as pre_gj
#import min_filter as min

def OTSU(arr):
    #hist = cv.calcHist([gray], [0], None, [256], [0, 256])  # 255*1的灰度直方图的数组
    gray_size = len(arr)  # 图像像素数
    k = 0  # 初始化灰度阈值
    best_k = 0  # 最佳阈值
    best_M = 0  # 衡量阈值性
    counts, bin_edges = np.histogram(arr, bins=int(max(arr)))
    p = []  # 灰度出现概率

    # for k in range(30,150):
    for i in range(len(counts)):
        p.insert(i, counts[i] / gray_size)  # 灰度集概率分布
    bb=int(max(arr))
    for k in range(0,int(max(arr))):
        u = 0  # 从1到k的累计出现概率的平均灰度级
        u_t = 0  # 从1到256的累计出现概率的平均灰度级
        σ2_0 = 0  # 类内方差
        σ2_1 = 0  # 类内方差
        σ2_t = 0  # 灰度级的总方差
        sum_0 = np.sum(counts[0:k + 1], axis=0)
        sum_1 = np.sum(counts[k + 1:bb], axis=0)

        w_0 = np.sum(p[0:k + 1:])
        w_1 = np.sum(p[k + 1:bb:])  # 各类的概率

        for i in range(k + 1):
            u = i * p[i] + u

        for i in range(len(counts)):
            u_t = i * p[i] + u_t

        u0 = u / w_0
        u1 = (u_t - u) / w_1  # 各类的平均灰度级

        for i in range(k + 1):
            σ2_0 = (p[i] / w_0) * np.square(i - u0) + σ2_0#背景
        for i in range(k + 1, bb):
            σ2_1 = (p[i] / w_1) * np.square(i - u1) + σ2_1  # 两类的类内方差,钢筋
        for i in range(bb):
            σ2_t = p[i] * np.square(i - u_t) + σ2_t  # 总方差

        σ2_w = w_0 * σ2_0 + w_1 * σ2_1  # 类内方差
        σ2_b = w_0 * w_1 * np.square(u1 - u0)  # 类间方差

        M = σ2_b / σ2_t  # 衡量阈值k的好坏
        if M > best_M:
            best_M = M;
            best_k = k;
    return best_M, best_k

from skimage import data, exposure, img_as_float

def get_image(path):
    # 获取图片
    img = cv.imread(path)
    gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)

    return img, gray


def smothing(gray):
    #均衡化
    gray = exposure.adjust_gamma(gray, 0.8)
    img1 = cv.equalizeHist(gray)
    img1 = cv.equalizeHist(img1)
    aussian = cv.GaussianBlur(img1, (5, 5), 1)
    dst = cv.medianBlur(aussian, 5)

    #dst = cv.medianBlur(dst, 7)
#    dst=min.max_min_value_filter(dst,ksize=2,mode=2)
    #dst=min.fillHole(dst)
    """锐化"""
    #kernel = np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]], np.float32)  # 锐化
    #dst = cv.filter2D(dst, -1, kernel=kernel)

    #blurred = cv2.GaussianBlur(gray, (3, 3), 0)

    return dst

def Thresh_and_blur(gray,gradient):
    #blurred = cv2.GaussianBlur(gradient, (5, 5), 0)
    #(_, thresh) = cv2.threshold(blurred, 0, 255, cv2.THRESH_BINARY)
    #otsu
    M, k = OTSU(gradient)
    print(M, k)
    ret, thresh= cv.threshold(gray, k, 255, cv.THRESH_BINARY)


    return thresh

def image_morphology(thresh):
    # 建立一个椭圆核函数
    kernel = cv.getStructuringElement(cv.MORPH_ELLIPSE, (3, 3))

    thresh = cv.dilate(thresh, kernel, iterations=1)

    closed = cv.erode(thresh, kernel, iterations=2)#腐蚀
    #closed = cv.morphologyEx(closed, cv.MORPH_ELLIPSE, kernel)

    opening = cv.morphologyEx(closed, cv.MORPH_OPEN, kernel)
    return closed

 

因为是第一个是代码片段是canny1.py  第二个代码片段为pic_process.py,调用的话先把两个代码片段分别复制到两个file里面，然后把图片路径改成自己的，运行就会自动保存到python所在路径。

参考文献：[1]李旭,王正勇,吴晓红,滕奇志.一种改进非极大值抑制的Canny边缘检测算法[J].成都信息工程学院学报,2011,26(05):564-569.