形态学处理

二值图形态学处理

二值图形态学基本算子

二值图形态学图像处理通常在目标图像中卷积式地移动一个结构元素并进行逻辑运算。

腐蚀
腐蚀运算的结果是结构元素$B$平移$x$后仍包含在目标图像$A$内的所有$x$组成的集合：

A ⊙ B = { x ∣ B + x ⊆ A } A\odot B=\{x|B+x\subseteq A\} A⊙B={x∣B+x⊆A}

也就是说用结构元素找出在目标图像内部可以完全放下该结构元素的值为 1 的区域。腐蚀具有缩小图像中的物体和消除图像中比结构元素小的成分的作用。

膨胀
膨胀运算的结果是结构元素$B$平移$x$后与目标图像$A$至少有 1 个非 0 元素相交时对应的结构元素的原点位置所组成的集合：
A ⊕ B = { x ∣ ( B + x ) ∩ A ≠ ∅ } A\oplus B=\{x|({B}+x)\cap A\neq\emptyset\} A⊕B={x∣(B+x)∩A=∅}

膨胀运算能扩张图像中的物体。

开
开运算是先腐蚀后膨胀：
$$A^{\circ }B=(A\odot B)\oplus B$$

能够消除小物体、断开狭窄细长的连接带和平滑图形边缘。

闭
闭运算是先膨胀后腐蚀：

$$A\cdot B=(A\oplus B)\odot B$$

能够填充物体内细小空洞、连接邻近物体和平滑图形内边缘。

击中击不中
对于目标图像$A$和结构元素$B$，如果$A\cap B\ne \varnothing$，那么称$B$击中A； 如果$A\cap B=\varnothing$， 那么称B击不中A。形态学击中击不中变换表达式为：

$$A⊛B=\left(A\odot B_1\right)\cap \left[A^c\odot B_2\right]$$

其中$A^c=1-A$ 。通常令$B_2=1-B_1$ 。击中击不中变换可以用来检测特定形状（结构元素的形状）所处位置。

import numpy as np
from matplotlib import pyplot as plt
import cv2

np.set_printoptions(threshold=np.inf)
plt.rcParams['font.sans-serif'] = ['SimHei']


def img_erode(bin_im, kernel):
    '''腐蚀'''
    kernel_w = kernel.shape[0]
    kernel_h = kernel.shape[1]

    erode_img = np.zeros(shape=bin_im.shape)
    for i in range(1, bin_im.shape[0] - kernel_w + 2):
        for j in range(1, bin_im.shape[1] - kernel_h + 2):
            a = bin_im[i - 1:i - 1 + kernel_w,
                j - 1:j - 1 + kernel_h]  # 找到每次迭代中对应的目标图像小矩阵
            if np.sum(a * kernel) == np.sum(kernel):  # 判定是否“完全重合”
                erode_img[i, j] = 1
    return erode_img


def img_dilate(bin_im, kernel):
    '''膨胀'''
    kernel_w = kernel.shape[0]
    kernel_h = kernel.shape[1]
    dilate_img = np.zeros(shape=bin_im.shape)
    for i in range(1, bin_im.shape[0] - kernel_w + 1 + 1):
        for j in range(1, bin_im.shape[1] - kernel_h + 1 + 1):
            a = bin_im[i - 1:i - 1 + kernel_w,
                j - 1:j - 1 + kernel_h]
            dilate_img[i, j] = np.max(a * kernel)  # 若“有重合”，则点乘后最大值为0
    return dilate_img


def img_open(bin_im, erope_k, dilate_k):
    '''开:先腐蚀再膨胀'''
    open_img = img_erode(bin_im, erope_k)
    open_img = img_dilate(open_img, dilate_k)
    return open_img


def img_close(bin_im, erope_k, dilate_k):
    '''闭:先膨胀再腐蚀'''
    close_img = img_dilate(bin_im, dilate_k)
    close_img = img_erode(close_img, erope_k)
    return close_img



def hit_and_miss(A, B1, B2):
    A_c = 1 - A
    A_c_B1 = img_erode(A, B1).astype(np.uint8)
    A_c_B2 = img_erode(A_c, B2).astype(np.uint8)
    return A_c_B1 & A_c_B2


def hit_and_miss_cv(A, B1, B2):
    A_c = 1 - A
    A_c_B1 = cv2.erode(A, B1).astype(np.uint8)
    A_c_B2 = cv2.erode(A_c, B2).astype(np.uint8)
    return A_c_B1 & A_c_B2


