图像距离变换不同 xyz 体素大小

Question

我想在二进制图像中找到每个体素到边界元素的最小距离，其中 z 体素大小与 xy 体素大小不同。也就是说，单个体素代表 225x110x110 (zyx) nm 的体积。

通常，我会用 scipy.ndimage.morphology.distance_transform_edt (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.morphology.distance_transform_edt.html) 做一些事情，但这给出了体素各向同性大小的假设：

dtrans_stack = np.zeros_like(segm_stack) # empty array to add to

### iterate over the t dimension and get distance transform
for t_iter in range(dtrans_stack.shape[0]):
    segm_ = segm_stack[t_iter, ...] # segmented image in single t
    neg_segm = np.ones_like(segm_) - segm_ # negative of the segmented image
    
    # get a ditance transform with isotropic voxel sizes
    dtrans_stack_iso = distance_transform_edt(segm_)
    dtrans_neg_stack_iso = -distance_transform_edt(neg_segm) # make distance in the segmented image negative
    dtrans_stack[t_iter, ...] = dtrans_stack_iso + dtrans_neg_stack_iso

我可以使用 scipy.spatial.distance.cdist (https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html) 蛮力做到这一点，但这需要很长时间，如果可以的话我宁愿避免

vox_multiplier = np.array([z_voxelsize, xy_voxelsize, xy_voxelsize]) # array of voxel sizes

## get a subset of coordinatess so I'm not wasting times in empty space
disk_size = 5 # size of disk for binary dilation
mip_tz = np.max(np.max(decon_stack, axis = 1), axis = 0)
thresh_li = threshold_li(mip_tz) # from from skimage.filters
mip_mask = mip_tz >= thresh_li
mip_mask = remove_small_objects(mip_mask) # from skimage.morphology
mip_dilated = binary_dilation(mip_mask, disk(disk_size)) # from skimage.morphology

# get the coordinates of the mask
coords = np.argwhere(mip_dilated == 1)
ycoords = coords[:, 0]
xcoords = coords[:, 1]

# get the lower and upper bounds of the xyz coordinates
ylb = np.min(ycoords)
yub = np.max(ycoords)
xlb = np.min(xcoords)
xub = np.max(xcoords)
zlb = 0
zub = zdims -1

# make zeros arrays of the proper size
dtrans_stack = np.zeros_like(segm_stack)
dtrans_stack_neg = np.zeros_like(segm_stack) # this will be the distance transform into the low inten area

for t_iter in range(dtrans_stack.shape[0]):
    segm_ = segm_stack[t_iter, ...]
    neg_segm_ = np.ones_like(segm_) - segm_ # negative of the segmented image
    
    # get the coordinats of segmented image and convert to nm 
    segm_coords = np.argwhere(segm_ == 1)
    segm_coords_nm = vox_multiplier * segm_coords
    neg_segm_coords = np.argwhere(neg_segm_ == 1)
    neg_segm_coords_nm = vox_multiplier * neg_segm_coords
    
    # make an empty arrays for the xy and z distance transforms
    dtrans_stack_x = np.zeros_like(segm_)
    dtrans_stack_y = np.zeros_like(segm_)
    dtrans_stack_z = np.zeros_like(segm_)
    dtrans_stack_neg_x = np.zeros_like(segm_)
    dtrans_stack_neg_y = np.zeros_like(segm_)
    dtrans_stack_neg_z = np.zeros_like(segm_)
    
    # iterate over the zyx and determine the minimum distance in nm from segmented image
    for z_iter in range(zlb, zub):
        for y_iter in range(ylb, yub):
            for x_iter in range(xlb, xub):
                coord_nm =  vox_multiplier* np.array([z_iter, y_iter, x_iter]) # change coords from pixel to nm
                coord_nm = coord_nm.reshape(1, 3) # reshape for distance calculateion
                dists_segm = distance.cdist(coord_nm, segm_coords_nm) # distance from the segmented image 
                dists_neg_segm = distance.cdist(coord_nm, neg_segm_coords_nm) # distance from the negative segmented image
                dtrans_stack[t_iter, z_iter, y_iter, x_iter] = np.min(dists_segm) # add minimum distance to distance transfrom stack
                dtrans_neg_stack[t_iter, z_iter, y_iter, x_iter] = np.min(dists_neg_segm)
 

这里是单个 zslice 分段图像的图像，如果这有助于清除问题 single z-slice of segmented image

Answer 1

通常，我会用scipy.ndimage.morphology.distance_transform_edt 做一些事情，但这给出了体素各向同性大小的假设：

它没有这样的事情！您正在寻找 sampling= 参数。来自latest version of the docs：

元素沿每个维度的间距。如果是一个序列，其长度必须等于输入等级；如果是单个数字，则用于所有轴。如果未指定，则隐含统一的网格间距。

如果您将像素视为小方块/立方体，“采样”或“间距”的措辞可能有点神秘，这可能就是您错过它的原因。在大多数情况下，最好将像素视为网格上的点样本，样本之间的间距固定。我推荐 Alvy Ray 的 a pixel is not a little square，以便更好地理解这个术语。