Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .Rbuildignore
Original file line number Diff line number Diff line change
Expand Up @@ -11,3 +11,4 @@
^.*\.Rproj$
^\.Rproj\.user$
^README\.Rmd$
^LICENSE\.md$
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Description: Fits polygon spatial masks to XY point coordinates such as
License: MIT + file LICENSE
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.3.2
RoxygenNote: 7.3.3
Depends:
R (>= 4.1.0)
Imports:
Expand Down
23 changes: 2 additions & 21 deletions LICENSE
Original file line number Diff line number Diff line change
@@ -1,21 +1,2 @@
MIT License

Copyright (c) 2026 Raredon Lab

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
YEAR: 2026
COPYRIGHT HOLDER: Raredon Laboratory
21 changes: 21 additions & 0 deletions LICENSE.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
# MIT License

Copyright (c) 2026 Raredon Laboratory

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
5 changes: 2 additions & 3 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,14 @@

export(fit_spatial_mask)
importFrom(parallel,mclapply)
importFrom(sf,"st_crs<-")
importFrom(sf,NA_crs_)
importFrom(sf,st_area)
importFrom(sf,st_as_sf)
importFrom(sf,st_bbox)
importFrom(sf,st_buffer)
importFrom(sf,st_convex_hull)
importFrom(sf,"st_crs<-")
importFrom(sf,st_covered_by)
importFrom(sf,st_difference)
importFrom(sf,st_distance)
importFrom(sf,st_geometry)
Expand All @@ -19,6 +20,4 @@ importFrom(sf,st_sfc)
importFrom(sf,st_simplify)
importFrom(sf,st_union)
importFrom(sf,st_within)
importFrom(stats,complete.cases)
importFrom(stats,filter)
importFrom(stats,quantile)
4 changes: 4 additions & 0 deletions R/TissueMask-package.R
Original file line number Diff line number Diff line change
Expand Up @@ -26,3 +26,7 @@
#'
#' @keywords internal
"_PACKAGE"

