Incorrect pixel data when creating DICOM Data using PyDicom

Viewed 15

Edit1: By changing the DICOM Bits stored etc. to 64, the image looks great when loaded back into Python now. However, my saved data still looks exactly as it did before when viewed in ImageJ. I suspect that it might be a bit issue with ImageJ for now.

Edit2: I changed temparray to type 'int16' and changed the dicom attributes back, and it looks good in ImageJ now. So in theory, I think I have solved my issue, but I'll leave this open in case anyone else comes across this as an issue in the future.

temparray = temparray.astype('int16')

I'm very new to python but I am currently trying to convert some text file data (specifically some ICRP data for a reference human) into a DICOM image to use as an "atlas" (i.e. a 3D array of integers where each integer corresponds to an organ or tissue within a human). The code itself was a bit of a headache, as the original ICRP dat file is very messy (different number of columns, extra characters etc.) but I think I have finally got there with the following:

filename = "/media/sf_VMachine_Shared_Path/ValidationData/ICRP110/AF/AF.dat"

#this code reads in the data as a direct line, removing all new line and tab characters
with open(filename, "r") as data1:
    fileimport = data1.read().replace('\n','').replace('\t','')

#split the import using spaces and remove empty elements (filter as empty elements return false)
fileimportarray = fileimport.split(" ")
fileimportarray = list(filter(None, fileimportarray))

#get matrix dimensions (different for male and female)
dimX = 299
dimY = 137
dimZ = 348

#create an empty 3D array the size of the final output data for the atlas and change the type to integers
atlasarray = np.array(fileimportarray).reshape(dimZ, dimY, dimX)
atlasarray = atlasarray.astype(int)

This code could probably be more succinct if I had more knowledge, but appears to work well when checked with Matplotlib.

I then tried to slice this data into 2D images and iterated through to generate the total 3D information using the following code (some of it from this very website). At this stage, I just want to generate a 3D dataset containing integer values in each voxel. Currently, I am using .tobytes for ds.PixelData but I have tried a few other options there without success.

#following code is edited from https://stackoverflow.com/questions/14350675/create-pydicom-file-from matplotlib import pyplot as plt
from pydicom.dataset import Dataset, FileDataset
from pydicom.uid import ExplicitVRLittleEndian
import pydicom._storage_sopclass_uids

print("Setting meta data...")
meta = pydicom.Dataset()
meta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.CTImageStorage
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.CTImageStorage
ds.PatientName = "Atlas"
ds.PatientID = "123456"
ds.Modality = "CT"
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=dimZ

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"0.3\0.3"
ds.PhotometricInterpretation = "MONOCHROME2"
ds.PixelRepresentation = 1

#loop through slices and make a file for each one
z=0
for z in range(dimZ-1):
    ds.Rows = dimY
    ds.Columns = dimX
    ds.InstanceNumber=z+1
    pydicom.dataset.validate_file_meta(ds.file_meta, enforce_standard=True)
    temparray=atlasarray[z,:,:]

    #test plots
    if z == 248:
        plt.imshow(temparray, interpolation='nearest')
        plt.show()

    #ds.PixelData = np.ascontiguousarray(temparray,dtype=int)
    ds.PixelData = temparray.tobytes()
    ds.save_as("output/atlas"+str(z+1)+".dcm")

print("DCM files saved in /home/gate/PycharmProjects/pythonProject/Learning")

#test the image by reimporting
from pydicom.data import get_testdata_files

testfile = "/home/gate/PycharmProjects/pythonProject/Learning/output/atlas248.dcm"
dataset = pydicom.dcmread(testfile, force=True)
plt.imshow(dataset.pixel_array, cmap=plt.cm.bone)
plt.show()

Looking at the data in an external program, such as ImageJ, results in incorrect images, so I thought I'd try generating a random slice within Python using Matplotlib and then show the same slice when reimported (in case it was an ImageJ issue with the header or something) but I get the same issue. Slice 248 is just a random slice with decent anatomy to make it easy to see.

This image below is the Matplotlib version of slice 248 and is exactly how I expect the image to look. matplotlib image of array slice 248

and what I am seeing is more akin to the following enter image description here

Before I checked the images with Matplotlib, I thought it might have been an array issue, but I am not convinced anymore. Additionally, when running the code, I get a message about the length of the pixel data in the dataset indicates that it contains padding and 245k bytes will be removed from the end of the data.

I'd be grateful for any help or pointers on this. I'm new to python (and not particularly au fait with the intricacies of DICOM either), so feel free to trash my code if needed. Thank you

1 Answers

See edits at top of page. Left question in case other people end up here via Google.

Related