main.R

####
# Test steps on single image
####

# global imports
import::from("cluster", "ellipsoidhull", .character_only = TRUE)
import::from(
  "imager",
  c("grabRect", "grayscale", "grow", "load.image", "shrink",
    "split_connected", "%inr%"),
  .character_only = TRUE
)
import::from("magrittr", c("%>%", "%$%"), .character_only = TRUE)
import::from("mclust", "densityMclust", .character_only = TRUE)
import::from("Momocs", c("import_jpg", "Out"), .character_only = TRUE)
import::from("purrr", c("discard", "keep", "map"), .character_only = TRUE)

# local imports
import::from(
  file.path("scripts", "image.R"),
  c("binarize", "correct_illumination", "invert_grayscale"),
  .character_only = TRUE
)
import::from(
  file.path("scripts", "pixset.R"),
  c("combine", "intersect_borders", "combine_bis", "rotation_angle"),
  .character_only = TRUE
)
import::from(
  file.path("scripts", "utils.R"),
  c("ehull", "cm_to_pixel"),
  .character_only = TRUE
)


# Load image -------------------------------------------------------------------

fpath <- file.path("data", "Example", "DSC_6177.JPG")
img <- load.image(fpath)


# Remove tray edges manually ---------------------------------------------------

# try to isolate tray edges
img_reduc <- grabRect(img, output = "im")


# Binarize ---------------------------------------------------------------------

# convert to grayscale
img_gr <- grayscale(img_reduc)
plot(img_gr, main = "grayscale")

# correct illumination and invert grayscale
img_gr_inv <- img_gr %>%
  correct_illumination() %>%
  invert_grayscale()
plot(img_gr_inv, main = "grayscale inverted")

# binarize (alternative is to use the imager::threshold function)
img_pix <- binarize(img_gr_inv, quantile = 0.95)
plot(img_pix, main = "binarized")


# Extract connected components and ruler ---------------------------------------

set_pix <- img_pix %>%
  # small dilatation to improve components detection
  grow(2L) %>%
  # extract components
  split_connected(., high_connectivity = TRUE)

# there are 3 kind of components: small disconnected bits (antenna, legs),
# medium (geriss), one large (ruler)
# fit mixture of gaussian distribution with 2 comps
set_size <- set_pix %>% lapply(sum) %>% unlist()
gm <- densityMclust(log10(set_size), G = 2L, plot = FALSE)
thresholds <- qnorm(
  c(5e-3, 1L - 5e-3),
  gm$parameters$mean[2L],
  sqrt(rev(gm$parameters$variance$sigmasq)[1L])
)
plot(gm, what = "density", data = log10(set_size), breaks = 200L)
abline(v = thresholds, col = "purple", lty = 2L)

# split components then shrink back
set_pix_small <- set_pix %>%
  keep(~ sum(.) %inr% c(2L, 10L**thresholds[1L])) %>%
  map(shrink, 2L)
set_pix_med <- set_pix %>%
  keep(~ (sum(.) %inr% 10L**thresholds)) %>%
  map(shrink, 2L)
ruler <- set_pix %>%
  keep((~ sum(.) > 10L**thresholds[2L])) %>%
  map(shrink, 2L)
rm(set_pix)

# remove pixsets that intersect the image borders
set_pix_med <- set_pix_med %>%
  discard(~ intersect_borders(., img_gr))

# plot first medium component
set_pix_med[[1L]] %>%
  imager::autocrop() %>%
  imager::pad(50L, "xy") %>%
  plot(main = "first medium component")


# Re-attach small component to medium ones -------------------------------------

set_pix_med <- set_pix_med %>%
  map(~ combine(., set_pix_small))
rm(set_pix_small)

set_pix_med[[1L]] %>%
  imager::autocrop() %>%
  imager::pad(50L, "xy") %>%
  plot(main = "re-attached first component")


# Get conversion factor from ruler  --------------------------------------------

# get the conversion factor cm to pixel
conv_factor <- cm_to_pixel(ruler[[1L]])

ruler_ct <- ruler[[1L]] %>%
  imager::contours(nlevels = 1L)
plot(
  ruler[[1L]], main = "ruler", xlim = c(1200L, 1800L), ylim = c(1500L, 1000L)
)
ruler_ct[[1L]] %$% {points(x, y, col = "red")}
ruler_ct[[3L]] %$% {points(x, y, col = "green")}