if __name__ == "__main__":
    img = cv2.imread('proj03-images/wirebond-mask.tif', 0)
    ret, bin_img = cv2.threshold(img, 127, 1, cv2.THRESH_BINARY)

    kernel1 = np.ones(shape=(3, 3))  # 结构元素，可以为非规则形状
    kernel2 = np.ones(shape=(3, 3))  # 同上


    erode_im = img_erode(bin_img, kernel1)
    dilate_im = img_dilate(bin_img, kernel1)
    open_im = img_open(bin_img, kernel1, kernel2)
    close_im = img_close(bin_img, kernel1, kernel2)

    # erode_im = cv2.erode(img, kernel2)
    # dilate_im = cv2.dilate(img, kernel2)
    # open_im = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel=kernel1)
    # close_im = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel=kernel1)


    kernel1 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]).astype(np.uint8)
    kernel2 = 1 - kernel1
    img_hit_and_miss = hit_and_miss(bin_img, kernel1, kernel2)
    #img_hit_and_miss = hit_and_miss_cv(bin_img, kernel1, kernel2)

    print(np.sum(img_hit_and_miss))


    plt.subplot(231)
    plt.title("(a)二值图像", y=-0.2)
    plt.imshow(bin_img, cmap="gray")
    plt.subplot(232)
    plt.title("(b)腐蚀运算", y=-0.2)
    plt.imshow(erode_im, cmap="gray")
    plt.subplot(233)
    plt.title("(c)膨胀运算", y=-0.2)
    plt.imshow(dilate_im, cmap="gray")
    plt.subplot(234)
    plt.title("(d)开运算", y=-0.2)
    plt.imshow(open_im, cmap="gray")
    plt.subplot(235)
    plt.title("(e)闭运算", y=-0.2)
    plt.imshow(close_im, cmap="gray")
    plt.subplot(236)
    plt.title("(f)击中-击不中", y=-0.2)
    plt.imshow(img_hit_and_miss, cmap="gray")
    plt.show()

二值图连通分量提取、区域标记

将原图中的一个非零点不断膨胀，然后与原图相交，得到连通分量，之后将连通分量从原图去掉；循环以上操作获取多个连通分量，直到原图全为 0：

$$X_{\mathrm{k}}=\left(\begin{array}{c}{X}_{k-1}\oplus B\end{array}\right)\cap A,k=1,2,3,\dots$$

$X_0$为标记点，可指定或随机选取。

import cv2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']


def extract_connected_components(src):
 # 提取连通分量
    num_pixels = []  # 保存连通分量中的像素个数
    img_pixels = []  # 保存对应连通分量的图像

    canvas = np.zeros(src.shape, dtype=np.uint8)  # 将每次提取的连通分量放到这里
    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 5))  # 3×3结构元做膨胀运算

    x, y = np.where(src > 0)  # 获取图像中值为1的点
    canvas[x[0]][y[0]] = 1  # 以第一个点为起始点

    while src.any():  # src 矩阵内部所有元素为零的时候，为False
        count = len(x)  # 求每个连通分量的像素点数
        tmp = canvas  # tmp暂存上一步的迭代操作

        src_dilate = cv2.dilate(canvas, kernel)  # 膨胀
        canvas = cv2.bitwise_and(src, src_dilate)  # 和原图求交集

        if (tmp == canvas).all():  # 判断两次迭代操作是否相等。相等的时候，一次提取结束；否则，继续迭代
            img_pixels.append(canvas)  # 将连通分量图片保存到列表
            src[canvas > 0] = 0  # 将该次的连通分量去掉
            canvas = np.zeros(src.shape, np.uint8)  # 将每次提取的连通分量放到这里

            x, y = np.where(src > 0)  # 获取图像中的前景像素点
            num_pixels.append(count - len(x))  # 将上次的点数 - 下次的点数 = 提取连通分量的像素点个数

            if len(x) != 0:  # src 还有前景像素的时候
                canvas[x[0]][y[0]] = 1  # 取起始点

    return num_pixels, img_pixels  # 返回每个连通分量个数和连通分量图像



def mark_connected_components(img,num):
    output = np.zeros((img.shape[0], img.shape[1], 3), np.uint8)
    for i in range(0, num):
        mask = dst[i] == 1
        output[:, :, 0][mask] = np.random.randint(0, 255)
        output[:, :, 1][mask] = np.random.randint(0, 255)
        output[:, :, 2][mask] = np.random.randint(0, 255)

    return output
img = cv2.imread("proj03-images/Chickenfilet with bones.tif", flags=0)


ret, img_bin = cv2.threshold(img, 200, 1, cv2.THRESH_BINARY)  # 二值化处理
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 5))  # 3×3结构元
img_erode = cv2.erode(img_bin, kernel)

num, dst = extract_connected_components(img_erode.copy())  # 传入拷贝图像
print(len(num))  # 打印连通分量的像素点个数

cv2.imshow('img', img_erode * 255)
cv2.waitKey()

output=mark_connected_components(img_erode,len(num))
cv2.imshow('output', output)
cv2.waitKey(0)


