【IP101】图像“瘦身魔法“详解：从基础细化到Zhang-Suen、Hilditch算法与中轴变换的完整代码实现-EW帮帮网

图像细化算法详解 🎨

欢迎来到图像处理的"瘦身工作室"！在这里，我们将学习如何像一位细心的雕刻家一样，将图像中的目标对象"瘦身"为单像素宽度的骨架。这个过程就像把一条粗线条逐渐"削瘦"成一条细线，同时保持其拓扑结构不变。让我们一起探索这个神奇的"瘦身魔法"吧！🚀

📚 目录

1. 算法原理
2. 应用场景
3. 基本细化算法
4. Hilditch细化算法
5. Zhang-Suen细化算法
6. 骨架提取
7. 中轴变换
8. 优化建议

1. 算法原理

图像细化的核心思想是通过迭代删除目标边缘像素来获得骨架。这个过程需要保证：

原则	说明	重要性
保持连通性	不能破坏目标的连接关系	⭐⭐⭐⭐⭐
适度细化	不能过度腐蚀导致目标消失	⭐⭐⭐⭐
中心定位	骨架应该位于物体的中心位置	⭐⭐⭐⭐

就像雕刻家精心雕琢木头一样，我们需要小心翼翼地"削掉"边缘像素，直到得到理想的骨架。🎯

2. 应用场景

图像细化算法在多个领域都有重要应用：🌟

应用领域	具体应用	技术要点
📝 字符识别	手写字符细化	特征提取和识别
👆 指纹识别	指纹骨架提取	指纹匹配和识别
🛣️ 道路提取	道路网络提取	地图制作和导航
🏥 医学图像	血管网络分析	疾病诊断辅助
🎯 模式识别	目标形状简化	特征匹配

3. 基本细化算法

3.1 基本原理

最基本的细化算法采用迭代的方式,每次迭代删除满足特定条件的边界点。判断一个点是否可以删除通常需要考虑其8邻域像素的分布情况。

💡 数学小贴士：像素邻域编号
$\begin{matrix} P_9 & P_2 & P_3 \\ P_8 & P_1 & P_4 \\ P_7 & P_6 & P_5 \end{matrix}$

每次迭代必须满足以下条件：

条件类型	具体条件	作用
边界点条件	当前点是边界点	确保只处理边缘
连通性条件	2 ≤ B(P1) ≤ 6	保持连通性
连续性条件	A(P1) = 1	避免过度细化
删除条件	P2 × P4 × P6 = 0 且 P4 × P6 × P8 = 0	保持结构完整

3.2 C++实现

void basic_thinning(const Mat& src, Mat& dst) {
    CV_Assert(!src.empty() && src.type() == CV_8UC1);

    src.copyTo(dst);
    bool has_changed;

    do {
        has_changed = false;
        Mat tmp = dst.clone();

        #pragma omp parallel for collapse(2)
        for (int y = 1; y < dst.rows - 1; y++) {
            for (int x = 1; x < dst.cols - 1; x++) {
                if (tmp.at<uchar>(y, x) == 0) continue;

                // 检查是否为边界点
                if (!is_boundary(tmp, y, x)) continue;

                // 计算P2到P9的值
                int p2 = tmp.at<uchar>(y-1, x) > 0;
                int p3 = tmp.at<uchar>(y-1, x+1) > 0;
                int p4 = tmp.at<uchar>(y, x+1) > 0;
                int p5 = tmp.at<uchar>(y+1, x+1) > 0;
                int p6 = tmp.at<uchar>(y+1, x) > 0;
                int p7 = tmp.at<uchar>(y+1, x-1) > 0;
                int p8 = tmp.at<uchar>(y, x-1) > 0;
                int p9 = tmp.at<uchar>(y-1, x-1) > 0;

                // 条件1：2 <= B(P1) <= 6
                int B = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;
                if (B < 2 || B > 6) continue;

                // 条件2：A(P1) = 1
                int A = count_transitions(tmp, y, x);
                if (A != 1) continue;

                // 条件3和4
                if ((p2 * p4 * p6 == 0) && (p4 * p6 * p8 == 0)) {
                    dst.at<uchar>(y, x) = 0;
                    has_changed = true;
                }
            }
        }
    } while (has_changed);
}

3.3 Python实现

