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.

and what I am seeing is more akin to the following

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