R/fit_cie_model.R

In rcaiman: CAnopy IMage ANalysis

#' Fit CIE sky model
#'
#' @description
#' Fit the CIE sky model to data sampled from a canopy photograph using
#' general-purpose optimization.
#'
#' @details
#' The method is based on \insertCite{Lang2010;textual}{rcaiman}.
#' For best results, the input data should show a linear relation between
#' digital numbers and the amount of light reaching the sensor. See
#' [extract_radiometry()] and [read_caim_raw()] for details.
#' As a compromise solution, [invert_gamma_correction()] can be used.
#'
#' @param rr list, typically the output of [extract_rr()].
#'   If generated by other means, it must contain:
#'   \describe{
#'     \item{`zenith_dn`}{numeric vector of length one.}
#'     \item{`sky_points`}{`data.frame` with columns `a` (azimuth, deg),
#'       `z` (zenith, deg), and `rr` (relative radiance).}
#'   }
#' @param sun_angles named numeric vector of length two, with components
#'   `z` and `a` in degrees, e.g., `c(z = 49.3, a = 123.1)`.
#'   See [estimate_sun_angles()] for details.
#' @param custom_sky_coef numeric vector of length five, or numeric matrix
#'   with five columns. Custom starting coefficients for optimization.
#'   If not provided, coefficients are initialized from standard skies.
#' @param std_sky_no numeric vector. Standard sky numbers as in [cie_table].
#'   If not provided, all are used.
#' @param general_sky_type character vector of length one. Must be `"Overcast"`,
#'   `"Clear"`, or `"Partly cloudy"`. See column `general_sky_type` in
#'   [cie_table] for details. If not provided, all sky types are used.
#' @param method character vector. Optimization methods passed to
#'   [stats::optim()]. See that function for supported names.
#'
#' @return List with the following components:
#' \describe{
#'   \item{`rr`}{The input `rr` with an added `pred` column in
#'     `sky_points`, containing predicted values.}
#'   \item{`opt_result`}{List returned by [stats::optim()].}
#'   \item{`coef`}{Numeric vector of length five. CIE model coefficients.}
#'   \item{`sun_angles`}{Numeric vector of length two. Sun zenith and azimuth
#'     (degrees).}
#'   \item{`method`}{Character string. Optimization method used.}
#'   \item{`start`}{Numeric vector of length five. Starting parameters.}
#'   \item{`metric`}{Numeric value. Mean squared deviation as in
#'     \insertCite{Gauch2003;textual}{rcaiman}.}
#' }
#'
#' @references \insertAllCited{}
#'
#' @export
#'
#' @section Background:
#' This function is based on \insertCite{Lang2010;textual}{rcaiman}. In theory,
#' the best result would be obtained with data showing a linear relation between
#' digital numbers and the amount of light reaching the sensor. See
#' [extract_radiometry()] and [read_caim_raw()] for further details. As a
#' compromise solution, [invert_gamma_correction()] can be used.
#'
#' @section Digitizing sky points with ImageJ:
#' The [point selection tool of ‘ImageJ’
#' software](https://imagej.net/ij/docs/guide/146-19.html#sec:Multi-point-Tool)
#' can be used to manually digitize points and create a CSV file from which to
#' read coordinates (see *Examples*). After digitizing the points on the image,
#' this is a recommended workflow: 1. Use the dropdown menu *Analyze > Measure*
#' to open the Results window. 2. Use *File > Save As...* to obtain the CSV
#' file.
#'
#' Use this code to create the input `sky_points` from ImageJ data:
#' ````
#' sky_points <- read.csv(path)
#' sky_points <- sky_points[c("Y", "X")]
#' colnames(sky_points) <- c("row", "col")
#' ````
#'
#' @section Digitizing sky points with QGIS:
#' To use the [QGIS software](https://qgis.org/) to manually digitize
#' points, drag and drop the image in an empty project,
#' create an new vector layer, digitize points manually, save the editions, and
#' close the project.
#'
#' To create the new vector layer, this is a recommended workflow:
#' 1. Fo to the dropdown menu *Layer > Create Layer > New Geopackage Layer...*
#' 2. Choose "point" in the Geometry type dropdown list
#' 3. Make sure the CRS is EPSG:7589.
#' 4. Click on the *Toogle Editing* icon
#' 5. Click on the *Add Points Feature* icon.
#'
#' Use this code to create the input `sky_points` from QGIS data:
#' ````
#' sky_points <- terra::vect(path)
#' sky_points <- terra::extract(caim, sky_points, cells = TRUE)
#' sky_points <- terra::rowColFromCell(caim, sky_points$cell) %>% as.data.frame()
#' colnames(sky_points) <- c("row", "col")
#' ````
#'
#' @examples
#' \dontrun{
#' caim <- read_caim()
#' z <- zenith_image(ncol(caim), lens())
#' a <- azimuth_image(z)
#'
#' # Manual method following Lang et al. (2010)
#' path <- system.file("external/sky_points.csv",
#'                     package = "rcaiman")
#' sky_points <- read.csv(path)
#' sky_points <- sky_points[c("Y", "X")]
#' colnames(sky_points) <- c("row", "col")
#' head(sky_points)
#' plot(caim$Blue)
#' points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)
#'
#' # x11()
#' # plot(caim$Blue)
#' # sun_angles <- click(c(z, a), 1) %>% as.numeric()
#' sun_angles <- c(z = 49.5, a = 27.4) #taken with above lines then hardcoded
#'
#' sun_row_col <- row_col_from_zenith_azimuth(z, a,
#'                                            sun_angles["z"],
#'                                            sun_angles["a"])
#' points(sun_row_col[2], nrow(caim) - sun_row_col[1],
#'        col = "yellow", pch = 8, cex = 3)
#'
#' rr <- extract_rr(caim$Blue, z, a, sky_points)
#'
#' set.seed(7)
#' model <- fit_cie_model(rr, sun_angles,
#'                            general_sky_type = "Clear")
#'
#' plot(model$rr$sky_points$rr, model$rr$sky_points$pred)
#' abline(0,1)
#' lm(model$rr$sky_points$pred~model$rr$sky_points$rr) %>% summary()
#'
#' sky <- cie_image(z, a, model$sun_angles, model$coef) * model$rr$zenith_dn
#'
#' plot(sky)
#' ratio <- caim$Blue/sky
#' plot(ratio)
#' plot(ratio > 1.05)
#' plot(ratio > 1.15)
#' }
fit_cie_model <- function(rr, sun_angles,
                          custom_sky_coef = NULL,
                          std_sky_no = NULL,
                          general_sky_type = NULL ,
                          method = c("Nelder-Mead", "BFGS", "CG", "SANN")) {

  # Validate inputs -----------------------------------------------------------
  # rr should be a list with a data frame component named `sky_points`
  if (!is.list(rr) || !("sky_points" %in% names(rr))) {
    stop("`rr` must be a list containing a component named `sky_points`.")
  }
  if (!is.data.frame(rr$sky_points)) {
    stop("`rr$sky_points` must be a data frame.")
  }
  # rr$sky_points should have the expected columns
  required_cols <- c("a", "z", "rr")
  if (!all(required_cols %in% names(rr$sky_points))) {
    stop(sprintf("`rr$sky_points` must contain columns %s.",
                 paste(sprintf('"%s"', required_cols), collapse = ", ")))
  }
  # sun_angles must be a named numeric vector of length 2 with names z and a
  if (!is.numeric(sun_angles) || length(sun_angles) != 2 ||
      !identical(names(sun_angles), c("z", "a"))) {
    stop("`sun_angles` must be a named numeric vector of length two with names 'z' and 'a' in that order.")
  }
  .assert_choice(general_sky_type, c("Overcast", "Partly cloudy", "Clear"),
                 allow_null = TRUE)
  if (!is.null(custom_sky_coef)) {
    if (is.vector(custom_sky_coef)) {
      if (!(is.numeric(custom_sky_coef) && length(custom_sky_coef) == 5)) {
        stop("`custom_sky_coef`, if vector, must be numeric of length five.")
      }
    } else if (is.matrix(custom_sky_coef)) {
      if (!(is.numeric(custom_sky_coef) && ncol(custom_sky_coef) == 5)) {
        stop("`custom_sky_coef`, if matrix, must be numeric with five columns.")
      }
    } else stop("`custom_sky_coef` must be NULL, vector of matrix.")
  }
  if (!is.null(std_sky_no)) {
    .check_vector(std_sky_no, "integerish", sign = "positive")
    if (max(std_sky_no) > 15) {
      stop("`std_sky_no` cannot be greater than 15.")
    }
  }

  # Manage the set of start parameter according to user choice
  skies <- rcaiman::cie_table
  if (!is.null(custom_sky_coef)) {
    if (is.vector(custom_sky_coef)) {
      skies <- skies[1,]
      skies[1, 1:5] <- custom_sky_coef
    } else {
      custom_sky_coef <- as.matrix(custom_sky_coef)
      stopifnot(is.numeric(custom_sky_coef))
      skies <- skies[1:nrow(custom_sky_coef),]
      skies[,1:5] <- custom_sky_coef
    }
    skies$general_sky_type <- "custom"
    skies$description <- "custom"
  }
  if (!is.null(std_sky_no)) {
    skies <- skies[std_sky_no,]
  } else if (!is.null(general_sky_type)) {
    indices <- skies$general_sky_type == general_sky_type
    skies <- skies[indices,]
  }

  # Retrieve coordinates and convert to radians
  AzP <- .degree2radian(rr$sky_points$a)
  Zp <- .degree2radian(rr$sky_points$z)

  vault <- expand.grid(seq(pi / 18, pi / 2, pi / 18),
                       seq(pi / 18, 2 * pi, pi / 18))
  vault <- rbind(matrix(c(0, 0), ncol = 2), as.matrix(vault))

  AzS <- .degree2radian(sun_angles["a"])
  Zs <- .degree2radian(sun_angles["z"])

  # When the c or e parameters are negative, the sun can be darker than the
  # sky, which is unrealistic. So, the code avoids that possibility by using
  # abs(). Dark suns can also be seen with a low values of a. It will depend
  # on the indicatrix function, but anything below -1 might be a problem.
  # Positive values of b create unrealistic values at the horizon. Positive
  # values of d produce negative values, which are unrealistc
  fcost <- function(params) {
    .a <- params[1]
    .b <- params[2]
    .c <- params[3]
    .d <- params[4]
    .e <- params[5]

    w <- 1e10
    penalty_param <- w * max(0, .b) + w * max(0, -.e)

    vault_values <- .cie_sky_model(vault[,2],
                                   vault[,1], AzS, Zs, .a, .b, .c, .d, .e)
    vault_values[vault_values > 0] <- 0

    x     <- .cie_sky_model(AzP, Zp, AzS, Zs, .a, .b, .c, .d, .e)
    x_sun <- .cie_sky_model(AzS, Zs, AzS, Zs, .a, .b, .c, .d, .e)
    penalty_behavior <- w * abs(.c) * max(0, max(x) - x_sun) +
      w * sum(vault_values)^2

    x <- .cie_sky_model(AzP, Zp, AzS, Zs, .a, .b, .c, .d, .e)
    residuals <- x - rr$sky_points$rr

    loss_value <- sqrt(mean(residuals^2)) / mean(rr$sky_points$rr)

    loss_value + penalty_behavior + penalty_param
  }
  fcost <- compiler::cmpfun(fcost)
  y <- rr$sky_points$rr

  # Try all start parameters (brute force approach)
  .optim_multi <- function(i, method) {
    start_params <- skies[i, 1:5] %>% as.numeric()
    .optim_single <- function(method) {
      fit <- tryCatch(
        stats::optim(par = start_params,
                     fn = fcost,
                     method = method),
        error = function(e) {
          warning("`optim` failed unexpectedly.")
          list(par = start_params, convergence = NULL)
        }
      )

      coef <- fit$par
      names(coef) <- c("a", "b", "c", "d", "e")
      pred <- .cie_sky_model(AzP, Zp, AzS, Zs,
                             .a = coef[1],
                             .b = coef[2],
                             .c = coef[3],
                             .d = coef[4],
                             .e = coef[5])

      #10.2134/agronj2003.1442
      x <- pred
      reg <- tryCatch(lm(x~y),
                      error = function(e) NULL,
                      warning = function(w) NULL)

      if (is.null(reg)) {
        MSE <- 1e10
      } else if (!is.na(summary(reg) %>% .$r.squared %>% suppressWarnings())){
        m <- stats::coef(reg)[2]
        r_squared <- summary(reg) %>% .$r.squared
        SB <- (mean(x) - mean(y))^2
        NU <- (1 - m)^2 * mean(x^2)
        LC <- (1 - r_squared) * mean(y^2)
        MSE <- SB + NU + LC
        MSE
      } else {
        MSE <- 1e10
      }

      rr$sky_points$pred <- x

      list(rr = rr,
           opt_result = fit,
           coef = coef,
           sun_angles = sun_angles,
           method = method,
           start = start_params,
           metric = unname(MSE))
    }
    lapply(method, .optim_single)
  }

  models <- Map(.optim_multi, 1:nrow(skies), method) %>% suppressWarnings()

  # Choose the best result
  ## best per parameter
  models <- lapply(models,
                   function(l) {
                        metrics <- lapply(seq_along(l),
                                          function(i) l[[i]]$metric) %>%
                                            unlist()
                        l[[which.min(metrics)]]
                   })
  ## best
  metrics <- lapply(seq_along(models), function(i) models[[i]]$metric) %>% unlist
  i <- which.min(metrics)
  models[[i]]
}