# Split superimposed components  -----------------------------------------------

# TODO create function that split connected pixset using morphological opening
# illustration with two elements

# morphological opening (erosion o dilation) to isolate abdomen
pixset_l <- set_pix_med[[4L]] %>%
  imager::mopening_square(6L) %>%
  imager::split_connected(high_connectivity = TRUE) %>%
  keep(~ (sum(.) > 1e2L))
plot(
  pixset_l %>% imager::parany() %>%
    imager::autocrop() %>% imager::pad(50L, "xy")
)

# further split apart the remaining parts
pixset_part <- (set_pix_med[[4L]] - pixset_l %>% imager::parany()) %>%
  imager::split_connected(high_connectivity = TRUE)
plot(
  (set_pix_med[[4L]] - pixset_l %>% imager::parany()) %>%
    imager::autocrop() %>% imager::pad(50L, "xy")
)

# try to re-attach parts to each component
pixset_l <- pixset_l %>%
  map(~ combine_bis(., pixset_part, padding = 25L))
rm(pixset_part)
plot(pixset_l[[1L]] %>% imager::autocrop() %>% imager::pad(50L, "xy"))
plot(pixset_l[[2L]] %>% imager::autocrop() %>% imager::pad(50L, "xy"))

# TODO better alternative: fit an ellipsoid hull around each abdomen (also
# increase both semi-axis length) then check which part satisfy the ellipsoid
# equation -> re-attach back if so


# Rotate along x-axis ----------------------------------------------------------

# get the rotation angle along x-axis
theta <- rotation_angle(set_pix_med[[1L]], size = 6L)
eh <- ehull(set_pix_med[[1L]], 6L)

set_pix_med[[1L]] %>%
  plot(main = "ellipsoid hull", xlim = c(800L, 950L), ylim = c(150L, 1L))
set_pix_med[[1L]] %>%
  imager::mopening_square(6L) %>% imager::as.pixset() %>% imager::where() %$%
  {points(x, y, cex = 0.25, col = "green")}
lines(predict(eh), col = "red")

# rotated pixset
set_pix_med[[1L]] %>%
  imager::imrotate(theta * 180L / pi, eh$loc[1L], eh$loc[2L], 1L) %>%
  plot(main = "rotated", xlim = c(800L, 950L), ylim = c(150L, 1L))


# TODO: Need to further orientate the head towards the y-axis (±pi/2).


# Export component and ruler to JPG (for Momocs) -------------------------------

# create temp directory
tmp_dir <- file.path("data", "tmp")
dir.create(tmp_dir, showWarnings = FALSE)
# export component
for (idx in seq_along(set_pix_med)) {
  set_pix_med[[idx]] %>%
    imager::autocrop(.) %>%
    imager::pad(50L, "xy") %>%
    {max(.) - .} %>%
    imager::as.cimg(.) %>%
    imager::save.image(
      file.path(tmp_dir, paste0("comp", "_", idx, ".jpg")),
      quality = 1.0
    )
}
# export ruler
ruler %>%
  .[[1L]] %>%
  imager::autocrop(.) %>%
  imager::pad(20L, "xy") %>%
  {max(.) - .} %>%
  imager::as.cimg(.) %>%
  imager::save.image(file.path(tmp_dir, "ruler.jpg"), quality = 1.0)


# Measurements through Momocs --------------------------------------------------

# import jpg file
coo <- import_jpg(jpg.paths = file.path(tmp_dir, "comp_1.jpg"))
# find contour
ct_bis <- Out(coo)
# plot contour
ct_bis[1L] %>% Momocs::coo_plot()

# TODO: some useful functions to investigate:
#  * Momocs::coo_oscillo : shape analysis (Fourier elliptic and the likes)
#
#  * Momocs::coo_intersect_segment, Momocs::coo_intersect_angle : find points of
#    a contour that intersect with a line defined by a segment (or an angle)
#
#  * Momocs::coo_right, Momocs::coo_left, Momocs::coo_top, Momocs::coo_down :
#    to keep only the right (left, top or down) part of the contour
#
#  * Momocs::coo_untiltx : correct rotational biases appearing after applying
#    sliding methods.
#
# There are more interesting functions, see the doc:
# https://cran.r-project.org/web/packages/Momocs/Momocs.pdf