####
# Utility functions for imager::pixset objects.
####

import::here("dplyr", "inner_join", .character_only = TRUE)
import::here(
  "imager",
  c(
    "as.cimg", "as.pixset", "get.stencil", "imshift", "is.pixset", "parany",
    "px.borders", "where"
  ),
  .character_only = TRUE
)
import::here("magrittr", "%>%", .character_only = TRUE)

# local imports
import::here("utils.R", c("ehull", "square_stencil"), .character_only = TRUE)


#' Get centroid of a pixset.
#'
#' @param pixset an imager::pixset
#'
#' @returns an integer vector
get_centroid <- function(pixset) {
  centroid <- pixset %>%
    where() %>%
    colMeans() %>%
    as.integer()
  return(centroid)
}


#' Center a pixset.
#'
#' @param pixset an imager::pixset
#'
#' @returns an imager::pixset
px_center <- function(pixset) {
  centroid <- get_centroid(pixset)
  delta <- dim(pixset)[1L:2L] / 2L - centroid
  out <- pixset %>%
    as.cimg() %>%
    imshift(delta_x = delta[1L], delta_y = delta[2L]) %>%
    as.pixset()
  return(out)
}


#' Intersect borders
#'
#' Test whether a pixset intersect an image borders.
#'
#' @param pixset an imager::pixset
#' @param img an imager::cimg
#'
#' @returns a vector of boolean
intersect_borders <- function(pixset, img) {
  # extract the image borders (as pixels coordinates)
  borders <- img %>%
    px.borders() %>%
    where()
  # count the number of pixels intersecting the borders
  pix_on_borders <- suppressMessages(
    pixset %>%
      where() %>%
      inner_join(borders) %>%
      nrow()
  )
  return(pix_on_borders > 0L)
}


#' Combine pixset
#'
#' Combine pixsets into a reference pixset according to a shared neighbourhood.
#' The neighbourhood is defined as a squared stencil centered around each
#' candidate pixset.
#'
#' @param pixset_ref an imager::pixset, the reference.
#' @param pixset_list a list of imager::pixset, the list of pixsets to combine
#' with the reference.
#' @param padding an integer, the stencil length
#'
#' @returns a vector of boolean
combine <- function(pixset_ref, pixset_list, padding = 50L) {
  # convert candidat to list if there is only one
  if (!is.list(pixset_list) & is.pixset(pixset_list)) {
    pixset_list <- list(pixset_list)
  }
  # iterate over parts
  for (idx in seq_along(pixset_list)) {
    # get centroid
    centroid <- get_centroid(pixset_list[[idx]])
    # define stencil for neighborhood
    stencil <- square_stencil(centroid, padding, pixset_ref)
    # check overlap
    overlap <- pixset_ref %>%
      get.stencil(stencil, x = centroid[1L], y = centroid[2L]) %>%
      sum()
    # merge if overlap
    if (overlap > 0L) {
      pixset_ref <- parany(list(pixset_ref, pixset_list[[idx]]))
    }
  }
  return(pixset_ref)
}

combine_bis <- function(pixset_ref, pixset_list, padding = 25L) {
  # convert candidat to list if there is only one
  if (!is.list(pixset_list) & is.pixset(pixset_list)) {
    pixset_list <- list(pixset_list)
  }
  # get centroid
  centroid <- get_centroid(pixset_ref)
  # define stencil for neighborhood
  stencil <- square_stencil(centroid, padding, pixset_ref)
  # iterate over parts
  for (idx in seq_along(pixset_list)) {
    # check overlap
    overlap <- pixset_list[[idx]] %>%
      get.stencil(stencil, x = centroid[1L], y = centroid[2L]) %>%
      sum()
    # merge if overlap
    if (overlap > 0L) {
      pixset_ref <- parany(list(pixset_ref, pixset_list[[idx]]))
    }
  }
  return(pixset_ref)
}


#' Rotation angle
#'
#' Get the rotation angle to align the abdomen along the x-axis. The function
#' fits an ellipsoid hull around the pixset to derive the rotation angle.
#'
#' @param pixset an imager::pixset
#' @param size an integer, the morphological opening factor (optional)
#'
#' @returns a real
rotation_angle <- function(pixset, size = 6L) {
  # compute ellipsoid hull
  eh <- ehull(pixset, size)
  # get semi-axis lengths
  l_term <- (eh$cov[1L, 1L] + eh$cov[2L, 2L]) / 2L
  r_term <- sum(
    c((eh$cov[1L, 1L] - eh$cov[2L, 2L]) / 2L, eh$cov[1L, 2L]) ** 2L
  )
  lambda_1 <- l_term + sqrt(r_term)
  lambda_2 <- l_term - sqrt(r_term)
  axis_l <- c(sqrt(lambda_1), sqrt(lambda_2))
  # get rotation angle (from x-axis)
  if (eh$cov[1L, 2L] != 0L) {
    theta <- pi - atan2(lambda_1 - eh$cov[1L, 1L], eh$cov[1L, 2L])
  } else {
    theta <- ifelse(eh$cov[1L, 1L] >= eh$cov[2L, 2L], 0L, pi / 2L)
  }
  return(theta)
}