数字图像处理 第九章 形态学图像处理

一、预备知识

    集合B的反射$\hat{B}$，定义为$\hat{B}$$\left\{\begin{array}{c}w\mid w=-b,b\in B\end{array}\right\}$，即关于原集合原点对称。
    集合A平移到点$z=(z_1,z_2)$，表示为$(A{)}_z$，定义为$(A{)}_z=\left\{\begin{array}{c}c\mid c=a+z,a\in A\end{array}\right\}$。
    如下图所示：

二、腐蚀和膨胀

2.1腐蚀

    作为$Z^2$中的集合A和B，表示为$A\ominus B$的B对A的腐蚀义为
$$A\ominus B=\left\{\begin{array}{c}z\mid (B{)}_2\subseteq A\end{array}\right\}$$
表面上，该式指出B对A的腐蚀是一个用$z$平移的B包含在A中的所有的点$z$的集合。其表达的等价形式为：
$$A\ominus B=\left\{\begin{array}{c}z\mid (B{)}_2\cap A^c=\varnothing \end{array}\right\}$$
其中$A^c是A$补集，$\varnothing$是空集。

2.2膨胀

    A和B是$Z^2$中的集合，表示$A\oplus B$的B对A 的膨胀定义为
$$A\oplus B=\left\{\begin{array}{c}z\mid (\hat{B}{)}_z\cap A\ne \varnothing \end{array}\right\}$$
这个公式是以B关于它的原点的映像，并且以z对映像进行平移为基础的。B对A的膨胀是所有位移z的集合，这样，$\hat{B}$和A至少有一个元素是重叠的。其表达的等价形式为：
$$A\oplus B=\left\{\begin{array}{c}z\mid [(\hat{B}{)}_z\cap A]\subseteq A\end{array}\right\}$$

Python实现图像的腐蚀与膨胀

代码：

import cv2 as cv
import numpy as np

# 如果图像过大处理速度慢，可采用以下函数裁剪图像大小
def caijian(img):
    h=img.shape[0]
    w=img.shape[1]
    print("原图像大小为:",h,w)
    img1=cv.resize(img, (300, 300), interpolation=cv.INTER_CUBIC)
    h1=img1.shape[0]
    w1=img1.shape[1]
    print("裁剪后图像大小为:",h1,w1)
    return img1

def rgb2gray(img):
    h=img.shape[0]
    w=img.shape[1]
    img1=np.zeros((h,w),np.uint8)
    for i in range(h):
        for j in range(w):
            img1[i,j]=0.144*img[i,j,0]+0.587*img[i,j,1]+0.299*img[i,j,2]
    return img1

# 二值化
def otsu(img):
    h=img.shape[0]
    w=img.shape[1]
    m=h*w   # 图像像素点总和
    otsuimg=np.zeros((h,w),np.uint8)
    threshold_max=threshold=0   # 定义临时阈值和最终阈值
    histogram=np.zeros(256,np.int32)   # 初始化各灰度级个数统计参数
    probability=np.zeros(256,np.float32)   # 初始化各灰度级占图像中的分布的统计参数
    for i in range (h):
        for j in range (w):
            s=img[i,j]
            histogram[s]+=1   # 统计灰度级中每个像素在整幅图像中的个数
    for k in range (256):
        probability[k]=histogram[k]/m   # 统计每个灰度级占图像中的分布
    for i in range (255):
        w0 = w1 = 0   # 定义前景像素点和背景像素点灰度级占图像中的分布
        fgs = bgs = 0   # 定义前景像素点灰度级总和and背景像素点灰度级总和
        for j in range (256):
            if j<=i:   # 当前i为分割阈值
                w0+=probability[j]   # 前景像素点占整幅图像的比例累加
                fgs+=j*probability[j]
            else:
                w1+=probability[j]   # 背景像素点占整幅图像的比例累加
                bgs+=j*probability[j]
        u0=fgs/w0   # 前景像素点的平均灰度
        u1=bgs/w1   # 背景像素点的平均灰度
        g=w0*w1*(u0-u1)**2   # 类间方差
        if g>=threshold_max:
            threshold_max=g
            threshold=i
    print(threshold)
    for i in range (h):
        for j in range (w):
            if img[i,j]>threshold:
                otsuimg[i,j]=255
            else:
                otsuimg[i,j]=0
    return otsuimg

