Create pydicom file from numpy array

Question

I'm trying to create a mew dicom image from a standard-sized (512 x 512 or 256 x 256) numpy array. It seems like this should be straightforward, and I've adapted my code from http://code.google.com/p/pydicom/source/browse/source/dicom/examples/write_new.py, which appears to execute the same process, but when I save the file, I can't view it either in RadiAnt or MicroDicom.

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

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 3

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 4

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