当前位置：首页 > news >正文

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

news 2025/5/18 14:31:59

图像细化算法详解 🎨

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

📚 目录

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) <= 6int B = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;if (B < 2 || B > 6) continue;// 条件2：A(P1) = 1int A = count_transitions(tmp, y, x);if (A != 1) continue;// 条件3和4if ((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() // 255changing = Truedef is_boundary(img, y, x):if img[y, x] == 0:return Falsefor dy in [-1, 0, 1]:for dx in [-1, 0, 1]:if dy == 0 and dx == 0:continueny, nx = y + dy, x + dxif 0 <= ny < img.shape[0] and 0 <= nx < img.shape[1]:if img[ny, nx] == 0:return Truereturn Falsedef count_transitions(img, y, x):values = [img[y-1, x],   # P2img[y-1, x+1], # P3img[y, x+1],   # P4img[y+1, x+1], # P5img[y+1, x],   # P6img[y+1, x-1], # P7img[y, x-1],   # P8img[y-1, x-1], # P9img[y-1, x]    # P2]count = 0for i in range(len(values)-1):if values[i] == 0 and values[i+1] == 1:count += 1return countwhile changing:changing = Falsetemp = 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:continueif 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) <= 6B = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9if B < 2 or B > 6:continue# 条件2：A(P1) = 1A = count_transitions(temp, y, x)if A != 1:continue# 条件3和4if (p2 * p4 * p6 == 0) and (p4 * p6 * p8 == 0):skeleton[y, x] = 0changing = True# 转换回0-255格式result = skeleton.astype(np.uint8) * 255return 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() // 255changing = Truedef count_nonzero_neighbors(img, y, x):count = 0for dy in [-1, 0, 1]:for dx in [-1, 0, 1]:if dy == 0 and dx == 0:continueny, nx = y + dy, x + dxif 0 <= ny < img.shape[0] and 0 <= nx < img.shape[1]:if img[ny, nx] > 0:count += 1return countdef count_transitions(img, y, x):values = [img[y-1, x],   # P2img[y-1, x+1], # P3img[y, x+1],   # P4img[y+1, x+1], # P5img[y+1, x],   # P6img[y+1, x-1], # P7img[y, x-1],   # P8img[y-1, x-1], # P9img[y-1, x]    # P2]count = 0for i in range(len(values)-1):if values[i] == 0 and values[i+1] == 1:count += 1return countwhile changing:changing = Falsetemp = 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:continueA = count_transitions(temp, y, x)if A != 1:continue# 计算连通性conn = 0if temp[y-1, x] > 0 and temp[y-1, x+1] > 0: conn += 1if temp[y-1, x+1] > 0 and temp[y, x+1] > 0: conn += 1if temp[y, x+1] > 0 and temp[y+1, x+1] > 0: conn += 1if temp[y+1, x+1] > 0 and temp[y+1, x] > 0: conn += 1if temp[y+1, x] > 0 and temp[y+1, x-1] > 0: conn += 1if temp[y+1, x-1] > 0 and temp[y, x-1] > 0: conn += 1if temp[y, x-1] > 0 and temp[y-1, x-1] > 0: conn += 1if temp[y-1, x-1] > 0 and temp[y-1, x] > 0: conn += 1if conn == 1:skeleton[y, x] = 0changing = True# 转换回0-255格式result = skeleton.astype(np.uint8) * 255return 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条件1if (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条件2if (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() // 255def zhang_suen_iteration(img, iter_type):changing = Falserows, 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 = 0for 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 += 1B = sum([p2,p3,p4,p5,p6,p7,p8,p9])m1 = p2 * p4 * p6 if iter_type == 0 else p2 * p4 * p8m2 = p4 * p6 * p8 if iter_type == 0 else p2 * p6 * p8if (A == 1 and B >= 2 and B <= 6 and m1 == 0 and m2 == 0):markers[i,j] = 1changing = Trueimg[markers == 1] = 0return img, changing# 迭代细化changing = Truewhile 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) * 255return 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:breakreturn 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 forfor (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