基于ITK Python包装器的μCT自动化处理：3D厚度图计算问题求助

我之前在处理μCT数据的中轴提取时也遇到过原生RegionalMaximaImageFilter的类似问题——默认的邻域设置和判断逻辑在3D场景下很容易丢掉细小的分支。下面是我摸索出来的两种解决方案，一种是纯ITK Python接口的自定义实现，另一种是结合NumPy加速的版本，都能更好地保留中轴分支：

核心思路 #

原生过滤器的问题通常出在邻域范围（默认可能只用6邻域，忽略了斜向的像素）和全局判断逻辑（没有严格限定只在物体内部比较）。我们需要自定义逻辑：

• 针对3D场景使用26邻域（所有相邻像素，包括斜向）
• 仅在二值化后的物体区域内判断局部最大值
• 允许灵活调整判断阈值（应对浮点距离变换的精度问题）

方案1：纯ITK Python自定义过滤器 #

这个方案完全基于ITK的迭代器接口，不需要额外依赖，适合需要严格遵循ITK pipeline的场景：

import itk

def custom_regional_maxima_3d(input_distance_image, input_binary_image):
    # 获取图像类型与基础信息
    ImageType = type(input_distance_image)
    region = input_distance_image.GetLargestPossibleRegion()
    size = region.GetSize()
    
    # 初始化输出图像，默认填充0
    output_image = ImageType.New()
    output_image.SetRegions(region)
    output_image.CopyInformation(input_distance_image)
    output_image.Allocate()
    output_image.FillBuffer(0)
    
    # 生成3D全邻域（26个方向）的偏移量，跳过中心像素
    offsets = []
    for dx in (-1, 0, 1):
        for dy in (-1, 0, 1):
            for dz in (-1, 0, 1):
                if dx == 0 and dy == 0 and dz == 0:
                    continue
                offsets.append(itk.Index[3]([dx, dy, dz]))
    
    # 创建迭代器遍历输入图像
    dist_iterator = itk.ImageRegionIteratorWithIndex[ImageType](input_distance_image, region)
    bin_iterator = itk.ImageRegionIteratorWithIndex[type(input_binary_image)](input_binary_image, region)
    
    while not dist_iterator.IsAtEnd():
        current_idx = dist_iterator.GetIndex()
        # 仅处理物体内部的像素（二值图中为1的区域）
        if bin_iterator.Get() == 1:
            current_val = dist_iterator.Get()
            is_local_max = True
            
            # 遍历所有邻域像素
            for offset in offsets:
                neighbor_idx = current_idx + offset
                # 检查邻域索引是否在图像范围内
                if (0 <= neighbor_idx[0] < size[0] and
                    0 <= neighbor_idx[1] < size[1] and
                    0 <= neighbor_idx[2] < size[2]):
                    # 切换到邻域位置，判断是否属于物体区域
                    bin_iterator.SetIndex(neighbor_idx)
                    if bin_iterator.Get() == 1:
                        neighbor_val = input_distance_image.GetPixel(neighbor_idx)
                        # 如果当前像素小于邻域像素，说明不是局部最大值
                        if current_val < neighbor_val - 1e-6:  # 加小阈值处理浮点精度
                            is_local_max = False
                            break
            
            # 标记局部最大值像素
            if is_local_max:
                output_image.SetPixel(current_idx, 1)
        
        # 移动迭代器到下一个像素
        dist_iterator.Next()
        bin_iterator.Next()
    
    return output_image

使用示例 #

假设你已经完成了二值化，得到binary_image（itk.Image[itk.UC, 3]类型），先计算距离变换，再用自定义函数提取中轴：

# 计算物体内部的距离变换（SignedDanielsson适合3D场景）
distance_filter = itk.SignedDanielssonDistanceMapImageFilter.New(Input=binary_image)
distance_filter.SetInsideIsPositive(True)
distance_filter.Update()
distance_image = distance_filter.GetOutput()

# 提取中轴（局部最大值）
skeleton_image = custom_regional_maxima_3d(distance_image, binary_image)

方案2：NumPy加速版本 #

如果你的μCT数据体积很大，纯ITK迭代器的速度会比较慢。结合NumPy和SciPy的滑动窗口函数能大幅提升效率：

import itk
import numpy as np
from scipy.ndimage import maximum_filter

def custom_regional_maxima_3d_numpy(distance_image, binary_image):
    # 将ITK图像转换为NumPy数组
    dist_np = itk.GetArrayFromImage(distance_image)
    bin_np = itk.GetArrayFromImage(binary_image)
    
    # 只保留物体区域的掩码
    object_mask = bin_np == 1
    
    # 用3x3x3窗口计算邻域最大值，背景区域设为负无穷避免干扰
    max_neighbor = maximum_filter(dist_np, size=3, mode='constant', cval=-np.inf)
    
    # 判断当前像素是否等于邻域最大值，且属于物体区域
    skeleton_np = np.where((object_mask) & (np.isclose(dist_np, max_neighbor, atol=1e-6)), 1, 0)
    
    # 将结果转回ITK图像，保留原图像的空间信息
    skeleton_image = itk.GetImageFromArray(skeleton_np)
    skeleton_image.CopyInformation(distance_image)
    
    return skeleton_image

关键优势 #

• 速度比纯ITK迭代器快10-100倍，适合处理GB级别的μCT数据
• np.isclose函数能灵活处理浮点距离变换的精度误差

额外提示 #

• 如果需要保留更细的分支，可以调整邻域大小（比如用5x5x5窗口，但会增加计算量）
• 若距离变换的输出是整数类型，可以去掉浮点精度阈值，直接用current_val >= neighbor_val判断

内容的提问来源于stack exchange，提问作者Thomas Janvier