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

如何解决尝试裁剪图像并使用 pydicom 保存 dicom，

尝试加载具有 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 ...

非常感谢任何帮助！

解决方法

能够通过使用 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 会生成全黑图像。我觉得这应该是最简单的方法，但上述方法最终对我有用，所以我很高兴能够弄清楚。