# # OpenCV
# num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(img_erode, connectivity=8)
#
#
# print('各连通量的点数：',np.sort(stats[1:,-1]))
#
# # 不同的连通域赋予不同的颜色
# output = np.zeros((img.shape[0], img.shape[1], 3), np.uint8)
# for i in range(1, num_labels):
#     mask = labels == i
#     output[:, :, 0][mask] = np.random.randint(0, 255)
#     output[:, :, 1][mask] = np.random.randint(0, 255)
#     output[:, :, 2][mask] = np.random.randint(0, 255)
# cv2.imshow('oginal', output)
# cv2.waitKey()
# cv2.destroyAllWindows()

对提取到的连通分量用不同颜色进行区域标记，得到最终结果：

原图	结果

二值图细化算法

使用击中击不中获取与结构元素相同形状的图像边缘像素，然后将其删除，通过迭代多
次最终得到单像素宽带骨架。对于结构元素组：
{ B } = { B 1 , B 2 , B 3 , … , B n } \{B\}=\left\{B_{1},B_{2},B_{3},\ldots,B_{n}\right\} {B}={B1​,B2​,B3​,…,Bn​}

对目标图像$A$使用击中-击不中细化：

$$A=A\otimes B_i=A-(A⊛B_i),i=1,2,...,n$$

如果在迭代过程中，$⊛$操作后$A$没有变化，则结束细化，输出细化结果。

import cv2

import numpy as np
def thin(img):
    hit1 = np.array([[0, 0, 0], [1, 1, 1], [1, 1, 1]])
    hit2 = np.array([[1, 0, 0], [1, 1, 0], [1, 1, 1]])
    hit3 = np.array([[0, 0, 0], [1, 1, 1], [1, 1, 1]])
    hit4 = np.array([[1, 1, 0], [1, 1, 0], [1, 1, 0]])
    hit5 = np.array([[1, 1, 1], [1, 1, 1], [0, 0, 0]])
    hit6 = np.array([[1, 1, 1], [0, 1, 1], [0, 0, 1]])
    hit7 = np.array([[0, 1, 1], [0, 1, 1], [0, 1, 1]])
    hit8 = np.array([[1, 1, 0], [1, 0, 0], [0, 0, 0]])
    kernels = [hit1, hit2, hit3, hit4, hit5, hit6, hit7, hit8]
    dst = img.copy()
    while True:
        tmp = dst
        for kernel in kernels:
            dst = dst - hit_and_miss(dst, kernel,  1-kernel)
        if (dst == tmp).all():
            return dst
        print(np.sum(dst != tmp))
img2 = cv2.imread("UTK.tif", 0)


ret2, img_bin2 = cv2.threshold(img2, 100, 1, cv2.THRESH_BINARY)  # 二值化处理



img_thin = thin(img_bin2)
cv2.imshow('img', img_bin2 * 255)
cv2.imshow('img_thin', img_thin * 255)
cv2.waitKey(0)

原图	结果

灰度图形态学处理

灰度图形态学基本算子

腐蚀
在灰度图像中，用结构元素$b(x,y)$对输入图像$f(x,y)$进行灰度腐蚀运算可表示为：

( f ⊙ b ) ( s , t ) = min ⁡ { f ( s + x , t + y ) ∣ ( s + x , t + y ) ∈ D f ; ( x , y ) ∈ D b } (f\odot b)(s,t)=\min\left\{f(s+x,t+y)|(s+x,t+y)\in D_{f};(x,y)\in D_{b}\right\} (f⊙b)(s,t)=min{f(s+x,t+y)∣(s+x,t+y)∈Df​;(x,y)∈Db​}

式中，$D_f$和$D_b$分别是$f(x,y)$和$b(x,y)$的定义域。求某点的腐蚀运算就是计算该点结构元素覆盖范围内各点最小值作为该点的腐蚀结果。 经腐蚀运算后，图像暗的区域范围将扩大。

膨胀
用结构元素$b(x,y)$对输入图像$f(x,y)$进行灰度膨胀运算可表示为：

( f ⊕ b ) ( s , t ) = max ⁡ { f ( s − x , t − y ) ∣ ( s − x , t − y ) ∈ D f ; ( x , y ) ∈ D b } (f\oplus b)(s,t)=\max\left\{f(s-x,t-y)|(s-x,t-y)\in D_f;(x,y)\in D_b\right\} (f⊕b)(s,t)=max{f(s−x,t−y)∣(s−x,t−y)∈Df​;(x,y)∈Db​}

求某点的膨胀运算就是计算该点局部范围内各点的最大值作为该点的膨胀结果。经膨胀运算后，图像亮的区域范围将扩大。

开
开运算是先腐蚀后膨胀：
$$f^{\text{॰}}b=(f\odot b)\oplus b$$

开运算操作可以消除相对于结构元素尺寸较小的亮细节。

