zhang-suen算法提取裂缝图骨架

一、算法介绍

zhang-suen是一种应用较广的细化算法，也是骨架提取算法常用的算法。

初始条件： 输入是一个二值化图像，零点表示背景，非零值表示前景。算法的目标是细化前景物体的边缘。

迭代条件：Zhang-Suen 算法包含两个迭代步骤，通常会交替执行直到不再发生任何改变。这两个步骤是：

第一步，算法会遍历图像的每个像素，对于前景色将会检查是否满足一下条件。

• 当前像素为白色
• 当前像素P的 8 邻域上下、左右、左上、左下、右上、右下中至少有一个黑色像素
• 当前像素的 8 邻域中，白色像素数量在2到6之间

如果以上条件全部满足，那么将当前像素标记为黑色

第二步中，算法再次遍历图像的每个像素，对于前景色将会检查其是否满足以下条件：

• 当前像素为白色
• 当前像素的 8 邻域中，至少有一个是黑色
• 当前像素的 8 邻域中，白色像素数量在2到6之间
• 当前像素的 8 邻域的黑色像素不是连通的（即不是一个连通分量）

如果以上条件全部满足，那么将当前像素标记为黑色（背景）。

迭代终止条件： 迭代步骤1和步骤2将交替执行，直到不再有像素可以被删除。如果经过一轮迭代后，没有像素被删除，那么算法终止，否则继续下一轮迭代。

结果： 经过多轮迭代后，Zhang-Suen 算法将细化图像中白色目标对象的边缘，最终得到一个具有一像素宽度的骨架。

具体的实现可能因编程语言和图像处理库的不同而有所不同，但上述描述了 Zhang-Suen 算法的基本原理和步骤。

如果想要选择手写算法，可以参考下面两篇

激光条纹中心提取——Zhang-Suen法python

视觉组学习内容：Zhang-Suen骨架提取算法

二、skimage实现骨架提取

原先我也使用过手写的，但大多还是有逻辑错误，运行超时等情况，后面调研发现在skimage的morphology模块中实现了zhang-suen算法，我们需要保证我们的图像是二值图，这里用的是Otsu阈值方法来自动选择阈值，然后将图像二值化，将图像分成前景和背景。

from skimage.filters import threshold_otsu,median
from skimage.morphology import skeletonize,dilation,disk
from skimage import io, morphology

import matplotlib.pyplot as plt
plt.switch_backend('TkAgg')
# 使用Otsu阈值方法进行二值化处理

path = r"D:\PythonProject\RoadCrack\dimension2_data\num\001.png"

image = io.imread(path, as_gray=True)
thresh = threshold_otsu(image)
binary = image > thresh
skeleton = skeletonize(binary)
io.imshow(skeleton)
# io.imsave('output.png', skeleton)
io.show()

原图：

运行效果： 

可以看见图片上面还有一点瑕疵，可以使用一些图像预处理，帮助减小图像中的毛刺或噪声，同时改善图像的质量。

根据我对我自己数据集的测试，采用膨胀、中值滤波和闭运算等算法可以有效的减小图像中的瑕疵。

我使用的是裂缝森林数据集，你可以尝试其他的一些形态学方法，尽量早到适用于自己数据集的组合，下面是我经过测试获得的组合，大家可以借鉴。

from skimage.filters import threshold_otsu,median
from skimage.morphology import skeletonize,dilation,disk
from skimage import io, morphology

import matplotlib.pyplot as plt
plt.switch_backend('TkAgg')
# 使用Otsu阈值方法进行二值化处理

path = r"D:\PythonProject\RoadCrack\dimension2_data\num\001.png"

image = io.imread(path, as_gray=True)
thresh = threshold_otsu(image)
binary = image > thresh

binary = dilation(binary, disk(3))
binary = median(binary, selem=morphology.disk(5))
binary = dilation(binary, disk(2))
binary = median(binary, selem=morphology.disk(5))
# 添加闭运算
selem = morphology.disk(3)
binary = morphology.closing(binary, selem)

skeleton = skeletonize(binary)

io.imshow(skeleton)
# io.imsave('output.png', skeleton)
io.show()

三、数据集问题

你可以从网上可以找到开源数据集进行测试，我这里采用的是裂缝森林的数据集，你可以从Kaggle上找到crackforest | Kaggle。如果你下载的数据集是mat后缀的，可以使用下面这个脚本进行转换。

# Mat2png.py

from os.path import isdir
from scipy import io
import os, sys
import numpy as np
from PIL import Image

if __name__ == '__main__':
    file_path = './groundTruth/'
    png_img_dir = './groundTruthPngImg/'
    if not isdir(png_img_dir):
        os.makedirs(png_img_dir)
    image_path_lists = os.listdir(file_path)
    images_path = []
    for index in range(len(image_path_lists)):
        image_file = os.path.join(file_path, image_path_lists[index])
        # print(image_file)#./CrackForest-dataset-master/groundTruth/001.mat
        images_path.append(image_file)
        image_mat = io.loadmat(image_file)
        segmentation_image = image_mat['groundTruth']['Segmentation'][0]
        segmentation_image_array = np.array(segmentation_image[0])
        image = Image.fromarray((segmentation_image_array - 1) * 255)
        png_image_path = os.path.join(png_img_dir, "%s.png" % image_path_lists[index][0:3])
        image.save(png_image_path)

转换后的图片可能还存在下面的问题，我获得的裂缝森林数据集不是标准的二值图，如果不是细心，真没发现它的标签存在问题

 远看没有任何的问题，如果你使用图片查看器或ps放大查看就会看到下面的情况：

这种图片是有问题的，所以必须要进行二值化操作。当然这里也可以转灰度后使用阈值法：

std = 127.5
mask[mask > std] = 255
mask[mask < std] = 0
mask = mask.astype("uint8")