Create pydicom file from numpy array

Question

I'm trying to create a new dicom image from a standard-sized (512 x 512 or 256 x 256) numpy array.

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 update :)

None of these answers worked for me. This is what I ended up with to save a valid monochrome 16bpp MR slice which is correctly displayed at least in Slicer, Radiant and 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")

Note the following:

• Going through FileDataset constructor as PyDicom docs suggest was failing to create a valid header for me
• validate_file_meta will create some missing elements in header for you (version)
• You need to specify endianness and explicit/implicit VR twice :/
• This method will allow you to create a valid volume as well as long as you update ImagePositionPatient and InstanceNumber for each slice accordingly
• Make sure your numpy array is cast to data format that has same number of bits as your BitsStored

Answer 2

Here is a functional version of the code I needed to write. It will write a 16-bit grayscale DICOM image from a given 2D array of pixels. According to the DICOM standard, the UIDs should be unique for each image and series, which this code doesn't worry about, because I don't know what the UIDs actually do. If anyone else does, I'll be happy to add it in.

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')

Answer 3

I was able to reduce @Corvin's great answer even more. Here is a minimalist code example allowing one to save a (dummy) 3D numpy array to a valid DICOM image that can be opened with Amide:

#!/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)

As one might observe, a lot of fields are missing if the output image.dcm file is passed to dciodvfy. The filling of these fields is left to the reader ;)

Answer 4

