Open jakubMitura14 opened 3 years ago
This package can load the Dicom files, but any further processing will require additional packages.
The shared data appears to have folders with multiple series. These are not supported by the package, but this can be remedied by defining a few supporting functions
using DICOM
function load_dicom(dir)
dcms = dcmdir_parse(dir)
loaded_dcms = Dict()
# 'dcms' could contain data for different series, so we have to filter by series
unique_series = unique([dcm.SeriesInstanceUID for dcm in dcms])
for (idx, series) in enumerate(unique_series)
dcms_in_series = filter(dcm -> dcm.SeriesInstanceUID == series, dcms)
pixeldata = extract_pixeldata(dcms_in_series)
loaded_dcms[idx] = (; pixeldata = pixeldata, dcms = dcms_in_series)
end
return loaded_dcms
end
function extract_pixeldata(dcm_array)
if length(dcm_array) == 1
return only(dcm_array).PixelData
else
return cat([dcm.PixelData for dcm in dcm_array]...; dims = 3)
end
end
The dicom files inside a folder can then be loaded by using load_dicom
defined above, e.g.
dcms = load_dicom("./sample1/DICOM/20110611/03240000")
and the pixeldata can be accessed by:
dcms[1].pixeldata # Empty
dcms[2].pixeldata # Empty
dcms[3].pixeldata # Scout
dcms[4].pixeldata # CT
dcms[5].pixeldata # PET
In the above example, the Dicom folder contained five series which is why dcms
had the keys 1, 2, 3, 4, 5
, with each number representing an individual series.
It looks like the CT and PET scans are already registered, but the CT scan dimensions are 512x512x425
while the PET scan is 200x200x68
. In case it is helpful, here are a few steps that could be done to bring the PET and CT data into the same space by using information from the DICOM header along with the Interpolations package.
(Disclaimer: I can't guarantee that the code will correctly do what it claims to do)
Assuming that the CT data will be downsampled to match the PET data, the following functions can be employed:
using Interpolations
# This ensures that first slice has lowest position, and last slice has highest position
# This assumption is required for `get_gridpoints()` to return correct slice locations
function sort_slices(scan)
slicelocations = [dcm.SliceLocation for dcm in scan.dcms]
if issorted(slicelocations)
return scan
elseif issorted(slicelocations; rev = true)
pixeldata = reverse(scan.pixeldata; dims = 3)
dcms = reverse(scan.dcms)
return (; pixeldata, dcms)
end
error("Could not sort scan")
end
function get_gridpoints(scan)
@assert ndims(scan.pixeldata) == 3
dcm = first(scan.dcms)
# [!] I might have gotten the x and y backwards below
# but this might not matter if nx = ny
(ox, oy, oz) = dcm.ImagePositionPatient
(dx, dy) = dcm.PixelSpacing
dz = dcm.SliceThickness
(nx, ny, nz) = size(scan.pixeldata)
x = ox:dx:ox+dx*(nx-1)
y = oy:dy:oy+dy*(ny-1)
z = oz:dz:oz+dz*(nz-1)
return (x, y, z)
end
function interpolate_to(input, target)
targetgrid = get_gridpoints(target)
inputgrid = get_gridpoints(input)
itp = LinearInterpolation(inputgrid, input.pixeldata, extrapolation_bc = Line())
return itp((targetgrid)...)
end
Wrapping it all together:
dcms = load_dicom("./sample1/DICOM/20110611/03240000")
ct = sort_slices(dcms[4])
pet = sort_slices(dcms[5])
interp_ct = interpolate_to(ct, pet)
# The following arrays should have the same dimensions and can be used for further processing
pet_voxels = pet.pixeldata
ct_voxels = interp_ct
Thank You For Help !!
Hello I have some anonymized Dicom fileset of FDG PET/CT (CT data + data from PET - nuclear medicine) I would like to it into tensor/ tensors for further work and processing . Is it possible using this package or some other julia package ? If not , and You will point me in a direction how to achieve it I will try to contribute to create such functionality.
The File in a link below https://drive.google.com/drive/folders/1jxJr2IdhozgSXKLutHt-eLfI9Jv2ihSZ?usp=sharing