# 图像腐蚀
def etch(img):
    h=img.shape[0]
    w=img.shape[1]
    img1=np.zeros((h,w),np.uint8)
    for i in range (1,h-1):
        for j in range (1,w-1):
            min=img[i,j]
            for k in range (i-1,i+2):
                for l in range (j-1,j+2):
                    if k<0|k>=h-1|l<0|l>=w-1:
                        continue
                    if img[k,l]<min:
                        min=img[k,l]
            img1[i,j]=min
    return img1

# 图像膨胀
def expand(img):
    h=img.shape[0]
    w=img.shape[1]
    img1=np.zeros((h,w),np.uint8)
    for i in range (1,h-1):
        for j in range (1,w-1):
            max=img[i,j]
            for k in range (i-1,i+2):
                for l in range (j-1,j+2):
                    if k<0|k>=h-1|l<0|l>=w-1:
                        continue
                    if img[k,l]>max:
                       max=img[k,l]
            img1[i,j]=max
    return img1

image = cv.imread("JMU.jpg")   # 输入图像，根据实际情况而改变路径
caijianimage=caijian(image)   # 裁剪图像大小
grayimage = rgb2gray(caijianimage)   # 转灰
otsuimage = otsu(grayimage)   # 二值化
etchimage=etch(otsuimage)   # 腐蚀
expandimage=expand(otsuimage)   # 膨胀
cv.imshow("caijianimage",caijianimage)   # 输出裁剪图像相当于原图
cv.imshow("grayimage",grayimage)   # 输出灰度图
cv.imshow("otsuimage",otsuimage)   # 输出二值化图像
cv.imshow("etchimage",etchimage)   # 输出腐蚀图像
cv.imshow("expandimage",expandimage)   # 输出膨胀图像
cv.waitKey(0)
cv.destroyAllWindows()

（代码参考链接：http://t.csdn.cn/F37Ez）
结果：

    从左至右共五幅图，分别是原图，灰度图，二值化图，腐蚀图，膨胀图。
分析：
    腐蚀缩小或细化了二值图像中的物体，我们可以将腐蚀看成是形态学滤波操作，这种操作将小于结构元的图像细节从图像中去除了。
    与腐蚀不同，腐蚀是一种收缩或细化操作，膨胀则会增长或粗化二值图像中的物体。

2.3对偶性

    膨胀和腐蚀彼此关于集合求补运算和反射运算是对偶的，即
$$(A\ominus B{)}^c=A^c\oplus \hat{B}$$
和
$$(A\oplus B{)}^c=A^c\ominus \hat{B}$$
    上述式子指出，B对A的腐蚀是$\hat{B}$对$A^c$的膨胀的补。

三、开操作和闭操作

    结构元B对集合A的开操作，表示为$A\circ B$，其定义为$$A\circ B=(A\ominus B)\oplus B$$
因此，B对A的开操作就是B对A的腐蚀，紧接着用B对结果进行膨胀。
    类似地，用结构元B对集合A的闭操作，表示为$A\bullet B$,定义如下$$A\bullet B=(A\oplus B)\ominus B$$
上式说明，B对集合A的闭操作就是简单地用B对A膨胀，紧接着用B对结果进行腐蚀。

Python实现图像的开闭操作

代码：

import cv2 as cv

# 图像的开闭操作
def open_image(image):
    gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
    ret, binary = cv.threshold(gray, 0, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)
    cv.imshow("binaryzation", binary)
    kernel = cv.getStructuringElement(cv.MORPH_RECT, (5, 5))
    binary = cv.morphologyEx(binary, cv.MORPH_OPEN, kernel)
    cv.imshow("open", binary)

