Skip to content
Snippets Groups Projects
Commit 06b3674b authored by Gilquin's avatar Gilquin
Browse files

init: add R script files

parent efd4d9ee
No related branches found
No related tags found
No related merge requests found
renv.lock 0 → 100644
This diff is collapsed.
library/
local/
cellar/
lock/
python/
sandbox/
staging/
This diff is collapsed.
{
"bioconductor.version": null,
"external.libraries": [],
"ignored.packages": [],
"package.dependency.fields": [
"Imports",
"Depends",
"LinkingTo"
],
"ppm.enabled": null,
"ppm.ignored.urls": [],
"r.version": null,
"snapshot.type": "implicit",
"use.cache": true,
"vcs.ignore.cellar": true,
"vcs.ignore.library": true,
"vcs.ignore.local": true,
"vcs.manage.ignores": true
}
####
# Plotting functions.
####
# global imports
import::here("dplyr", "filter", .character_only = TRUE)
import::here(
"ggplot2",
c(
"aes", "facet_wrap", "geom_histogram", "geom_raster", "ggplot",
"scale_x_continuous", "scale_y_continuous", "scale_fill_viridis_c"
),
.character_only = TRUE
)
import::here("magrittr", "%>%", .character_only = TRUE)
# local import
import::here("image.R", "img_to_df", .character_only = TRUE)
#' Channel histogram
#'
#' Histogram ggplot of each color channel of an image.
#'
#' @param img an imager::cimg
#' @param cc a character vector, subset of ("R", "G", "B").
#'
#' @returns a ggplot2 graph
channel_hist <- function(img, cc = c("R", "G", "B")) {
if (all(cc %in% c("R", "G", "B"))) {
gg <- img %>%
img_to_df() %>%
filter(channel %in% cc) %>%
ggplot(aes(value, col = channel)) +
geom_histogram(bins = 30L) +
facet_wrap(~channel)
} else {
stop("Invalid 'cc' argument.")
}
return(gg)
}
#' Raster channel
#'
#' Raster ggplot of an image channel.
#'
#' @param img an imager::cimg
#' @param cc character, specify the image channel, one of "R", "G", "B".
#'
#' @returns a ggplot2 graph
raster_channel <- function(img, cc = NULL) {
if (cc %in% c("R", "G", "B")) {
gg <- img %>%
img_to_df() %>%
filter(channel == cc) %>%
ggplot(aes(x, y, fill = value)) +
geom_raster() +
scale_x_continuous(expand = c(0L, 0L)) +
scale_y_continuous(expand = c(0L, 0L), trans = scales::reverse_trans()) +
scale_fill_viridis_c(direction = 1L, option = "plasma")
} else {
stop("Invalid 'cc' argument.")
}
return(gg)
}
####
# Utility functions for imager::cimg objects.
####
# global imports
import::here("dplyr", c("mutate", "sample_n"), .character_only = TRUE)
import::here(
"imager",
c("grayscale", "HSVtoRGB", "imappend", "imsplit", "map_il", "RGBtoHSV"),
.character_only = TRUE
)
import::here("magrittr", "%>%", .character_only = TRUE)
import::here("mclust", "densityMclust", .character_only = TRUE)
import::here("purrr", "modify_at", .character_only = TRUE)
# local imports
import::here("utils.R", "sample_histogram", .character_only = TRUE)
#' DataFrame conversion
#'
#' Convert image to dataframe and expose color channel.
#'
#' @param img an imager::cimg
#'
#' @returns a data.frame
img_to_df <- function(img) {
out <- img %>%
as.data.frame() %>%
mutate(channel = factor(cc, labels = c("R", "G", "B")))
return(out)
}
#' Histogram equalization
#'
#' Flatten histogram by replacing the pixel value of an image by their rank.
#'
#' @param img an imager::cimg
#'
#' @returns an imager::cimg
hist_eq <- function(img) {
return(as.cimg(ecdf(img)(img), dim = dim(img)))
}
#' Enhance contrast
#'
#' Enhance the contrasts of an image by running an histogram equalization
#' separately on each channel and combining the results.
#'
#' @param img an imager::cimg
#'
#' @returns an imager::cimg
enhance_contrast <- function(img) {
out <- img %>%
imsplit("cc") %>%
map_il(hist_eq) %>%
imappend("cc")
return(out)
}
#' Reduce saturation
#'
#' Reduce the saturation of an image through HSV conversion.
#'
#' @param img an imager::cimg
#' @param ratio an integer, how much to divide the saturation by.
#'
#' @returns an imager::cimg
desaturation <- function(img, ratio = 2L) {
out <- img %>%
RGBtoHSV() %>%
imsplit("cc") %>%
modify_at(2L, ~ . / ratio) %>%
imappend("cc") %>%
HSVtoRGB()
return(out)
}
#' Correct illumination
#'
#' Correct a gray-scaled image illumination by fitting a linear model and
#' removing the spatial trend.
#'
#' @param img an imager::cimg
#' @param nsamples an integer, pixel subsampling value.
#'
#' @returns an imager::cimg object
correct_illumination <- function(img, nsamples = 1e4L) {
# convert to grayscale if needed
if (rev(dim(img))[1L] > 1L) {
img <- grayscale(img)
}
# linear regression trend
trend <- img %>%
as.data.frame() %>%
sample_n(nsamples) %>%
lm(value ~ x * y, data = .) %>%
predict(img)
out <- img - trend
return(out)
}
#' Invert grayscale image
#'
#' @param img an imager::cimg
#'
#' @returns an imager::cimg
invert_grayscale <- function(img) {
# convert to grayscale if needed
if (rev(dim(img))[1L] > 1L) {
img <- grayscale(img)
}
out <- max(img) - img
return(out)
}
#' Binarize
#'
#' Binarize a grayscale image.
#'
#' @param img an imager::cimg
#' @param quantile a real, the quantile level used for thresholding.
#' @param ... mclust::densityMclust parameters
#'
#' @returns an imager::cimg
binarize <- function(img, quantile = 0.95, ...) {
# convert to grayscale if needed and invert
if (rev(dim(img))[1L] > 1L) {
stop("A grayscale image is expected.")
}
# sample
sample <- sample_histogram(img)
# fit Gaussian mixture
gm <- densityMclust(sample, G = 1L, plot = FALSE, ...)
# threshold based on 95% quantile
threshold <- qnorm(
quantile, gm$parameters$mean[1L], sqrt(gm$parameters$variance$sigmasq[1L])
)
out <- img > threshold
return(out)
}
####
# 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
pixset %>%
imager::imrotate(theta * 180L / pi, ehull$loc[1L], ehull$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
####
# 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)
}
####
# 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)
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment