####
# Additional utility functions.
####

# global imports
import::here("cluster", "ellipsoidhull", .character_only = TRUE)
import::here(
  "imager",
  c(
    "as.pixset", "contours", "is.cimg", "is.pixset", "%inr%",
    "mopening_square", "where"
  ),
  .character_only = TRUE
)
import::here("mclust", "densityMclust", .character_only = TRUE)
import::here("purrr", c("keep", "map"), .character_only = TRUE)


#' Sample pixel values from histogram
#'
#' @param obj an array-like object (img, pixset)
#' @param ratio a real, the sampling ratio.
#'
#' @returns an integer vector
sample_histogram <- function(obj, ratio = 0.02) {
  n_samples <- min(1e4L, prod(dim(obj)[1L:2L]))
  n_breaks <- as.integer(n_samples * ratio)
  hist <- hist(obj, breaks = n_breaks, plot = FALSE)
  bins <- with(
    hist, sample(length(mids), n_samples, p = density, replace = TRUE)
  )
  sample <- runif(
    length(bins), hist$breaks[bins], hist$breaks[bins + 1L]
  )
  return(sample)
}


#' Square stencil
#'
#' Define a square stencil within the boundaries of an image cropping points
#' that are out of bounds.
#'
#' @param centroid an integer vector, the stencil centroid.
#' @param pad an integer, the padding around the stencil centroid.
#' @param obj an imager::cimg or an imager::pixset
#'
#' @returns a data.frame
square_stencil <- function(centroid, pad, obj) {

  # get image dimensions
  if (is.cimg(obj) | is.pixset(obj)) {
    dims <- as.integer(dim(obj)[1L:2L])
  } else {
    stop("Invalid 'obj' argument.")
  }
  # check bounds
  bounds <- c(
    ifelse(centroid - pad < 1L, -centroid + 1L, -pad),
    ifelse(pad + centroid > dims, dims - centroid, pad)
  )
  # define stencil
  stencil <- expand.grid(
    dx = seq(bounds[1L], bounds[3L]),
    dy = seq(bounds[2L], bounds[4L])
  )
  return(stencil)
}


#' Centimeter to pixel
#'
#' Get the centimeter to pixel conversion factor from a ruler defined as
#' a pixset.
#'
#' @param ruler an imager::pixset
#' @param quantile a real, the quantile level used for thresholding.
#'
#' @returns a data.frame
cm_to_pixel <- function(ruler, quantile = 5e-3) {

  # get ruler contours size
  ct_size <- ruler %>%
    contours(nlevels = 1L) %>%
    map(~ length(.$x)) %>%
    unlist()
  # estimate the distribution by a Gaussian
  gm <- densityMclust(log10(ct_size), G = 1L, plot = FALSE)
  # threshold to discriminate square contours
  thresholds <- qnorm(
    c(quantile, 1L - quantile),
    gm$parameters$mean[1L],
    sqrt(gm$parameters$variance$sigmasq[1L])
  )
  # get the conversion factor cm to pixel
  conv_factor <- ct_size %>%
    keep(~ (. %inr% 10L ** thresholds)) %>%
    {. / 4L} %>%
    mean()
  return(conv_factor)
}


#' Ellipsoid hull
#'
#' Fit an ellipsoid hull for a pixset with optional morphological opening.
#'
#' @param pixset an imager::pixset
#' @param size an integer, the morphological opening factor.
#'
#' @returns a data.frame
ehull <- function(pixset, size = 6L) {
  out <- pixset %>%
    mopening_square(size) %>%
    as.pixset() %>%
    where() %>%
    as.matrix() %>%
    ellipsoidhull()
  return(out)
}