def basic_thinning(img_path):
    """
    使用基本细化算法进行图像细化
    """
    # 读取图像
    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise ValueError(f"无法读取图像: {img_path}")

    # 二值化
    _, binary = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)

    # 转换为0和1格式
    skeleton = binary.copy() // 255
    changing = True

    def is_boundary(img, y, x):
        if img[y, x] == 0:
            return False
        for dy in [-1, 0, 1]:
            for dx in [-1, 0, 1]:
                if dy == 0 and dx == 0:
                    continue
                ny, nx = y + dy, x + dx
                if 0 <= ny < img.shape[0] and 0 <= nx < img.shape[1]:
                    if img[ny, nx] == 0:
                        return True
        return False

    def count_transitions(img, y, x):
        values = [
            img[y-1, x],   # P2
            img[y-1, x+1], # P3
            img[y, x+1],   # P4
            img[y+1, x+1], # P5
            img[y+1, x],   # P6
            img[y+1, x-1], # P7
            img[y, x-1],   # P8
            img[y-1, x-1], # P9
            img[y-1, x]    # P2
        ]
        count = 0
        for i in range(len(values)-1):
            if values[i] == 0 and values[i+1] == 1:
                count += 1
        return count

    while changing:
        changing = False
        temp = skeleton.copy()

        for y in range(1, skeleton.shape[0]-1):
            for x in range(1, skeleton.shape[1]-1):
                if temp[y, x] == 0:
                    continue

                if not is_boundary(temp, y, x):
                    continue

                # 计算P2到P9的值
                p2, p3, p4, p5, p6, p7, p8, p9 = (
                    temp[y-1, x], temp[y-1, x+1], temp[y, x+1], temp[y+1, x+1],
                    temp[y+1, x], temp[y+1, x-1], temp[y, x-1], temp[y-1, x-1]
                )

                # 条件1：2 <= B(P1) <= 6
                B = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
                if B < 2 or B > 6:
                    continue

                # 条件2：A(P1) = 1
                A = count_transitions(temp, y, x)
                if A != 1:
                    continue

                # 条件3和4
                if (p2 * p4 * p6 == 0) and (p4 * p6 * p8 == 0):
                    skeleton[y, x] = 0
                    changing = True

    # 转换回0-255格式
    result = skeleton.astype(np.uint8) * 255
    return result

4. Hilditch细化算法

4.1 基本原理

Hilditch细化算法是一种改进的细化算法，根据以下条件判断是否删除当前点：

连通性条件：2 ≤ B(P1) ≤ 6
连续性条件：A(P1) = 1
端点条件：P2 + P4 + P6 + P8 ≥ 1
删除条件：P2 × P4 × P6 = 0 且 P4 × P6 × P8 = 0

4.2 C++实现

void hilditch_thinning(const Mat& src, Mat& dst) {
    CV_Assert(!src.empty() && src.type() == CV_8UC1);

    src.copyTo(dst);
    bool has_changed;

    do {
        has_changed = false;
        Mat tmp = dst.clone();

        #pragma omp parallel for collapse(2)
        for (int y = 1; y < dst.rows - 1; y++) {
            for (int x = 1; x < dst.cols - 1; x++) {
                if (tmp.at<uchar>(y, x) == 0) continue;

                // 计算Hilditch算法条件
                int B = count_nonzero_neighbors(tmp, y, x);
                if (B < 2 || B > 6) continue;

                int A = count_transitions(tmp, y, x);
                if (A != 1) continue;

                // 计算连通性
                int conn = 0;
                if (tmp.at<uchar>(y-1, x) > 0 && tmp.at<uchar>(y-1, x+1) > 0) conn++;
                if (tmp.at<uchar>(y-1, x+1) > 0 && tmp.at<uchar>(y, x+1) > 0) conn++;
                if (tmp.at<uchar>(y, x+1) > 0 && tmp.at<uchar>(y+1, x+1) > 0) conn++;
                if (tmp.at<uchar>(y+1, x+1) > 0 && tmp.at<uchar>(y+1, x) > 0) conn++;
                if (tmp.at<uchar>(y+1, x) > 0 && tmp.at<uchar>(y+1, x-1) > 0) conn++;
                if (tmp.at<uchar>(y+1, x-1) > 0 && tmp.at<uchar>(y, x-1) > 0) conn++;
                if (tmp.at<uchar>(y, x-1) > 0 && tmp.at<uchar>(y-1, x-1) > 0) conn++;
                if (tmp.at<uchar>(y-1, x-1) > 0 && tmp.at<uchar>(y-1, x) > 0) conn++;

                if (conn == 1) {
                    dst.at<uchar>(y, x) = 0;
                    has_changed = true;
                }
            }
        }
    } while (has_changed);
}

4.3 Python实现

