Dale Black (Jan 25 2022 at 03:07): Dale Black (Jan 25 2022 at 03:07):Is anyone here proficient with medical image processing that could help me debug what might be going on with this DICOM loading function that I am porting from Python to Julia? I am seeing strange behavior where the 3D DICOM images are loading in Julia in almost the exact same way, but the slices seem to be off by some factor. I am not sure where to go from here
Python function:
import os, os.path
import pydicom
def dcm_reader(dcm_path):
dcm_files = []
for (dirpathpythonames, filenames) in os.walk(dcm_path,topdown=False):
for filename in filenames:
try:
if not filename == 'DIRFILE':
dcm_file = str(os.path.join(dirpath, filename))
pydicom.read_file(dcm_file, stop_before_pixels = True)
dcm_files.append(dcm_file)
except:
pass
read_RefDs = True
while read_RefDs:
for index in range(len(dcm_files)):
try:
RefDs = pydicom.read_file(dcm_files[index], stop_before_pixels = False)
read_RefDs = False
break
except:
pass
slice_thick_ori = RefDs.SliceThickness
ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(dcm_files))
dcm_array = np.zeros([ConstPixelDims[0],ConstPixelDims[1],len(dcm_files)],\
dtype=RefDs.pixel_array.dtype)
instances = []
for filenameDCM in dcm_files:
try:
ds = pydicom.read_file(filenameDCM, stop_before_pixels = True)
instances.append(int(ds.InstanceNumber))
except:
pass
instances.sort()
index = 0
for filenameDCM in dcm_files:
try:
ds = pydicom.read_file(filenameDCM)
dcm_array[:,:,instances.index(ds.InstanceNumber)] = ds.pixel_array
if ds.InstanceNumber in instances[:2]:
if ds.InstanceNumber == instances[0]:
loc_1 = ds.SliceLocation
else:
loc_2 = ds.SliceLocation
index += 1
except:
pass
try:
RefDs.SliceThickness = abs(loc_1 - loc_2)
except:
pass
dcm_array = dcm_array * RefDs.RescaleSlope + RefDs.RescaleIntercept
return RefDs, dcm_array, slice_thick_ori
Julia function:
using DICOM
function dcm_reader(dcm_path)
dcm_files = []
for (dirpath, dirnames, filenames) in walkdir(dcm_path, topdown=false)
for filename in filenames
try
if (filename == "DIRFILE") == false
dcm_file = string(dirpath, "/", filename)
dcm_parse(dcm_file)
push!(dcm_files, dcm_file)
end
catch
nothing
end
end
end
read_RefDs = true
local RefDs
while read_RefDs
for index in range(1, length(dcm_files))
try
RefDs = dcm_parse(dcm_files[index])
read_RefDs = false
break
catch
nothing
end
end
end
header = RefDs.meta
slice_thick_ori = header[(0x0018, 0x0050)]
rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])
ConstPixelDims = (rows, cols, length(dcm_files))
dcm_array = zeros(ConstPixelDims...)
instances = []
for filenameDCM in dcm_files
try
ds = dcm_parse(filenameDCM)
head = ds.meta
InstanceNumber = head[(0x0020, 0x0013)]
push!(instances, InstanceNumber)
catch
nothing
end
end
sort(instances)
index = 0
for filenameDCM in dcm_files
try
ds = dcm_parse(filenameDCM)
head = ds.meta
InstanceNumber = head[(0x0020, 0x0013)]
index = findall(x -> x==InstanceNumber, instances)
pixel_array = head[(0x7fe0, 0x0010)]
dcm_array[:, :, index] = pixel_array
if InstanceNumber in instances[1:3]
if InstanceNumber == instances[1]
SliceLocation = head[(0x0020, 0x1041)]
loc_1 = SliceLocation
else
SliceLocation = head[(0x0020, 0x1041)]
loc_2 = SliceLocation
end
end
index += 1
catch
nothing
end
end
try
SliceThickness = head[(0x0018, 0x0050)]
SliceThickness = abs(loc_1 - loc_2)
catch
nothing
end
RescaleSlope = header[(0x0028, 0x1053)]
RescaleIntercept = header[(0x0028, 0x1052)]
dcm_array = dcm_array .* RescaleSlope .+ RescaleIntercept
return RefDs.meta, dcm_array, slice_thick_ori
end
This produces nearly identical 3D arrays with a slight difference in the slice number it seems like?
Python:
dcm_array.max()
# returns : 1275.0
dcm_array.min()
# returns: -2048.0
Julia:
maximum(dcm_array)
# returns: 1275.0
minimum(dcm_array)
# returns: -2048.0
These slices correspond to one another when comparing the Python array vs the Julia array (so it seems like it's off by a factor of 3):
image.png
image.png
But, within the same arrays, these two slices also correspond to one another (so now it seems like it's off by a factor of 5):
image.png
image.png
Zachary P Christensen (Jan 25 2022 at 04:06):Maybe take a look at https://github.com/JuliaNeuroscience/NIfTI.jl/blob/master/examples/dicom2nifti.jl
Dale Black (Jan 25 2022 at 04:51):I have gone through that before and used it for many of my other projects, so thanks for putting that out there btw! For this specifically, I need the DICOMs to load in a way that is identical to the way the python function loads it because there are a ton of downstream functions that expect this. I am more confused about why this simple port of the python function is producing such surprising results? Do you know if dcm_parse from DICOM.jl loads the pixels differently than pydicom does?
Dale Black (Jan 25 2022 at 04:53):I can't think of what else would be responsible for this shift, other than some underlying difference in DICOM.jl vs pydicom, because the rest of the function seems so straightforward.
Zachary P Christensen (Jan 25 2022 at 05:06):Is there any possibility that walking the directory is different?
Dale Black (Jan 25 2022 at 15:46):I can check, but I don't think that was the problem because the files that are output in a downstream function make sense. It might have something to with how I am using InstanceNumber in the Julia code to order the slices compared to the Python code. I will investigate those two possibilities today
Zachary P Christensen (Jan 25 2022 at 16:00):I'm pretty terrible at converting Python code so it's hard for me to diagnose some of this. although that does seem like the most likely place for things getting mixed up.
Dale Black (Jan 25 2022 at 20:23):Wow it was actually super simple.. In one of the lines I had sort(instances) when I needed sort!(instances)
Last updated: Oct 02 2023 at 04:34 UTC