1

我正在开发一个预处理函数，该函数将 DICOM 文件作为输入并返回 3D np.array（图像堆栈）。问题是我需要保持ImagePositionPatient[2]输出数组中已处理图像之间的关联和相对位置。

例如，如果一个切片ImagePositionPatient[2] == 5映射到返回堆栈中位置 3 的已处理切片，我需要返回另一个数组，该数组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

4

1 回答 1

1

我将尝试至少对答案给出一些提示。正如评论中所讨论的那样，由于需要插值，调整大小可能会删除原始位置的处理数据 - 所以最后你必须想出一个解决方案，或者通过将调整大小目标更改为实际的倍数分辨率，或者改为返回插值位置。

基本思想是让您的位置数组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 img, z_out

顺便说一句：如果您的图像具有横向（轴向）方向，您的代码只会按预期工作，否则您必须从Image Position Patientand计算 z 位置数组Image Orientation Patient，而不仅仅是使用图像位置的 z 分量。

于 2020-06-07T11:02:43.787 回答