change usage of original kernel;

allow custom kernels

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"

src/QuantitativeSusceptibilityMappingTGV.jl

@@ -7,6 +7,6 @@ include("tgv_helper.jl")
 include("laplacian.jl")
 include("oblique_stencil.jl")
 
-export qsm_tgv, get_laplace_phase3, get_laplace_phase_del, get_laplace_phase_romeo
+export qsm_tgv, get_laplace_phase3, get_laplace_phase_del, get_laplace_phase_romeo, stencil
 
 end

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

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