'Remove Section' GUI implementation (#54)

* fixed typos and some formatting in documentation

* temporary solution to mac window resize bug

* defined outline for delete_range function

* added delete_range! to the exported functions

* delte_range! basic implementation

* added 'ignore selection' button to GUI

* reset button resets ignored regions

* bug fix

* update fitted curve when removing section of f

* delete_range! updates residuals

* delete selection button changes

* comments and docs change

* docs change

* removed temporary bug fix

* added ker_struct to inv_out_1D type

* f scales when filtering region

* updated filter parameter in 1D inversion results and residuals calculation

* 1D gui buttons upate

* changing how residuals are calculated and plotte in 1D case

* filter_selection only changes filter, changes to inv_out_1D

* fit and residuals updating moved from filter_selection! to just before plotting

* added more csv example files

* fixed PFG plot and filtered selections case

* rename filter function

---------

Co-authored-by: aris <aris.mav@hotmail.com>

Author: louob20

…ple_data/csv_IR/Chesire_sandstone_IR.csv …_data/csv_files/Chesire_sandstone_IR.csvexample_data/csv_IR/Chesire_sandstone_IR.csv renamed to example_data/csv_files/Chesire_sandstone_IR.csv example_data/csv_IR/Chesire_sandstone_IR.csv renamed to example_data/csv_files/Chesire_sandstone_IR.csv

@@ -0,0 +1,32 @@
+0.0007,0.346025
+0.0014,0.220525
+0.0021,0.147435
+0.0028,0.111799
+0.0035,0.0880444
+0.0042,0.0739969
+0.0049,0.0549486
+0.0056,0.0584712
+0.0063,0.0473388
+0.007,0.0434845
+0.0077,0.0281661
+0.0084,0.0363114
+0.0091,0.0402291
+0.0098,0.0319523
+0.0105,0.0415405
+0.0112,0.0320954
+0.0119,0.0252967
+0.0126,0.0283395
+0.0133,0.0177866
+0.014,0.0195451
+0.0147,0.0153591
+0.0154,0.0215488
+0.0161,0.0234149
+0.0168,0.0254711
+0.0175,0.0252673
+0.0182,0.0198021
+0.0189,0.0196933
+0.0196,0.0243988
+0.0203,0.0293147
+0.021,0.0200053
+0.0217,0.0167435
+0.0224,0.0238447

example_data/csv_files/graphene_CPMG.csv

@@ -1,4 +1,4 @@
-These are some IR (inversion recovery) data points corresponding to Cheshire sandstone, for testing the invert function.
+These are some IR (inversion recovery) data points corresponding to Cheshire sandstone, and some PFG (pulsed field gradient) data corresponding to water, for testing the invert function.
 
 You can use them as
 

example_data/csv_IR/info.md example_data/csv_files/info.mdexample_data/csv_IR/info.md renamed to example_data/csv_files/info.md

@@ -0,0 +1,16 @@
+0.0168598,8.93818
+0.101159,6.9632
+0.185458,5.40806
+0.269757,4.23223
+0.354056,3.30936
+0.438355,2.59608
+0.522654,2.02982
+0.606953,1.59765
+0.691252,1.25206
+0.775551,0.986033
+0.85985,0.777756
+0.944149,0.610387
+1.02845,0.4785
+1.11275,0.37976
+1.19705,0.299643
+1.28135,0.235117

example_data/csv_files/water_PFG.csv

@@ -102,20 +102,20 @@ function Makie.plot(res_mat::AbstractVecOrMat{NMRInversions.inv_out_1D}; selecti
     if res.seq in [NMRInversions.IR]
         ax1 = Axis(fig[1:5, 1:5], xlabel="time (s)", ylabel="Signal (a.u.)")
         ax2 = Axis(fig[1:5, 6:10], xlabel="T₁ (s)", xscale=log10)
-        ax3 = Axis(fig[6:10,1:5], xlabel= "index", ylabel="Residuals (a.u.)")
+        ax3 = Axis(fig[6:10,1:5], xlabel= "time (s)", ylabel="Residuals (a.u.)")
 
     elseif res.seq in [NMRInversions.CPMG]
         ax1 = Axis(fig[1:5, 1:5], xlabel="time (s)", ylabel="Signal (a.u.)")
         ax2 = Axis(fig[1:5, 6:10], xlabel="T₂ (s)", xscale=log10)
-        ax3 = Axis(fig[6:10,1:5], xlabel= "index", ylabel="Residuals (a.u.)")
+        ax3 = Axis(fig[6:10,1:5], xlabel= "time (s)", ylabel="Residuals (a.u.)")
 
     elseif res.seq in [NMRInversions.PFG]
         ax1 = Axis(fig[1:5, 1:5], xlabel="b factor (s/m² e-9)", ylabel="Signal (a.u.)")
         ax2 = Axis(fig[1:5, 6:10], xlabel="D (m²/s)", xscale=log10)
-        ax3 = Axis(fig[6:10,1:5], xlabel= "index", ylabel="Residuals (a.u.)")
+        ax3 = Axis(fig[6:10,1:5], xlabel= "b factor (s/m² e-9)", ylabel="Residuals (a.u.)")
     end
 
-    #=linkxaxes!(ax1, ax3)=#
+    linkxaxes!(ax1, ax3)
 
     for r in res_mat
         draw_on_axes(ax1, ax2, ax3, r, selections = selections)
@@ -157,9 +157,10 @@ function Makie.plot(res::NMRInversions.inv_out_1D)
     vlines!(fig.content[2], int_low, color=:red)
     vlines!(fig.content[2], int_high, color=:red)
 
-    button_label = Button(fig[8,6:10], label = "Save selection")
-    button_reset = Button(fig[9,6:10], label = "Reset selections")
-    button_save = Button(fig[10,6:10], label = "Save and exit")
+    button_label = Button(fig[8,6:10], label = "Label current selection")
+    button_filter = Button(fig[9,6:10], label = "Filter-out current selection")
+    button_reset = Button(fig[10,6:8], label = "Reset selections")
+    button_save = Button(fig[10,8:10], label = "Save and exit")
 
     on(button_label.clicks) do _
         push!(res.selections, slider.interval[])
@@ -171,8 +172,19 @@ function Makie.plot(res::NMRInversions.inv_out_1D)
         vlines!(fig.content[2], int_high, color=:red)
     end
 
+    on(button_filter.clicks) do _
+        empty!(fig.content[1])
+        empty!(fig.content[2])
+        empty!(fig.content[3])
+        filter_range!(res, slider.interval[])
+        draw_on_axes(fig.content[1], fig.content[2], fig.content[3], res, selections = true)
+        vlines!(fig.content[2], int_low, color=:red)
+        vlines!(fig.content[2], int_high, color=:red)
+    end
+
     on(button_reset.clicks) do _
         empty!(res.selections)
+        res.filter = ones(length(res.X))
         empty!(fig.content[1])
         empty!(fig.content[2])
         empty!(fig.content[3])
@@ -187,18 +199,26 @@ function Makie.plot(res::NMRInversions.inv_out_1D)
         save(savedir, f)
     end
 
+
     return fig
 end
 
 
 function draw_on_axes(ax1, ax2, ax3, res::NMRInversions.inv_out_1D; selections = false)
 
-    c = length(ax2.scene.plots)
+    # apply filter to f and update fitted curve and residuals
+    f_prime = res.filter .* res.f
+    yfit_prime = create_kernel(
+        res.seq, res.xfit, (res.seq == PFG ? res.X .* 1e9 : res.X)) * f_prime
+    r_prime = create_kernel(
+        res.seq, res.x, (res.seq == PFG ? res.X .* 1e9 : res.X)) * f_prime - res.y
 
+    c = length(ax2.scene.plots)
     scatter!(ax1, res.x, real.(res.y), colormap=:tab10, colorrange=(1, 10), color=c)
-    lines!(ax1, res.xfit, res.yfit, colormap=:tab10, colorrange=(1, 10), color=c)
-    lines!(ax2, res.X, res.f, colormap=:tab10, colorrange=(1, 10), color=c)
-    lines!(ax3, res.r, colormap=:tab10, colorrange=(1, 10), color=c)
+    lines!(ax1, res.xfit, yfit_prime, colormap=:tab10, colorrange=(1, 10), color=c)
+    lines!(ax2, res.X, f_prime, colormap=:tab10, colorrange=(1, 10), color=c)
+    lines!(ax3, res.x, r_prime, colormap=:tab10, colorrange=(1, 10), color=c)
+    scatter!(ax3, res.x, r_prime, colormap=:tab10, colorrange=(1, 10), color=c)
 
 
     if selections
@@ -230,11 +250,11 @@ function draw_on_axes(ax1, ax2, ax3, res::NMRInversions.inv_out_1D; selections =
                 ax2, 
                 res.X[s[1]:s[2]],
                 zeros(length(res.f[s[1]:s[2]])),
-                res.f[s[1]:s[2]],
+                (res.f .* res.filter)[s[1]:s[2]],
                 color = clr, alpha = 0.5
             )
 
-            height = (1 - (i-1) * 0.1) * maximum(res.f) 
+            height = (1 - (i-1) * 0.1) * maximum(res.f .* res.filter) 
             text!(
                 ax2, res.X[2], height, 
                 text = Xlabel[1]*"$(round(wa[i], sigdigits = 2))"*Xlabel[2], 
@@ -622,8 +642,11 @@ function Makie.plot(res::NMRInversions.inv_out_2D)
         end
 
     end ## BUTTON CLICKS
-
+    
+    
     gui
+    
+
 end
 
 end

ext/gui_ext.jl

@@ -59,25 +59,26 @@ function invert(seq::Type{<:pulse_sequence1D}, x::AbstractArray, y::Vector;
     if isa(alpha, Real)
 
         α = alpha
-        f, r = solve_regularization(ker_struct.K, ker_struct.g, α, solver)
+        f, _ = solve_regularization(ker_struct.K, ker_struct.g, α, solver)
 
     else
 
-        f, r, α = find_alpha(ker_struct, solver, alpha)
+        f, _, α = find_alpha(ker_struct, solver, alpha)
 
     end
 
     x_fit = exp10.(range(log10(1e-8), log10(1.1 * x[end]), 512))
     y_fit = create_kernel(seq, x_fit, X) * f
 
     isreal(y) ? SNR = NaN : SNR = calc_snr(y)
+    
+    r = create_kernel(seq, x, X) * f - y
 
     if seq == PFG
         X .= X ./ 1e9
     end
 
-
-    return inv_out_1D(seq, x, y, x_fit, y_fit, X, f, r, SNR, α, [],"")
+    return inv_out_1D(seq, x, y, x_fit, y_fit, X, f, r, SNR, α, [], ones(length(X)),"")
 
 end
 

src/inversion_functions.jl

@@ -247,3 +247,19 @@ function weighted_averages(r::inv_out_2D)
 
     return wa_T1, wa_T2
 end
+
+
+export filter_range!
+"""
+    filter_range!(res::inv_res_1D , range)
+Apply selected range to the filter, scaling it to keep the integral of `f` constant.
+"""
+function filter_range!(res::inv_out_1D, range)
+
+    integral = sum(res.f)
+    res.filter[range[1]:range[2]] .= 0
+    new_integral = sum(res.f .* res.filter)
+    scale = new_integral != 0 ? integral / new_integral : 0
+    res.filter .= res.filter .* scale
+
+end

src/misc.jl

@@ -341,7 +341,7 @@ end
 
 export inv_out_1D
 """
-    inv_out_1D(seq, x, y, xfit, yfit, X, f, r, SNR, α, selections, title)
+    inv_out_1D(seq, x, y, xfit, yfit, X, f, r, SNR, α, selections, filter, title)
 
 Output of the invert function for 1D pulse sequences.
 A structure containing the following fields:
@@ -355,7 +355,8 @@ A structure containing the following fields:
 - `r`, the residuals.
 - `SNR`, the signal-to-noise ratio.
 - `α`, the regularization parameter.
-- `selections`, a vector of tuples whose elements are indicate the first and last index of the selected peaks. 
+- `selections`, a vector of tuples whose elements indicate the first and last index of the selected peaks.
+- `filter`, a vector representing the mask used to filter and scale the data.
 - `title`, a title describing the data.
 
 """
@@ -371,6 +372,7 @@ mutable struct inv_out_1D
     SNR::Real
     alpha::Real
     selections::Vector{Tuple}
+    filter::Vector
     title::String
 end