When I try to convert a nifti file to a dicom series with ImageSeriesWriter, I encounter weird behavior when inspecting the results with ITK-SNAP. For a few cases the intensity range shifts from only positive double values in the nifti image to containing negative integer values in the dicom series. (datatypes are also defined as double and unsigned int for the image templates)
The only steps I take are reading the image, flipping it along y, and writing the dicom files.
Below is an image of the layer inspector, where you can see the same dataset, one time in nifti (on the right) and one time in dicom (on the left).
Additionally, here is an image from the datasets themselves. Again, nifti is on the right and dicom on the left. Here the intensity cut off in the bladder is obvious, but why does it occur when I donât do any processing on it? Do I need to set the data types for reading and writing differently?
That definitely looks like overflow. It would be explained if you used 16-bit signed int for output. But rescale intensity filter should help. Maybe DICOM does not like types other than 8-bit or 16-bit ints? Can you try rescaling to 16-bit unsigned and writing that to DICOM?
These are PET images and the ânativeâ pixel data type for PET images is float. When you want to store these back in DICOM you need to use the rescale slope (â0028|1053â) and rescale intercept (â0028|1052â) to select the number of digits after the decimal point you want to keep. SimpleITK code illustrating this is shown here. I suspect that if you follow a similar approach in ITK youâll get the expected results.
with the in- and output images being double as pixel type and four-dimensional. Is this wrong? I mean in theory, I just want to rescale the intensity range, so the double type should not be the issue I suppose?
@zivy, thank you for pointing me to that script, I will test that one. For the DICOM series, I have already created a header accounting for the Bit-related tags:
Not sure why that isnât working as you are exercising the same functionality as the SimpleITK example (itâs just one wrapping layer over ITK/GDCM). One thing to clarify, the DICOM dictionary based approach should be applied to the original data, no need to apply the RescaleFilter.
Code snippet doesnât show how the slope and intercept are set in the dictionary. When looking at the specific line reported in the error, it appears we shouldnât have reached it as the parameter there is int and there is another method that deals with float. So possibly something in the way the values were set in the dictionary?
That is what I am wondering about as well, it should be the same as in SimpleITK. And thanks, I was aware of that and just to confirm that I am not using the rescaling together with the dictionary approach.
Not sure why the rescale slope and intercept are ignored when set as float.
Before we dig any deeper, can you try and slightly modify the SimpleITK script (read the nifti image instead of the array creation) and confirm that the script does what you want. Once that is done, we can try to work backwards and see where the settings differ.
What happens if you set slope to 1000? The meaning might be inverse of what you seem to be trying.
If you do itk.WriteImage<itk::Image<unsigned short, 3>>(rescale->GetOutput(), "debug.nrrd"); or itk.WriteImage<itk::Image<double, 3>>(rescale->GetOutput(), "debug.nrrd");? It is better to visualize in Slicer, as ITK-SNAP converts intensities to short internally.
I modified the script to read the nifti 4D volume and removed the reading, because I just needed the output. In case it helps, below the script:
import SimpleITK as sitk
import sys
import time
import os
nifti_path = "path to nifti file"
output_path = "path to output folder"
def writeSlices(series_tag_values, new_img, out_dir, slice, volume):
image_slice = new_img[:, :, slice, volume]
# Tags shared by the series.
list(map(lambda tag_value: image_slice.SetMetaData(tag_value[0], tag_value[1]), series_tag_values))
# Slice specific tags.
# Instance Creation Date
image_slice.SetMetaData("0008|0012", time.strftime("%Y%m%d"))
# Instance Creation Time
image_slice.SetMetaData("0008|0013", time.strftime("%H%M%S"))
# Setting the type to CT so that the slice location is preserved and
# the thickness is carried over.
image_slice.SetMetaData("0008|0060", "PT")
# (0020, 0032) image position patient determines the 3D spacing between
# slices.
# Image Position (Patient)
image_slice.SetMetaData("0020|0032", '\\'.join(map(str, new_img.TransformIndexToPhysicalPoint((0, 0, slice, volume)))))
# Instance Number
image_slice.SetMetaData("0020,0013", str(slice))
# Write to the output directory and add the extension dcm, to force
# writing in DICOM format.
image_name = str(volume) + '_' + str(slice) + '.dcm'
print(image_name)
writer.SetFileName(os.path.join(out_dir, image_name))
writer.Execute(image_slice)
if __name__ == '__main__':
# Read nifti file
niftiReader = sitk.ImageFileReader()
niftiReader.SetFileName(nifti_path)
img_new = niftiReader.Execute()
number_of_volumes = img_new.GetSize()[3]
number_of_slices = img_new.GetSize()[2]
print(str(number_of_slices)+"/"+str(number_of_volumes))
writer = sitk.ImageFileWriter()
# Use the study/series/frame of reference information given in the meta-data
# dictionary and not the automatically generated information from the file IO
writer.KeepOriginalImageUIDOn()
modification_time = time.strftime("%H%M%S")
modification_date = time.strftime("%Y%m%d")
# Copy some of the tags and add the relevant tags indicating the change.
# For the series instance UID (0020|000e), each of the components is a number,
# cannot start with zero, and separated by a '.' We create a unique series ID
# using the date and time. Tags of interest:
direction = img_new.GetDirection()
series_tag_values = [
("0008|0031", modification_time), # Series Time
("0008|0021", modification_date), # Series Date
("0008|0008", "DERIVED\\SECONDARY"), # Image Type
("0020|000e", "1.2.826.0.1.3680043.2.1125."
+ modification_date + ".1" + modification_time), # Series Instance UID
("0020|0037", '\\'.join(map(str, (direction[0], direction[3], direction[6],
direction[1], direction[4],
direction[7])))), # Image Orientation
# (Patient)
("0008|103e", "Created-SimpleITK") # Series Description
]
rescale_slope = 0.001 # keep three digits after the decimal point
series_tag_values = series_tag_values + [
('0028|1053', str(rescale_slope)), # rescale slope
('0028|1052', '0'), # rescale intercept
('0028|0100', '16'), # bits allocated
('0028|0101', '16'), # bits stored
('0028|0102', '15'), # high bit
('0028|0103', '1')] # pixel representation
# Write slices to output directory
print(range(number_of_slices))
print(range(number_of_volumes))
for v in range(number_of_volumes):
for s in range(number_of_slices):
writeSlices(series_tag_values, img_new, output_path, s, v)
sys.exit(0)
The results are not really⌠useable I would say, but I can imagine that this is my fault due to potential errors I made when editing the script:
Good, so it is clear that the images themselves are cropped when it comes to their intensity.
I got it working now with the rescaler as well, which confirms that additionally. The result is shown below:
The speckles are gone, and the image looks good now. The issue I face now is that the intensity values are altered (obviously) and not the same anymore between the two datasets for the same cursor position. There is no way to retain the intensities and do the conversion without the speckle issue, I suppose?
Also, when I do the following:
std::cout << "Input image intensity range is " << rescaler->GetInputMinimum() << " to " << rescaler->GetInputMaximum() << std::endl;
I get as output: Input image intensity range is 3.40282e+38 to 0, which is kinda unexpected.
This will not work properly before rescaler->Update(), which might be called implicitly via writer->Update().
It is probably possible, but you must not use rescaling. Going via setting DICOM tags rescale slope (and rescale intercept) is your best bet. Whoever wrote this page in 2008 thought that writing float pixels with GDCM was impossible. I am not sure whether that has changed.
But looking at the DICOM standard, it should be possible. You might have to use DCMTK instead of GDCM.
As far as I know, writing float with GDCM is still not possible. At least for me, it crashes every time I try to set the pixel type to float.
I tried to build ITK with DCMTK, but I get the same errors as in here. Is there any workaround, building DCMTK by itself, or applying the âpatchâ manually?
Edit: It was possible to first build DCMTK and afterwords set the flag in ITK to use system DCMTK.
Are there any examples on writing DICOM with DCMTK as ImageIO within and Series Writer? Because I tried replacing the GDCM IO by DCMTK IO and now no files get written.
@zivy, sorry to bump you again directly, but I also found this post about the same issue, where the intercept and slope wonât get updated. Maybe this issue wasnât resolved till now?
SimpleITK/ITK can write float pixels to DICOM via GDCM using float values for slope, intercept. The DicomSeriesFromArray example illustrates this:
python DicomSeriesFromArray.py output_dir float64
Can you share the original nifti file? If yes, Iâll try and see what settings enable saving it as DICOM with float values, what modifications need to be done to your current SimpleITK script. Once we figure that out youâll be able to set your ITK code accordingly.
The rescale slope/intercept in the original script were not appropriate for your data. Please try the following script (does not retain the exact values, but nearly, youâll have to play with this if you want an exact copy):
import SimpleITK as sitk
import sys
import time
import os
def writeSlices(series_tag_values, new_img, out_dir, i):
image_slice = new_img[:, :, i]
# Tags shared by the series.
list(map(lambda tag_value: image_slice.SetMetaData(tag_value[0],
tag_value[1]),
series_tag_values))
# Slice specific tags.
# Instance Creation Date
image_slice.SetMetaData("0008|0012", time.strftime("%Y%m%d"))
# Instance Creation Time
image_slice.SetMetaData("0008|0013", time.strftime("%H%M%S"))
# Setting the type to CT so that the slice location is preserved and
# the thickness is carried over.
image_slice.SetMetaData("0008|0060", "PT")
# (0020, 0032) image position patient determines the 3D spacing between
# slices.
# Image Position (Patient)
image_slice.SetMetaData("0020|0032", '\\'.join(
map(str, new_img.TransformIndexToPhysicalPoint((0, 0, i)))))
# Instance Number
image_slice.SetMetaData("0020,0013", str(i))
# Write to the output directory and add the extension dcm, to force
# writing in DICOM format.
writer.SetFileName(os.path.join(out_dir, str(i) + '.dcm'))
writer.Execute(image_slice)
if len(sys.argv) < 3:
print("Usage: python " + __file__ + " <input_file> <output_directory>")
sys.exit(1)
new_img = sitk.ReadImage(sys.argv[1])
# Write the 3D image as a series
# IMPORTANT: There are many DICOM tags that need to be updated when you modify
# an original image. This is a delicate opration and requires
# knowledge of the DICOM standard. This example only modifies some.
# For a more complete list of tags that need to be modified see:
# http://gdcm.sourceforge.net/wiki/index.php/Writing_DICOM
# If it is critical for your work to generate valid DICOM files,
# It is recommended to use David Clunie's Dicom3tools to validate
# the files:
# http://www.dclunie.com/dicom3tools.html
writer = sitk.ImageFileWriter()
# Use the study/series/frame of reference information given in the meta-data
# dictionary and not the automatically generated information from the file IO
writer.KeepOriginalImageUIDOn()
modification_time = time.strftime("%H%M%S")
modification_date = time.strftime("%Y%m%d")
# Copy some of the tags and add the relevant tags indicating the change.
# For the series instance UID (0020|000e), each of the components is a number,
# cannot start with zero, and separated by a '.' We create a unique series ID
# using the date and time. Tags of interest:
direction = new_img.GetDirection()
series_tag_values = [
("0008|0031", modification_time), # Series Time
("0008|0021", modification_date), # Series Date
("0008|0008", "DERIVED\\SECONDARY"), # Image Type
("0020|000e", "1.2.826.0.1.3680043.2.1125."
+ modification_date + ".1" + modification_time), # Series Instance UID
("0020|0037", '\\'.join(map(str, (direction[0], direction[3], direction[6],
direction[1], direction[4],
direction[7])))), # Image Orientation
# (Patient)
("0008|103e", "Created-SimpleITK") # Series Description
]
# If we want to write floating point values, we need to use the rescale
# slope, "0028|1053", to select the number of digits we want to keep. We
# also need to specify additional pixel storage and representation
# information.
arr_view = sitk.GetArrayViewFromImage(new_img)
i_min = arr_view.min()
i_max = arr_view.max()
# the bulk intensity data range is [-2^15,+2^15], stored using 2's complement
# we rescale the actual data to fit in this discrete range.
rescale_slope = (i_min-i_max)/65536
rescale_intercept = i_min
series_tag_values = series_tag_values + [
('0028|1053', str(rescale_slope)), # rescale slope
('0028|1052', str(rescale_intercept)), # rescale intercept
('0028|0100', '16'), # bits allocated
('0028|0101', '16'), # bits stored,
('0028|0102', '15'), # high bit
('0028|0103', '1')] # pixel representation
# Write slices to output directory
list(map(lambda i: writeSlices(series_tag_values, new_img, sys.argv[2], i),
range(new_img.GetDepth())))
sys.exit(0)
Thank you @zivy, I just tested your script with the slight modification to just take the first volume for a quick run and the results look very good already. Really thank you very much for that - very helpful.
The remaining question would be how to apply the rescale and intercept to the written series within ITK for C++? Because here I still face the issue that the slope and intercept do not get applied.
Happy to hear that youâre closer to solving the problem.
So, I donât think there is much you need to do beyond a correct configuration of the meta-data information. Going over the C++ snippets above I already identified one difference. The C++ code has PixelRepresentation as 0 and the SimpleITK code has it as 1.
My advice is that you run the SimpleITK code and print out the contents of the meta-data dictionary and do the same for the C++ code. Compare the two and youâll see what you need to set to get the C++ code to work.