def close_image(image):
    gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
    ret, binary = cv.threshold(gray, 0, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)
    kernel = cv.getStructuringElement(cv.MORPH_RECT, (5, 5))
    binary = cv.morphologyEx(binary, cv.MORPH_CLOSE, kernel)
    cv.imshow("close", binary)

src = cv.imread("JMU.jpg")
cv.imshow("original", src)
open_image(src)
close_image(src)
cv.waitKey(0)
cv.destroyAllWindows()

（代码参考：http://t.csdn.cn/uUDGq）
结果：

分析：
    开操作一般会平滑物体的轮廓、断开较窄的的狭颈并消除细的突出物。闭操作同样也会平滑轮廓的一部分，但与开操作相反，它通常会弥合较窄的间断和细长的沟壑，消除小的孔洞，填补轮廓线中的断裂。

四、一些基本的形态学算法

4.1边界提取

    表示为$\beta (A)$的集合的边界可以通过先用B对A腐蚀，而后执行A和腐蚀的结果之间的集合之差得到，即
$$\beta (A)=A-(A\ominus B)$$
其中B是一个适当的结构元。

4.2孔洞填充

    孔洞填充定义为
$$X_k=(X_{k-1}\oplus B)\cap A^c,k=1,2,3,...$$
如果$X_k=X_{k-1}$，则算法在迭代的第k步结束。Xk 和A的并集包含被填充的集合和它的边界。
    条件膨胀：如果对上述公式的左部不加限制，则 上述公式的膨胀将填充整个区域。利用$A^c$的交集将 结果限制在感兴趣区域内，实现条件膨胀。

4.3连通分量的提取

    令Y表示一个包含于集合A中的连通分量，并假设Y 中的一个点p是已知的。用下列迭代式生成Y的所有元素：
$$X_k=(X_{k-1}\oplus B)\cap A,k=1,2,3,...$$
如果$X_k=X_{k-1}$，迭代过程结束，$X_k$包含输入图像中的所有的连通分量。

4.4凸壳

    如果连接集合A内任意两个点的直线段都在A的内部，则A是凸形的。集合S的凸壳H是包含S的最小凸集合，H-S称为S的凸缺。
    求取集合A的凸壳C(A)的简单形态学算法：
    令$B^i$表示4个结构元素,i=1,2,3,4
    $X_k^i=(X_{k1}⊛B^i\cup A)$,i =1,2,3,4  k=1,2,3,4,…
其中$X_0^i=A$。令$D^i=X_k^i$，当该过程收敛时，则A的凸壳为
$$C(A)=\bigcup _{i=1}^4{D}^i$$
    先对A用$B^1$运用击中或击不中变换，反复使用，当不 再发生变化时，执行与A的并集运算，用$D^1$表示结果。上述过程用$B^2$重复，直到不发生变化，最后得到的4个D的并集组成了A的凸壳。

4.5细化

    细化过程根据击中或击不中变换定义
$$A\otimes B=A-(A⊛B)=A\cap (A⊛B{)}^c$$
    定义结构元素序列为
$$\left\{\begin{array}{c}B\end{array}\right\}=\left\{\begin{array}{c}{B}^1,B^2,B^3,....B^n\end{array}\right\}$$
$B_i是B_{i-1}$旋转后的形式，如在$B_4$中旋转$90^{\circ }$
    用结构元素序列定义细化为
$$A\otimes \left\{\begin{array}{c}B\end{array}\right\}=((...((A\otimes B^1)\otimes B^2)...)\otimes B^n)$$
即连续使用$B^1,B^2,...,B^n$对A细化。

4.6粗化

    粗化和细化在形态学上是对偶过程，定义为
$$A\odot B=A\cup (A⊛B)$$
    用结构元素序列定义粗化为
$$A\odot \left\{\begin{array}{c}B\end{array}\right\}=((...((A\odot B^1)\odot B^2)...)\odot B^n)$$
    粗化可以通过细化算法求补集实现： 先对所讨论集合的背景进行细化，然后对结果求补集。