From 8d4d7710e19f36b282fce5b4be9ae831844dbac0 Mon Sep 17 00:00:00 2001 From: aris Date: Wed, 6 Nov 2024 17:24:30 +0000 Subject: [PATCH] added interractivity for 1D inversion plots --- docs/src/functions.md | 4 +- ext/gui_ext.jl | 168 +++++++++++++++++++++++++++++-------- src/inversion_functions.jl | 3 +- src/misc.jl | 17 +++- src/types.jl | 6 +- 5 files changed, 156 insertions(+), 42 deletions(-) diff --git a/docs/src/functions.md b/docs/src/functions.md index b2c46c1..363a849 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -127,8 +127,8 @@ We basically take the `plot()` function offered by GLMakie and extend it to type For example, for 1D inversions, we have: ```@docs -plot(::NMRInversions.inv_out_1D) -plot!(::Union{Makie.Figure,Makie.GridPosition}, ::NMRInversions.inv_out_1D ) +plot(res::inv_out_1D) +plot(res_mat::VecOrMat{inv_out_1D};kwargs...) ``` diff --git a/ext/gui_ext.jl b/ext/gui_ext.jl index de97cff..8bf1c24 100644 --- a/ext/gui_ext.jl +++ b/ext/gui_ext.jl @@ -73,65 +73,165 @@ end ## Plots for 1D inversions """ - plot(res::NMRInversions.inv_out_1D...; kwargs...) + plot(res_mat::VecOrMat{inv_out_1D};kwargs...) -Plot the results contained in a `inv_out_1D` structure. -This function can take any number of `inv_out_1D` structures as input. +Plot the results contained in a vector of `inv_out_1D` structures. + +The keyword (optional) arguments are: + +- `selections` : Whether to highlight the selections in the plots (default is `false`). """ -function Makie.plot(res::NMRInversions.inv_out_1D...; kwargs...) +function Makie.plot(res_mat::AbstractVecOrMat{NMRInversions.inv_out_1D}; selections = false) - f = Figure(size=(800, 400)) - plot!(f, res...; kwargs...) + fig = Figure(size=(700, 600)) - return f + res = res_mat[1] + # Make axes + 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= "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= "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= "b factor (s/m² e-9)", ylabel="Residuals (a.u.)") + end + + linkxaxes!(ax1, ax3) + for r in res_mat + draw_on_axes(ax1, ax2,ax3, r, selections = selections) + end + + return fig end """ - plot!(fig, res...; title) + plot(res::inv_out_1D) -Plot the results contained in a `inv_out_1D` structure on a figure or a grid position object. -This function can take any number of `inv_out_1D` structures as input. +Run the GUI to plot the 1D inversion results and select peaks you want to label. +""" +function Makie.plot(res::NMRInversions.inv_out_1D) -The arguments are: + fig = plot([res], selections = true) -- `fig` : The figure or grid position object. -- `res` : One or more `inv_out_1D` structures containing the fit results. + slider = IntervalSlider(fig[6,6:10], range = [1:length(res.X)...], startvalues = (1, length(res.X))) -Keyword (optional) arguments: -- `title` : Title of the plot. -""" -function Makie.plot!(fig::Union{Makie.Figure,Makie.GridPosition}, res::NMRInversions.inv_out_1D...; - title="") + int_low= lift(slider.interval) do i + return res.X[i[1]] + end - # Make axes - if res[1].seq in [NMRInversions.IR] - ax1 = Axis(fig[:, 1], xlabel="time (s)", ylabel="Signal (a.u.)") - ax2 = Axis(fig[:, 2], xlabel="T₁ (s)", xscale=log10) + int_high= lift(slider.interval) do i + return res.X[i[2]] + end - elseif res[1].seq in [NMRInversions.CPMG] - ax1 = Axis(fig[:, 1], xlabel="time (s)", ylabel="Signal (a.u.)") - ax2 = Axis(fig[:, 2], xlabel="T₂ (s)", xscale=log10) + weighted_average = lift(slider.interval) do i + return res.f[i[1]:i[2]]' * res.X[i[1]:i[2]] / sum(res.f[i[1]:i[2]]) + end - elseif res[1].seq in [NMRInversions.PFG] - ax1 = Axis(fig[:, 1], xlabel="b factor (s/m² e-9)", ylabel="Signal (a.u.)") - ax2 = Axis(fig[:, 2], xlabel="D (m²/s)", xscale=log10) + labeltext = lift(slider.interval) do i + return "Selected range from $(round(res.X[i[1]], sigdigits=2)) to $(round(res.X[i[2]], sigdigits=2)) + \n Weighted average = $(round(weighted_average[], sigdigits=2))" + end + + Label(fig[7,6:10], labeltext, tellwidth=false) + 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") + + on(button_label.clicks) do _ + push!(res.selections, slider.interval[]) + empty!(fig.content[1]) + empty!(fig.content[2]) + empty!(fig.content[3]) + 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 - for r in res - draw_on_axes(ax1, ax2, r) + on(button_reset.clicks) do _ + empty!(res.selections) + empty!(fig.content[1]) + empty!(fig.content[2]) + empty!(fig.content[3]) + draw_on_axes(fig.content[1], fig.content[2], fig.content[3], res) + vlines!(fig.content[2], int_low, color=:red) + vlines!(fig.content[2], int_high, color=:red) end + on(button_save.clicks) do _ + savedir = NMRInversions.save_file(res.title, filterlist = "png") + f = plot([res], selections = true) + save(savedir, f) + end + + return fig end -function draw_on_axes(ax1, ax2, res::NMRInversions.inv_out_1D) - scatter!(ax1, res.x, res.y) - lines!(ax1, res.xfit, res.yfit) - lines!(ax2, res.X, res.f) +function draw_on_axes(ax1, ax2, ax3, res::NMRInversions.inv_out_1D; selections = false) + + c = length(ax2.scene.plots) + + scatter!(ax1, res.x, 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) + scatter!(ax3, res.x, res.r, colormap=:tab10, colorrange=(1, 10), color=c) + lines!(ax3, res.x, res.r, colormap=:tab10, colorrange=(1, 10), color=c) + if selections + + if res.seq in [NMRInversions.IR] + Xlabel = [" = ", " (s)"] + elseif res.seq in [NMRInversions.CPMG] + Xlabel = [" = ", " (s)"] + elseif res.seq in [NMRInversions.PFG] + Xlabel = [" = ", " (m²/s)"] + end + wa = weighted_averages(res) + for (i,s) in enumerate(res.selections) + clr = only( + Makie.numbers_to_colors( + [i], + Makie.to_colormap(:tab10), + identity, + Vec2(1, 10), + Makie.automatic, + Makie.automatic, + RGBAf(0, 0, 0, 0) + ) + ) + + band!( + ax2, + res.X[s[1]:s[2]], + zeros(length(res.f[s[1]:s[2]])), + res.f[s[1]:s[2]], + color = clr, alpha = 0.5 + ) + + height = (1 - (i-1) * 0.1) * maximum(res.f) + text!( + ax2, res.X[2], height, + text = Xlabel[1]*"$(round(wa[i], sigdigits = 2))"*Xlabel[2], + align = (:left, :top), color = clr + ) + + #=vlines!(ax2, res.X[s[1]], color = clr, alpha = 0.8)=# + #=vlines!(ax2, res.X[s[2]],color = clr, alpha = 0.8)=# + end + end end diff --git a/src/inversion_functions.jl b/src/inversion_functions.jl index 0393f50..9c9d390 100644 --- a/src/inversion_functions.jl +++ b/src/inversion_functions.jl @@ -82,9 +82,8 @@ function invert(seq::Type{<:pulse_sequence1D}, x::AbstractArray, y::Vector; X .= X ./ 1e9 end - weighted_average = f'X / sum(f) - return inv_out_1D(seq, x, y, x_fit, y_fit, X, f, r, SNR, α, weighted_average,"") + return inv_out_1D(seq, x, y, x_fit, y_fit, X, f, r, SNR, α, [],"") end diff --git a/src/misc.jl b/src/misc.jl index 4fcfbad..69ca9e6 100644 --- a/src/misc.jl +++ b/src/misc.jl @@ -244,4 +244,19 @@ function compress_data(t_direct::AbstractVector, G::AbstractMatrix, bins::Int=64 usv1 = svd(sqrt(W0) * K1) #paper (13) # surface(G, camera=(110, 25), xscale=:log10) -end \ No newline at end of file +end + + +export weighted_averages +""" + weighted_averages(r::inv_out_1D) +Return a vector with the weighted averages for the selections in `r.selections`. +""" +function weighted_averages(r::inv_out_1D) + + wa = Vector(undef, length(r.selections)) + for (i,s) in enumerate(r.selections) + wa[i] = r.f[s[1]:s[2]]' * r.X[s[1]:s[2]] / sum(r.f[s[1]:s[2]]) + end + return wa +end diff --git a/src/types.jl b/src/types.jl index 53765d1..117da31 100644 --- a/src/types.jl +++ b/src/types.jl @@ -207,7 +207,7 @@ end export inv_out_1D """ - inv_out_1D(seq, x, y, xfit, yfit, X, f, r, SNR, α, wa) + inv_out_1D(seq, x, y, xfit, yfit, X, f, r, SNR, α, selections, title) Output of the invert function for 1D pulse sequences. A structure containing the following fields: @@ -221,7 +221,7 @@ A structure containing the following fields: - `r`, the residuals. - `SNR`, the signal-to-noise ratio. - `α`, the regularization parameter. -- `wa`, the weighted average of the inversion results (e.g. the mean relaxation time or diffusion coefficient). +- `selections`, a vector of tuples whose elements are indicate the first and last index of the selected peaks. - `title`, a title describing the data. """ @@ -236,7 +236,7 @@ mutable struct inv_out_1D r::Vector SNR::Real alpha::Real - wa::Real + selections::Vector{Tuple} title::String end