python - 尝试裁剪图像并使用 pydicom 保存 dicom，

Question

尝试加载具有 JPEG2000 压缩的胸部 X 射线 DICOM 文件，提取像素阵列，裁剪它，然后另存为新的 DICOM 文件。在 Windows10 和 MacOS 机器上试过这个，但得到类似的错误。在conda环境下运行Python 3.6.13、GDCM 2.8.0、OpenJpeg 2.3.1、Pillow 8.1.2(先安装OpenJPEG和GDCM，再安装Pillow和Pydicom)。

我的初始代码:

file_list = [f.path for f in os.scandir(basepath)]
ds = pydicom.dcmread(file_list[0])
arr = ds.pixel_array
arr = arr[500:1500,500:1500]
ds.Rows = arr.shape[0]
ds.Columns = arr.shape[1]
ds.PixelData = arr.tobytes()
outputpath = os.path.join(basepath, "test.dcm")
ds.save_as(outputpath)

后续错误:ValueError: With tag (7fe0, 0010) got exception: (7FE0,0010) Pixel Data has an undefined length indicated that it's compressed,但数据未按要求封装。有关更多信息，请参见 pydicom.encaps.encapsulate()

然后我尝试将 ds.PixelData 行修改为 ds.PixelData = pydicom.encaps.encapsulate([arr.tobytes()]) 以创建 .dcm没有错误，但是当我打开 .dcm 进行查看时，它没有显示任何图像(全黑)。

我的下一次尝试是看看我是否需要以某种方式压缩回 JPEG2000，所以我尝试:

arr = Image.fromarray(arr)
output = io.BytesIO()
arr.save(output, format='JPEG2000')

但随后出现错误:OSError: encoder jpeg2k not available。我也试过 format='JPEG' 但它告诉我 OSError: cannot write mode I;16 as JPEG ...

非常感谢任何帮助!

Answer 1

能够通过使用 imagecodecs 库和 jpeg2k_encode 函数实现此功能。一个潜在的陷阱是您需要 .copy() 数组以满足函数的 C 连续要求，如果需要，您可以通过运行 arr_crop.flag 来确认这一点。这是最适合我的更新代码:

import os
import numpy as np
import matplotlib.pyplot as plt
import pydicom
from pydicom.encaps import encapsulate
from pydicom.uid import JPEG2000
from imagecodecs import jpeg2k_encode

file_list = [f.path for f in os.scandir(basepath)]
ds = pydicom.dcmread(file_list[0])
arr = ds.pixel_array
#Need to copy() to meet jpeg2k_encodes C contiguous requirement
arr_crop = arr[500:1500,500:1500].copy() 
# jpeg2k_encode to perform JPEG2000 compression
arr_jpeg2k = jpeg2k_encode(arr_crop)
# convert from bytearray to bytes before saving to PixelData
arr_jpeg2k = bytes(arr_jpeg2k)
ds.Rows = arr_crop.shape[0]
ds.Columns = arr_crop.shape[1]
ds.PixelData = encapsulate([arr_jpeg2k])
outputpath = os.path.join(basepath, "test.dcm")
ds.save_as(outputpath)

我最后还使用了 interactivecrop 包来相对快速地获得我需要的裁剪指数(一个提示，以防 future 的人在 jupyter 中尝试这个)。如果有帮助，这里有一段代码(在上面运行之前运行):

from interactivecrop.interactivecrop import main as crop
file_names = [os.path.split(f)[1].split(".")[0] for f in file_list]
image_list = []
for x in file_list:
    ds = pydicom.dcmread(x)
    arr = ds.pixel_array
    image_list.append(arr)
crop(image_list, file_names, optimize=True)
#After cropping all images, will print a dictionary
#copied and pasted this dictionary to a new cell as crop_dict
#use the below function to convert the output to actual indices
def convert_crop_to_index(fname, crop_dict):
    x = [crop_dict[fname][1], crop_dict[fname][1] + crop_dict[fname][3]]
    y = [crop_dict[fname][0], crop_dict[fname][0] + crop_dict[fname][2]]
    return x, y
arr_crop = arr[x[0]:x[1],y[0]:y[1]].copy()

一直无法弄清楚为什么 ds.decompress() 并保存解压后的 dicom 会生成全黑图像。我觉得这应该是最简单的方法，但上面的方法最终对我有用，所以我很高兴能够弄清楚。