From 45792bcf96f35773acf342f65928e2a0136b4e73 Mon Sep 17 00:00:00 2001 From: Korbinian Eckstein Date: Fri, 15 Mar 2024 14:01:07 +0100 Subject: [PATCH] change usage of original kernel; allow custom kernels --- Project.toml | 2 +- src/tgv.jl | 14 ++--- test.jl | 143 +++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 8 deletions(-) create mode 100644 test.jl diff --git a/Project.toml b/Project.toml index e4f9a8a..3ff781a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantitativeSusceptibilityMappingTGV" uuid = "bd393529-335a-4aed-902f-5de61cc7ff49" authors = ["Korbinian Eckstein korbinian90@gmail.com"] -version = "0.4.2" +version = "0.5.0" [deps] FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" diff --git a/src/tgv.jl b/src/tgv.jl index 3bb435e..a14bb95 100644 --- a/src/tgv.jl +++ b/src/tgv.jl @@ -7,7 +7,7 @@ Example: chi = qsm_tgv(randn(20,20,20), trues(20,20,20), [1,1,2]; TE=0.025, fieldstrength=3) """ -function qsm_tgv(phase, mask, res; TE, B0_dir=[0, 0, 1], fieldstrength=3, regularization=2, alpha=get_default_alpha(regularization), step_size=3, iterations=get_default_iterations(res, step_size), erosions=3, dedimensionalize=false, correct_laplacian=true, laplacian=get_laplace_phase_del, type=Float32, gpu=nothing, nblocks=32, orig_kernel=false) +function qsm_tgv(phase, mask, res; TE, B0_dir=[0, 0, 1], fieldstrength=3, regularization=2, alpha=get_default_alpha(regularization), step_size=3, iterations=get_default_iterations(res, step_size), erosions=3, dedimensionalize=false, correct_laplacian=true, laplacian=get_laplace_phase_del, type=Float32, gpu=nothing, nblocks=32, kernel=nothing) device, cu = select_device(gpu) phase, res, alpha, fieldstrength, mask = adjust_types(type, phase, res, alpha, fieldstrength, mask) @@ -28,7 +28,7 @@ function qsm_tgv(phase, mask, res; TE, B0_dir=[0, 0, 1], fieldstrength=3, regula laplace_phi0 .-= mean(laplace_phi0[mask0]) # mean correction to avoid background artefact end - alphainv, tau, sigma, resinv, laplace_kernel, dipole_kernel = set_parameters(alpha, res, B0_dir, cu; orig_kernel) + alphainv, tau, sigma, resinv, laplace_kernel, dipole_kernel = set_parameters(alpha, res, B0_dir, cu; kernel) laplace_phi0, mask, mask0 = cu(laplace_phi0), cu(mask), cu(mask0) # send to device @@ -142,17 +142,17 @@ function dipole_kernel_orig(res) return dipole_kernel end -function set_parameters(alpha, res, B0_dir, cu; orig_kernel=false) - if orig_kernel isa AbstractArray - dipole_kernel = orig_kernel - +function set_parameters(alpha, res, B0_dir, cu; kernel=nothing) + if kernel isa AbstractArray + dipole_kernel = kernel + grad_norm_sqr = 4 * sum(res .^ -2) grad_norm = sqrt(grad_norm_sqr) wave_norm = sum(abs.(dipole_kernel)) norm_matrix = [0 grad_norm 1; 0 0 grad_norm; grad_norm_sqr wave_norm 0] F = svd(norm_matrix) norm_sqr = first(F.S)^2 - elseif orig_kernel + elseif kernel == :original dipole_kernel = dipole_kernel_orig(res) grad_norm_sqr = 4 * sum(res .^ -2) norm_sqr = 2 * grad_norm_sqr^2 + 1 diff --git a/test.jl b/test.jl new file mode 100644 index 0000000..f2ef7f7 --- /dev/null +++ b/test.jl @@ -0,0 +1,143 @@ +using DICOM, QuantitativeSusceptibilityMappingTGV, CUDA, MriResearchTools, NPZ + + +fn_phase = "/home/korbi/data/MRI/data/josef_streaking/P141_3depi.dicom.non-interpolated/EnhancedMRImage_EPI3D_sag_0.67mm_iso_Pha_00001.dcm" +fn_mag = "/home/korbi/data/MRI/data/josef_streaking/P141_3depi.dicom.non-interpolated/EnhancedMRImage_EPI3D_sag_0.67mm_iso_Mag_00003.dcm" +fn_mask = "/home/korbi/data/MRI/data/josef_streaking/P141_3depi.dicom.non-interpolated/wip_QSM_20240222T074855_Mask.npy" + +phase = dcm_parse(fn_phase)[tag"PixelData"] .* 2π ./ 4096 .- π +mag = dcm_parse(fn_mag)[tag"PixelData"] +mask = npzread(fn_mask) .!= 0 + +savenii(mask, "mask", ".") + +voxel_size = [0.67, 0.67, 0.67] # in [mm] +TE = 0.035 +B0 = 2.8949 +B0_dir = [1, 0, 0] + +chi = qsm_tgv(phase, mask, voxel_size; TE=TE, fieldstrength=B0, B0_dir=B0_dir, gpu=CUDA) +savenii(chi, "chi", ".") + + + + + +## + + +phase = dcm_parse("/home/korbi/data/MRI/Siemens_QSM/data/EnhancedMRImage_3depi_0.65iso_te27_ef17_00001.dcm")[tag"PixelData"] .* 2π ./ 4096 .- π +mag = dcm_parse("/home/korbi/data/MRI/Siemens_QSM/data/EnhancedMRImage_3depi_0.65iso_te27_ef17_00002.dcm")[tag"PixelData"] +mask = dcm_parse("/home/korbi/data/MRI/Siemens_QSM/data/EnhancedMRImage_3depi_0.65iso_te27_ef17_Qsm_00002.dcm")[tag"PixelData"] .!= 0x0800 + +minp, maxp = extrema(phase) + +savenii(mask, "mask", ".") + +voxel_size = [0.325521, 0.325521, 0.65] # in [mm] +TE = 0.027 # in [s] +fieldstrength = 2.89497 # in [T] +B0_dir = [0, 0, 1] + +chi_orig = qsm_tgv(phase, mask, voxel_size; iterations=1000:1000:8000, TE=TE, fieldstrength=fieldstrength, gpu=CUDA, B0_dir=B0_dir, orig_kernel=true); + +for (i, chi) in enumerate(chi_orig) + savenii(chi, "chi_orig_$(i*1000)", ".") +end + +# chi_orig = qsm_tgv(phase, mask, voxel_size; iterations=500, TE=TE, fieldstrength=fieldstrength, gpu=CUDA, B0_dir=B0_dir); +# savenii(chi_orig, "chi_orig_500", ".") + +## First time install NPZ for reading npy files +# import Pkg +# Pkg.add("NPZ") + + +using QuantitativeSusceptibilityMappingTGV, MriResearchTools, NPZ, CUDA + +mask = npzread("/home/korbi/data/MRI/data/josef_ice_run/wip_iSWI_20220913T080705_Mask_2.npy") .!= 0; # convert to boolean +phase = npzread("/home/korbi/data/MRI/data/josef_ice_run/wip_iSWI_20220913T080705_Pha_2.npy") + +voxel_size = [1.0, 1.0, 1.0] # in [mm] +TE = 0.021 # in [s] +fieldstrength = 3 # in [T] +B0_dir = [0, 0, 1] + +chi = qsm_tgv(phase, mask, voxel_size; TE=TE, fieldstrength=fieldstrength, gpu=CUDA, B0_dir=B0_dir); +chi_orig = qsm_tgv(phase, mask, voxel_size; TE=TE, fieldstrength=fieldstrength, gpu=CUDA, B0_dir=B0_dir, orig_kernel=true); + +savenii(chi, "chi", ".") +savenii(chi_orig, "chi_orig", ".") +savenii((chi .- chi_orig) ./ chi_orig, "chi_relative", ".") # relative difference + +# Change regularization strength (1-4) +# chi = qsm_tgv(phase, mask, res; TE, fieldstrength, regularization=1); + +## + +using MriResearchTools, QuantitativeSusceptibilityMappingTGV, CUDA + +phase = readphase("/home/korbi/data/MRI/data/ashley_qsm_test/bids-3depi/sub-1/ses-064mm-epi/anat/sub-1_ses-064mm-epi_run-1_part-phase_T2starw.nii") +mask = niread("/home/korbi/data/MRI/data/ashley_qsm_test/result/mask.nii") .!= 0 + +voxel_size = [1.0, 1.0, 1.0] # in [mm] +TE = 0.021 # in [s] +fieldstrength = 3 # in [T] +B0_dir = [0, 0, 1] + +chi_orig = qsm_tgv(phase, mask, voxel_size; TE, fieldstrength, orig_kernel=true, gpu=CUDA); +chi_st = qsm_tgv(phase, mask, voxel_size; TE, fieldstrength, orig_kernel=:test, gpu=CUDA); + + +savenii(chi_orig, "chi_orig", ".") +savenii(chi_st, "chi_st", ".") + + +savenii((chi_st .- chi_orig) ./ abs.(chi_orig), "chi_st_relative", ".") # relative differenceswapdims(phase, 1, 3) +savenii(chi_st .- chi_orig, "chi_diff", ".") + +## Rotation + +# Rotate 90 degrees around y-axis +phase_x = permutedims(phase, [3, 2, 1]) +mask_x = permutedims(mask, [3, 2, 1]) +chi_x = qsm_tgv(phase_x, mask_x, voxel_size; TE, fieldstrength, orig_kernel=dipole_kernel_orig_x(voxel_size), gpu=CUDA); +savenii(chi_x, "chi_x", ".") + +## + +function dipole_kernel_orig_z(res) + dipole_kernel = zeros(Float32, 3, 3, 3) + dipole_kernel[3, 2, 2] = 1/3 / res[1]^2 + dipole_kernel[2, 3, 2] = 1/3 / res[2]^2 + dipole_kernel[2, 2, 3] = -2/3 / res[3]^2 + dipole_kernel[2, 2, 1] = -2/3 / res[3]^2 + dipole_kernel[2, 1, 2] = 1/3 / res[2]^2 + dipole_kernel[1, 2, 2] = 1/3 / res[1]^2 + dipole_kernel[2, 2, 2] = -sum(dipole_kernel) + return dipole_kernel +end + +function dipole_kernel_orig_x(res) + dipole_kernel = zeros(Float32, 3, 3, 3) + dipole_kernel[3, 2, 2] = -2/3 / res[1]^2 + dipole_kernel[2, 3, 2] = 1/3 / res[2]^2 + dipole_kernel[2, 2, 3] = 1/3 / res[3]^2 + dipole_kernel[2, 2, 1] = 1/3 / res[3]^2 + dipole_kernel[2, 1, 2] = 1/3 / res[2]^2 + dipole_kernel[1, 2, 2] = -2/3 / res[1]^2 + dipole_kernel[2, 2, 2] = -sum(dipole_kernel) + return dipole_kernel +end + +function dipole_kernel_orig_y(res) + dipole_kernel = zeros(Float32, 3, 3, 3) + dipole_kernel[3, 2, 2] = 1/3 / res[1]^2 + dipole_kernel[2, 3, 2] = -2/3 / res[2]^2 + dipole_kernel[2, 2, 3] = 1/3 / res[3]^2 + dipole_kernel[2, 2, 1] = 1/3 / res[3]^2 + dipole_kernel[2, 1, 2] = -2/3 / res[2]^2 + dipole_kernel[1, 2, 2] = 1/3 / res[1]^2 + dipole_kernel[2, 2, 2] = -sum(dipole_kernel) + return dipole_kernel +end \ No newline at end of file