Nothing
#' @rdname subtr_baseline
#' @title Automated background subtraction for spectral data
#'
#' @description
#' This baseline correction routine iteratively finds the baseline of a spectrum
#' using a polynomial fitting or accepts a manual baseline.
#'
#' @details
#' This is a translation of Michael Stephen Chen's MATLAB code written for the
#' \code{imodpolyfit} routine from Zhao et al. 2007.
#'
#' @param x a list object of class \code{OpenSpecy}.
#' @param type one of \code{"polynomial"} or \code{"manual"} depending on
#' whether you want spectra to be corrected with a manual baseline or with
#' polynomial baseline fitting.
#' @param degree the degree of the polynomial. Must be less than the number of
#' unique points when raw is \code{FALSE}. Typically a good fit can be
#' found with a 8th order polynomial.
#' @param raw if \code{TRUE}, use raw and not orthogonal polynomials.
#' @param baseline an \code{OpenSpecy} object containing the baseline data to be
#' subtracted.
#' @param make_rel logical; if \code{TRUE} spectra are automatically normalized
#' with \code{\link{make_rel}()}.
#' @param \ldots further arguments passed to \code{\link[stats]{poly}()}.
#'
#' @return
#' \code{subtr_baseline()} returns a data frame containing two columns named
#' \code{"wavenumber"} and \code{"intensity"}.
#'
#' @examples
#' data("raman_hdpe")
#'
#' subtr_baseline(raman_hdpe)
#'
#' @author
#' Win Cowger, Zacharias Steinmetz
#'
#' @seealso
#' \code{\link[stats]{poly}()};
#' \code{\link{smooth_intens}()}
#'
#' @references
#' Chen MS (2020). Michaelstchen/ModPolyFit. \emph{MATLAB}.
#' Retrieved from \url{https://github.com/michaelstchen/modPolyFit}
#' (Original work published July 28, 2015)
#'
#' Zhao J, Lui H, McLean DI, Zeng H (2007). “Automated Autofluorescence
#' Background Subtraction Algorithm for Biomedical Raman Spectroscopy.”
#' \emph{Applied Spectroscopy}, \strong{61}(11), 1225–1232.
#' \doi{10.1366/000370207782597003}.
#'
#' @importFrom stats terms model.frame sd lm poly approx
#' @importFrom data.table .SD
#' @export
subtr_baseline <- function(x, ...) {
UseMethod("subtr_baseline")
}
#' @rdname subtr_baseline
#'
#' @export
subtr_baseline.default <- function(x, ...) {
stop("object 'x' needs to be of class 'OpenSpecy'", call. = F)
}
#' @rdname subtr_baseline
#'
#' @export
subtr_baseline.OpenSpecy <- function(x, type = "polynomial",
degree = 8, raw = FALSE,
baseline, make_rel = TRUE, ...) {
if(type == "polynomial") {
sbg <- x$spectra[, lapply(.SD, .subtr_bl_poly, x = x$wavenumber,
degree = degree, raw = raw, ...)]
} else if(type == "manual") {
if(!is_OpenSpecy(baseline))
stop("'baseline' needs to be of class 'OpenSpecy'", call. = F)
sbg <- x$spectra[, lapply(.SD, .subtr_bl_manual, x = x$wavenumber,
bl_y = baseline$spectra[[1]],
bl_x = baseline$wavenumber)]
} else stop("'type' must be either 'polynomial' or 'manual'", call. = F)
if (make_rel) x$spectra <- sbg[, lapply(.SD, make_rel)] else x$spectra <- sbg
return(x)
}
.subtr_bl_poly <- function(y, x, degree, raw, ...) {
xout <- x
yin <- y
dev_prev <- 0 # standard deviation residuals for the last iteration of polyfit;
# set initially to 0
first_iter <- TRUE
criteria_met <- FALSE
while (!criteria_met) {
# Predict the intensity using the polynomial of specified length
paramVector <- lm(y ~ stats::poly(x, degree = degree, raw = raw, ...))
residual <- paramVector$residuals
mod_poly <- paramVector$fitted.values
dev_curr <- sd(residual)
# Remove peaks
if (first_iter) {
peaks <- c()
for (i in 1:length(y)) {
if (y[i] > mod_poly[i] + dev_curr) {
peaks <- c(peaks,i)
}
}
y <- y[-peaks]
mod_poly <- mod_poly[-peaks]
x <- x[-peaks]
first_iter <- FALSE
}
# Replace data with lower value if polynomial is lower
for (j in 1:length(y)) {
if (mod_poly[j] + dev_curr > y[j]) {
y[j] <- y[j]
} else {
y[j] <- mod_poly[j]
}
}
# Test criteria
criteria_met <- abs((dev_curr - dev_prev) / dev_curr) <= 0.05
# Approximate the intensity back to the original wavelengths, allows below
# the peak to be interpolated
if(criteria_met) {
ysbg <- approx(x, y, xout = xout, rule = 2, method = "linear",
ties = mean)[2] |>
unlist() |>
unname()
return(yin - ysbg)
}
# Update previous residual metric
dev_prev <- dev_curr
}
}
.subtr_bl_manual <- function(y, x, bl_y, bl_x, ...) {
y - approx(bl_x, bl_y, xout = x, rule = 2, method = "linear", ties = mean)$y
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.