Merge pull request #560 from LucSteinbuch/new_crosstalk_detection

New: functions underlaying paper "A novel tool to access crosstalk in luminescence detection" Anna-Maartje de Boer, Luc Steinbuch

Thanks @LucSteinbuch for this valuable contribution and all the effort you have undertaken to implement it. @mcol For aligning the code with the package structure and the review. I will make a few additional changes after the merge, but they are all of technical nature only.

NAMESPACE

@@ -94,6 +94,7 @@ export(analyse_baSAR)
 export(analyse_pIRIRSequence)
 export(analyse_portableOSL)
 export(apply_CosmicRayRemoval)
+export(apply_Crosstalk)
 export(apply_EfficiencyCorrection)
 export(bin)
 export(bin_RLum.Data)
@@ -114,6 +115,7 @@ export(calc_IEU)
 export(calc_Lamothe2003)
 export(calc_MaxDose)
 export(calc_MinDose)
+export(calc_MoransI)
 export(calc_OSLLxTxDecomposed)
 export(calc_OSLLxTxRatio)
 export(calc_SourceDoseRate)
@@ -184,6 +186,7 @@ export(plot_FilterCombinations)
 export(plot_GrowthCurve)
 export(plot_Histogram)
 export(plot_KDE)
+export(plot_MoranScatterplot)
 export(plot_NRt)
 export(plot_OSLAgeSummary)
 export(plot_RLum)
@@ -195,6 +198,7 @@ export(plot_RLum.Results)
 export(plot_ROI)
 export(plot_RadialPlot)
 export(plot_Risoe.BINfileData)
+export(plot_SingleGrainDisc)
 export(plot_ViolinPlot)
 export(read_BIN2R)
 export(read_Daybreak2R)
@@ -302,6 +306,7 @@ importFrom(stats,coef)
 importFrom(stats,complete.cases)
 importFrom(stats,confint)
 importFrom(stats,density)
+importFrom(stats,dist)
 importFrom(stats,dnorm)
 importFrom(stats,fitted)
 importFrom(stats,formula)

NEWS.Rmd

@@ -10,6 +10,10 @@ header-includes:
 
 ## New functions
 
