make B0 combination flexible and add get_B0_snr function

Project.toml

@@ -1,7 +1,7 @@
 name = "MriResearchTools"
 uuid = "557dad86-b9bd-4533-8525-1a39684d020f"
 authors = ["Korbinian Eckstein"]
-version = "3.1.6"
+version = "3.2.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

README.md

@@ -46,7 +46,7 @@ savenii(unwrapped, "unwrapped", outputfolder, header(phase))
 **Function Reference:** https://korbinian90.github.io/MriResearchTools.jl/dev
 
 [ROMEO](https://github.com/korbinian90/ROMEO.jl) 3D/4D Phase Unwrapping  
-`romeo` `unwrap` `unwrap_individual` `romeovoxelquality` `mask_from_voxelquality`
+`romeo` `unwrap` `unwrap_individual` `romeovoxelquality` `mask_from_voxelquality` `calculateB0_unwrapped` `get_B0_snr`
 
 Laplacian unwrapping  
 `laplacianunwrap`

src/MriResearchTools.jl

@@ -47,7 +47,7 @@ export  readphase, readmag, niread, write_emptynii,
         brain_mask,
         robustrescale,
         #combine_echoes,
-        calculateB0_unwrapped,
+        calculateB0_unwrapped, get_B0_snr,
         romeovoxelquality,
         getHIP,
         laplacianunwrap, laplacianunwrap!,

src/romeofunctions.jl

@@ -10,12 +10,35 @@ The phase offsets have to be removed prior.
 
 See also [`mcpc3ds`](@ref)
 """
-function calculateB0_unwrapped(unwrapped_phase, mag, TEs)
+function calculateB0_unwrapped(unwrapped_phase, mag, TEs, type=:phase_snr)
     dims = 4
     TEs = to_dim(TEs, 4)
-    B0 = (1000 / 2π) * sum(unwrapped_phase .* mag .* mag .* TEs; dims) ./ sum(mag .* mag .* TEs.^2; dims) |> I -> dropdims(I; dims)
+    weight = get_B0_phase_weighting(mag, TEs, type)
+    B0 = (1000 / 2π) * sum(unwrapped_phase .* weight; dims) ./ sum(weight .* TEs; dims) |> I -> dropdims(I; dims)
     B0[.!isfinite.(B0)] .= 0
     return B0
 end
 
+function get_B0_phase_weighting(mag, TEs, type)
+    if type == :phase_snr
+        mag .* mag .* TEs
+    elseif type == :average
+        to_dim(ones(length(TEs)), 4)
+    elseif type == :TEs
+        TEs
+    elseif type == :mag
+        mag
+    elseif type == :magTEs
+        mag .* TEs
+    elseif type == :simulated_mag
+        mag = to_dim(exp.(-TEs / 20), 4)
+        mag .* mag .* TEs
+    end
+end
+
+function get_B0_snr(mag, TEs, type=:phase_snr)
+    weight = get_B0_phase_weighting(mag, to_dim(TEs, 4), type)
+    sum(mag .* weight; dims=4) ./ sum(weight; dims=4)
+end
+
 romeovoxelquality = voxelquality