DicomDir Pet CT

[jakubMitura14]

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

[notZaki]

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.

[notZaki]

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

[jakubMitura14]

Thank You For Help !!