Add invivo demo as a test

.gitignore

@@ -7,3 +7,4 @@ Manifest.toml
 /dev/
 /docs/build/
 /docs/site/
+test/demo_data

Project.toml

@@ -4,6 +4,8 @@ authors = ["Zaki A"]
 version = "0.2.0"
 
 [deps]
+CancerImagingArchive = "1c555f4c-8d6b-45c3-bbc1-c64989f538fe"
+DICOM = "a26e6606-dd52-5f6a-a97f-4f611373d757"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 NIfTI = "a3a9e032-41b5-5fc4-967a-a6b7a19844d3"

src/Perfusion.jl

@@ -1,7 +1,11 @@
 module Perfusion
 
 include("utils.jl")
-export @extract, interquartile_mean, percent_error
+export @extract, interquartile_mean, percent_error, apply_mask, crop
+
+using DICOM
+include("load_dicom.jl")
+export load_vfa_dicom, load_dce_dicom
 
 using LinearAlgebra: norm
 using LsqFit
@@ -37,4 +41,7 @@ export fit_rrm_nls, fit_rrm_lls, fit_crrm_nls, fit_crrm_lls
 export fit_errm_lls, fit_cerrm_lls
 export fit_rrift, fit_rrift_with_crrm, fit_rrift_with_cerrm
 
+using CancerImagingArchive
+include("invivo_demo.jl")
+
 end # module

src/invivo_demo.jl

