opencv第四次作业（图像腐蚀，直线检测，数学形态应用）

在这里插入代码片
```import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import scipy.signal as signal
import cv2 as cv
import random
import easygui as g
import imutils
import time
import math

def Canny(img):
    # Gray scale
    def BGR2GRAY(img):
        b = img[:, :, 0].copy()
        g = img[:, :, 1].copy()
        r = img[:, :, 2].copy()

        # Gray scale
        out = 0.2126 * r + 0.7152 * g + 0.0722 * b
        out = out.astype(np.uint8)

        return out

    # Gaussian filter for grayscale
    def gaussian_filter(img, K_size=3, sigma=1.3):

        if len(img.shape) == 3:
            H, W, C = img.shape
            gray = False
        else:
            img = np.expand_dims(img, axis=-1)
            H, W, C = img.shape
            gray = True

        ## Zero padding
        pad = K_size // 2
        out = np.zeros([H + pad * 2, W + pad * 2, C], dtype=np.float)
        out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float)

        ## prepare Kernel
        K = np.zeros((K_size, K_size), dtype=np.float)
        for x in range(-pad, -pad + K_size):
            for y in range(-pad, -pad + K_size):
                K[y + pad, x + pad] = np.exp(- (x ** 2 + y ** 2) / (2 * sigma * sigma))
        # K /= (sigma * np.sqrt(2 * np.pi))
        K /= (2 * np.pi * sigma * sigma)
        K /= K.sum()           #得到了我们自己的高斯模板 k

        tmp = out.copy()

        # filtering
        for y in range(H):
            for x in range(W):
                for c in range(C):
                    out[pad + y, pad + x, c] = np.sum(K * tmp[y: y + K_size, x: x + K_size, c])

        out = np.clip(out, 0, 255)
        out = out[pad: pad + H, pad: pad + W]    #原图经过pading操作后比原图大pading，我们这里裁减掉padding，得到图片原始的真实大小
        out = out.astype(np.uint8)

        if gray:
            out = out[..., 0]      #如果是灰度图，我们将其降维成两维，因为之前将图像的维度上升到了三维

        return out

    # sobel filter
    def sobel_filter(img, K_size=3):
        if len(img.shape) == 3:
            H, W, C = img.shape
        else:
            H, W = img.shape

        # Zero padding
        pad = K_size // 2
        out = np.zeros((H + pad * 2, W + pad * 2), dtype=np.float)
        out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float)
        tmp = out.copy()

        out_v = out.copy()
        out_h = out.copy()

        ## Sobel vertical
        Kv = [[1., 2., 1.], [0., 0., 0.], [-1., -2., -1.]]
        ## Sobel horizontal
        Kh = [[1., 0., -1.], [2., 0., -2.], [1., 0., -1.]]

        # filtering
        for y in range(H):
            for x in range(W):
                out_v[pad + y, pad + x] = np.sum(Kv * (tmp[y: y + K_size, x: x + K_size]))
                out_h[pad + y, pad + x] = np.sum(Kh * (tmp[y: y + K_size, x: x + K_size]))

        out_v = np.clip(out_v, 0, 255)
        out_h = np.clip(out_h, 0, 255)

        out_v = out_v[pad: pad + H, pad: pad + W]
        out_v = out_v.astype(np.uint8)
        out_h = out_h[pad: pad + H, pad: pad + W]
        out_h = out_h.astype(np.uint8)

        return out_v, out_h                 #输出水平方向与竖直方向的sober算子的输出矩阵

    def get_edge_angle(fx, fy):
        # get edge strength
        edge = np.sqrt(np.power(fx.astype(np.float32), 2) + np.power(fy.astype(np.float32), 2))
        edge = np.clip(edge, 0, 255)

        fx = np.maximum(fx, 1e-10)     #maximum输出最大值，这里为了防止我们输出一个零的数值，因此这里特地把0设置为0.00000001
        # fx[np.abs(fx) <= 1e-5] = 1e-5

        # get edge angle
        angle = np.arctan(fy / fx)     #这里如果x=0,我们可以认为垂直的导数无限大，因此我们这里除以0.0000001没毛病

        return edge, angle      #返回边界图的矩阵以及相应的角度矩阵

    def angle_quantization(angle):
        angle = angle / np.pi * 180             #将角度转化弧度制
        angle[angle < -22.5] = 180 + angle[angle < -22.5]
        _angle = np.zeros_like(angle, dtype=np.uint8)
        _angle[np.where(angle <= 22.5)] = 0                             #这里我们设置四个阈值，使得角度判断更正为对0，,45，90,135的判断
        _angle[np.where((angle > 22.5) & (angle <= 67.5))] = 45
        _angle[np.where((angle > 67.5) & (angle <= 112.5))] = 90
        _angle[np.where((angle > 112.5) & (angle <= 157.5))] = 135          #以45度作为一次划分

        return _angle

    def non_maximum_suppression(angle, edge):           #非极大值抑制
        H, W = angle.shape
        _edge = edge.copy()

        for y in range(H):
            for x in range(W):
                if angle[y, x] == 0:
                    dx1, dy1, dx2, dy2 = -1, 0, 1, 0                #因为无论哪个方向都有两种定义方式，因此设置两种的梯度方向信息
                elif angle[y, x] == 45:
                    dx1, dy1, dx2, dy2 = 1, 1, -1, -1
                elif angle[y, x] == 90:
                    dx1, dy1, dx2, dy2 = 0, -1, 0, 1
                elif angle[y, x] == 135:
                    dx1, dy1, dx2, dy2 = -1, 1, 1, -1
                if x == 0:
                    dx1 = max(dx1, 0)
                    dx2 = max(dx2, 0)               #这些步骤的作用是避免我的图像边界点加上梯度信息后超出了边界范围，让最小值为0
                if x == W - 1:
                    dx1 = min(dx1, 0)
                    dx2 = min(dx2, 0)
                if y == 0:
                    dy1 = max(dy1, 0)
                    dy2 = max(dy2, 0)
                if y == H - 1:
                    dy1 = min(dy1, 0)
                    dy2 = min(dy2, 0)
                if max(max(edge[y, x], edge[y + dy1, x + dx1]), edge[y + dy2, x + dx2]) != edge[y, x]:      #根据定义的点的信息来判断当前点与周围点是否为极大值的点，如果不是就舍弃掉。
                    _edge[y, x] = 0

        return _edge

    def hysterisis(edge, HT=100, LT=30):
        H, W = edge.shape

        # Histeresis threshold
        edge[edge >= HT] = 255
        edge[edge <= LT] = 0                #进行我们自己的二值化处理

        _edge = np.zeros((H + 2, W + 2), dtype=np.float32)          #因为我们自己设置的图片是经过了3*3的padding，因此每个人边界平均长了2
        _edge[1: H + 1, 1: W + 1] = edge

        ## 8 - Nearest neighbor
        nn = np.array(((1., 1., 1.), (1., 0., 1.), (1., 1., 1.)), dtype=np.float32)

        for y in range(1, H + 2):
            for x in range(1, W + 2):
                if _edge[y, x] < LT or _edge[y, x] > HT:
                    continue
                if np.max(_edge[y - 1:y + 2, x - 1:x + 2] * nn) >= HT:
                    _edge[y, x] = 255                       #膨胀一下
                else:
                    _edge[y, x] = 0

        edge = _edge[1:H + 1, 1:W + 1]              #裁掉多余的部分

        return edge

    # grayscale
    gray = BGR2GRAY(img)            #灰度化

    # gaussian filtering
    gaussian = gaussian_filter(gray, K_size=5, sigma=1.4)           #进行高斯处理，高斯的模板为5*5，sigma为1.4（平滑）

    # sobel filtering
    fy, fx = sobel_filter(gaussian, K_size=3)                       #进行sobel处理，返回水平方向与垂直方向两个矩阵（提取边缘信息）

    # get edge strength, angle
    edge, angle = get_edge_angle(fx, fy)                #根据水平和竖直方向的矩阵返回边界和角度，此时的角度还为弧度制

    # angle quantization
    angle = angle_quantization(angle)                      #将角度进行转化，转化为多少多少度的形式

    # non maximum suppression
    edge = non_maximum_suppression(angle, edge)             #进行非极大值抑制，将每个像素点根据自身以及自己的梯度信息将周围邻域进行比较， 如果不是局部最大值，舍弃掉
    #相当于保留最强的边界信息
    # hysterisis threshold
    out = hysterisis(edge, 100, 30)                         #二值化

    return out


# 霍夫变换实现检测图像中的20条直线
def Hough_Line(edge, img):
    ## Voting
    def voting(edge):
        H, W = edge.shape

        drho = 1
        dtheta = 1

        # get rho max length
        rho_max = np.ceil(np.sqrt(H ** 2 + W ** 2)).astype(np.int)      #ceil为向上取整

        # hough table
        hough = np.zeros((rho_max, 180), dtype=np.int)    #建立hough映射表

        # get index of edge
        # ind[0] 是 符合条件的纵坐标，ind[1]是符合条件的横坐标
        ind = np.where(edge == 255)   #找到符合条件的、满足非极大值抑制的实际点

        ## hough transformation
        # zip函数返回元组
        for y, x in zip(ind[0], ind[1]):
            for theta in range(0, 180, dtheta):
                # get polar coordinat4s
                t = np.pi / 180 * theta                     #将角度在转化为弧度制，方便下面公式的计算
                rho = int(x * np.cos(t) + y * np.sin(t))    #计算hough空间的纵轴的值

                # vote
                hough[rho, theta] += 1                      #对于符合条件的hough空间进行统计，满足条件个数加一

        out = hough.astype(np.uint8)

        return out

    # non maximum suppression
    def non_maximum_suppression(hough):
        rho_max, _ = hough.shape

        ## non maximum suppression
        for y in range(rho_max):
            for x in range(180):
                # get 8 nearest neighbor
                x1 = max(x - 1, 0)
                x2 = min(x + 2, 180)
                y1 = max(y - 1, 0)
                y2 = min(y + 2, rho_max - 1)
                if np.max(hough[y1:y2, x1:x2]) == hough[y, x] and hough[y, x] != 0:             #在hough空间中也进行一次非最大值抑制
                    pass
                    # hough[y,x] = 255
                else:
                    hough[y, x] = 0                     #注意:hough空间中存储的是[y,x]的个数

        return hough

    def inverse_hough(hough, img):
        H, W, _ = img.shape
        rho_max, _ = hough.shape

        out = img.copy()

        # get x, y index of hough table
        # np.ravel 将多维数组降为1维
        # argsort  将数组元素从小到大排序，返回索引
        # [::-1]   反序->从大到小
        # [:20]    前20个
        ind_x = np.argsort(hough.ravel())[::-1][:20]                #取hough空间中最大的20个的索引作为我们的目标
        ind_y = ind_x.copy()
        thetas = ind_x % 180                #因为原数据为（最大可能边长*180），因此我们对180取余数即可得到我们的角度值的大小
        rhos = ind_y // 180                 #因为原数据为（最大可能边长*180）。我们除以180，向下取整，即可得到我们的hough的rhos
        # each theta and rho
        for theta, rho in zip(thetas, rhos):            #通过zip得到每个点的角度和长度
            # theta[radian] -> angle[degree]
            t = np.pi / 180. * theta                    #转化为弧度制
            line1 = 0
            line2 = 0
            arctan1 = 0
            arectan2 = 0
            x1=[]
            y1=[]
            arctan1_y = 0
            line1_end = 0
            # hough -> (x,y)
            for x in range(W):
                if np.sin(t) != 0:
                    y = - (np.cos(t) / np.sin(t)) * x + (rho) / np.sin(t)
                    y = int(y)
                    if y >= H or y < 0:
                        continue
                    out[y, x] = [0, 255, 255]
                    line1 +=1
                    if len(x1) <= 5:
                        x1.append((x,y))
            if (x1[1][0]-x1[0][0])==0 or (x1[1][1]-x1[0][1])==0:
                line1_end = line1
            else:
                arctan1=(x1[1][1]-x1[0][1])/(x1[1][0]-x1[0][0])
                arctan1_y = arctan1 * line1
                line1_end = math.sqrt(arctan1_y**2+line1**2)
            print('定x的线的长度为 %d' % line1_end)
            for y in range(H):
                if np.cos(t) != 0:
                    x = - (np.sin(t) / np.cos(t)) * y + (rho) / np.cos(t)
                    x = int(x)
                    if x >= W or x < 0:
                        continue
                    out[y, x] = [0, 0, 255]
                    if len(y1) <= 5:
                        y1.append((x,y))
                    line2 +=1
            if (y1[1][0] - y1[0][0]) ==0 or (y1[1][1] - y1[0][1])==0:
                line2_end = line2
            else:
                arctan2 = (y1[1][1] - y1[0][1]) / (y1[1][0] - y1[0][0])
                arectan2_y = arectan2 * line2
                line2_end = math.sqrt(arectan2_y**2+line2**2)
            print('定y的线的长度为 %d' % line2_end)
        out = out.astype(np.uint8)

        return out

    # voting
    hough = voting(edge)

    # non maximum suppression
    hough = non_maximum_suppression(hough)

    # inverse hough
    out = inverse_hough(hough, img)

    return out

def otsu_binarization(img, th=128):
    H, W = img.shape
    out = img.copy()

    max_sigma = 0
    max_t = 0

    # determine threshold
    for _t in range(1, 255):
        v0 = out[np.where(out < _t)]
        m0 = np.mean(v0) if len(v0) > 0 else 0.
        w0 = len(v0) / (H * W)
        v1 = out[np.where(out >= _t)]
        m1 = np.mean(v1) if len(v1) > 0 else 0.
        w1 = len(v1) / (H * W)
        sigma = w0 * w1 * ((m0 - m1) ** 2)
        if sigma > max_sigma:
            max_sigma = sigma
            max_t = _t

    # Binarization
    print("threshold >>", max_t)
    th = max_t
    out[out < th] = 0
    out[out >= th] = 255

    return out


# Morphology Dilate
def Morphology_Dilate(img, Dil_time=1):
    H, W = img.shape
    out = img.copy()

    # kernel
    # MF = np.array(((0, 1, 0),
    #                (1, 0, 1),
    #                (0, 1, 0)), dtype=np.int)
    MF = np.ones((11, 11))
    # each erode
    for i in range(Dil_time):
        tmp = np.pad(out, (5, 5), 'edge')
        # erode
        for y in range(5, H):
            for x in range(5, W):
                if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) >255:
                    out[y, x] = 255
    out_gap = out - img
    out_gap = out_gap.astype(np.uint8)                  #膨胀图减去二值图
    cv.imshow('pengzhang - erzhitu',out_gap)
    cv.imshow('dilate',out)

def Morphology_Dilate1(img, Dil_time=1):
    H, W = img.shape
    out = img.copy()

    # kernel
    # MF = np.array(((0, 1, 0),
    #                (1, 0, 1),
    #                (0, 1, 0)), dtype=np.int)
    MF = np.ones((11, 11))
    # each erode
    for i in range(Dil_time):
        tmp = np.pad(out, (5, 5), 'edge')
        # erode
        for y in range(5, H):
            for x in range(5, W):
                if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) > 255:
                    out[y, x] = 255
    return out



# Morphology Erode
def Morphology_Erode(img, Erode_time=1):
    H, W = img.shape
    out = img.copy()

    # kernel
    # MF = np.array(((0, 1, 0),
    #                (1, 0, 1),
    #                (0, 1, 0)), dtype=np.int)
    MF = np.ones((11,11))
    # each erode
    for i in range(Erode_time):
        tmp = np.pad(out, (5, 5), 'edge')
        # erode
        for y in range(5, H):
            for x in range(5, W):
                if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) < 255 * 121:
                    out[y, x] = 0
    out_gap = img-out
    out_gap = out_gap.astype(np.uint8)                              #原图减去复试图
    cv.imshow('erode',out)
    cv.imshow('yuan tu - fushi ',out_gap)
    cv.waitKey(0)

def Morphology_Erode1(img, Erode_time=1):
    H, W = img.shape
    out = img.copy()

    # kernel
    # MF = np.array(((0, 1, 0),
    #                (1, 0, 1),
    #                (0, 1, 0)), dtype=np.int)
    MF = np.ones((11, 11))
    # each erode
    for i in range(Erode_time):
        tmp = np.pad(out, (5, 5), 'edge')
        # erode
        for y in range(5, H):
            for x in range(5, W):
                if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) < 255 * 121:
                    out[y, x] = 0
    return out

def hough_fun(img):
    edge = Canny(img)
                # Hough
    out = Hough_Line(edge, img)

    out = out.astype(np.uint8)
    cv.imshow('hough',out)
    cv.waitKey(0)

def open_and_close(img, Erode_time=1):
    img = img[:,:,0]*0.11 + img[:,:,1] *0.59+img[:,:,2] *0.3
    img= img.astype(np.uint8)
    img = otsu_binarization(img)
    H, W = img.shape
    out = img.copy()
    out_open = Morphology_Erode1(out)     #开运算先腐蚀后膨胀
    out_open= Morphology_Dilate1(out_open)

    out_close = Morphology_Dilate1(out)
    out_close = Morphology_Erode1(out_close)
    #闭运算先膨胀后腐蚀
    cv.imshow('open',out_open)
    cv.imshow('close',out_close)
    cv.waitKey(0)


if __name__ == '__main__':
    msg = "请输入您想要完成的任务（建议您第一步先打开图片）"
    title = '第四次作业'
    choice = ('打开图片', '退出')
    a = g.buttonbox(msg=msg, title=title, choices=choice)
    if a == '打开图片':
        filename = g.fileopenbox(msg="请打开一个jpg文件")
        img = cv.imread(filename)
        msg1 = "选择您想要实现的功能"
        title1 = '第四次作业'
        choice1 = ('图像腐蚀', '直线检测', '数学形态应用','显示原图','显示灰度图','重新选择图片','退出')
        q = 1
        while q:
            b = g.buttonbox(msg=msg1, title=title1, choices=choice1)
            if b == '图像腐蚀':
               img_gray = img[:, :, 0] * 0.11 + img[:, :, 1] * 0.59 + img[:, :, 2] * 0.3
               img_gray = img_gray.astype(np.uint8)
               ostu=otsu_binarization(img_gray)
               cv.imshow('binary image',ostu)
               Morphology_Dilate(ostu)
               Morphology_Erode(ostu)

            elif b == '直线检测':
                hough_fun(img)

            elif b == '数学形态应用':
               open_and_close(img)
            elif b == '显示原图':
                cv.imshow('原图',img)
            elif b == '显示灰度图':
                img_gray1 = img[:, :, 0] * 0.11 + img[:, :, 1] * 0.59 + img[:, :, 2] * 0.3
                img_gray1 = img_gray1.astype(np.uint8)
                cv.imshow('原图',img_gray1)
            elif b == '重新选择图片':
                filename = g.fileopenbox(msg="请打开一个jpg文件")
                img = cv.imread(filename)
            else:
                q = 0