added weighted average functions

docs/src/functions.md

@@ -151,3 +151,14 @@ simple plot by running `plot([results])`. The latter can be used with
 matrices of `inv_out_2D` structures, and the layout of the resulting plot will
 reflect the layout of the matrix.
 
+# Miscellaneous functions
+
+Once you have selected some peaks in your inversion results through the GUI,
+you might want to extract the weighted averages of these selected peaks,
+to get the underlying relaxation times or diffusion coefficients.
+The following functions do the job:
+
+```@docs
+weighted_averages(r::inv_out_1D)
+weighted_averages(r::inv_out_2D)
+```

docs/src/tutorial.md

@@ -26,8 +26,9 @@ data = import_spinsolve()
 !!! info
     Since we called the `import_spinsolve` function without an argument, 
     it'll open a file dialog for us to select the files we want to import.
-    (note that `import_spinsolve` requires two files, the `aqcu.par` file,
-    plus the file containing the experiment data.)
+    (note that `import_spinsolve` requires two files, the `aqcu.par` file
+    containing aqcuistion parameters, plus the `.dat`file containing the 
+    experiment data.)
 
 Now we have the data imported, the inversion can be performed using a single line of code!
 
@@ -39,15 +40,18 @@ results = invert(data)
     The `results` variable above is an `inv_out_1D` or `inv_out_2D` structure, 
     which contains all the relevant information produced by the `inversion` function.
     To access that information, we can look at the fields of the structure using the dot notation.
-    The field names can be shown by using the REPL help mode (typing ? at the julia> prompt), 
-    and typing the variable's name (in this case ?results). 
+    The field names contained in the structure can be shown by using the REPL help mode 
+    (typing ? at the julia> prompt), and typing the variable's name (in this case, `?results`). 
     Alternatively, running `@doc results` will also give you the same answers.
+
 The results can easily be visualised through the GLMakie extension of the package.
 
 ```julia
 using GLMakie
 plot(results)
 ```
+This will open a GUI with tools to interactively extract some information from the inversion results,
+by selecting regions and labelling them accordingly.
 
 !!! info
     The `plot` function of GLMakie is modified by this package 
@@ -56,16 +60,14 @@ plot(results)
     on how your plots look, it's best to create them from scratch 
     using all the tools available in GLMakie.
 
-Or, if the plot is all you need, one line of code is enough:
-
+The process above can also be achieved by a single line of code:
 ```julia
 using NMRInversions, GLMakie
 plot(invert(import_spinsolve()))
 ```
 
 Note that the workflow above can work for both 1D and 2D inversions!
 
-
 # Using the expfit function
 
 In a similar way, we can perform various exponential 

ext/gui_ext.jl

@@ -199,7 +199,9 @@ function draw_on_axes(ax1, ax2, ax3, res::NMRInversions.inv_out_1D; selections =
         elseif res.seq in [NMRInversions.PFG]
             Xlabel = ["<D> = ", " (m²/s)"]
         end
+
         wa = weighted_averages(res)
+
         for (i,s) in enumerate(res.selections)
             clr = only(
                 Makie.numbers_to_colors(

src/misc.jl

@@ -250,13 +250,51 @@ end
 export weighted_averages
 """
     weighted_averages(r::inv_out_1D)
-Return a vector with the weighted averages for the selections in `r.selections`.
+Return a vector with the weighted averages for the selections in the input structure.
 """
 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]]) 
+        println("The weighted average of peak $(i) is: $(round(wa[i], sigdigits=2))")
     end
     return wa
 end
+
+
+"""
+    weighted_averages(r::inv_out_2D)
+Return two vectors with the weighted averages for the selections in the input structure, one for each dimension.
+"""
+function weighted_averages(r::inv_out_2D)
+
+    wa_T1 = Vector(undef, length(r.selections))
+    wa_T2 = Vector(undef, length(r.selections))
+
+    z = r.F' .* r.filter'
+
+    x = range(0, 1, size(z, 1))
+    y = range(0, 1, size(z, 2))
+    points = [[i, j] for i in x, j in y]
+    mask = zeros(size(points))
+
+    for (i,s) in enumerate(r.selections)
+
+        mask .= [PolygonOps.inpolygon(p, s; in=1, on=1, out=0) for p in points]
+        spo = mask .* z
+
+        indir_dist = vec(sum(spo, dims=2))
+        dir_dist = vec(sum(spo, dims=1))
+
+        wa_T1[i] = indir_dist' *  r.X_indir / sum(spo)
+        wa_T2[i] = dir_dist' *  r.X_dir / sum(spo)
+
+        println("The weighted averages of selection $(collect('a':'z')[i]) are:")
+        println("<T₁> = $(round(wa_T1[i], sigdigits=2)), <T₂> = $(round(wa_T2[i], sigdigits=2)), "*
+                "<T₁>/<T₂> = $(round(wa_T1[i]/wa_T2[i], sigdigits=2))")
+        println()
+    end
+
+    return wa_T1, wa_T2
+end