diff --git a/.Rbuildignore b/.Rbuildignore
index 8a8b2e1..e733708 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -11,3 +11,4 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 ^README\.Rmd$
+^LICENSE\.md$
diff --git a/DESCRIPTION b/DESCRIPTION
index f2266b3..ec46a3d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -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:
diff --git a/LICENSE b/LICENSE
index 41f0f84..a46a28d 100644
--- a/LICENSE
+++ b/LICENSE
@@ -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
diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000..3fb6ede
--- /dev/null
+++ b/LICENSE.md
@@ -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.
diff --git a/NAMESPACE b/NAMESPACE
index c924fcc..e860c9b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -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)
@@ -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)
diff --git a/R/TissueMask-package.R b/R/TissueMask-package.R
index 9ea3b00..45e8cb0 100644
--- a/R/TissueMask-package.R
+++ b/R/TissueMask-package.R
@@ -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"))
diff --git a/R/fit_spatial_mask.R b/R/fit_spatial_mask.R
index 3269ffb..cbb2060 100644
--- a/R/fit_spatial_mask.R
+++ b/R/fit_spatial_mask.R
@@ -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.
@@ -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.
@@ -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
@@ -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
@@ -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. Douglas–Peucker 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
@@ -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,
 
@@ -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,
@@ -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
@@ -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]
@@ -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)
@@ -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)
   }