@@ -0,0 +1,74 @@
+const gbm_study_uid = "1.3.6.1.4.1.14519.5.2.1.4591.4001.304604545029494418165835320551"
+
+function download_invivo_studies(study=gbm_study_uid; destination::AbstractString, overwrite = false)
+    make_folder(destination)
+    (vfa_folder, dce_folder) = download_vfa_and_dce(study; destination, overwrite)
+    return (vfa_folder, dce_folder)
+end
+
+function download_vfa_and_dce(study_id; destination, overwrite = false)
+    vfa_folder = joinpath(destination, study_id, "vfa")
+    dce_folder = joinpath(destination, study_id, "dce")
+
+    if isdir(vfa_folder) && isdir(dce_folder)
+        if !overwrite
+            return (vfa_folder, dce_folder)
+        end
+    end
+
+    gbm_series = tcia_series(study = study_id)
+    vfa_series = find_vfa_series(gbm_series)
+    dce_series = find_dce_series(gbm_series)
+
+    if !isdir(vfa_folder) || overwrite == true
+        download_series(vfa_series, destination = vfa_folder)
+    end
+
+    if !isdir(dce_folder) || overwrite == true
+        download_series(dce_series, destination = dce_folder)
+    end
+    return (vfa_folder, dce_folder)
+end
+
+find_vfa_series(series_dataframe) = find_in_description("MAP", series_dataframe)
+find_dce_series(series_dataframe) = find_in_description("DYN", series_dataframe)
+
+function find_in_description(word_to_find::AbstractString, series_dataframe)
+    descriptions = series_dataframe.SeriesDescription
+    found_indices = findall(occursin.(word_to_find, descriptions))
+    if length(found_indices) < 1
+        errant_id = series_dataframe.PatientID[1]
+        errant_study = series_dataframe.StudyInstanceUID[1]
+        error("No single series in $errant_id found with $word_to_find in their description.
+              The full study UID is $errant_study\n")
+    elseif length(found_indices) > 1
+        errant_id = series_dataframe.PatientID[1]
+        @warn "Multiple series in $errant_id found containing $word_to_find in their description.
+              An arbitrary series will be used among the found series"
+        found_index = found_indices[end]
+    else
+        found_index = found_indices[1]
+    end
+    found_series = series_dataframe.SeriesInstanceUID[found_index]
+    return found_series
+end
+
+function download_series(series_id::AbstractString; destination)
+    make_folder(destination; remove_existing = true)
+    zip_file = joinpath(destination, "downloaded.zip")
+    tcia_images(series = series_id, file = zip_file)
+    unzip_command = `unzip -o $zip_file -d $destination`
+    run(unzip_command)
+    rm(zip_file)
+    return destination
+end
+
+function make_folder(desired_folder::AbstractString; remove_existing = false)
+    if !isdir(desired_folder)
+        mkpath(desired_folder)
+    elseif remove_existing
+        rm(desired_folder, recursive = true)
+        mkpath(desired_folder)
+    end
+    return desired_folder
+end

src/load_dicom.jl

@@ -0,0 +1,71 @@
+function compute_concentration(; vfa_folder, dce_folder, r1=3.3/1000, num_baseline::Int=3)
+    vfa = load_vfa_dicom(folder = vfa_folder)
+    dce = load_dce_dicom(folder = dce_folder)
+
+    relaxation_maps = fit_relaxation(:despot; vfa...).estimates
+    concentration = signal_to_concentration(dce.signal, R10=1.0./relaxation_maps.T1, angle=dce.angle, TR=dce.TR, r1=r1, BAF=num_baseline)
+    return (ct = concentration, t = dce.timepoints, relaxation_maps...)
+end
+
+function load_vfa_dicom(; folder)
+    dicom_data = load_dicom_from_folder(folder)
+    number_of_images = length(dicom_data)
+    unique_flip_angles = unique(lookup.(dicom_data, "Flip Angle"))
+    number_of_flip_angles = length(unique_flip_angles)
+    number_of_slices = Int(number_of_images / number_of_flip_angles)
+
+    dummy_image = lookup(dicom_data[1], "Pixel Data")
+    image_size = size(dummy_image)
+    signal_data = zeros(image_size..., number_of_images)
+    flip_angles = zeros(number_of_images)
+    for dicom in dicom_data
+        instance = lookup(dicom, "Instance Number")
+        signal_data[:,:,instance] = lookup(dicom, "Pixel Data")
+        flip_angles[instance] = lookup(dicom, "Flip Angle")
+    end
+    signal_data = reshape(signal_data, (image_size..., number_of_slices, number_of_flip_angles))
+    flip_angles = reshape(flip_angles, (number_of_slices, number_of_flip_angles))[1,:]
+    @. flip_angles = deg2rad(flip_angles)
+    TR = lookup(dicom_data[1], "Repetition Time")
+    return (signal = signal_data, angles = flip_angles, TR = TR)
+end
+
+function load_dce_dicom(; folder, num_slices::Int = 16)
+    dicom_data = load_dicom_from_folder(folder)
+    number_of_images = length(dicom_data)
+    number_of_timepoints = Int(number_of_images / num_slices)
+
+    dummy_image = lookup(dicom_data[1], "Pixel Data")
+    image_size = size(dummy_image)
+    signal_data = zeros(image_size..., number_of_images)
+    timepoints = zeros(number_of_images)
+    for dicom in dicom_data
+        instance = lookup(dicom, "Instance Number")
+        signal_data[:,:,instance] = lookup(dicom, "Pixel Data")
+        timepoints[instance] = lookup(dicom, "Trigger Time")
+    end
+    signal_data = reshape(signal_data, (image_size..., num_slices, number_of_timepoints))
+    timepoints = reshape(timepoints, (num_slices, number_of_timepoints))[1,:] ./ 1000 ./ 60
+    TR = lookup(dicom_data[1], "Repetition Time")
+    flip_angle = deg2rad(lookup(dicom_data[1], "Flip Angle"))
+    return (signal = signal_data, timepoints = vec(timepoints), TR = TR, angle = flip_angle)
+end
+
+function load_dicom_from_folder(folder)
+    dicom_files = get_dicom_files_in_folder(folder = folder)
+    number_of_files = length(dicom_files)
+    dummy_data = dcm_parse(dicom_files[1])
+    dicom_data = Vector{typeof(dummy_data)}(undef, number_of_files)
+    for (index, file) in enumerate(dicom_files)
+        dicom_data[index] = dcm_parse(file)
+    end
+    return dicom_data
+end
+
+function get_dicom_files_in_folder(; folder)
+    files = readdir(folder)
+    dicom_files = joinpath.(folder, files[is_dicom.(files)])
+    return dicom_files
+end
+
+is_dicom(file::AbstractString) = splitext(file)[2] == ".dcm"

src/utils.jl

@@ -88,3 +88,24 @@ function unify_size(element, desired_size)
     end
 end
 
+function apply_mask(; data, mask)
+    @assert length(size(data)) >= length(size(mask))
+    mask_indices = findall(mask)
+    return data[mask_indices, :]
+end
+
+function crop(data::AbstractArray; mask=nothing)
+    if isnothing(mask)
+        mask = @. !isnan(data) & (data > 0)
+    end
+    xlim = findall(vec(sum(mask, dims = [2,3])) .> 0)
+    ylim = findall(vec(sum(mask, dims = [1,3])) .> 0)
+    return data[xlim, ylim, :]
+end
+
+function crop(data::NamedTuple; mask=nothing)
+    for key in keys(data)
+        data[key] = crop(data[key]; mask = mask)
+    end
+    return data
+end

test/runtests.jl

@@ -356,3 +356,16 @@ end
     @test std(estimates_constrained.rel_ve) < std(estimates.rel_ve)
     @test std(estimates_constrained.kep) < std(estimates.kep)
 end
+
+@testset "In-vivo demo" begin
+    destination = "./demo_data"
+    (vfa_folder, dce_folder) = Perfusion.download_invivo_studies(; destination)
+    vfa = load_vfa_dicom(folder = vfa_folder)
+    dce = load_dce_dicom(folder = dce_folder)
+
+    relaxation_maps = fit_relaxation(:despot; vfa...).estimates
+
+    bolus_arrival_frame = 3
+    r1 = 3.3/1000
+    concentration = signal_to_concentration(dce.signal; R10=1.0./relaxation_maps.T1, angle=dce.angle, TR=dce.TR, r1=r1, BAF=bolus_arrival_frame)
+end