def hilditch_thinning(img_path):
    """
    使用Hilditch算法进行图像细化
    """
    # 读取图像
    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise ValueError(f"无法读取图像: {img_path}")

    # 二值化
    _, binary = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)

    # 转换为0和1格式
    skeleton = binary.copy() // 255
    changing = True

    def count_nonzero_neighbors(img, y, x):
        count = 0
        for dy in [-1, 0, 1]:
            for dx in [-1, 0, 1]:
                if dy == 0 and dx == 0:
                    continue
                ny, nx = y + dy, x + dx
                if 0 <= ny < img.shape[0] and 0 <= nx < img.shape[1]:
                    if img[ny, nx] > 0:
                        count += 1
        return count

    def count_transitions(img, y, x):
        values = [
            img[y-1, x],   # P2
            img[y-1, x+1], # P3
            img[y, x+1],   # P4
            img[y+1, x+1], # P5
            img[y+1, x],   # P6
            img[y+1, x-1], # P7
            img[y, x-1],   # P8
            img[y-1, x-1], # P9
            img[y-1, x]    # P2
        ]
        count = 0
        for i in range(len(values)-1):
            if values[i] == 0 and values[i+1] == 1:
                count += 1
        return count

    while changing:
        changing = False
        temp = skeleton.copy()

        for y in range(1, skeleton.shape[0]-1):
            for x in range(1, skeleton.shape[1]-1):
                if temp[y, x] == 0:
                    continue

                # 计算Hilditch算法条件
                B = count_nonzero_neighbors(temp, y, x)
                if B < 2 or B > 6:
                    continue

                A = count_transitions(temp, y, x)
                if A != 1:
                    continue

                # 计算连通性
                conn = 0
                if temp[y-1, x] > 0 and temp[y-1, x+1] > 0: conn += 1
                if temp[y-1, x+1] > 0 and temp[y, x+1] > 0: conn += 1
                if temp[y, x+1] > 0 and temp[y+1, x+1] > 0: conn += 1
                if temp[y+1, x+1] > 0 and temp[y+1, x] > 0: conn += 1
                if temp[y+1, x] > 0 and temp[y+1, x-1] > 0: conn += 1
                if temp[y+1, x-1] > 0 and temp[y, x-1] > 0: conn += 1
                if temp[y, x-1] > 0 and temp[y-1, x-1] > 0: conn += 1
                if temp[y-1, x-1] > 0 and temp[y-1, x] > 0: conn += 1

                if conn == 1:
                    skeleton[y, x] = 0
                    changing = True

    # 转换回0-255格式
    result = skeleton.astype(np.uint8) * 255
    return result

5. Zhang-Suen细化算法

5.1 基本原理

Zhang-Suen细化算法是一种改进的细化算法，通过两次迭代来细化图像：

第一次迭代条件：

2 ≤ B(P1) ≤ 6
A(P1) = 1
P2 × P4 × P6 = 0
P4 × P6 × P8 = 0

第二次迭代条件：

2 ≤ B(P1) ≤ 6
A(P1) = 1
P2 × P4 × P8 = 0
P2 × P6 × P8 = 0

5.2 C++实现

void zhang_suen_thinning(const Mat& src, Mat& dst) {
    CV_Assert(!src.empty() && src.type() == CV_8UC1);

    src.copyTo(dst);
    bool has_changed;

    do {
        has_changed = false;

        // 第一次迭代
        Mat tmp = dst.clone();
        #pragma omp parallel for collapse(2)
        for (int y = 1; y < dst.rows - 1; y++) {
            for (int x = 1; x < dst.cols - 1; x++) {
                if (tmp.at<uchar>(y, x) == 0) continue;

                int B = count_nonzero_neighbors(tmp, y, x);
                if (B < 2 || B > 6) continue;

                int A = count_transitions(tmp, y, x);
                if (A != 1) continue;

                // Zhang-Suen条件1
                if (tmp.at<uchar>(y-1, x) * tmp.at<uchar>(y, x+1) * tmp.at<uchar>(y+1, x) == 0 &&
                    tmp.at<uchar>(y, x+1) * tmp.at<uchar>(y+1, x) * tmp.at<uchar>(y, x-1) == 0) {
                    dst.at<uchar>(y, x) = 0;
                    has_changed = true;
                }
            }
        }

        // 第二次迭代
        tmp = dst.clone();
        #pragma omp parallel for collapse(2)
        for (int y = 1; y < dst.rows - 1; y++) {
            for (int x = 1; x < dst.cols - 1; x++) {
                if (tmp.at<uchar>(y, x) == 0) continue;

                int B = count_nonzero_neighbors(tmp, y, x);
                if (B < 2 || B > 6) continue;

                int A = count_transitions(tmp, y, x);
                if (A != 1) continue;

                // Zhang-Suen条件2
                if (tmp.at<uchar>(y-1, x) * tmp.at<uchar>(y, x+1) * tmp.at<uchar>(y, x-1) == 0 &&
                    tmp.at<uchar>(y-1, x) * tmp.at<uchar>(y+1, x) * tmp.at<uchar>(y, x-1) == 0) {
                    dst.at<uchar>(y, x) = 0;
                    has_changed = true;
                }
            }
        }
    } while (has_changed);
}

