write 2D array to dicom

[Elhacene Matene]

Hello,

I try to write dicom from numpy array (2D), i used this code :

import dicom, dicom.UID
from dicom.dataset import Dataset, FileDataset
import numpy as np
import struct
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.
    H1 = [1, 0, 0, 0, 1, 0, 0, 0, 1]
    H2 = [0.044023399999999997, 0.041792500000000003, 0.22416422222222224]
    file_meta = Dataset()
    file_meta.MediaStorageSOPClassUID = '1.2.840.10008.5.1.4.1.1.2' #CT image
    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.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 = '1.2.840.10008.5.1.4.1.1.2'
    ds.SecondaryCaptureDeviceManufctur = 'Python 2.7.3'

    ## These are the necessary imaging components of the FileDataset object.

    # Image Orientation from CT
    ds.ImageOrientationPatient = [str(H1[0]), str(H1[1]), str(H1[2]), str(H1[3]), str(H1[4]), str(H1[5]), str(H1[6]), str(H1[7]), str(H1[8])]
        #the upper left coordinates
    ds.ImagePositionPatient = [str(H2[0]), str(H2[1]), str(H2[2])]
    ds.Columns = pixel_array.shape[0]
    ds.Rows = pixel_array.shape[1]
    ds.Widh = pixel_array.shape[0]
    ds.Height = pixel_array.shape[1]
    ds.SliceThickness = str(0.5)
    ds.PixelSpacing = [str(0.5), str(0.5)]
    ds.SliceLocation = str(0.22416422222222224)
    ds.ImagerPixelSpacing = [str(0.5), str(0.5)]
    ds.PixelRepresentation = 0
    ds.WindowCenter = 0
    ds.WindowWidth = 1000
    ds.RescaleIntercept = -1024
    ds.RescaleSlope = 1
    ds.RescaleType = 'HU'
    ds.PhotometricInterpretation = "MONOCHROME2"
    ds.ColorType = 'grayscale'
    pixel_array = pixel_array.transpose(1,0)
#    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')




the problem, i get a dicom with a different shape of data.
Is there a probleme of data ordering in numpy ?

thank you for your help

Best regards

[Darcy Mason]

I ran your code, and there are several problems.  First is that I get a int64 type for the array coming from __main__, when I suspect you want uint16 (I see the code to convert to 16 was there and commented out).

Second, when I read the file back in and try to look at ds.pixel_array, I get exceptions because BitsAllocated and (then SamplesPerPixel) are not stored.  The SamplesPerPixel exception went away when I set BitsAllocated (correctly in this case) to 64. Having done all that, I read the array back and the shape is (256, 256) as expected.

Also I'm not sure why there is a transpose operation in the write_dicom function, although the Rows and Columns are set as if there is one.

Can you be more specific about what shapes you are trying to write and what you get back?

-Darcy