canny 算法 python实现， 双边滤波--自适应阈值改进--形态学操作

# -*- coding: utf-8 -*-
import numpy as np
import cv2
import os
import csv
#高斯滤波
def smooth(image, sigma=1.4, length=5):
    # Compute gaussian filter
    k = length // 2
    gaussian = np.zeros([length, length])
    for i in range(length):
        for j in range(length):
            gaussian[i, j] = np.exp(-((i - k) ** 2 + (j - k) ** 2) / (2 * sigma ** 2))
    gaussian /= 2 * np.pi * sigma ** 2
    gaussian = gaussian / np.sum(gaussian)

    W, H = image.shape
    new_image = np.zeros([W - k * 2, H - k * 2])

    for i in range(W - 2 * k):
        for j in range(H - 2 * k):
            new_image[i, j] = np.sum(image[i:i + length, j:j + length] * gaussian)

    new_image = np.uint8(new_image)

    return new_image


# 双边滤波
def bilateral_filter(image, d=3, sigmaColor=20, sigmaSpace=40):
    # 应用双边滤波
    filtered_image = cv2.bilateralFilter(image, d, sigmaColor, sigmaSpace)
    return filtered_image

#梯度
def get_gradient_and_direction(image):
    Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

    W, H = image.shape
    gradients = np.zeros([W - 2, H - 2])
    direction = np.zeros([W - 2, H - 2])

    for i in range(W - 2):
        for j in range(H - 2):
            dx = np.sum(image[i:i + 3, j:j + 3] * Gx)
            dy = np.sum(image[i:i + 3, j:j + 3] * Gy)
            gradients[i, j] = np.sqrt(dx ** 2 + dy ** 2)
            if dx == 0:
                direction[i, j] = np.pi / 2
            else:
                direction[i, j] = np.arctan(dy / dx)

    gradients = np.uint8(gradients)

    return gradients, direction


# 自适应阈值
def otsu_thresholds(image, high_ratio=0.5):
    hight_value, otsu_thresh = cv2.threshold(image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    high_thresh = hight_value
    low_thresh = int(high_ratio * high_thresh)
    return low_thresh, high_thresh


#非极大抑制
def NMS(gradients, direction):
    W, H = gradients.shape
    nms = np.copy(gradients[1:-1, 1:-1])

    for i in range(1, W - 1):
        for j in range(1, H - 1):
            theta = direction[i, j]
            weight = np.tan(theta)
            if theta > np.pi / 4:
                d1 = [0, 1]
                d2 = [1, 1]
                weight = 1 / weight
            elif theta >= 0:
                d1 = [1, 0]
                d2 = [1, 1]
            elif theta >= - np.pi / 4:
                d1 = [1, 0]
                d2 = [1, -1]
                weight *= -1
            else:
                d1 = [0, -1]
                d2 = [1, -1]
                weight = -1 / weight

            g1 = gradients[i + d1[0], j + d1[1]]
            g2 = gradients[i + d2[0], j + d2[1]]
            g3 = gradients[i - d1[0], j - d1[1]]
            g4 = gradients[i - d2[0], j - d2[1]]

            grade_count1 = g1 * weight + g2 * (1 - weight)
            grade_count2 = g3 * weight + g4 * (1 - weight)

            if grade_count1 > gradients[i, j] or grade_count2 > gradients[i, j]:
                nms[i - 1, j - 1] = 0

    return nms

#高低阈值
def double_threshold(nms, threshold1, threshold2):
    visited = np.zeros_like(nms)
    output_image = nms.copy()
    W, H = output_image.shape

    def dfs(i, j):
        if i >= W or i < 0 or j >= H or j < 0 or visited[i, j] == 1:
            return
        visited[i, j] = 1
        if output_image[i, j] > threshold1:
            output_image[i, j] = 255
            dfs(i - 1, j - 1)
            dfs(i - 1, j)
            dfs(i - 1, j + 1)
            dfs(i, j - 1)
            dfs(i, j + 1)
            dfs(i + 1, j - 1)
            dfs(i + 1, j)
            dfs(i + 1, j + 1)
        else:
            output_image[i, j] = 0

    for w in range(W):
        for h in range(H):
            if visited[w, h] == 1:
                continue
            if output_image[w, h] >= threshold2:
                dfs(w, h)
            elif output_image[w, h] <= threshold1:
                output_image[w, h] = 0
                visited[w, h] = 1

    for w in range(W):
        for h in range(H):
            if visited[w, h] == 0:
                output_image[w, h] = 0
    return output_image

#形态学-膨胀-闭操作
def dilate_close(image, kernel_dilaterd=(3,3), kernel_close=(3,3)):
    image = cv2.dilate(image,  np.ones(kernel_dilaterd, np.uint8), iterations=1)  #膨胀
    closed= cv2.morphologyEx(image, cv2.MORPH_CLOSE, np.ones(kernel_close, np.uint8), iterations=1)  #闭操作
    return closed

#改进的Canny  （传灰度图）
def canny_binary_filer_otus(image):
    smoothed_image = smooth(image)   #高斯滤波
    #smoothed_image = bilateral_filter(image, d=9, sigmaColor=35, sigmaSpace=75)  # 双边滤波
    gradients, direction = get_gradient_and_direction(smoothed_image)  # 梯度计算
    nms = NMS(gradients, direction)  # 非极大抑制
    low_thresh, high_thresh = otsu_thresholds(smoothed_image, 0.5)
    output_image = double_threshold(nms, low_thresh, high_thresh)  # 双阈值
    output_image = dilate_close(output_image, kernel_dilaterd=(3, 3), kernel_close=(3, 3))  # 形态学(膨胀-闭操作)
    return output_image