如何在预处理函数后保持图像和坐标（z 轴）之间的关联？
How can I keep the association between image and coordinate (z-axis) after a preprocess function?
我正在研究一个预处理函数，它将 DICOM 文件作为输入，returns 是 3D np.array（图像堆栈）。问题是我需要保持 ImagePositionPatient[2] 与输出数组中已处理图像的相对位置之间的关联。
例如，如果 ImagePositionPatient[2] == 5 的切片映射到 returned 堆栈中位置 3 的已处理切片，我需要 return 另一个具有 5在第三个位置，所有原始切片都相同。对于在通过插值或填充处理期间创建的切片，数组应包含一个占位符值，如 -99999。
我把我的代码贴在这里。
编辑：新的简化版本
def lung_segmentation(patient_dir):
"""
Load the dicom files of a patient, build a 3D image of the scan, normalize it to (1mm x 1mm x 1mm) and segment
the lungs
:param patient_dir: directory of dcm files
:return: a numpy array of size (384, 288, 384)
"""
""" LOAD THE IMAGE """
# Initialize image and get dcm files
dcm_list = glob(patient_dir + '/*.dcm')
img = np.zeros((len(dcm_list), 512, 512), dtype='float32') # inizializza un
# vettore di len(..) di matrici di 0 e di ampiezza 512x512
z = []
# For each dcm file, get the corresponding slice, normalize HU values, and store the Z position of the slice
for i, f in enumerate(dcm_list):
dcm = dicom.read_file(f)
img[i] = float(dcm.RescaleSlope) * dcm.pixel_array.astype('float32') + float(dcm.RescaleIntercept)
z.append(dcm.ImagePositionPatient[-1])
# Get spacing and reorder slices
spacing = list(map(float, dcm.PixelSpacing)) + [np.median(np.diff(np.sort(z)))]
print("LO SPACING e: "+str(spacing))
# spacing = list(map(lambda dcm, z: dcm.PixelSpacing + [np.median(np.diff(np.sort(z)))]))
img = img[np.argsort(z)]
""" NORMALIZE HU AND RESOLUTION """
# Clip and normalize
img = np.clip(img, -1024, 4000) # clippa con minimo a 1024 e max a 4k
img = (img + 1024.) / (4000 + 1024.)
# Rescale 1mm x 1mm x 1mm
new_shape = map(lambda x, y: int(x * y), img.shape, spacing[::-1])
old_shape = img.shape
img = resize(img, new_shape, preserve_range=True)
print('nuova shape calcolata'+ str(img.shape)+' con calcolo eseguito su img_shape: '+str(old_shape)+' * '+str(spacing[::-1]))
lungmask = np.zeros(img.shape) # WE NEED LUNGMASK FOR CODE BELOW
lungmask[int(img.shape[0]/2 - img.shape[0]/4) : int(img.shape[0]/2 + img.shape[0]/4),
int(img.shape[1]/2 - img.shape[1]/4) : int(img.shape[1]/2 + img.shape[1]/4),
int(img.shape[2]/2 - img.shape[2]/4) : int(img.shape[2]/2 + img.shape[2]/4)] = 1
# I set to value = 1 some pixel for executing code below, free to change
""" CENTER AND PAD TO GET SHAPE (384, 288, 384) """
# Center the image
sum_x = np.sum(lungmask, axis=(0, 1))
sum_y = np.sum(lungmask, axis=(0, 2))
sum_z = np.sum(lungmask, axis=(1, 2))
mx = np.nonzero(sum_x)[0][0]
Mx = len(sum_x) - np.nonzero(sum_x[::-1])[0][0]
my = np.nonzero(sum_y)[0][0]
My = len(sum_y) - np.nonzero(sum_y[::-1])[0][0]
mz = np.nonzero(sum_z)[0][0]
Mz = len(sum_z) - np.nonzero(sum_z[::-1])[0][0]
img = img * lungmask
img = img[mz:Mz, my:My, mx:Mx]
# Pad the image to (384, 288, 384)
nz, nr, nc = img.shape
pad1 = int((384 - nz) / 2)
pad2 = 384 - nz - pad1
pad3 = int((288 - nr) / 2)
pad4 = 288 - nr - pad3
pad5 = int((384 - nc) / 2)
pad6 = 384 - nc - pad5
# Crop images too big
if pad1 < 0:
img = img[:, -pad1:384 - pad2]
pad1 = pad2 = 0
if img.shape[0] == 383:
pad1 = 1
if pad3 < 0:
img = img[:, :, -pad3:288 - pad4]
pad3 = pad4 = 0
if img.shape[1] == 287:
pad3 = 1
if pad5 < 0:
img = img[:, :, -pad5:384 - pad6]
pad5 = pad6 = 0
if img.shape[2] == 383:
pad5 = 1
# Pad
img = np.pad(img, pad_width=((pad1 - 4, pad2 + 4), (pad3, pad4), (pad5, pad6)), mode='constant')
# The -4 / +4 is here for "historical" reasons, but it can be removed
return img
调整大小方法等的参考库是skimage
我会尝试至少给出一些答案的提示。正如评论中所讨论的那样，由于需要插值，调整大小可能会删除原始位置的已处理数据 - 所以最后你必须想出一个解决方案，或者通过将调整大小目标更改为实际的倍数分辨率，或通过 return 插值位置代替。
基本思想是让您的位置数组 z 进行变换，就像图像在 z 方向上一样。所以对于处理中改变处理图像的z位置的每个操作，都必须对z进行类似的操作。
假设您有 5 个切片，切片距离为 3mm：
>>> z
[0, 6, 3, 12, 9]
我们可以从中创建一个 numpy 数组以便于处理：
z_out = np.array(y)
这对应于未处理的 img 列表。
现在您对图像列表进行排序，因此您还必须对 z_out:
进行排序
img = img[np.argsort(z)]
z_out = np.sort(z_out)
>>> z_out
[0, 3, 6, 9, 12]
接下来，调整图像大小，引入插值切片。
我将在这里假设调整大小已完成，以便在调整大小时切片距离是目标分辨率的倍数。在这种情况下，您要计算插值切片的数量，并用相应的占位符值填充新的位置数组：
slice_distance = int((max(z) - min(z)) / (len(z) - 1))
nr_interpolated = slice_distance - 1 # you may adapt this to your algorithm
index_diff = np.arange(len(z) - 1) # used to adapt the insertion index
for i in range(nr_interpolated):
index = index_diff * (i + 1) + 1 # insertion index for placeholders
z_out = np.insert(z_out, index, -99999) # insert placeholder for interpolated positions
这为您提供了 z 数组，其中填充了图像数组中出现插值切片的占位符值：
>>> z_out
[0, -99999, -999999, 3, -99999, -999999, 6, -99999, -999999, 9, -99999, -999999, 12]
然后你必须在z方向上做与图像相同的填充：
img = np.pad(img, pad_width=((pad1 - 4, pad2 + 4), (pad3, pad4), (pad5, pad6)), mode='constant')
# use 'minimum' so that the placeholder is used
z_out = np.pad(z_out, pad_width=(pad1 - 4, pad2 + 4), mode='minimum')
为简单起见，假设填充值为 1 和 3，这将得到：
>>> z_out
[-99999, 0, -99999, -999999, 3, -99999, -999999, 6, -99999, -999999, 9, -99999, -999999, 12, -99999, -999999, -99999]
如果你在z方向上有更多的变换，你必须对z_out做相应的改变。如果你完成了，你可以return你的位置列表连同图像列表：
return img, z_out
顺便说一句：如果您的图像具有横向（轴向）方向，您的代码只会按预期工作，否则您必须从 [=23 计算 z 位置数组=] 和 Image Orientation Patient，而不是仅仅使用图像位置的 z 分量。
我正在研究一个预处理函数，它将 DICOM 文件作为输入，returns 是 3D np.array（图像堆栈）。问题是我需要保持 ImagePositionPatient[2] 与输出数组中已处理图像的相对位置之间的关联。
例如，如果 ImagePositionPatient[2] == 5 的切片映射到 returned 堆栈中位置 3 的已处理切片，我需要 return 另一个具有 5在第三个位置，所有原始切片都相同。对于在通过插值或填充处理期间创建的切片，数组应包含一个占位符值，如 -99999。
我把我的代码贴在这里。 编辑：新的简化版本
def lung_segmentation(patient_dir):
"""
Load the dicom files of a patient, build a 3D image of the scan, normalize it to (1mm x 1mm x 1mm) and segment
the lungs
:param patient_dir: directory of dcm files
:return: a numpy array of size (384, 288, 384)
"""
""" LOAD THE IMAGE """
# Initialize image and get dcm files
dcm_list = glob(patient_dir + '/*.dcm')
img = np.zeros((len(dcm_list), 512, 512), dtype='float32') # inizializza un
# vettore di len(..) di matrici di 0 e di ampiezza 512x512
z = []
# For each dcm file, get the corresponding slice, normalize HU values, and store the Z position of the slice
for i, f in enumerate(dcm_list):
dcm = dicom.read_file(f)
img[i] = float(dcm.RescaleSlope) * dcm.pixel_array.astype('float32') + float(dcm.RescaleIntercept)
z.append(dcm.ImagePositionPatient[-1])
# Get spacing and reorder slices
spacing = list(map(float, dcm.PixelSpacing)) + [np.median(np.diff(np.sort(z)))]
print("LO SPACING e: "+str(spacing))
# spacing = list(map(lambda dcm, z: dcm.PixelSpacing + [np.median(np.diff(np.sort(z)))]))
img = img[np.argsort(z)]
""" NORMALIZE HU AND RESOLUTION """
# Clip and normalize
img = np.clip(img, -1024, 4000) # clippa con minimo a 1024 e max a 4k
img = (img + 1024.) / (4000 + 1024.)
# Rescale 1mm x 1mm x 1mm
new_shape = map(lambda x, y: int(x * y), img.shape, spacing[::-1])
old_shape = img.shape
img = resize(img, new_shape, preserve_range=True)
print('nuova shape calcolata'+ str(img.shape)+' con calcolo eseguito su img_shape: '+str(old_shape)+' * '+str(spacing[::-1]))
lungmask = np.zeros(img.shape) # WE NEED LUNGMASK FOR CODE BELOW
lungmask[int(img.shape[0]/2 - img.shape[0]/4) : int(img.shape[0]/2 + img.shape[0]/4),
int(img.shape[1]/2 - img.shape[1]/4) : int(img.shape[1]/2 + img.shape[1]/4),
int(img.shape[2]/2 - img.shape[2]/4) : int(img.shape[2]/2 + img.shape[2]/4)] = 1
# I set to value = 1 some pixel for executing code below, free to change
""" CENTER AND PAD TO GET SHAPE (384, 288, 384) """
# Center the image
sum_x = np.sum(lungmask, axis=(0, 1))
sum_y = np.sum(lungmask, axis=(0, 2))
sum_z = np.sum(lungmask, axis=(1, 2))
mx = np.nonzero(sum_x)[0][0]
Mx = len(sum_x) - np.nonzero(sum_x[::-1])[0][0]
my = np.nonzero(sum_y)[0][0]
My = len(sum_y) - np.nonzero(sum_y[::-1])[0][0]
mz = np.nonzero(sum_z)[0][0]
Mz = len(sum_z) - np.nonzero(sum_z[::-1])[0][0]
img = img * lungmask
img = img[mz:Mz, my:My, mx:Mx]
# Pad the image to (384, 288, 384)
nz, nr, nc = img.shape
pad1 = int((384 - nz) / 2)
pad2 = 384 - nz - pad1
pad3 = int((288 - nr) / 2)
pad4 = 288 - nr - pad3
pad5 = int((384 - nc) / 2)
pad6 = 384 - nc - pad5
# Crop images too big
if pad1 < 0:
img = img[:, -pad1:384 - pad2]
pad1 = pad2 = 0
if img.shape[0] == 383:
pad1 = 1
if pad3 < 0:
img = img[:, :, -pad3:288 - pad4]
pad3 = pad4 = 0
if img.shape[1] == 287:
pad3 = 1
if pad5 < 0:
img = img[:, :, -pad5:384 - pad6]
pad5 = pad6 = 0
if img.shape[2] == 383:
pad5 = 1
# Pad
img = np.pad(img, pad_width=((pad1 - 4, pad2 + 4), (pad3, pad4), (pad5, pad6)), mode='constant')
# The -4 / +4 is here for "historical" reasons, but it can be removed
return img
调整大小方法等的参考库是skimage
我会尝试至少给出一些答案的提示。正如评论中所讨论的那样，由于需要插值，调整大小可能会删除原始位置的已处理数据 - 所以最后你必须想出一个解决方案，或者通过将调整大小目标更改为实际的倍数分辨率，或通过 return 插值位置代替。
基本思想是让您的位置数组 z 进行变换，就像图像在 z 方向上一样。所以对于处理中改变处理图像的z位置的每个操作，都必须对z进行类似的操作。
假设您有 5 个切片，切片距离为 3mm：
>>> z
[0, 6, 3, 12, 9]
我们可以从中创建一个 numpy 数组以便于处理：
z_out = np.array(y)
这对应于未处理的 img 列表。
现在您对图像列表进行排序，因此您还必须对 z_out:
img = img[np.argsort(z)]
z_out = np.sort(z_out)
>>> z_out
[0, 3, 6, 9, 12]
接下来，调整图像大小，引入插值切片。
我将在这里假设调整大小已完成，以便在调整大小时切片距离是目标分辨率的倍数。在这种情况下，您要计算插值切片的数量，并用相应的占位符值填充新的位置数组：
slice_distance = int((max(z) - min(z)) / (len(z) - 1))
nr_interpolated = slice_distance - 1 # you may adapt this to your algorithm
index_diff = np.arange(len(z) - 1) # used to adapt the insertion index
for i in range(nr_interpolated):
index = index_diff * (i + 1) + 1 # insertion index for placeholders
z_out = np.insert(z_out, index, -99999) # insert placeholder for interpolated positions
这为您提供了 z 数组，其中填充了图像数组中出现插值切片的占位符值：
>>> z_out
[0, -99999, -999999, 3, -99999, -999999, 6, -99999, -999999, 9, -99999, -999999, 12]
然后你必须在z方向上做与图像相同的填充：
img = np.pad(img, pad_width=((pad1 - 4, pad2 + 4), (pad3, pad4), (pad5, pad6)), mode='constant')
# use 'minimum' so that the placeholder is used
z_out = np.pad(z_out, pad_width=(pad1 - 4, pad2 + 4), mode='minimum')
为简单起见，假设填充值为 1 和 3，这将得到：
>>> z_out
[-99999, 0, -99999, -999999, 3, -99999, -999999, 6, -99999, -999999, 9, -99999, -999999, 12, -99999, -999999, -99999]
如果你在z方向上有更多的变换，你必须对z_out做相应的改变。如果你完成了，你可以return你的位置列表连同图像列表：
return img, z_out
顺便说一句：如果您的图像具有横向（轴向）方向，您的代码只会按预期工作，否则您必须从 [=23 计算 z 位置数组=] 和 Image Orientation Patient，而不是仅仅使用图像位置的 z 分量。