尝试使用 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（在安装 Pillow 和 Pydicom 之前先安装 OpenJPEG 和 GDCM）。

我的初始代码：

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 indicating that it's compressed, but the data isn't encapsulated as required. See pydicom.encaps.encapsulate() for more information

然后我尝试将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 包来相对快速地获得我需要的作物指数（一个提示，以防未来人们在 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 会生成全黑图像。我觉得这应该是最简单的方法，但上述方法最终对我有用，所以我很高兴能够弄清楚。