# Suppress R CMD check NOTE: "no visible binding for global variable 'x'/'y'"
# These are ggplot2 aesthetic mappings in fit_spatial_mask(plot = TRUE).
utils::globalVariables(c("x", "y"))
88 changes: 52 additions & 36 deletions R/fit_spatial_mask.R
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#' Fit a Spatial Boundary Mask to Point Coordinate Data
#'
#' @description
#' Fits a polygon (or multi-polygon) spatial mask an [sf::sfc] geometry
#' object to a set of XY point coordinates. Suitable for cell centroids,
#' Fits a polygon (or multi-polygon) spatial mask (an [sf::sfc] geometry
#' object) to a set of XY point coordinates. Suitable for cell centroids,
#' single-molecule transcript locations, or any spatial point process data.
#' The mask faithfully captures tissue shape including internal voids (vessel
#' lumens, necrotic cores, tissue tears) and disconnected fragments.
Expand All @@ -20,9 +20,9 @@
#' 1. Bin point coordinates onto a regular grid.
#' 2. Convolve the occupancy grid with a 2-D Gaussian of width `raster_sigma`
#' (in coordinate units).
#' 3. Threshold at `raster_threshold × max`. Cells above threshold are "inside".
#' 3. Threshold at `raster_threshold * max`. Cells above threshold are "inside".
#' 4. Convert each "inside" cell to an axis-aligned rectangle and dissolve via
#' GEOS union. Holes and islands emerge from the geometry naturally no
#' GEOS union. Holes and islands emerge from the geometry naturally -- no
#' ring-winding logic required.
#' 5. Apply a morphological close (buffer out then in) to smooth staircase
#' boundaries.
Expand All @@ -31,15 +31,15 @@
#' names differ, the first two columns are used and renamed. Rows containing
#' `NA` in `x` or `y` are silently dropped.
#' @param method Character scalar. Masking algorithm. One of:
#' * `"raster"` *(default)* Gaussian smoothing on a rasterised grid,
#' * `"raster"` *(default)* -- Gaussian smoothing on a rasterised grid,
#' binary threshold, then GEOS union of "on" cells. Supports holes and
#' islands. Recommended for most datasets.
#' * `"kde"` 2-D kernel density estimate via [MASS::kde2d()] contoured
#' * `"kde"` -- 2-D kernel density estimate via [MASS::kde2d()] contoured
#' with [isoband::isobands()]. Supports holes and islands. Requires
#' packages **MASS** and **isoband**.
#' * `"concave"` Concave hull via [concaveman::concaveman()]. No holes or
#' * `"concave"` -- Concave hull via [concaveman::concaveman()]. No holes or
#' islands. Requires package **concaveman**.
#' * `"convex"` Convex hull via [sf::st_convex_hull()]. No holes or
#' * `"convex"` -- Convex hull via [sf::st_convex_hull()]. No holes or
#' islands. No additional packages required.
#' @param concavity Numeric. Relative concavity for `method = "concave"`.
#' Lower values produce tighter, more concave boundaries; higher values
Expand All @@ -62,11 +62,11 @@
#' @param raster_sigma Numeric. Standard deviation of the Gaussian kernel in
#' **coordinate units** (the same units as `x` and `y`). This is the key
#' tuning parameter:
#' * Larger `raster_sigma` holes fill in, islands merge, boundary smooths.
#' * Smaller `raster_sigma` holes and fine structure are preserved.
#' * `NULL` (default) auto-set to 3% of the domain width.
#' * Larger `raster_sigma` -> holes fill in, islands merge, boundary smooths.
#' * Smaller `raster_sigma` -> holes and fine structure are preserved.
#' * `NULL` (default) -> auto-set to 3% of the domain width.
#' @param raster_threshold Numeric in `(0, 1)`. A grid cell is "inside" when
#' its smoothed value exceeds `raster_threshold × max(smoothed_grid)`.
#' its smoothed value exceeds `raster_threshold * max(smoothed_grid)`.
#' Higher values shrink the mask; lower values grow it. Default `0.15`.
#' @param raster_min_pts Integer. Minimum number of points that must fall in
#' a grid cell before it is eligible for the convolution. Cells with fewer
Expand All @@ -75,8 +75,8 @@
#' [sf::st_buffer()]. Zero means no expansion. Default `0`.
#' @param smooth_mask Logical. If `TRUE`, the mask boundary is simplified with
#' [sf::st_simplify()] after fitting. Default `FALSE`.
#' @param smooth_tolerance Numeric. DouglasPeucker tolerance passed to
#' [sf::st_simplify()]. `NULL` auto-set to 1% of the bounding-box
#' @param smooth_tolerance Numeric. Douglas-Peucker tolerance passed to
#' [sf::st_simplify()]. `NULL` -> auto-set to 1% of the bounding-box
#' diagonal. Default `NULL`.
#' @param n_cores Integer. Number of cores for parallel row-smoothing in
#' [parallel::mclapply()] (Unix/macOS only; ignored on Windows). Default
Expand Down Expand Up @@ -105,27 +105,27 @@
#' # Raster method (recommended default)
#' mask <- fit_spatial_mask(coords, verbose = FALSE)
#'
#' # Convex hull instant, no additional packages
#' # Convex hull -- instant, no additional packages
#' mask_cv <- fit_spatial_mask(coords, method = "convex", verbose = FALSE)
#'
#' # Concave hull requires 'concaveman'
#' # Concave hull -- requires 'concaveman'
#' \donttest{
#' mask_cc <- fit_spatial_mask(coords, method = "concave",
#' concavity = 1.5, verbose = FALSE)
#' }
#'
#' # Verify all points are inside the mask
#' # Verify all points are covered by the mask (boundary counts as covered)
#' library(sf)
#' pts <- sf::st_as_sf(coords, coords = c("x", "y"), crs = sf::NA_crs_)
#' all(sf::st_within(pts, sf::st_union(mask), sparse = FALSE)[, 1])
#' all(sf::st_covered_by(pts, sf::st_union(mask), sparse = FALSE)[, 1])
#'
#' @seealso [sf::st_sfc()], [sf::st_union()], [sf::st_buffer()]
#'
#' @export
#' @importFrom sf st_as_sf st_convex_hull st_union st_sfc st_polygon
#' @importFrom sf st_buffer st_distance st_simplify st_bbox st_area
#' @importFrom sf st_crs<- st_make_valid st_geometry st_sf st_within
#' @importFrom sf st_difference NA_crs_
#' @importFrom sf "st_crs<-" st_make_valid st_geometry st_sf st_within
#' @importFrom sf st_difference NA_crs_ st_covered_by
fit_spatial_mask <- function(
coords,

Expand Down Expand Up @@ -178,7 +178,7 @@ fit_spatial_mask <- function(
if (n_pts < 3L)
stop("Need at least 3 coordinate pairs to fit a mask (got ", n_pts, ").")
if (n_pts > 50000L && method == "kde" && verbose)
message("Large n = ", n_pts, " with method='kde' consider method='raster'.")
message("Large n = ", n_pts, " with method='kde' -- consider method='raster'.")

# ---- 2. Fit mask -----------------------------------------------------------
mask_poly <- switch(method,
Expand Down Expand Up @@ -235,20 +235,20 @@ fit_spatial_mask <- function(
.build_sf_with_holes(bands, crs)
},

# -- Raster: Gaussian sum threshold cell union -------------------------
# -- Raster: Gaussian sum -> threshold -> cell union -----------------------
#
# Algorithm:
# 1. Bin points onto a regular grid (raster_resolution × raster_resolution).
# 1. Bin points onto a regular grid (raster_resolution x raster_resolution).
# 2. Convolve the binary occupancy grid with a 2-D Gaussian of width
# raster_sigma (coordinate units). Equivalent to placing an isotropic
# Gaussian at every point and summing them all.
# 3. Threshold at raster_threshold × max(field). Cells above threshold
# 3. Threshold at raster_threshold * max(field). Cells above threshold
# are "inside"; cells below are "outside".
# 4. Convert each "inside" cell to an axis-aligned rectangle and dissolve
# via GEOS union. Holes and islands emerge from the geometry naturally
# no ring-winding logic required.
# -- no ring-winding logic required.
# 5. Smooth the staircase boundary with a morphological close
# (buffer out then in by ~0.6 × cell diagonal).
# (buffer out then in by ~0.6 * cell diagonal).
"raster" = {

# Grid cell edges
Expand All @@ -275,16 +275,29 @@ fit_spatial_mask <- function(
sigma_cu <- if (is.null(raster_sigma)) 0.03 * domain_w else raster_sigma
sigma_gc <- sigma_cu / hx # grid-cell units for convolution kernel

# Clamp sigma_gc so the Gaussian kernel fits within the grid.
# stats::filter() requires filter length <= time series length.
# Kernel length = 2 * ceiling(3 * sigma_gc) + 1 <= raster_resolution
# => sigma_gc <= floor((raster_resolution - 1) / 2) / 3
max_sigma_gc <- floor((raster_resolution - 1L) / 2L) / 3
if (sigma_gc > max_sigma_gc) {
if (verbose)
message(" raster_sigma clamped: sigma_gc = ", round(sigma_gc, 1),
" exceeds grid limit (", round(max_sigma_gc, 1), " gc).",
" Increase raster_resolution for finer control at large sigma.")
sigma_gc <- max_sigma_gc
}

if (verbose)
message(" raster_sigma = ", round(sigma_cu, 4),
" coord units (", round(sigma_gc, 1), " grid cells)")

# Convolve binary occupancy with Gaussian sum of point Gaussians
# Convolve binary occupancy with Gaussian -- sum of point Gaussians
indicator <- (cm >= raster_min_pts) * 1.0
sm <- .gaussian_smooth_2d(indicator, sigma = sigma_gc,
n_cores = n_cores)

# Threshold binary
# Threshold -> binary
thresh <- raster_threshold * max(sm, na.rm = TRUE)
binary <- sm >= thresh
on_cells <- which(binary, arr.ind = TRUE) # [row = yi, col = xi]
Expand All @@ -307,11 +320,11 @@ fit_spatial_mask <- function(
)))
})

# GEOS dissolve topology (holes, islands) emerges automatically
# GEOS dissolve -- topology (holes, islands) emerges automatically
raw <- sf::st_union(sf::st_sfc(rects, crs = crs))

# Morphological close: smooth staircase boundary.
# Buffer out then in by ~0.6 × cell diagonal.
# Buffer out then in by ~0.6 * cell diagonal.
cell_diag <- sqrt(hx^2 + hy^2)
raw <- sf::st_buffer(raw, dist = cell_diag * 0.6)
raw <- sf::st_buffer(raw, dist = -cell_diag * 0.6)
Expand All @@ -329,14 +342,17 @@ fit_spatial_mask <- function(
mask_poly <- sf::st_buffer(mask_poly, dist = buffer_dist)

# ---- 4. Containment guarantee ----------------------------------------------
# Use st_covered_by (not st_within): boundary points count as covered.
# st_within is strict interior containment and returns FALSE for hull vertices,
# which would trigger a zero-distance corrective buffer that does nothing.
pts_sfc <- sf::st_geometry(
sf::st_as_sf(coords, coords = c("x", "y"), crs = crs))
contained <- sf::st_within(pts_sfc, sf::st_union(mask_poly), sparse = FALSE)
n_out <- sum(!contained[, 1L])
mask_union <- sf::st_union(mask_poly)
contained <- sf::st_covered_by(pts_sfc, mask_union, sparse = FALSE)
n_out <- sum(!contained[, 1L])
if (n_out > 0L) {
warning(n_out, " point(s) outside mask — applying corrective buffer.")
dists <- sf::st_distance(pts_sfc[!contained[, 1L]],
sf::st_union(mask_poly))
warning(n_out, " point(s) outside mask - applying corrective buffer.")
dists <- sf::st_distance(pts_sfc[!contained[, 1L]], mask_union)
mask_poly <- sf::st_buffer(mask_poly,
dist = max(as.numeric(dists)) * 1.05)
}
Expand Down Expand Up @@ -387,7 +403,7 @@ fit_spatial_mask <- function(
ggplot2::coord_sf(datum = NA) +
ggplot2::labs(
title = paste0("TissueMask | method = '", method, "'"),
subtitle = paste0(n_pts, " points \u00b7 area = ",
subtitle = paste0(n_pts, " points | area = ",
round(suppressWarnings(
as.numeric(sf::st_area(mask_poly))), 2)),
x = "X", y = "Y") +
Expand Down
38 changes: 38 additions & 0 deletions man/TissueMask-package.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

23 changes: 23 additions & 0 deletions man/dot-build_sf_with_holes.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

18 changes: 18 additions & 0 deletions man/dot-gaussian_kernel_1d.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading