从 numpy 数组创建 pydicom 文件

Question

我正在尝试从标准大小（512 x 512 或 256 x 256）的 numpy 数组创建新的 dicom 图像。

import dicom, dicom.UID
from dicom.dataset import Dataset, FileDataset

def write_dicom(pixel_array,filename):
    
    file_meta = Dataset()
    ds = FileDataset(filename, {},file_meta = file_meta,preamble="\0"*128)
    ds.PixelData = pixel_array.tostring()
    ds.save_as(filename)
    return

if __name__ == "__main__":
    import numpy as np
    pixel_array = np.tile(np.arange(256).reshape(16,16), (16,16)) * 4
    write_dicom(pixel_array,'pretty.dcm')

Answer 1

2020 年更新 :)

这些答案都不适合我。这就是我最终保存有效的单色 16bpp MR 切片的结果，该切片至少在 Slicer、Radiant 和 MicroDicom 中正确显示：

import pydicom
from pydicom.dataset import Dataset, FileDataset
from pydicom.uid import ExplicitVRLittleEndian
import pydicom._storage_sopclass_uids

image2d = image2d.astype(np.uint16)

print("Setting file meta information...")
# Populate required values for file meta information

meta = pydicom.Dataset()
meta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.MRImageStorage
meta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid()
meta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian  

ds = Dataset()
ds.file_meta = meta

ds.is_little_endian = True
ds.is_implicit_VR = False

ds.SOPClassUID = pydicom._storage_sopclass_uids.MRImageStorage
ds.PatientName = "Test^Firstname"
ds.PatientID = "123456"

ds.Modality = "MR"
ds.SeriesInstanceUID = pydicom.uid.generate_uid()
ds.StudyInstanceUID = pydicom.uid.generate_uid()
ds.FrameOfReferenceUID = pydicom.uid.generate_uid()

ds.BitsStored = 16
ds.BitsAllocated = 16
ds.SamplesPerPixel = 1
ds.HighBit = 15

ds.ImagesInAcquisition = "1"

ds.Rows = image2d.shape[0]
ds.Columns = image2d.shape[1]
ds.InstanceNumber = 1

ds.ImagePositionPatient = r"0\0\1"
ds.ImageOrientationPatient = r"1\0\0\0\-1\0"
ds.ImageType = r"ORIGINAL\PRIMARY\AXIAL"

ds.RescaleIntercept = "0"
ds.RescaleSlope = "1"
ds.PixelSpacing = r"1\1"
ds.PhotometricInterpretation = "MONOCHROME2"
ds.PixelRepresentation = 1

pydicom.dataset.validate_file_meta(ds.file_meta, enforce_standard=True)

print("Setting pixel data...")
ds.PixelData = image2d.tobytes()

ds.save_as(r"out.dcm")

注意以下几点：

• 按照 PyDicom 文档的建议，通过 FileDataset 构造函数未能为我创建有效的标头
• validate_file_meta 会为您在标题中创建一些缺失的元素（版本）
• 您需要指定字节顺序和显式/隐式 VR 两次：/
• 只要您相应地更新每个切片的 ImagePositionPatient 和 InstanceNumber，此方法将允许您创建有效的卷
• 确保将您的 numpy 数组转换为与您的 BitsStored 具有相同位数的数据格式

Answer 2

我能够进一步减少@Corvin 的好答案。这是一个极简代码示例，允许将（虚拟）3D numpy 数组保存到可以使用Amide 打开的有效 DICOM 图像：

#!/usr/bin/python3
import numpy
import pydicom
import pydicom._storage_sopclass_uids

# dummy image
image = numpy.random.randint(2**16, size=(512, 512, 512), dtype=numpy.uint16)

# metadata
fileMeta = pydicom.Dataset()
fileMeta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.CTImageStorage
fileMeta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid()
fileMeta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian

# dataset
ds = pydicom.Dataset()
ds.file_meta = fileMeta

ds.Rows = image.shape[0]
ds.Columns = image.shape[1]
ds.NumberOfFrames = image.shape[2]

ds.PixelSpacing = [1, 1] # in mm
ds.SliceThickness = 1 # in mm

ds.BitsAllocated = 16
ds.PixelRepresentation = 1
ds.PixelData = image.tobytes()

# save
ds.save_as('image.dcm', write_like_original=False)

正如人们可能观察到的，如果将输出image.dcm 文件传递​​给dciodvfy，则会丢失很多字段。这些字段的填写留给读者；）

Answer 3

DICOM 是一种非常复杂的格式。有很多方言，兼容性是一个运气问题。你也可以试试nibabel，也许它的方言对 RadiAnt 或 MicroDicom 更有吸引力。

一般来说，我建议尽可能使用 Nifti 格式。它的标准更加简洁，不兼容的情况很少见。 nibabel 也支持这个。

Answer 4

为需要它的人提供一个工作配置和一个问题。 问题在另一个线程Create a Dicom from multiple jpg images 对我有用的是没有压缩的灰度。每次压缩尝试都失败了，我不知道为什么：

# Populate required values for file meta information
meta = pydicom.Dataset()
meta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian
meta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.MRImageStorage
meta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid()

# build dataset
ds = Dataset()
ds.file_meta = meta
ds.fix_meta_info()

# unknown options
ds.is_little_endian = True
ds.is_implicit_VR = False
ds.SOPClassUID = pydicom._storage_sopclass_uids.MRImageStorage
ds.SeriesInstanceUID = pydicom.uid.generate_uid()
ds.StudyInstanceUID = pydicom.uid.generate_uid()
ds.FrameOfReferenceUID = pydicom.uid.generate_uid()
ds.BitsStored = 16
ds.BitsAllocated = 16
ds.SamplesPerPixel = 1
ds.HighBit = 15
ds.ImagesInAcquisition = "1"
ds.InstanceNumber = 1
ds.ImagePositionPatient = r"0\0\1"
ds.ImageOrientationPatient = r"1\0\0\0\-1\0"
ds.ImageType = r"ORIGINAL\PRIMARY\AXIAL"
ds.RescaleIntercept = "0"
ds.RescaleSlope = "1"
ds.PixelRepresentation = 1

# Case options
ds.PatientName = "Anonymous"
ds.PatientID = "123456"
ds.Modality = "MR"
ds.StudyDate = '20200225'
ds.ContentDate = '20200225'

# convert image to grayscale
img = Image.open(filename).convert('L')
img.save(filename)

# open image, decode and ensure_even stream
with open(filename, 'rb') as f:
    arr = decode(f)

def ensure_even(stream):
    # Very important for some viewers
    if len(stream) % 2:
        return stream + b"\x00"
    return stream

# required for pixel handler
ds.BitsStored = 8
ds.BitsAllocated = 8
ds.HighBit = 7
ds.PixelRepresentation = 0

# grayscale without compression WORKS
ds.PhotometricInterpretation = "MONOCHROME2"
ds.SamplesPerPixel = 1  # 1 color = 1 sample per pixel
ds.file_meta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian
ds.PixelData = ensure_even(arr.tobytes())

# JPEGBaseline compressed DOES NOT WORK
# ds.PixelData = encapsulate([ensure_even(arr.tobytes())])
# ds.PhotometricInterpretation = "YBR_FULL"
# ds.SamplesPerPixel = 3  # 3 colors = 3 sampleperpixel
# ds.file_meta.TransferSyntaxUID = pydicom.uid.JPEGBaseline
# ds.compress(pydicom.uid.JPEGBaseline)

# JPEGExtended compressed DOES NOT WORK
# ds.PixelData = encapsulate([ensure_even(arr.tobytes())])
# ds.PhotometricInterpretation = "YBR_FULL_422"
# ds.SamplesPerPixel = 3  # 3 colors = 3 sampleperpixel
# ds.file_meta.TransferSyntaxUID = pydicom.uid.JPEGExtended
# ds.compress(pydicom.uid.JPEGExtended)

