等离子体细丝图像分析代码报错修复及测量需求

需求概述 #

• 处理多张顶部被遮挡截断的等离子体细丝照片，实现：
  • 测量细丝的长度与最大厚度
  • 对遮挡处的截断部分进行插值连接，还原完整形态
• 基于开源代码修改后运行报错，需修复并完善功能

报错原因分析 #

核心报错来自代码最后可视化部分：

img1[v[0], v[1]] = [0, 0, 255]

• 原代码用cv2.IMREAD_GRAYSCALE读取单通道灰度图，却赋值三通道BGR颜色值，类型不匹配导致赋值失败
• 合并路径时未处理重复点，干扰后续计算与可视化

修复后完整代码 #

# -*- coding: utf-8 -*-
"""
等离子体细丝长度、厚度测量及遮挡插值修复
"""

import cv2
import numpy as np
from skimage import morphology, graph, measure
from skan import Skeleton
from scipy.interpolate import interp1d

MAX_JUNCTION = 4    # 最大连接点尺寸
MAX_ANGLE = 80      # 连接点允许的最大角度偏差
DELTA = 3           # 端点到内部点的距离，用于估算端点方向
MIN_PATH_LENGTH = 100  # 过滤短路径阈值

def angle(v1, v2):
    """计算两个向量之间的夹角（度数）"""
    rad = np.arctan2(v2[1], v2[0]) - np.arctan2(v1[1], v1[0])
    return np.abs((np.rad2deg(rad) % 360) - 180)

def calculate_fiber_thickness(binary_img, skeleton):
    """计算细丝的最大厚度（像素单位）"""
    # 距离变换：计算每个前景点到最近背景点的距离（半径）
    distance = cv2.distanceTransform(binary_img, cv2.DIST_L2, 5)
    # 提取骨架上的距离值，厚度为半径的2倍
    skeleton_thickness = distance[skeleton > 0]
    return np.max(skeleton_thickness) * 2 if len(skeleton_thickness) > 0 else 0

def interpolate_blocked_region(path, blocked_mask):
    """对遮挡区域的截断路径进行三次样条插值补全"""
    # 按y坐标排序（适配垂直方向的细丝）
    path_sorted = sorted(path, key=lambda p: p[0])
    y_coords = np.array([p[0] for p in path_sorted])
    x_coords = np.array([p[1] for p in path_sorted])
    
    # 筛选未被遮挡的点（mask=0表示未遮挡）
    unblocked_idx = [i for i, (y, x) in enumerate(path_sorted) if blocked_mask[y, x] == 0]
    if len(unblocked_idx) < 2:
        return path_sorted  # 点太少无法插值，返回原路径
    
    # 三次样条插值补全
    f = interp1d(y_coords[unblocked_idx], x_coords[unblocked_idx], kind='cubic', fill_value="extrapolate")
    full_y = np.arange(np.min(y_coords), np.max(y_coords)+1)
    full_x = f(full_y).astype(int)
    
    return list(zip(full_y, full_x))

# ---------------------- 主流程 ----------------------
# 读取彩色图像（避免单通道赋值三通道颜色的错误）
img = cv2.imread('DSC_7987.jpg')
gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 高斯模糊+Otsu二值化，反转图像使细丝为白色前景
blur = cv2.GaussianBlur(gray_img,(5,5),0)
ret, thresh = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
thresh = cv2.bitwise_not(thresh)

# 骨架提取与去噪
skeleton = morphology.skeletonize(thresh, method='lee')
skeleton = morphology.remove_small_objects(skeleton.astype(bool), MIN_PATH_LENGTH, connectivity=2)
skeleton = skeleton.astype(np.uint8) * 255

# 分割骨架为独立路径
g = Skeleton(skeleton)
lengths = np.array(g.path_lengths())
paths = [list(np.array(g.path_coordinates(i)).astype(int)) for i in range(g.n_paths) if lengths[i] > MIN_PATH_LENGTH]

# 提取端点及对应方向向量
endpoints = []
for i, p in enumerate(paths):
    if len(p) > DELTA:
        # 起点+方向向量（内部点指向起点）
        start_point = p[0]
        start_dir = np.subtract(p[DELTA], start_point)
        endpoints.append([start_point, start_dir, i])
        # 终点+方向向量（内部点指向终点）
        end_point = p[-1]
        end_dir = np.subtract(p[-1 - DELTA], end_point)
        endpoints.append([end_point, end_dir, i])

# 匹配可合并的截断端点对
angles = []
costs = np.where(skeleton > 0, 1, 255)  # 路径代价矩阵：骨架区域代价低

