小胶质细胞骨架分支重连与背景去除技术问题咨询

1. 优化骨架分支重连效果 #

自定义重连函数失效通常是因为只判断距离，忽略了端点方向一致性和局部骨架连续性。可以结合骨架端点的梯度方向匹配+距离阈值实现精准重连：

• 步骤1：提取骨架端点。用skimage.morphology.skeletonize生成骨架后，通过邻域像素计数筛选端点（端点的3x3邻域内仅包含1个其他骨架像素）。
• 步骤2：计算端点对的距离与方向向量。方向向量通过端点邻域的相邻骨架像素位置推导，确保分支方向对齐。
• 步骤3：仅连接满足以下条件的端点对：
  • 欧氏距离小于10-20像素（根据图像分辨率调整）
  • 方向向量夹角小于30°（余弦相似度大于0.866）
• 代码示例：

import numpy as np
from skimage.morphology import skeletonize
from scipy.ndimage import convolve

# 假设binary_img是二值化后的图像
skeleton = skeletonize(binary_img)

# 提取骨架端点
def get_endpoints(skeleton):
    kernel = np.ones((3,3), dtype=np.uint8)
    neighbor_count = convolve(skeleton.astype(np.uint8), kernel, mode='same')
    # 骨架像素自身占1个计数，邻域内仅1个其他骨架像素时为端点
    endpoints = (skeleton == 1) & (neighbor_count == 2)
    return np.argwhere(endpoints)

endpoints = get_endpoints(skeleton)

# 遍历端点对并连接符合条件的分支
for i in range(len(endpoints)):
    for j in range(i+1, len(endpoints)):
        y1, x1 = endpoints[i]
        y2, x2 = endpoints[j]
        distance = np.sqrt((x2-x1)**2 + (y2-y1)**2)
        if distance > 15:
            continue
        
        # 计算端点方向向量
        def get_dir(y, x):
            neighbors = np.argwhere(skeleton[y-1:y+2, x-1:x+2] == 1)
            # 排除自身位置（邻域中心）
            dirs = neighbors - [1,1]
            return dirs[np.any(dirs != [0,0], axis=1)][0]
        
        dir1 = get_dir(y1, x1)
        dir2 = get_dir(y2, x2)
        cos_angle = np.dot(dir1, dir2) / (np.linalg.norm(dir1)*np.linalg.norm(dir2))
        
        if cos_angle > np.cos(np.radians(30)):
            # 用Bresenham算法绘制连接线
            dx, dy = abs(x2-x1), abs(y2-y1)
            sx = 1 if x1 < x2 else -1
            sy = 1 if y1 < y2 else -1
            err = dx - dy
            while x1 != x2 or y1 != y2:
                skeleton[y1, x1] = 1
                e2 = 2 * err
                if e2 > -dy:
                    err -= dy
                    x1 += sx
                if e2 < dx:
                    err += dx
                    y1 += sy
            skeleton[y2, x2] = 1

2. 去除未连接胞体的背景丝状物 #

核心思路是先定位胞体，再保留仅与胞体连通的结构：

• 步骤1：提取胞体区域。对二值图像做10x10圆盘的形态学闭运算，填充细分支后，用面积阈值（如>500像素）筛选出胞体掩码。
• 步骤2：形态学重建。以胞体掩码为种子，原始二值图像为掩膜，通过膨胀重建仅保留与胞体连通的区域。
• 步骤3：过滤骨架，仅保留重建区域内的骨架像素。
• 代码示例：

from skimage.morphology import closing, disk, reconstruction
from skimage.measure import label, regionprops

# 提取胞体掩码
closed = closing(binary_img, disk(10))
labeled_closed = label(closed)
cell_body_mask = np.zeros_like(closed)
for prop in regionprops(labeled_closed):
    if prop.area > 500:
        cell_body_mask[labeled_closed == prop.label] = 1

# 重建仅连接胞体的区域
connected_regions = reconstruction(cell_body_mask, binary_img, method='dilation')

# 生成过滤后的骨架
clean_skeleton = skeletonize(connected_regions)

3. 识别并保留连接大胞体的结构 #

通过连通分量分析+大胞体属性筛选实现：

• 步骤1：筛选大胞体。对胞体掩码做连通标记，保留面积>800像素的区域作为大胞体掩码。
• 步骤2：标记二值图像的所有连通分量，仅保留与大胞体掩码有重叠的分量。
• 步骤3：基于筛选后的分量生成最终骨架。
• 代码示例：

# 筛选大胞体
labeled_cells = label(cell_body_mask)
large_cell_mask = np.zeros_like(cell_body_mask)
for prop in regionprops(labeled_cells):
    if prop.area > 800:
        large_cell_mask[labeled_cells == prop.label] = 1

# 保留与大胞体连通的分量
labeled_binary = label(binary_img)
final_mask = np.zeros_like(binary_img)
for prop in regionprops(labeled_binary):
    overlap = np.sum(large_cell_mask[labeled_binary == prop.label])
    if overlap > 0:
        final_mask[labeled_binary == prop.label] = 1

# 最终骨架
final_skeleton = skeletonize(final_mask)

内容的提问来源于stack exchange，提问作者Beto Santiago