@@ -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") +
diff --git a/man/TissueMask-package.Rd b/man/TissueMask-package.Rd
new file mode 100644
index 0000000..4cd7179
--- /dev/null
+++ b/man/TissueMask-package.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/TissueMask-package.R
+\docType{package}
+\name{TissueMask-package}
+\alias{TissueMask}
+\alias{TissueMask-package}
+\title{TissueMask: Spatial Boundary Mask Fitting for Point Coordinate Data}
+\description{
+TissueMask provides \code{\link[=fit_spatial_mask]{fit_spatial_mask()}}, which fits polygon spatial masks
+to XY point coordinate data such as cell centroids or transcript locations
+from spatial transcriptomics and spatial proteomics experiments.
+
+Four masking methods are available, supporting geometries of increasing
+topological complexity — from simple convex hulls to multi-polygon masks
+with internal voids and disconnected islands:\tabular{llll}{
+   Method \tab Holes/Islands \tab Speed \tab Recommended for \cr
+   \code{"convex"} \tab No \tab Instant \tab Sanity checks \cr
+   \code{"concave"} \tab No \tab Fast \tab Simple non-convex shapes \cr
+   \code{"kde"} \tab Yes \tab Moderate \tab Dense, smooth datasets \cr
+   \code{"raster"} \tab Yes \tab Fast \tab Most use cases (default) \cr
+}
+}
+\section{Part of TissueSuite}{
+
+TissueMask is the first component of the \strong{TissueSuite} package family
+developed at the Raredon Laboratory, Yale School of Medicine. The companion
+package \strong{TissueField} uses TissueMask output to estimate steady-state
+molecular concentration fields via a diffusion–clearance PDE.
+}
+
+\seealso{
+\code{\link[=fit_spatial_mask]{fit_spatial_mask()}}
+}
+\author{
+\strong{Maintainer}: Raredon Laboratory \email{raredonlab@yale.edu}
+
+}
+\keyword{internal}
diff --git a/man/dot-build_sf_with_holes.Rd b/man/dot-build_sf_with_holes.Rd
new file mode 100644
index 0000000..3ac2ae6
--- /dev/null
+++ b/man/dot-build_sf_with_holes.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.build_sf_with_holes}
+\alias{.build_sf_with_holes}
+\title{Assemble sf geometry with holes from isoband output}
+\usage{
+.build_sf_with_holes(bands, crs)
+}
+\arguments{
+\item{bands}{List of band objects from \code{\link[isoband:isobands]{isoband::isobands()}}, each with
+\verb{$x}, \verb{$y}, and optionally \verb{$id} components.}
+
+\item{crs}{CRS specification for the output \link[sf:sfc]{sf::sfc} object.}
+}
+\value{
+An \link[sf:sfc]{sf::sfc} geometry (single or multi-polygon with holes).
+}
+\description{
+Used only by \code{method = "kde"}. Takes a list of isoband band objects,
+extracts rings, sorts by area, and assembles exterior polygons with
+interior holes using GEOS set operations.
+}
+\keyword{internal}
diff --git a/man/dot-gaussian_kernel_1d.Rd b/man/dot-gaussian_kernel_1d.Rd
new file mode 100644
index 0000000..aa9fc58
--- /dev/null
+++ b/man/dot-gaussian_kernel_1d.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.gaussian_kernel_1d}
+\alias{.gaussian_kernel_1d}
+\title{Normalised 1-D Gaussian kernel}
+\usage{
+.gaussian_kernel_1d(sigma)
+}
+\arguments{
+\item{sigma}{Standard deviation in grid-cell units.}
+}
+\value{
+Numeric vector of kernel weights (sums to 1).
+}
+\description{
+Normalised 1-D Gaussian kernel
+}
+\keyword{internal}
diff --git a/man/dot-gaussian_smooth_2d.Rd b/man/dot-gaussian_smooth_2d.Rd
new file mode 100644
index 0000000..30baf36
--- /dev/null
+++ b/man/dot-gaussian_smooth_2d.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.gaussian_smooth_2d}
+\alias{.gaussian_smooth_2d}
+\title{Separable 2-D Gaussian convolution}
+\usage{
+.gaussian_smooth_2d(mat, sigma, n_cores = 1L)
+}
+\arguments{
+\item{mat}{Numeric matrix to smooth.}
+
+\item{sigma}{Gaussian standard deviation in grid-cell units.}
+
+\item{n_cores}{Integer. Number of cores for parallel row smoothing
+(Unix only). Default \code{1L}.}
+}
+\value{
+Smoothed numeric matrix.
+}
+\description{
+Applies a separable Gaussian blur to a matrix by convolving rows then
+columns with a 1-D Gaussian kernel. \code{sigma} is in \strong{grid-cell units}.
+For parallel row processing on Unix, set \code{n_cores > 1}.
+}
+\keyword{internal}
diff --git a/man/dot-signed_area.Rd b/man/dot-signed_area.Rd
new file mode 100644
index 0000000..85a1192
--- /dev/null
+++ b/man/dot-signed_area.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.signed_area}
+\alias{.signed_area}
+\title{Signed area of a ring matrix}
+\usage{
+.signed_area(ring_mat)
+}
+\arguments{
+\item{ring_mat}{Numeric matrix with columns for x and y coordinates.
+Must have at least 3 rows.}
+}
+\value{
+Numeric scalar (signed area).
+}
+\description{
+Computes the signed area of a closed ring via the shoelace formula.
+Positive = counter-clockwise (exterior), negative = clockwise (hole).
+}
+\keyword{internal}
diff --git a/man/dot-smooth_cols.Rd b/man/dot-smooth_cols.Rd
new file mode 100644
index 0000000..860de4b
--- /dev/null
+++ b/man/dot-smooth_cols.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.smooth_cols}
+\alias{.smooth_cols}
+\title{Column-wise 1-D convolution}
+\usage{
+.smooth_cols(mat, k)
+}
+\arguments{
+\item{mat}{Numeric matrix.}
+
+\item{k}{Filter kernel vector.}
+}
+\value{
+Matrix with same dimensions as \code{mat}, edge NA values set to 0.
+}
+\description{
+Column-wise 1-D convolution
+}
+\keyword{internal}
diff --git a/man/dot-smooth_rows.Rd b/man/dot-smooth_rows.Rd
new file mode 100644
index 0000000..4e75e18
--- /dev/null
+++ b/man/dot-smooth_rows.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.smooth_rows}
+\alias{.smooth_rows}
+\title{Row-wise 1-D convolution}
+\usage{
+.smooth_rows(mat, k)
+}
+\arguments{
+\item{mat}{Numeric matrix.}
+
+\item{k}{Filter kernel vector.}
+}
+\value{
+Matrix with same dimensions as \code{mat}, edge NA values set to 0.
+}
+\description{
+Row-wise 1-D convolution
+}
+\keyword{internal}
diff --git a/man/fit_spatial_mask.Rd b/man/fit_spatial_mask.Rd
new file mode 100644
index 0000000..a054abd
--- /dev/null
+++ b/man/fit_spatial_mask.Rd
@@ -0,0 +1,175 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fit_spatial_mask.R
+\name{fit_spatial_mask}
+\alias{fit_spatial_mask}
+\title{Fit a Spatial Boundary Mask to Point Coordinate Data}
+\usage{
+fit_spatial_mask(
+  coords,
+  method = "raster",
+  concavity = 2,
+  length_threshold = 0,
+  kde_bandwidth = NULL,
+  kde_threshold = 0.05,
+  kde_resolution = 256L,
+  raster_resolution = 256L,
+  raster_sigma = NULL,
+  raster_threshold = 0.15,
+  raster_min_pts = 1L,
+  buffer_dist = 0,
+  smooth_mask = FALSE,
+  smooth_tolerance = NULL,
+  n_cores = 1L,
+  crs = NA,
+  plot = FALSE,
+  verbose = TRUE
+)
+}
+\arguments{
+\item{coords}{A data frame or matrix with columns \code{x} and \code{y}. If column
+names differ, the first two columns are used and renamed. Rows containing
+\code{NA} in \code{x} or \code{y} are silently dropped.}
+
+\item{method}{Character scalar. Masking algorithm. One of:
+\itemize{
+\item \code{"raster"} \emph{(default)} -- Gaussian smoothing on a rasterised grid,
+binary threshold, then GEOS union of "on" cells. Supports holes and
+islands. Recommended for most datasets.
+\item \code{"kde"} -- 2-D kernel density estimate via \code{\link[MASS:kde2d]{MASS::kde2d()}} contoured
+with \code{\link[isoband:isobands]{isoband::isobands()}}. Supports holes and islands. Requires
+packages \strong{MASS} and \strong{isoband}.
+\item \code{"concave"} -- Concave hull via \code{\link[concaveman:concaveman]{concaveman::concaveman()}}. No holes or
+islands. Requires package \strong{concaveman}.
+\item \code{"convex"} -- Convex hull via \code{\link[sf:geos_unary]{sf::st_convex_hull()}}. No holes or
+islands. No additional packages required.
+}}
+
+\item{concavity}{Numeric. Relative concavity for \code{method = "concave"}.
+Lower values produce tighter, more concave boundaries; higher values
+approach the convex hull. Default \code{2}.}
+
+\item{length_threshold}{Numeric. Minimum edge length for
+\code{method = "concave"}. Edges shorter than this are not included.
+Default \code{0}.}
+
+\item{kde_bandwidth}{Numeric scalar or length-2 vector giving the KDE
+bandwidth(s) in the x and y directions. \code{NULL} uses
+\code{\link[MASS:bandwidth.nrd]{MASS::bandwidth.nrd()}} applied independently to each axis. Only used
+when \code{method = "kde"}. Default \code{NULL}.}
+
+\item{kde_threshold}{Numeric in \verb{(0, 1)}. KDE probability quantile used
+as the lower contour level. Smaller values produce larger masks.
+Default \code{0.05}.}
+
+\item{kde_resolution}{Integer. Grid resolution for the KDE density
+estimate. Larger values give finer boundaries at the cost of memory.
+Default \code{256L}.}
+
+\item{raster_resolution}{Integer. Number of grid cells along each axis for
+\code{method = "raster"}. Default \code{256L}.}
+
+\item{raster_sigma}{Numeric. Standard deviation of the Gaussian kernel in
+\strong{coordinate units} (the same units as \code{x} and \code{y}). This is the key
+tuning parameter:
+\itemize{
+\item Larger \code{raster_sigma} -> holes fill in, islands merge, boundary smooths.
+\item Smaller \code{raster_sigma} -> holes and fine structure are preserved.
+\item \code{NULL} (default) -> auto-set to 3\% of the domain width.
+}}
+
+\item{raster_threshold}{Numeric in \verb{(0, 1)}. A grid cell is "inside" when
+its smoothed value exceeds \code{raster_threshold * max(smoothed_grid)}.
+Higher values shrink the mask; lower values grow it. Default \code{0.15}.}
+
+\item{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
+points are treated as empty. Default \code{1L}.}
+
+\item{buffer_dist}{Numeric. Distance to expand the mask after fitting via
+\code{\link[sf:geos_unary]{sf::st_buffer()}}. Zero means no expansion. Default \code{0}.}
+
+\item{smooth_mask}{Logical. If \code{TRUE}, the mask boundary is simplified with
+\code{\link[sf:geos_unary]{sf::st_simplify()}} after fitting. Default \code{FALSE}.}
+
+\item{smooth_tolerance}{Numeric. Douglas-Peucker tolerance passed to
+\code{\link[sf:geos_unary]{sf::st_simplify()}}. \code{NULL} -> auto-set to 1\% of the bounding-box
+diagonal. Default \code{NULL}.}
+
+\item{n_cores}{Integer. Number of cores for parallel row-smoothing in
+\code{\link[parallel:mclapply]{parallel::mclapply()}} (Unix/macOS only; ignored on Windows). Default
+\code{1L}.}
+
+\item{crs}{CRS specification accepted by \code{\link[sf:st_crs]{sf::st_crs()}}: an EPSG integer,
+a WKT string, a PROJ string, or \code{NA} for no CRS. Default \code{NA}.}
+
+\item{plot}{Logical. If \code{TRUE}, a quick \link[ggplot2:ggplot2-package]{ggplot2::ggplot2} scatter + polygon plot is
+printed. Requires \strong{ggplot2} to be installed. Default \code{FALSE}.}
+
+\item{verbose}{Logical. If \code{TRUE}, prints a summary table and intermediate
+messages (e.g., \code{raster_sigma} in grid-cell units). Default \code{TRUE}.}
+}
+\value{
+An \link[sf:sfc]{sf::sfc} geometry collection of class \code{sfc_POLYGON} or
+\code{sfc_MULTIPOLYGON}. The CRS is set to \code{crs}. All input points (after NA
+removal) are guaranteed to lie inside or on the boundary of the returned
+mask. If any points fall outside the raw fitted mask a corrective buffer
+is applied automatically and a warning is issued.
+}
+\description{
+Fits a polygon (or multi-polygon) spatial mask (an \link[sf:sfc]{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.
+
+Four algorithms are available:\tabular{llll}{
+   \code{method} \tab Topology \tab Speed \tab Best for \cr
+   \code{"raster"} \tab Holes + islands \tab Fast \tab Most use cases (default) \cr
+   \code{"kde"} \tab Holes + islands \tab Moderate \tab Dense, smooth datasets \cr
+   \code{"concave"} \tab No holes \tab Fast \tab Simple non-convex shapes \cr
+   \code{"convex"} \tab No holes \tab Instant \tab Quick checks \cr
+}
+
+
+\strong{Raster method overview} (\code{method = "raster"}, v3):
+\enumerate{
+\item Bin point coordinates onto a regular grid.
+\item Convolve the occupancy grid with a 2-D Gaussian of width \code{raster_sigma}
+(in coordinate units).
+\item Threshold at \code{raster_threshold * max}. Cells above threshold are "inside".
+\item 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.
+\item Apply a morphological close (buffer out then in) to smooth staircase
+boundaries.
+}
+}
+\examples{
+set.seed(42)
+# Two-cluster point cloud
+coords <- data.frame(
+  x = c(rnorm(200, 0, 3), rnorm(200, 10, 3)),
+  y = c(rnorm(200, 0, 3), rnorm(200, 10, 3))
+)
+
+# Raster method (recommended default)
+mask <- fit_spatial_mask(coords, verbose = FALSE)
+
+# Convex hull -- instant, no additional packages
+mask_cv <- fit_spatial_mask(coords, method = "convex", verbose = FALSE)
+
+# Concave hull -- requires 'concaveman'
+\donttest{
+mask_cc <- fit_spatial_mask(coords, method = "concave",
+                            concavity = 1.5, verbose = FALSE)
+}
+
+# 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_covered_by(pts, sf::st_union(mask), sparse = FALSE)[, 1])
+
+}
+\seealso{
+\code{\link[sf:sfc]{sf::st_sfc()}}, \code{\link[sf:geos_combine]{sf::st_union()}}, \code{\link[sf:geos_unary]{sf::st_buffer()}}
+}
diff --git a/tests/testthat/test-fit_spatial_mask.R b/tests/testthat/test-fit_spatial_mask.R
index a1e0830..bdf1979 100644
--- a/tests/testthat/test-fit_spatial_mask.R
+++ b/tests/testthat/test-fit_spatial_mask.R
@@ -84,8 +84,8 @@ test_that("all points are contained in the raster mask", {
   coords    <- make_two_clusters()
   mask      <- fit_spatial_mask(coords, method = "raster", verbose = FALSE)
   pts_sfc   <- sf::st_as_sf(coords, coords = c("x", "y"), crs = sf::NA_crs_)
-  contained <- sf::st_within(sf::st_geometry(pts_sfc),
-                              sf::st_union(mask), sparse = FALSE)
+  contained <- sf::st_covered_by(sf::st_geometry(pts_sfc),
+                                  sf::st_union(mask), sparse = FALSE)
   expect_true(all(contained[, 1L]))
 })
 
@@ -93,8 +93,8 @@ test_that("all points are contained in the convex mask", {
   coords    <- make_circle()
   mask      <- fit_spatial_mask(coords, method = "convex", verbose = FALSE)
   pts_sfc   <- sf::st_as_sf(coords, coords = c("x", "y"), crs = sf::NA_crs_)
-  contained <- sf::st_within(sf::st_geometry(pts_sfc),
-                              sf::st_union(mask), sparse = FALSE)
+  contained <- sf::st_covered_by(sf::st_geometry(pts_sfc),
+                                  sf::st_union(mask), sparse = FALSE)
   expect_true(all(contained[, 1L]))
 })
 
@@ -240,10 +240,13 @@ test_that("two-cluster cloud produces a multi-part or larger mask", {
 # Verbose output
 # ---------------------------------------------------------------------------
 
-test_that("verbose = FALSE suppresses all output", {
-  coords <- make_circle()
-  expect_silent(
-    fit_spatial_mask(coords, method = "convex", verbose = FALSE)
+test_that("verbose = FALSE suppresses messages and printed output", {
+  # Use raster on a dense cloud — raster does not produce warnings in normal use.
+  coords <- make_two_clusters()
+  expect_no_message(
+    suppressWarnings(
+      fit_spatial_mask(coords, method = "raster", verbose = FALSE)
+    )
   )
 })
 
diff --git a/vignettes/advanced-topology.Rmd b/vignettes/advanced-topology.Rmd
index bbf1ace..fb3db83 100644
--- a/vignettes/advanced-topology.Rmd
+++ b/vignettes/advanced-topology.Rmd
@@ -241,7 +241,7 @@ same units as `x` and `y`). This makes it directly interpretable:
 
 ```{r sigma-sweep, fig.height=5}
 # Donut cloud — clear hole with adjustable fill behaviour
-sigmas <- c(0.15, 0.3, 0.6, 1.2, 2.5, 5.0)
+sigmas <- c(0.15, 0.3, 0.6, 1.0, 1.4, 1.7)
 pal    <- c("#2c3e50", "#2980b9", "#27ae60",
             "#f39c12", "#e74c3c", "#8e44ad")
 
@@ -269,9 +269,10 @@ wrap_plots(sw_plots, ncol = 3) +
 ```
 
 **Interpretation:** The hole is visible at sigma = 0.15 – 0.6. Around
-sigma = 1.2 it begins to fill. By sigma = 5.0 the annulus has merged into
-a solid disc. Choose `raster_sigma` to be just larger than the typical
-inter-point gap, but smaller than the smallest void you want to preserve.
+sigma = 1.0–1.4 it begins to fill. By sigma = 1.7 (roughly 17% of the
+domain width) the annulus has merged into a solid disc. Choose `raster_sigma`
+to be just larger than the typical inter-point gap, but smaller than the
+smallest void you want to preserve.
 
 A useful heuristic: set `raster_sigma` to roughly 2–3× the typical
 nearest-neighbour distance in your point cloud.
diff --git a/vignettes/getting-started.Rmd b/vignettes/getting-started.Rmd
index 6e9c0b9..2207931 100644
--- a/vignettes/getting-started.Rmd
+++ b/vignettes/getting-started.Rmd
@@ -237,7 +237,7 @@ coords_gap <- rbind(
   data.frame(x = rnorm(300, 8, 1.5), y = rnorm(300, 0, 1.5))
 )
 
-sigmas <- c(0.5, 1.5, 3.5)
+sigmas <- c(0.5, 1.5, 3.0)
 pal2   <- c("#2980b9", "#27ae60", "#8e44ad")
 sw2 <- mapply(function(sig, col) {
   m <- fit_spatial_mask(coords_gap, method = "raster",
@@ -316,10 +316,11 @@ mask_test <- fit_spatial_mask(coords_test, method = "raster", verbose = FALSE)
 
 pts_sfc   <- sf::st_as_sf(coords_test, coords = c("x", "y"),
                            crs = sf::NA_crs_)
-contained <- sf::st_within(sf::st_geometry(pts_sfc),
-                            sf::st_union(mask_test), sparse = FALSE)
+# Use st_covered_by (not st_within): boundary points count as covered.
+contained <- sf::st_covered_by(sf::st_geometry(pts_sfc),
+                                sf::st_union(mask_test), sparse = FALSE)
 
-cat("Points inside mask:", sum(contained[, 1]),
+cat("Points inside or on mask boundary:", sum(contained[, 1]),
     "/", nrow(coords_test), "\n")
 ```