From 96bc1b26e79736d7657e4cbe0754543fe77a2d14 Mon Sep 17 00:00:00 2001
From: Brent Sterckx <brent.sterckx@kuleuven.be>
Date: Mon, 30 Jun 2025 11:02:24 +0200
Subject: [PATCH] geom_topo no extrapolation outside convex hull

geom_topo now only shows fill when inside of the convex hull defined by the channel locations, i.e. not extrapolation outside of the convex hull.
---
 DESCRIPTION            |  3 ++-
 R/get_scalpmap.R       | 48 +++++++++++++++++++++++++++++++++++-------
 R/ggplot2-extensions.R |  4 ++--
 3 files changed, 44 insertions(+), 11 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 9ffe1eaf..fd35f26f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -36,7 +36,8 @@ Imports:
     Rcpp,
     data.table,
     isoband,
-    bslib
+    bslib,
+    lpSolveAPI
 Depends:
   R (>= 4.0.0)
 RoxygenNote: 7.3.2
diff --git a/R/get_scalpmap.R b/R/get_scalpmap.R
index 2095b310..1c4c8cba 100644
--- a/R/get_scalpmap.R
+++ b/R/get_scalpmap.R
@@ -277,6 +277,19 @@ biharmonic <- function(data,
                  nrow = grid_res,
                  ncol = grid_res))
 
+  lprec <- lpSolveAPI::make.lp(nrow = 0, ncol = nrow(data))
+  lpSolveAPI::lp.control(lprec, sense="min", epslevel = 'tight')
+  lpSolveAPI::set.objfn(lprec, numeric(nrow(data)))
+  lpSolveAPI::add.constraint(lprec, data$x, "=", NaN)
+  lpSolveAPI::add.constraint(lprec, data$y, "=", NaN)
+  lpSolveAPI::add.constraint(lprec, rep(1, nrow(data)), "=", 1)
+
+  fun <- function(p) {
+    lpSolveAPI::set.constr.value(lprec, constraints = 1, p[1]) # Adjust x constraint value
+    lpSolveAPI::set.constr.value(lprec, constraints = 2, p[2]) # Adjust y constraint value
+    return(solve(lprec) == 0)
+  }
+
   xy_coords <- unique(data[, c("x", "y")])
   xy <- xy_coords[, 1, drop = TRUE] + xy_coords[, 2, drop = TRUE] * sqrt(as.complex(-1))
   d <- matrix(rep(xy,
@@ -327,23 +340,26 @@ biharmonic <- function(data,
                         values_to = "fill",
                         names_transform = list(y = as.numeric))
 
+  # Check whether each grid point is in the convex hull defined by the elctrode locations. If not, do not plot fill.
+  data <- data[apply(cbind(c(data$x), c(data$y)), 1, fun), ]
+
   if (identical(interp_limit,
                 "head")) {
-     if (is.null(r)) {
-       circ_scale <- max_elec * 1.01
-     } else {
-       circ_scale <- r * 1.01
-     }
+    if (is.null(r)) {
+      circ_scale <- max_elec * 1.01
+    } else {
+      circ_scale <- r * 1.01
+    }
 
-   } else {
+  } else {
 
-     # add 20% or 20 mm buffer past furthest electrode, whichever is smaller
+    # add 20% or 20 mm buffer past furthest electrode, whichever is smaller
     if (r < max_elec) {
       circ_scale <- min(max_elec * 1.20, max_elec + 20)
     } else {
       circ_scale <- min(r * 1.20, r + 20)
     }
-   }
+  }
   data[sqrt(data$x ^ 2 + data$y ^ 2) <= circ_scale, ]
 
 }
@@ -374,6 +390,22 @@ fit_gam_topo <- function(data,
                                data,
                                type = "response")
 
+  lprec <- lpSolveAPI::make.lp(nrow = 0, ncol = nrow(data))
+  lpSolveAPI::lp.control(lprec, sense="min", epslevel = 'tight')
+  lpSolveAPI::set.objfn(lprec, numeric(nrow(data)))
+  lpSolveAPI::add.constraint(lprec, data$x, "=", NaN)
+  lpSolveAPI::add.constraint(lprec, data$y, "=", NaN)
+  lpSolveAPI::add.constraint(lprec, rep(1, nrow(data)), "=", 1)
+
+  fun <- function(p) {
+    lpSolveAPI::set.constr.value(lprec, constraints = 1, p[1]) # Adjust x constraint value
+    lpSolveAPI::set.constr.value(lprec, constraints = 2, p[2]) # Adjust y constraint value
+    return(solve(lprec) == 0)
+  }
+
+  # Check whether each grid point is in the convex hull defined by the elctrode locations. If not, do not plot fill.
+  data <- data[apply(cbind(c(data$x), c(data$y)), 1, fun), ]
+
   if (identical(interp_limit,
                 "head")) {
     if (is.null(r)) {
diff --git a/R/ggplot2-extensions.R b/R/ggplot2-extensions.R
index cb6fd6fd..761eb170 100644
--- a/R/ggplot2-extensions.R
+++ b/R/ggplot2-extensions.R
@@ -72,7 +72,7 @@ stat_scalpmap <- function(mapping = NULL,
                           inherit.aes = TRUE,
                           grid_res = 200,
                           interpolate = FALSE,
-                          interp_limit = c("skirt", "head"),
+                          interp_limit = c("skirt", "head", "convex_hull"),
                           method = "biharmonic",
                           r = NULL,
                           ...) {
@@ -115,7 +115,7 @@ StatScalpmap <-
                                             r = NULL) {
 
                      interp_limit <- match.arg(interp_limit,
-                                               c("skirt", "head"))
+                                               c("skirt", "head", 'convex_hull'))
                      data <- aggregate(fill ~ x + y,
                                        data = data,
                                        FUN = mean)