# JPEG2000 compressed DOES NOT WORK
# ds.PhotometricInterpretation = "RGB"
# ds.SamplesPerPixel = 3  # 3 colors = 3 sampleperpixel
# ds.file_meta.TransferSyntaxUID = pydicom.uid.JPEG2000
# ds.PixelData = encapsulate([ensure_even(arr.tobytes())])
# ds.compress(pydicom.uid.JPEG2000)

# Image shape
ds['PixelData'].is_undefined_length = False
array_shape = arr.shape
ds.Rows = array_shape[0]
ds.Columns = array_shape[1]

# validate and save
pydicom.dataset.validate_file_meta(ds.file_meta, enforce_standard=True)
new_filename = filename.replace('.jpg', name + '.dcm')
ds.save_as(new_filename, write_like_original=False)

Answer 5

Corvin 的 2020 年更新几乎对我有用。 元数据仍未写入文件，因此在读取它时引发了以下异常：

pydicom.errors.InvalidDicomError: 文件缺少 DICOM 文件元信息标头或标头中缺少“DICM”前缀。

为了解决这个问题并将元数据写入 dicom 文件，我需要将 enforce_standard=True 添加到 the save_as() 调用：

ds.save_as(filename=out_filename, enforce_standard=True) 

Answer 6

上面的示例有效，但导致许多工具抱怨 DICOM，甚至使用 itk/SimpleITK 作为堆栈时根本无法读取它们。我发现从 numpy 制作 DICOM 的最佳方法是使用 SimpleITK 工具并逐个切片生成 DICOM。一个基本示例 (https://github.com/zivy/SimpleITK/blob/8e94451e4c0e90bcc6a1ffdd7bc3d56c81f58d80/Examples/DicomSeriesReadModifyWrite/DicomSeriesReadModifySeriesWrite.py) 展示了如何在堆栈中加载、执行转换然后重新保存文件，但这可以通过使用

import SimpleITK as sitk
filtered_image = sitk.GetImageFromArray(my_numpy_array)

最终输出图像中的标签数量非常多，因此手动创建所有标签非常繁琐。此外，SimpleITK 支持 8、16、32 位图像以及 RGB，因此比在 pydicom 中制作它们要容易得多。

(0008, 0008) Image Type                          CS: ['DERIVED', 'SECONDARY']
(0008, 0016) SOP Class UID                       UI: Secondary Capture Image Storage
(0008, 0018) SOP Instance UID                    UI: 1.2.826.0.1.3680043.2.1125.1.35596048796922805578234000521866725
(0008, 0020) Study Date                          DA: '20170803'
(0008, 0021) Series Date                         DA: '20170803'
(0008, 0023) Content Date                        DA: 0
(0008, 0030) Study Time                          TM: '080429.171808'
(0008, 0031) Series Time                         TM: '080429'
(0008, 0033) Content Time                        TM: 0
(0008, 0050) Accession Number                    SH: ''
(0008, 0060) Modality                            CS: 'OT'
(0008, 0064) Conversion Type                     CS: 'WSD'
(0008, 0090) Referring Physician's Name          PN: ''
(0010, 0010) Patient's Name                      PN: ''
(0010, 0020) Patient ID                          LO: ''
(0010, 0030) Patient's Birth Date                DA: ''
(0010, 0040) Patient's Sex                       CS: ''
(0018, 2010) Nominal Scanned Pixel Spacing       DS: ['1', '3']
(0020, 000d) Study Instance UID                  UI: 1.2.826.0.1.3680043.2.1125.1.33389357207068897066210100430826006
(0020, 000e) Series Instance UID                 UI: 1.2.826.0.1.3680043.2.1125.1.51488923827429438625199681257282809
(0020, 0010) Study ID                            SH: ''
(0020, 0011) Series Number                       IS: ''
(0020, 0013) Instance Number                     IS: ''
(0020, 0020) Patient Orientation                 CS: ''
(0020, 0052) Frame of Reference UID              UI: 1.2.826.0.1.3680043.2.1125.1.35696880630664441938326682384062489
(0028, 0002) Samples per Pixel                   US: 1
(0028, 0004) Photometric Interpretation          CS: 'MONOCHROME2'
(0028, 0010) Rows                                US: 40
(0028, 0011) Columns                             US: 50
(0028, 0100) Bits Allocated                      US: 32
(0028, 0101) Bits Stored                         US: 32
(0028, 0102) High Bit                            US: 31
(0028, 0103) Pixel Representation                US: 1
(0028, 1052) Rescale Intercept                   DS: "0"
(0028, 1053) Rescale Slope                       DS: "1"
(0028, 1054) Rescale Type                        LO: 'US'
(7fe0, 0010) Pixel Data                          OW: Array of 8000 bytes

Answer 7

这是我需要编写的代码的功能版本。它将从给定的 2D 像素阵列写入 16 位灰度 DICOM 图像。根据 DICOM 标准，每个图像和系列的 UID 应该是唯一的，这段代码并不担心，因为我不知道 UID 实际做什么。如果有其他人，我会很乐意添加它。

import dicom, dicom.UID
from dicom.dataset import Dataset, FileDataset
import numpy as np
import datetime, time

def write_dicom(pixel_array,filename):
    """
    INPUTS:
    pixel_array: 2D numpy ndarray.  If pixel_array is larger than 2D, errors.
    filename: string name for the output file.
    """

    ## This code block was taken from the output of a MATLAB secondary
    ## capture.  I do not know what the long dotted UIDs mean, but
    ## this code works.
    file_meta = Dataset()
    file_meta.MediaStorageSOPClassUID = 'Secondary Capture Image Storage'
    file_meta.MediaStorageSOPInstanceUID = '1.3.6.1.4.1.9590.100.1.1.111165684411017669021768385720736873780'
    file_meta.ImplementationClassUID = '1.3.6.1.4.1.9590.100.1.0.100.4.0'
    ds = FileDataset(filename, {},file_meta = file_meta,preamble="\0"*128)
    ds.Modality = 'WSD'
    ds.ContentDate = str(datetime.date.today()).replace('-','')
    ds.ContentTime = str(time.time()) #milliseconds since the epoch
    ds.StudyInstanceUID =  '1.3.6.1.4.1.9590.100.1.1.124313977412360175234271287472804872093'
    ds.SeriesInstanceUID = '1.3.6.1.4.1.9590.100.1.1.369231118011061003403421859172643143649'
    ds.SOPInstanceUID =    '1.3.6.1.4.1.9590.100.1.1.111165684411017669021768385720736873780'
    ds.SOPClassUID = 'Secondary Capture Image Storage'
    ds.SecondaryCaptureDeviceManufctur = 'Python 2.7.3'

    ## These are the necessary imaging components of the FileDataset object.
    ds.SamplesPerPixel = 1
    ds.PhotometricInterpretation = "MONOCHROME2"
    ds.PixelRepresentation = 0
    ds.HighBit = 15
    ds.BitsStored = 16
    ds.BitsAllocated = 16
    ds.SmallestImagePixelValue = '\\x00\\x00'
    ds.LargestImagePixelValue = '\\xff\\xff'
    ds.Columns = pixel_array.shape[0]
    ds.Rows = pixel_array.shape[1]
    if pixel_array.dtype != np.uint16:
        pixel_array = pixel_array.astype(np.uint16)
    ds.PixelData = pixel_array.tostring()

    ds.save_as(filename)
    return



if __name__ == "__main__":
#    pixel_array = np.arange(256*256).reshape(256,256)
#    pixel_array = np.tile(np.arange(256).reshape(16,16),(16,16))
    x = np.arange(16).reshape(16,1)
    pixel_array = (x + x.T) * 32
    pixel_array = np.tile(pixel_array,(16,16))
    write_dicom(pixel_array,'pretty.dcm')