问题描述
我想在 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,x_iter] = np.min(dists_segm) # add minimum distance to distance transfrom stack
dtrans_neg_stack[t_iter,x_iter] = np.min(dists_neg_segm)
如果这有助于清除问题，这里是单个 zslice 的分割图像的图像 single z-slice of segmented image
解决方法
通常，我会用 scipy.ndimage.morphology.distance_transform_edt 做一些事情，但这给出了体素的各向同性大小的假设：
它没有这样的事情！您正在寻找 sampling= 参数。来自latest version of the docs：
元素沿每个维度的间距。如果是序列，则长度必须等于输入秩；如果是单个数字，则用于所有轴。如果未指定，则隐含网格间距为 1。
如果您将像素视为小方块/立方体，那么“采样”或“间距”这两个词可能有点神秘，这可能就是您错过它的原因。在大多数情况下，最好将像素视为网格上的点样本，样本之间的间距固定。我推荐 Alvy Ray 的 a pixel is not a little square 以更好地理解该术语。