From 6056cda75228fff88b93a0a6d823074cedaf6fff Mon Sep 17 00:00:00 2001
From: korbinian90 <korbinian90@gmail.com>
Date: Tue, 3 Dec 2024 13:34:01 +0100
Subject: [PATCH] add option to use external QSM mask

---
 ext/ClearswiApp/argparse.jl |  4 ++++
 ext/ClearswiApp/caller.jl   |  8 +++++++-
 src/phase_processing.jl     | 13 +++++++------
 src/utility.jl              |  5 +++--
 test/utility_test.jl        |  1 +
 5 files changed, 22 insertions(+), 9 deletions(-)

diff --git a/ext/ClearswiApp/argparse.jl b/ext/ClearswiApp/argparse.jl
index 2c9abbf..f0b20ef 100644
--- a/ext/ClearswiApp/argparse.jl
+++ b/ext/ClearswiApp/argparse.jl
@@ -27,6 +27,10 @@ function getargs(args::AbstractVector, version)
             help = """When activated uses RTS QSM for phase weighting.
             """
             action = :store_true
+        "--qsm-mask"
+            help = """The mask used for QSM. Use a custom mask, if the qsm_mask.nii 
+                is not good for your data."""
+            default = nothing
         "--mag-combine"
             help = """SNR | average | echo <n> | SE <te>.
                 Magnitude combination algorithm. echo <n> selects a specific
diff --git a/ext/ClearswiApp/caller.jl b/ext/ClearswiApp/caller.jl
index 3243e5c..1c4b366 100644
--- a/ext/ClearswiApp/caller.jl
+++ b/ext/ClearswiApp/caller.jl
@@ -85,8 +85,14 @@ function CLEARSWI.clearswi_main(args; version="1.3.3")
     phase_scaling_strength = try parse(Int, settings["phase-scaling-strength"]) catch; parse(Float32, settings["phase-scaling-strength"]) end
     writesteps = settings["writesteps"]
     qsm = settings["qsm"]
+    
+    if !isnothing(settings["qsm-mask"])
+        qsm_mask = readmag(settings["qsm-mask"]) .!= 0
+    else
+        qsm_mask = nothing
+    end
 
-    options = Options(;mag_combine, mag_sens, mag_softplus, phase_unwrap, phase_hp_sigma, phase_scaling_type, phase_scaling_strength, writesteps, qsm)
+    options = Options(;mag_combine, mag_sens, mag_softplus, phase_unwrap, phase_hp_sigma, phase_scaling_type, phase_scaling_strength, writesteps, qsm, qsm_mask=qsm_mask)
 
     swi = calculateSWI(data, options)
     mip = createIntensityProjection(swi, minimum, parse(Int, settings["mip-slices"]))
diff --git a/src/phase_processing.jl b/src/phase_processing.jl
index 39eb0d3..280d073 100644
--- a/src/phase_processing.jl
+++ b/src/phase_processing.jl
@@ -1,8 +1,7 @@
 function getswiphase(data, options)
     mask = robustmask(view(data.mag,:,:,:,1))
     savenii(mask, "maskforphase", options.writesteps, data.header)
-    # TODO output not readable
-    combined = getcombinedphase(data, options, mask)
+    combined = getcombinedphase(data, options, mask, options.qsm_mask)
     swiphase = createphasemask!(combined, mask, options.phase_scaling_type, options.phase_scaling_strength)
     savenii(swiphase, "swiphase", options.writesteps, data.header)
     return swiphase
@@ -42,7 +41,7 @@ function createphasemask!(swiphase, mask, phase_scaling_type, phase_scaling_stre
     return swiphase
 end
 
-function getcombinedphase(data, options, mask)
+function getcombinedphase(data, options, mask, qsm_mask)
     phase = data.phase
     mag = data.mag
     TEs = to_dim(data.TEs, 4)
@@ -51,7 +50,7 @@ function getcombinedphase(data, options, mask)
 
     if options.qsm
         vsz = data.header.pixdim[2:4]
-        return qsm_contrast(phase, mag, TEs, σ, vsz, save)
+        return qsm_contrast(phase, mag, TEs, qsm_mask, σ, vsz, save)
 
     elseif options.phase_unwrap == :laplacian
         return laplacian_combine(phase, mag, TEs, mask, σ, save)
@@ -66,8 +65,10 @@ function getcombinedphase(data, options, mask)
     error("Unwrapping $(options.phase_unwrap) ($(typeof(options.phase_unwrap))) not defined!")
 end
 
-function qsm_contrast(phase, mag, TEs, σ, vsz, save)
-    mask = qsm_mask_filled(phase[:,:,:,1])
+function qsm_contrast(phase, mag, TEs, mask, σ, vsz, save)
+    if isnothing(mask)
+        mask = qsm_mask_filled(phase[:,:,:,1])
+    end
     save(mask, "qsm_mask")
     combined = qsm_average(phase, mag, mask, TEs, vsz, B0=3) # uses laplacian
     save(combined, "qsm_average_laplacian")
diff --git a/src/utility.jl b/src/utility.jl
index 89a33f0..3f6aab2 100644
--- a/src/utility.jl
+++ b/src/utility.jl
@@ -62,9 +62,10 @@ struct Options
     phase_scaling_strength::Real
     writesteps::Union{AbstractString, Nothing}
     qsm::Bool
+    qsm_mask::Union{AbstractArray, Nothing}
 end
-function Options(; mag_combine=:SNR, mag_sens=nothing, mag_softplus=true, phase_unwrap=:laplacian, phase_hp_sigma=[4,4,0], phase_scaling_type=:tanh, phase_scaling_strength=4, writesteps=nothing, qsm=false)
-    Options(mag_combine, mag_sens, mag_softplus, phase_unwrap, phase_hp_sigma, phase_scaling_type, phase_scaling_strength, writesteps, qsm)
+function Options(; mag_combine=:SNR, mag_sens=nothing, mag_softplus=true, phase_unwrap=:laplacian, phase_hp_sigma=[4,4,0], phase_scaling_type=:tanh, phase_scaling_strength=4, writesteps=nothing, qsm=false, qsm_mask=nothing)
+    Options(mag_combine, mag_sens, mag_softplus, phase_unwrap, phase_hp_sigma, phase_scaling_type, phase_scaling_strength, writesteps, qsm, qsm_mask)
 end
 
 """
diff --git a/test/utility_test.jl b/test/utility_test.jl
index e535965..f2984ce 100644
--- a/test/utility_test.jl
+++ b/test/utility_test.jl
@@ -25,6 +25,7 @@ phase_scaling_type: linear
 phase_scaling_strength: 4
 writesteps: tmp
 qsm: false
+qsm_mask: nothing
 CLEARSWI.jl github version-tag: $(pkgversion(CLEARSWI))
 """
 file = "$tmp_folder/settings_swi.txt"