R/spc.rubberband.R

In hyperSpec: Work with Hyperspectral Data, i.e. Spectra + Meta Information (Spatial, Time, Concentration, ...)

##' Rubberband baseline
##'
##' Baseline with support points determined from a convex hull of the spectrum.
##' 
##' Use \code{debuglevel >= 1} to obtain debug plots, either directly via function argument or by setting hyperSpec's \code{debuglevel} option.
##' @title Rubberband baseline correction
##' @param spc hyperSpec object
##' @param ... further parameters handed to \code{\link[stats]{smooth.spline}}
##' @param upper logical indicating whether the lower or upper part of the hull should be used
##' @param noise noise level to be taken into account
##' @param spline logical indicating whether the baseline should be an interpolating spline through
##' the support points or piecewise linear.
##' @return hyperSpec object containing the baselines
##' @rdname spc-rubberband
##' @author Claudia Beleites
##' @seealso \code{\link[hyperSpec]{spc.fit.poly}}, \code{\link[hyperSpec]{spc.fit.poly.below}}
##' 
##' \code{vignette ("baseline")}
##' 
##' \code{\link[hyperSpec]{hy.setOptions}}
##'
##' @note This function is still experimental
##' @export
##' @examples
##' plot (paracetamol [,, 175 ~ 1800])
##' bl <- spc.rubberband (paracetamol [,, 175 ~ 1800], noise = 300, df = 20)
##' plot (bl, add = TRUE, col = 2)
##' 
##' plot (paracetamol [,, 175 ~ 1800] - bl)

spc.rubberband <- function (spc, ..., upper = FALSE, noise = 0, spline = TRUE){
  spc <- orderwl (spc)

  if (upper) spc@data$spc <- -spc@data$spc
  
  spc@data$spc <- .rubberband (spc@wavelength, spc@data$spc, 
                               noise = noise, spline = spline, ...)

  if (upper) spc@data$spc <- -spc@data$spc

  spc
}

##' @importFrom grDevices chull
.rubberband <- function (x, y, noise, spline, ..., debuglevel = hy.getOption ("debuglevel")){
  for (s in seq_len (nrow (y))){
    use <- which (!is.na (y [s,]))
    
    pts <- chull (x [use], y [s,use])
    pts <- use [pts]
    
    if (debuglevel >= 1L){
    	plot (x, y [s, ], type = "l")
    	points (x [pts], y [s, pts], pch = 1, col = matlab.dark.palette (length (pts)))
    }

    ## `chull` returns points in cw order
    ## => points between ncol (y) and 1 are lower part of hull
    imax <- which.max (pts) - 1
    
    ## if necessary, rotate pts so that ncol (y) is at position 1
    if (imax > 0L)
    	pts <- c (pts [- seq_len (imax)], pts [seq_len (imax)])

    ## now keep only pts until column index 1
    pts <- pts [1 : which.min (pts)]
    
    ## check whether first and last point are minima, 
    ## if not remove them.
    ## If they are minima, 2nd and 2nd last point do not appear in pts
    ## last point: 
    if (pts [2] == pts [1] - 1) pts <- pts [-1]
    
    ## now sort ascending (anyways needed later on)
    pts <- rev (pts)
    
    ## fist point:
    if (pts [2] == pts [1] + 1) pts <- pts [-1]
    
    if (debuglevel >= 1L){
    	points (x [pts], y [s, pts], pch = 19, col = matlab.dark.palette (length (pts)), cex = 0.7)
    }

    tmp <- approx (x = x [pts], y = y [s, pts], xout= x, method="linear")$y
    
    if (spline){
      pts <- which (y [s,] <= tmp + noise)

      if (length (pts) > 3)
        tmp <- predict (smooth.spline (x[pts], y[s, pts], ...)$fit, x, 0)$y 
      else 
        tmp <- spline (x [pts], y [s, pts], xout = x)$y
        
    }
    
    y [s, ] <- tmp
    
  }
  
  y  
}

.test (spc.rubberband) <- function (){
  context ("spc.rubberband")
  
  ## use data that yields fairly stable baseline solution
  paracetamol <- paracetamol [,, 300 ~ 550]
  
  test_that("spectrum containing NA inside", {
    tmp <- paracetamol
    tmp [[,, 400]] <- NA
    
    coefs <- spc.rubberband (tmp)
    expect_equal(
      coefs [[,, !is.na (tmp)]],
      spc.rubberband(paracetamol [,, !is.na (tmp)]) [[]]
    )
    
    ## bug was: all coefficients were silently 0 
    expect_true (all (abs (coefs [[]]) > sqrt (.Machine$double.eps)))
  })
  
  test_that ("spectrum containing NA at first wavelength (issue #95)", {
    tmp <- paracetamol
    tmp [[,, 1, wl.index = TRUE]] <- NA

    coefs <- spc.rubberband (tmp)
    expect_equal(
      coefs [[,, !is.na (tmp)]],
      spc.rubberband(paracetamol [,, !is.na (tmp)]) [[]]
    )
  })

  test_that ("spectrum containing NA at end", {
    tmp <- paracetamol [1]
    tmp [[,, nwl (paracetamol), wl.index = TRUE]] <- NA
    
    coefs <- spc.rubberband (tmp)
    expect_equal(
      coefs [[,, !is.na (tmp)]],
      spc.rubberband(paracetamol [1,, !is.na (tmp)]) [[]]
    )
  })
  
}