5.3 Python实现

def zhang_suen_thinning(img_path):
    """
    使用Zhang-Suen算法进行图像细化
    """
    # 读取图像
    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise ValueError(f"无法读取图像: {img_path}")

    # 二值化
    _, binary = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)

    # 转换为0和1格式
    skeleton = binary.copy() // 255

    def zhang_suen_iteration(img, iter_type):
        changing = False
        rows, cols = img.shape

        # 创建标记数组
        markers = np.zeros_like(img)

        for i in range(1, rows-1):
            for j in range(1, cols-1):
                if img[i,j] == 1:
                    # 获取8邻域
                    p2,p3,p4,p5,p6,p7,p8,p9 = (img[i-1,j], img[i-1,j+1], img[i,j+1],
                                              img[i+1,j+1], img[i+1,j], img[i+1,j-1],
                                              img[i,j-1], img[i-1,j-1])

                    # 计算条件
                    A = 0
                    for k in range(len([p2,p3,p4,p5,p6,p7,p8,p9])-1):
                        if [p2,p3,p4,p5,p6,p7,p8,p9][k] == 0 and [p2,p3,p4,p5,p6,p7,p8,p9][k+1] == 1:
                            A += 1
                    B = sum([p2,p3,p4,p5,p6,p7,p8,p9])

                    m1 = p2 * p4 * p6 if iter_type == 0 else p2 * p4 * p8
                    m2 = p4 * p6 * p8 if iter_type == 0 else p2 * p6 * p8

                    if (A == 1 and B >= 2 and B <= 6 and m1 == 0 and m2 == 0):
                        markers[i,j] = 1
                        changing = True

        img[markers == 1] = 0
        return img, changing

    # 迭代细化
    changing = True
    while changing:
        skeleton, changing1 = zhang_suen_iteration(skeleton, 0)
        skeleton, changing2 = zhang_suen_iteration(skeleton, 1)
        changing = changing1 or changing2

    # 转换回0-255格式
    result = skeleton.astype(np.uint8) * 255
    return result

6. 骨架提取

6.1 基本原理

骨架提取是一种特殊的细化算法，它通过形态学操作或距离变换来提取图像的中心线。骨架具有以下特点：

保持原图像的拓扑结构
位于物体的中心位置
宽度为单像素
保持物体的连通性

6.2 C++实现

void skeleton_extraction(const Mat& src, Mat& dst) {
    CV_Assert(!src.empty() && src.type() == CV_8UC1);

    // 使用距离变换和局部最大值提取骨架
    Mat dist;
    distanceTransform(src, dist, DIST_L2, DIST_MASK_PRECISE);

    dst = Mat::zeros(src.size(), CV_8UC1);

    #pragma omp parallel for collapse(2)
    for (int y = 1; y < src.rows - 1; y++) {
        for (int x = 1; x < src.cols - 1; x++) {
            if (src.at<uchar>(y, x) == 0) continue;

            // 检查是否为局部最大值
            float center = dist.at<float>(y, x);
            bool is_local_max = true;

            for (int dy = -1; dy <= 1 && is_local_max; dy++) {
                for (int dx = -1; dx <= 1; dx++) {
                    if (dy == 0 && dx == 0) continue;
                    if (dist.at<float>(y+dy, x+dx) > center) {
                        is_local_max = false;
                        break;
                    }
                }
            }

            if (is_local_max) {
                dst.at<uchar>(y, x) = 255;
            }
        }
    }
}

6.3 Python实现