for i1 in range(len(endpoints)):
    for i2 in range(i1 + 1, len(endpoints)):
        e1, d1, p1 = endpoints[i1]
        e2, d2, p2 = endpoints[i2]
        if p1 != p2:
            # 计算端点间最短路径，判断是否属于同一连接点
            p_route, c_route = graph.route_through_array(costs, tuple(e1), tuple(e2))
            if c_route <= MAX_JUNCTION:
                # 计算方向夹角，判断是否为同一细丝的截断部分
                deg = angle(d1, d2)
                if deg <= MAX_ANGLE:
                    angles.append((deg, i1, i2, p_route))

# 按夹角从小到大合并路径（优先合并方向最一致的端点）
angles.sort(key=lambda a: a[0])

for deg, i1, i2, p_route in angles:
    e1, e2 = endpoints[i1], endpoints[i2]
    if e1 and e2:
        p1_idx, p2_idx = e1[2], e2[2]
        # 合并路径：路径1 + 连接路径 + 反转后的路径2（保证方向一致）
        merged_path = paths[p1_idx] + p_route + paths[p2_idx][::-1]
        # 去重，避免重复点干扰
        merged_path = list({tuple(p): p for p in merged_path}.values())
        paths[p1_idx] = merged_path
        # 更新其他端点的路径索引
        for i, e in enumerate(endpoints):
            if e and e[2] == p2_idx:
                e[2] = p1_idx
        # 标记已合并的路径和端点
        paths[p2_idx] = []
        endpoints[i1] = None
        endpoints[i2] = None

# ---------------------- 测量与插值处理 ----------------------
# 生成遮挡区域掩码（假设顶部白色区域为遮挡，可根据实际情况调整阈值）
blocked_mask = np.zeros_like(gray_img, dtype=np.uint8)
blocked_region = gray_img > 240  # 白色遮挡区域
blocked_mask[blocked_region] = 1

# 处理每条有效路径
for idx, path in enumerate(paths):
    if not path:
        continue
    
    # 1. 遮挡区域插值还原完整路径
    interpolated_path = interpolate_blocked_region(path, blocked_mask)
    
    # 2. 计算细丝长度（像素单位，可根据比例尺转换为物理单位）
    if len(interpolated_path) > len(path):
        # 插值后的路径重新计算欧氏距离总和
        points = np.array(interpolated_path)
        fiber_length = np.sum(np.sqrt(np.sum(np.diff(points, axis=0)**2, axis=1)))
    else:
        fiber_length = lengths[idx]
    
    # 3. 计算最大厚度
    max_thickness = calculate_fiber_thickness(thresh, skeleton)
    
    # 4. 可视化结果
    result_img = img.copy()
    # 绘制插值后的完整路径（红色线条）
    for i in range(len(interpolated_path)-1):
        cv2.line(result_img, tuple(interpolated_path[i][::-1]), tuple(interpolated_path[i+1][::-1]), (0,0,255), 2)
    # 标注测量结果
    cv2.putText(result_img, f"Length: {fiber_length:.2f} px", (20, 30), cv2.FONT_HERSHEY_SIMPLEX, 1, (0,255,0), 2)
    cv2.putText(result_img, f"Max Thickness: {max_thickness:.2f} px", (20, 70), cv2.FONT_HERSHEY_SIMPLEX, 1, (0,255,0), 2)
    
    cv2.imshow(f"Fiber {idx+1} - Result", result_img)
    cv2.waitKey(0)

cv2.destroyAllWindows()

关键修改与功能说明 #

1. 报错修复 #

• 将图像读取改为彩色图，后续转灰度处理，解决单通道图赋值三通道颜色的类型不匹配问题
• 合并路径时添加去重逻辑，避免重复点干扰可视化与计算

2. 遮挡区域插值 #

• 新增interpolate_blocked_region函数，使用三次样条插值补全遮挡区域的截断路径，还原细丝完整形态
• 通过遮挡掩码识别顶部白色遮挡区域，仅对该区域的截断部分进行插值

3. 长度与厚度测量 #

• 长度计算：利用skan的原生路径长度计算，合并/插值后的路径重新计算欧氏距离总和，保证精度
• 厚度计算：通过距离变换计算骨架上每个点对应的细丝半径，最大厚度为半径的2倍

4. 可视化优化 #

• 用红色线条绘制插值后的完整路径
• 在图像上直接标注长度和最大厚度的测量结果，直观展示分析数据

内容的提问来源于stack exchange，提问作者Jan K