+* `apply_Crosstalk()`, `calc_MoransI()`, `plot_SingleGrainDisc()` and
+`plot_MoranScatterplot()` were contributed by Anna-Maartje de Boer and Luc
+Steinbuch (#560).
+
 * `calc_EED_Model()` models incomplete and heterogeneous bleaching of
 mobile grains after Guibert et al. (2017). Along with the function, the
 new `ExampleData.MortarData` data set was added.

NEWS.md

@@ -5,10 +5,14 @@
 
 <!-- NEWS.md was auto-generated by NEWS.Rmd. Please DO NOT edit by hand!-->
 
-# Changes in version 0.9.26.9000-107 (2025-01-31)
+# Changes in version 0.9.26.9000-107 (2025-02-07)
 
 ## New functions
 
+- `apply_Crosstalk()`, `calc_MoransI()`, `plot_SingleGrainDisc()` and
+  `plot_MoranScatterplot()` were contributed by Anna-Maartje de Boer and
+  Luc Steinbuch (#560).
+
 - `calc_EED_Model()` models incomplete and heterogeneous bleaching of
   mobile grains after Guibert et al. (2017). Along with the function,
   the new `ExampleData.MortarData` data set was added.

R/Luminescence-package.R

@@ -110,7 +110,7 @@
 #'
 #' @importFrom graphics plot plot.default frame abline mtext text lines par layout lines arrows axTicks axis barplot box boxplot contour curve grconvertX grconvertY hist legend persp points polygon rug segments title grid close.screen screen split.screen
 #' @importFrom grDevices adjustcolor axisTicks colorRampPalette gray.colors rgb topo.colors xy.coords dev.off
-#' @importFrom stats formula approx as.formula complete.cases density dnorm glm integrate lm median na.exclude na.omit nls nls.control pchisq pnorm quantile rnorm runif sd smooth smooth.spline spline t.test uniroot var weighted.mean setNames coef confint predict update residuals fitted qf
+#' @importFrom stats formula approx as.formula complete.cases density dist dnorm glm integrate lm median na.exclude na.omit nls nls.control pchisq pnorm quantile rnorm runif sd smooth smooth.spline spline t.test uniroot var weighted.mean setNames coef confint predict update residuals fitted qf
 #' @importFrom parallel parLapply makeCluster stopCluster
 #' @importFrom httr GET accept_json status_code content
 #'

R/apply_Crosstalk.R

@@ -0,0 +1,105 @@
+#' @title Apply crosstalk
+#'
+#' @description Add crosstalk, evenly spread in rook (top-right-bottom-left)
+#' directions, to all grain whole locations on one measurement discs
+#' (=on position on a measurement wheel in a reader). An
+#' added crosstalk value of as example 0.2 means that 0.2 of the value of
+#' the central grain is added to each grain in the rook directions. This is an additive
+#' action: the central grain itself is not affected by this operation (but will on its
+#' turn increase because of crosstalk from its neighbours).
+#' This function is used for simulations: can this added crosstalk be detected?
+#'
+#' @details If an element in `object` is `NA`, it is internally set to 0, so it will not
+#' be added.
+#'
+#' @param object [RLum.Results-class] or [numeric] (**required**): containing
+#' a numerical vector of length 100, representing one or more measurement
+#' discs ("positions") in a reader.
+#' Each element in the vector represents one grain hole location on a disc.
+#'
+#' @param n_crosstalk [numeric] (*with default*): A single number quantifying the added
+#' crosstalk. Defaults for testing purposes to 0.2. Can be any number, even negative,
+#' but for realistic simulations we suggest something between 0 and 0.25.
+#'
+#' @keywords crosstalk, simulation
+#'
+#' @return A vector of size 100, with the value at each grain hole location including simulated crosstalk.
+#'
+#' @author Anna-Maartje de Boer, Luc Steinbuch, Wageningen University & Research, 2025
+#'
+#' @references
+#' de Boer, A-M., Steinbuch, L., Heuvelink, G.B.M., Wallinga, J., 2025.
+#' A novel tool to assess crosstalk in single-grain luminescence detection.
+#' Submitted.
+#'
+#' @examples
+#' ## Create artificial disc observation
+#' observations <- set_RLum(class = "RLum.Results",
+#'                          data = list(vn_values = rep(x = c(1,2), each = 50))
+#'                        )
+#' hist(get_RLum(object = observations))
+#'
+#' ## Add crosstalk (with default set to 0.2), and visualize the difference
+#' ## in the resulting histogram.
+#' observations_with_simulated_crosstalk <- apply_Crosstalk(observations)
+#' hist(observations_with_simulated_crosstalk)
+#'
+#' @md
+#' @export
+apply_Crosstalk <- function(object,
+                            n_crosstalk = 0.2
+) {
+  .set_function_name("apply_Crosstalk")
+  on.exit(.unset_function_name(), add = TRUE)
+
+  ## Validate input arguments -----------------------
+
+  .validate_class(object, c("RLum.Results", "numeric", "integer"))
+  ## To add: validation on `object`
+  #  - should contain a numerical vector of length 100
+
+  .validate_class(n_crosstalk, c("numeric"))
+  ## To add: validate on `n_crosstalk`
+  #  - should be a single numerical value
+
+  ## Set variables  -----------------------
+  if(is.numeric(object))
+  {
+    vn_values <- object
+  } else
+  {
+    vn_values <- get_RLum(object)
+  }
+
+
+  vb_na_s <- is.na(vn_values)
+
+  vn_values[vb_na_s] <- 0
+
+  ## Prepare multiplication matrix -----------------------
+  ## Note: the physical appearance of a measurement disc,
+  ## 10x10 grains, is hard coded here
+  df_disc_locations <- data.frame(location = 1:100,
+                                  x = rep(1:10, times = 10),
+                                  y = rep(1:10, each = 10))
+
+  # Calculate matrix with euclidean distances
+  mn_dist <- as.matrix(dist(df_disc_locations[,c("x", "y")]))
+
+  ## All distances equal to one are subject to crosstalk
+  mn_crosstalk <- ifelse(mn_dist == 1, yes = n_crosstalk, no = 0)
+
+  ## The diagonal -- all distances equal to zero -- are one
+  diag(mn_crosstalk) <- 1
+
+  ## Matrix calculations  -----------------------
+
+  vn_sig_with_crosstalk <- as.numeric(mn_crosstalk %*% vn_values)
+
+  ## Assign original NA's to output
+  vn_sig_with_crosstalk[vb_na_s] <- NA
+
+
+  return(vn_sig_with_crosstalk)
+
+}