From 26ca32b9d30fa1f3187f50fcac13032126d8187e Mon Sep 17 00:00:00 2001 From: "Win Cowger, PhD" Date: Thu, 26 Sep 2024 16:23:06 -0700 Subject: [PATCH] add spatial smoothing on objects --- NAMESPACE | 3 + R/spatial_smooth.R | 136 +++++++++++++++++++++++++++++++ man/spatial_smooth.Rd | 34 ++++++++ tests/testthat/test-adj_intens.R | 5 +- 4 files changed, 177 insertions(+), 1 deletion(-) create mode 100644 R/spatial_smooth.R create mode 100644 man/spatial_smooth.Rd diff --git a/NAMESPACE b/NAMESPACE index b8b9f9e5..44937bfe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -130,6 +130,7 @@ export(run_app) export(sample_spec) export(sig_noise) export(smooth_intens) +export(spatial_smooth) export(spec_res) export(split_spec) export(subtr_baseline) @@ -148,6 +149,8 @@ importFrom(data.table,is.data.table) importFrom(data.table,melt) importFrom(data.table,rbindlist) importFrom(data.table,setDT) +importFrom(data.table,setkey) +importFrom(data.table,setnames) importFrom(data.table,setorder) importFrom(data.table,transpose) importFrom(digest,digest) diff --git a/R/spatial_smooth.R b/R/spatial_smooth.R new file mode 100644 index 00000000..d20e3540 --- /dev/null +++ b/R/spatial_smooth.R @@ -0,0 +1,136 @@ +#' @title Spatial Smoothing of OpenSpecy Objects +#' +#' @description +#' Applies spatial smoothing to an `OpenSpecy` object using a Gaussian filter. +#' +#' @param x an `OpenSpecy` object. +#' @param sigma a numeric vector specifying the standard deviations for the +#' Gaussian kernel in the x and y dimensions, respectively. +#' @param \dots further arguments passed to or from other methods. +#' +#' @details +#' This function performs spatial smoothing on the spectral data in an `OpenSpecy` object. +#' It assumes that the spatial coordinates are provided in the `metadata` element of the object, +#' specifically in the `x` and `y` columns, and that there is a `col_id` column in `metadata` that +#' matches the column names in the `spectra` data.table. +#' +#' @return +#' An `OpenSpecy` object with smoothed spectra. +#' +#' @author +#' Win Cowger +#' +#' @seealso +#' \code{\link{as_OpenSpecy}()}, \code{\link[mmand]{gaussianSmooth}()} +#' +#' @importFrom data.table data.table melt dcast setkey setnames +#' @importFrom mmand gaussianSmooth +#' @export +spatial_smooth <- function(x, sigma = c(1, 1, 1), ...) { + # Check that x is an OpenSpecy object + if (!inherits(x, "OpenSpecy")) { + stop("x must be an OpenSpecy object.", call. = FALSE) + } + + # Check that metadata contains 'x', 'y', and 'col_id' + if (is.null(x$metadata) || !all(c("x", "y", "col_id") %in% names(x$metadata))) { + stop("The OpenSpecy object must have 'x', 'y', and 'col_id' columns in its metadata.", call. = FALSE) + } + + # Extract wavenumbers and spectra + wavenumbers <- x$wavenumber + spectra <- x$spectra # data.table with spectra columns + + # Extract coordinates and ensure they are numeric + coords <- x$metadata[, .(col_id, x, y)] + #coords[, `:=`(x = as.numeric(x), y = as.numeric(y))] + + # Ensure that 'col_id's in metadata match column names in 'spectra' + if (!all(coords$col_id %in% names(spectra))) { + stop("Not all 'col_id's in metadata match column names in spectra.", call. = FALSE) + } + + # Reshape spectra data.table to long format + dt_long <- data.table::melt( + spectra, + measure.vars = names(spectra), + variable.name = "col_id", + value.name = "intensity" + ) + + # Add wavenumbers and coordinates to dt_long + dt_long[, wavenumber := rep(wavenumbers, times = ncol(spectra))] + # Merge coordinates into dt_long using 'col_id' + dt_long <- merge(dt_long, coords, by = "col_id", all.x = TRUE) + + # Check for any missing coordinates + if (any(is.na(dt_long$x) | is.na(dt_long$y))) { + stop("Missing spatial coordinates for some spectra.", call. = FALSE) + } + + # Get unique x, y, and wavenumber values + x_vals <- sort(unique(dt_long$x)) + y_vals <- sort(unique(dt_long$y)) + z_vals <- wavenumbers # wavenumbers are already sorted + + nx <- length(x_vals) + ny <- length(y_vals) + nz <- length(z_vals) + + # Map x, y, and wavenumber to indices + dt_long[, x_idx := match(x, x_vals)] + dt_long[, y_idx := match(y, y_vals)] + dt_long[, z_idx := match(wavenumber, z_vals)] + + # Initialize empty array + arr <- array(NA_real_, dim = c(nx, ny, nz)) + + # Fill the array with intensity values + idx <- cbind(dt_long$x_idx, dt_long$y_idx, dt_long$z_idx) + arr[idx] <- dt_long$intensity + + # Apply Gaussian smoothing (sigma for x and y; 0 for wavenumber dimension) + arr_smoothed <- mmand::gaussianSmooth(arr, sigma = sigma) + + # Convert the smoothed array back to data.table + idx_smoothed <- which(!is.na(arr_smoothed), arr.ind = TRUE) + dt_smoothed <- data.table::data.table( + x_idx = idx_smoothed[, 1], + y_idx = idx_smoothed[, 2], + z_idx = idx_smoothed[, 3], + intensity = arr_smoothed[idx_smoothed] + ) + + # Map indices back to x, y, and wavenumber + dt_smoothed[, x := x_vals[x_idx]] + dt_smoothed[, y := y_vals[y_idx]] + dt_smoothed[, wavenumber := z_vals[z_idx]] + + # Map back to 'col_id' using coordinates + # Create a lookup table for (x, y) to col_id + coord_to_colid <- coords[, .(x, y, col_id)] + setkey(coord_to_colid, x, y) + + # Merge col_id into dt_smoothed + dt_smoothed <- merge(dt_smoothed, coord_to_colid, by = c("x", "y"), all.x = TRUE) + + # Reshape data to wide format + dt_wide <- data.table::dcast( + dt_smoothed, + wavenumber ~ col_id, + value.var = "intensity" + ) + + # Ensure the columns are in the same order as original spectra + spectra_smoothed <- dt_wide[, names(spectra), with = FALSE] + + # Create new OpenSpecy object with smoothed spectra + x_smoothed <- as_OpenSpecy( + x = dt_wide$wavenumber, + spectra = spectra_smoothed, + metadata = x$metadata, + session_id = TRUE + ) + + return(x_smoothed) +} \ No newline at end of file diff --git a/man/spatial_smooth.Rd b/man/spatial_smooth.Rd new file mode 100644 index 00000000..9a9452d6 --- /dev/null +++ b/man/spatial_smooth.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial_smooth.R +\name{spatial_smooth} +\alias{spatial_smooth} +\title{Spatial Smoothing of OpenSpecy Objects} +\usage{ +spatial_smooth(x, sigma = c(1, 1, 1), ...) +} +\arguments{ +\item{x}{an \code{OpenSpecy} object.} + +\item{sigma}{a numeric vector specifying the standard deviations for the +Gaussian kernel in the x and y dimensions, respectively.} + +\item{\dots}{further arguments passed to or from other methods.} +} +\value{ +An \code{OpenSpecy} object with smoothed spectra. +} +\description{ +Applies spatial smoothing to an \code{OpenSpecy} object using a Gaussian filter. +} +\details{ +This function performs spatial smoothing on the spectral data in an \code{OpenSpecy} object. +It assumes that the spatial coordinates are provided in the \code{metadata} element of the object, +specifically in the \code{x} and \code{y} columns, and that there is a \code{col_id} column in \code{metadata} that +matches the column names in the \code{spectra} data.table. +} +\seealso{ +\code{\link{as_OpenSpecy}()}, \code{\link[mmand]{gaussianSmooth}()} +} +\author{ +Win Cowger +} diff --git a/tests/testthat/test-adj_intens.R b/tests/testthat/test-adj_intens.R index 4ddf4dc8..f75559c6 100644 --- a/tests/testthat/test-adj_intens.R +++ b/tests/testthat/test-adj_intens.R @@ -4,7 +4,7 @@ raman_hdpe$spectra$intensity2 <- raman_hdpe$spectra$intensity * 2 raman_hdpe$metadata <- rbind(raman_hdpe$metadata, raman_hdpe$metadata) test_that("adj_intens() handles input errors correctly", { - adj_intens(1:1000) |> expect_error() + adj_intens("abc") |> expect_error() }) test_that("adj_intens() works as expected", { @@ -13,6 +13,9 @@ test_that("adj_intens() works as expected", { adj_intens(raman_hdpe, type = "transmittance") |> expect_silent() adj_intens(raman_hdpe, type = "transmission") |> expect_error() + adj_intens(raman_hdpe, type = "reflectance", log_exp = "log") |> expect_silent() + adj_intens(raman_hdpe, type = "transmittance", log_exp = "exp") |> expect_silent() + expect_equal( cor(raman_hdpe$spectra$intensity, adj$spectra$intensity), 1, ignore_attr = F)