闭
闭运算是先膨胀后腐蚀：

$$f\bullet b=(f\oplus b)\odot b$$

闭运算操作可以消除相对于结构元素尺寸较小的暗细节。

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']


def img_erode(image, kernel, num=2):
    iChange = image.copy()
    num = int((kernel.shape[1] - 1) / 2)
    image = cv2.copyMakeBorder(image, num, num, num, num, 0)
    w = iChange.shape[1]
    h = iChange.shape[0]
    for i in range(h):
        for j in range(w):
            a = image[i:i + 2 * num+1, j:j + 2 * num+1]  # - kernel
            iChange[i, j] = np.min(a)
    return iChange


def img_dilate(image, kernel, num=2):
    iChange = image.copy()
    num = int((kernel.shape[1] - 1) / 2)
    image = cv2.copyMakeBorder(image, num, num, num, num, 0)
    w = iChange.shape[1]
    h = iChange.shape[0]
    for i in range(h):
        for j in range(w):
            a = image[i:i + 2 * num+1, j:j + 2 * num+1]  # - kernel
            iChange[i, j] = np.max(a)
    return iChange


def img_open(img, erope_k, dilate_k):
    '''开:先腐蚀再膨胀'''
    open_img = img_erode(img, erope_k)
    open_img = img_dilate(open_img, dilate_k)
    return open_img


def img_close(img, erope_k, dilate_k):
    '''闭:先膨胀再腐蚀'''
    close_img = img_dilate(img, dilate_k)
    close_img = img_erode(close_img, erope_k)
    return close_img


if __name__ == '__main__':
    img = cv2.imread("ckt_board_section.tif", flags=0)


    kernel1 = np.ones(shape=(5, 5))  # 此处声明结构元素，可以自定义为非规则形状
    kernel2 = np.ones(shape=(5, 5))  # 同上

    erode_im = img_erode(img, kernel1)
    dilate_im = img_dilate(img, kernel1)
    open_im = img_open(img, kernel1, kernel2)
    close_im = img_close(img, kernel1, kernel2)

    # erode_im = cv2.erode(img, kernel2)
    # dilate_im = cv2.dilate(img, kernel2)
    # open_im1 = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel=kernel1)
    # close_im = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel=kernel1)

    plt.subplot(231)
    plt.title("(a)二值图像", y=-0.2)
    plt.imshow(img, cmap="gray")
    plt.subplot(232)
    plt.title("(b)腐蚀运算", y=-0.2)
    plt.imshow(erode_im, cmap="gray")
    plt.subplot(233)
    plt.title("(c)膨胀运算", y=-0.2)
    plt.imshow(dilate_im, cmap="gray")
    plt.subplot(234)
    plt.title("(d)开运算", y=-0.2)
    plt.imshow(open_im, cmap="gray")
    plt.subplot(235)
    plt.title("(e)闭运算", y=-0.2)
    plt.imshow(close_im, cmap="gray")
    plt.show()

灰度图形态学梯度

灰度图形态学梯度用膨胀运算结果减去腐蚀运算的结果，表达式为：

$$g=(f\oplus b)-(f\odot b)$$

可用于提取灰度图的边缘。

原图	结果

灰度图 tophat 算法

灰度图 tophat 算法用原图像减去开运算结果，表达式为：
$$h=f-f^{\circ }b$$

可用于提取出灰度图中较亮的细节。

原图	结果

import cv2
def grad(img, erope_k, dilate_k):
    '''形态学梯度'''
    dilate_img = img_dilate(img.copy(), dilate_k)
    erode_img = img_erode(img.copy(), erope_k)
    return dilate_img - erode_img


def tophat(img, erope_k, dilate_k):
    open_img = img_open(img, erope_k, dilate_k).astype(np.uint8)
    # open_img=cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel=dilate_k)
    return (img - open_img)


if __name__ == '__main__':
    img1 = cv2.imread("headCT-Vandy.tif", 0)
    img2 = cv2.imread("rice_image_with_intensity_gradient.tif", 0)
    kernel3 = np.ones(shape=(5, 5))  # 此处声明结构元素，可以自定义为非规则形状
    kernel4 = np.ones(shape=(21, 21))  # 同上

    # 实现
    # grad_im = grad(img1, kernel3, kernel3)
    # tophat_im = tophat(img2, kernel4, kernel4)

    # OpenCV
    grad_im = cv2.morphologyEx(img1, cv2.MORPH_GRADIENT, kernel3)
    tophat_im = cv2.morphologyEx(img2, cv2.MORPH_TOPHAT, kernel4)

    cv2.imshow('grad_im', grad_im)
    cv2.waitKey(0)
    cv2.imshow('img2', img2)
    cv2.imshow('tophat_im', tophat_im)
    cv2.waitKey(0)