The above example works but causes many tools to complain about the DICOMs and they cannot even be read at all using itk/SimpleITK as a stack. The best way I have found for making DICOMs from numpy is by using the SimpleITK tools and generating the DICOMs slice-by-slice. A basic example (https://github.com/zivy/SimpleITK/blob/8e94451e4c0e90bcc6a1ffdd7bc3d56c81f58d80/Examples/DicomSeriesReadModifyWrite/DicomSeriesReadModifySeriesWrite.py) shows how to load in a stack, perform a transformation and then resave the files, but this can easily be modified by using the

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

The number of tags ultimately in output image is quite large and so manually creating all of them is tedious. Additionally SimpleITK supports 8, 16, 32-bit images as well as RGB so it is much easier than making them in 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 5

Corvin's 2020 update almost worked for me. The meta was still not written to the file, so when reading it the following exception was raised:

pydicom.errors.InvalidDicomError: File is missing DICOM File Meta Information header or the 'DICM' prefix is missing from the header.

In order to fix this and write the meta into the dicom file, I needed to add enforce_standard=True to the save_as() call:

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

Answer 6

DICOM is a really complicated format. There are many dialects, and compatibilty is rather a question of luck. You could alternatively try nibabel, maybe its dialect is more appealing to RadiAnt or MicroDicom.

In general, I'd recommend using Nifti-format whenever possible. Its standard is much more concise, and incompatibilities are rare. nibabel also supports this.

Answer 7

One working config for those who need it and one question. Question is in another thread Create a Dicom from multiple jpg images What worked for me was greyscale without compression. Every attempt at compression fails miserably I don't know why:

# 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 8

For a 3D CT scan, you can use the following code

def vol_to_dicom_for_ct(path_img_ct, patient_name, patient_id, path_dicom):

    """
    Converts a .nrrd/.mha/.nifti file into its .dcm files

    Params
    ------
    path_img_ct: str, the path of the .nrrd/.mha/.nifti file
    patient_name: str
    patient_id: str
    path_dicom: str, the final output directory
    
    Note: Verify the output with dciodvfy
        - Ref 1: https://www.dclunie.com/dicom3tools/workinprogress/index.html
        - Ref 2: https://manpages.debian.org/unstable/dicom3tools/dciodvfy.1.en.html
    """

    try:
        
        import sys
        import copy
        import random
        import shutil
        import subprocess
        import numpy as np

        if Path(path_img_ct).exists():
            
            try:
                import pydicom
                import pydicom._storage_sopclass_uids
            except:
                subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--user', 'pydicom'])
                import pydicom  

            try:
                import SimpleITK as sitk
            except:
                subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--user', 'SimpleITK']) # 2.1.1
                import SimpleITK as sitk
            
            try:
                import matplotlib.pyplot as plt
            except:
                subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--user', 'matplotlib']) # 2.1.1
                import matplotlib.pyplot as plt

            # Step 0 - Create save directory
            if Path(path_dicom).exists():
                shutil.rmtree(path_dicom)
            Path(path_dicom).mkdir(exist_ok=True, parents=True)
            
            # Step 1 - Get volume params
            img_ct      = sitk.ReadImage(str(path_img_ct))
            img_spacing = tuple(img_ct.GetSpacing())
            img_origin  = tuple(img_ct.GetOrigin()) # --> dicom.ImagePositionPatient
            img_array   = sitk.GetArrayFromImage(img_ct).astype(np.int16) # [D,H,W]

            # Step 2 - Create dicom dataset
            ds                     = pydicom.dataset.Dataset()
            ds.FrameOfReferenceUID = pydicom.uid.generate_uid() # this will stay the same for all .dcm files of a volume
            
            # Step 2.1 - Modality details
            ds.SOPClassUID         = pydicom._storage_sopclass_uids.CTImageStorage
            ds.Modality            = 'CT'
            ds.ImageType           = ['ORIGINAL', 'PRIMARY', 'AXIAL']

            # Step 2.2 - Image Details
            ds.PixelSpacing               = [float(img_spacing[0]), float(img_spacing[1])]
            ds.SliceThickness             = str(img_spacing[-1])
            ds.Rows                       = img_array.shape[1]
            ds.Columns                    = img_array.shape[2]
            
            ds.PatientPosition            = 'HFS'
            ds.ImageOrientationPatient    = [1, 0, 0, 0, 1, 0]
            ds.PositionReferenceIndicator = 'SN'
            
            ds.SamplesPerPixel            = 1
            ds.PhotometricInterpretation  = 'MONOCHROME2'
            ds.BitsAllocated              = 16
            ds.BitsStored                 = 16
            ds.HighBit                    = 15
            ds.PixelRepresentation        = 1
            
            ds.RescaleIntercept           = "0.0"
            ds.RescaleSlope               = "1.0"
            ds.RescaleType                = 'HU'
            
            # Step 3.1 - Metadata
            fileMeta = pydicom.Dataset()
            fileMeta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.CTImageStorage
            fileMeta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid() # this will change for each .dcm file of a volume
            fileMeta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian
            ds.file_meta = fileMeta

            # Step 3.2 - Include study details
            ds.StudyInstanceUID    = pydicom.uid.generate_uid()
            ds.StudyDescription    = ''
            ds.StudyDate           = '19000101'                   # needed to create DICOMDIR
            ds.StudyID             = str(random.randint(0,1000))  # needed to create DICOMDIR

            # Step 3.3 - Include series details
            ds.SeriesInstanceUID   = pydicom.uid.generate_uid()
            ds.SeriesDescription   = ''
            ds.SeriesNumber        = str(random.randint(0,1000)) # needed to create DICOMDIR

            # Step 3.4 - Include patient details
            ds.PatientName      = patient_name
            ds.PatientID        = patient_id

            # Step 3.5 - Manufacturer details
            ds.Manufacturer           = 'MICCAI2015'
            ds.ReferringPhysicianName = 'Mody'                   # needed for identification in RayStation
            ds.ManufacturerModelName  = 'test_offsite'

            # Step 4 - Make slices
            for slice_id in range(img_array.shape[0]):
                
                # Step 4.1 - Slice identifier
                random_uuid = pydicom.uid.generate_uid()
                ds.file_meta.MediaStorageSOPInstanceUID = random_uuid
                ds.SOPInstanceUID = random_uuid
                ds.InstanceNumber = str(slice_id+1)
                
                vol_origin_tmp          = list(copy.deepcopy(img_origin))
                vol_origin_tmp[-1]      += img_spacing[-1]*slice_id
                ds.ImagePositionPatient = vol_origin_tmp

                # Step 4.2 - Slice data
                img_slice               = img_array[slice_id,:,:]
                # plt.imshow(img_slice); plt.savefig(str(Path(path_dicom, '{}.png'.format(slice_id)))); plt.close()
                ds.PixelData            = img_slice.tobytes()

                save_path               = Path(path_dicom).joinpath(str(ds.file_meta.MediaStorageSOPInstanceUID) + '.dcm')
                ds.save_as(str(save_path), write_like_original=False)

            return ds.StudyInstanceUID, ds.SeriesInstanceUID

        else:
            print (' - [ERROR][vol_to_dicom_for_ct()] Error in path: path_img_ct: ', path_img_ct)
            return None, None

    except:
        traceback.print_exc()