def skeleton_extraction(img_path):
    """
    使用形态学操作提取图像骨架
    """
    # 读取图像
    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise ValueError(f"无法读取图像: {img_path}")

    # 二值化
    _, binary = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)

    # 创建结构元素
    kernel = cv2.getStructuringElement(cv2.MORPH_CROSS, (3,3))

    # 初始化骨架图像
    skeleton = np.zeros_like(binary)

    # 迭代提取骨架
    while True:
        # 形态学开运算
        eroded = cv2.erode(binary, kernel)
        opened = cv2.dilate(eroded, kernel)

        # 提取骨架点
        temp = cv2.subtract(binary, opened)

        # 更新骨架和二值图像
        skeleton = cv2.bitwise_or(skeleton, temp)
        binary = eroded.copy()

        # 当图像为空时停止迭代
        if cv2.countNonZero(binary) == 0:
            break

    return skeleton

7. 中轴变换

7.1 基本原理

中轴变换(Medial Axis Transform, MAT)是一种将二维形状转换为骨架的技术，其特点是：

骨架上的每个点都是到边界的最远点
保持了物体的拓扑结构
可以用于形状分析和描述
常用于计算机视觉和模式识别

7.2 C++实现

void medial_axis_transform(const Mat& src, Mat& dst, Mat& dist_transform) {
    CV_Assert(!src.empty() && src.type() == CV_8UC1);

    // 计算距离变换
    distanceTransform(src, dist_transform, DIST_L2, DIST_MASK_PRECISE);

    // 提取中轴
    dst = Mat::zeros(src.size(), CV_8UC1);

    #pragma omp parallel for
    for (int y = 1; y < src.rows - 1; y++) {
        for (int x = 1; x < src.cols - 1; x++) {
            if (src.at<uchar>(y, x) == 0) continue;

            float center = dist_transform.at<float>(y, x);
            bool is_medial_axis = false;

            // 检查梯度方向
            for (int dy = -1; dy <= 1; dy++) {
                for (int dx = -1; dx <= 1; dx++) {
                    if (dy == 0 && dx == 0) continue;
                    float neighbor = dist_transform.at<float>(y+dy, x+dx);

                    // 如果在相反方向上有相同的距离值，则为中轴点
                    if (abs(center - neighbor) < 1e-5) {
                        int opposite_y = y - dy;
                        int opposite_x = x - dx;
                        if (opposite_y >= 0 && opposite_y < src.rows &&
                            opposite_x >= 0 && opposite_x < src.cols) {
                            float opposite = dist_transform.at<float>(opposite_y, opposite_x);
                            if (abs(center - opposite) < 1e-5) {
                                is_medial_axis = true;
                                break;
                            }
                        }
                    }
                }
                if (is_medial_axis) break;
            }

            if (is_medial_axis) {
                dst.at<uchar>(y, x) = 255;
            }
        }
    }
}

7.3 Python实现

def medial_axis_transform(img_path):
    """
    计算图像的中轴变换
    """
    # 读取图像
    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise ValueError(f"无法读取图像: {img_path}")

    # 二值化
    _, binary = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)

    # 计算距离变换
    dist_transform = cv2.distanceTransform(binary, cv2.DIST_L2, 5)

    # 归一化距离变换结果
    cv2.normalize(dist_transform, dist_transform, 0, 255, cv2.NORM_MINMAX)

    # 提取局部最大值作为中轴点
    kernel = np.ones((3,3), dtype=np.uint8)
    dilated = cv2.dilate(dist_transform, kernel)
    medial_axis = (dist_transform == dilated) & (dist_transform > 20)

    # 转换为uint8类型
    result = medial_axis.astype(np.uint8) * 255

    # 转换为彩色图像
    result = cv2.cvtColor(result, cv2.COLOR_GRAY2BGR)

    return result

8. 优化建议

要提高细化算法的效果，可以考虑以下优化方向：🔧

8.1 预处理优化

使用中值滤波去除噪声
进行适当的二值化处理
填充小孔洞

8.2 算法选择

对于简单图像，使用基本细化算法
对于复杂图像，使用Zhang-Suen算法
需要保持拓扑结构时，选择Hilditch算法

8.3 后处理优化

去除毛刺
平滑骨架
修复断点

8.4 并行化处理

使用GPU加速
图像分块处理
多线程优化

🎯 总结

图像细化是一个既优雅又实用的算法，就像一位细心的雕刻家，它能够将复杂的图像简化为最本质的骨架结构。掌握这个算法，就像拥有了一把"瘦身魔法棒"，能够帮助我们更好地理解和分析图像！🎨✨

记住，好的细化效果需要：

合适的预处理
正确的算法选择
细致的参数调优
适当的后处理

让我们一起探索图像细化的奥秘，创造更多精彩的应用！🚀

【IP101】图像“瘦身魔法“详解：从基础细化到Zhang-Suen、Hilditch算法与中轴变换的完整代码实现