# R/stack_thresh.R In autothresholdr: An R Port of the 'ImageJ' Plugin 'Auto Threshold'

#### Documented in mean_stack_thresh

#' Threshold every image frame in an image stack based on their mean.
#'
#' An [ijtiff_img][ijtiff::ijtiff_img] is a 4-dimensional array indexed by
#' img[y, x, channel, frame]. For each channel (which consists of a stack of
#' frames), this function finds a threshold based on the sum all of the frames,
#' uses this to create a mask and then applies this mask to every frame in the
#' stack (so for a given pillar in the image stack, either all the pixels
#' therein are thresholded away or all are untouched, where pillar x,y of
#' channel ch is img[y, x, ch, ]).
#'
#' It's called mean_stack_thresh() and not sum_stack_thresh() because its
#' easier for people to visualize the mean of an image series than to visualize
#' the sum, but for the sake of this procedure, both are equivalent, except for
#' the fact that the thresholding routine invoked inside this function prefers
#' integers, which we get by using a sum but not by using a mean.
#'
#' \itemize{\item Values greater than or equal to the found threshold
#' \emph{pass} the thresholding and values less than the threshold \emph{fail}
#' the thresholding.
#'
#' \item{For ignore_white = TRUE, if the maximum value in the array is one of
#' 2^8-1, 2^16-1 or 2^32-1, then those max values are ignored.
#' That's because they're the white values in 8, 16 and 32-bit images
#' respectively (and these are the common image bit sizes to work with). This
#' guesswork has to be done because R does not know how many bits the image
#' was on disk. This guess is very unlikely to be wrong, and if it is, the
#' consequences are negligible anyway. If you're very concerned, then just
#' specify the white value as an integer in this ignore_white argument.}
#'
#' \item{If you have set ignore_black = TRUE and/or ignore_white = TRUE but
#' you are still getting error/warning messages telling you to try them, then
#' your chosen method is not working for the given array, so you should try a
#' different method.}
#'
#' \item For a given array, if all values are less than 2^8, saturated value
#' is 2^8 - 1, otherwise, saturated value is 2^16 - 1. }
#'
#' @param img A 4-dimensional array in the style of an
#'   [ijtiff_img][ijtiff::ijtiff_img] (indexed by img[y, x, channel, frame])
#'   or a 3-dimensional array which is a single channel of an
#'   [ijtiff_img][ijtiff::ijtiff_img] (indexed by img[y, x, frame]).
#' @param method The name of the thresholding method you wish to use. The
#'   available methods are "IJDefault", "Huang", "Huang2", "Intermodes",
#'   "IsoData", "Li", "MaxEntropy", "Mean", "MinErrorI", "Minimum",
#'   "Moments", "Otsu", "Percentile", "RenyiEntropy", "Shanbhag",
#'   "Triangle" and "Yen". Partial matching is performed i.e. method = "h"
#'   is enough to get you "Huang" and method = "in" is enough to get you
#'   "Intermodes". To perform \emph{manual} thresholding (where you set the
#'   threshold yourself), supply the threshold here as a number e.g. method =
#'   3.8 (so note that this would \emph{not} select the third method in the
#'   above list of methods). This manual threshold will then be used to
#'   threshold the sum stack to create a 2D mask and then this mask will be
#'   applied to all frames in the stack. If you want a different method for each
#'   channel, specify this parameter as a vector or list, one element per
#'   channel.
#' @param ignore_black Ignore black pixels/elements (zeros) when performing the
#'   thresholding?
#' @param ignore_white Ignore white pixels when performing the thresholding? If
#'   set to TRUE, the function makes a good guess as to what the white
#'   (saturated) value would be (see 'Details'). If this is set to a number, all
#'   pixels with value greater than or equal to that number are ignored.
#' @param fail When using auto_thresh_apply_mask(), to what value do you wish
#'   to set the pixels which fail to exceed the threshold? fail = 'saturate'
#'   sets them to saturated value (see 'Details'). fail = 'zero' sets them to
#'   zero. You can also specify directly here a natural number (must be between
#'   0 and 2^16 - 1) to use.
#' @param ignore_na This should be TRUE if NAs in int_arr should be
#'   ignored or FALSE if you want the presence of NAs in int_arr to throw
#'   an error.
#'
#' @return An object of class [stack_threshed_img] which is the thresholded
#'   image (an array in the style of an [ijtiff_img][ijtiff::ijtiff_img]).
#'   Pillars not exceeding the threshold are set to the fail value (default
#'   NA).
#'
#' @examples
#' \donttest{
#'   package = "autothresholdr"
#' ))
#' ijtiff::display(img[, , 1, 1])
#' ijtiff::display(img[, , 1, 1])
#' }
#' @export
mean_stack_thresh <- function(img, method, fail = NA,
ignore_black = FALSE, ignore_white = FALSE,
ignore_na = FALSE) {
checkmate::assert_array(img, min.d = 3, max.d = 4)
checkmate::assert_numeric(img, lower = 0)
if (all(is.na(img))) {
custom_stop(
"img cannot be all NAs.",
"Every element of your img is NA."
)
}
if (max(img, na.rm = TRUE) > .Machine$integer.max) { custom_stop( "All elements of img must be in the integer range.", " Your maximum element is {format(max(img, na.rm = TRUE), scientific = FALSE)}, which is greater than the maximum allowed value for an int, {format(.Machine$integer.max, scientific = FALSE)}.
"
)
}
if (!isTRUE(all.equal(floor(img), img, check.attributes = FALSE))) {
img[is.na(img)] <- 0
FALSE,
as.vector(img) %>% {
. != floor(.)
}
)
custom_stop(
"img must be an array of integers.",
"
Element {format(bad_index, scientific = FALSE)} of img is
{format(img[bad_index], scientific = FALSE)}, which is not an integer.
"
)
}
if (length(dim(img)) == 3) {
dim(img) %<>% {
c(.[1:2], 1, .[3])
}
}
d <- dim(img)
n_ch <- dim(img)[3]
out <- array(as.vector(img), dim = d)
if (length(method) == 1) method %<>% rep(n_ch)
thresh <- as.list(seq_len(n_ch))
for (i in seq_len(n_ch)) {
if (length(unique(as.vector(img[, , i, ]))) == 1) {
custom_stop(
"Constant arrays cannot bbe thresholded.",
"The array given for thresholding is constant. ",
"All values are equal to {format(img[[1]], scientific = FALSE)}."
)
}
fail <- translate_fail(img, fail)
if (is.numeric(method[[i]])) {
thresh[[i]] <- method[[i]]
} else {
sum_stack <- sum_pillars(structure(img[, , i, ], dim = dim(img)[-3]))
scaling_factor <- 1 # we do this in case the sum stack has big elements
max32int <- 2^31 - 1
mx <- max(sum_stack)
if (mx > max32int) {
scaling_factor <- max32int / mx
sum_stack <- round(sum_stack * scaling_factor)
}
thresh[[i]] <- auto_thresh(sum_stack, method[[i]],
ignore_black = ignore_black,
ignore_white = ignore_white,
ignore_na = ignore_na
) %T>% {
thresh_atts <- attributes(.)
. <- . / (scaling_factor * d[4])
attributes(.) <- thresh_atts
}
}
if (inherits(thresh[[i]], "integer") &&
!isTRUE(checkmate::check_integerish(as.vector(thresh[[i]])))) {
class(thresh[[i]]) %<>% setdiff("integer")
}
mean_stack <- mean_pillars(structure(img[, , i, ], dim = dim(img)[-3]))
out[, , i, ][set_indices] <- fail
}
if (length(thresh) == 1) thresh <- thresh[[1]]
stack_threshed_img(
img = out, thresh = thresh, fail_value = fail,
stack_thresh_method = "mean"
)
}

#' Threshold every image frame in a stack based on their median.
#'
#' An [ijtiff_img][ijtiff::ijtiff_img] is a 4-dimensional array indexed by
#' img[y, x, channel, frame]. For each channel (which consists of a stack of
#' frames), this function finds a threshold based on all of the frames, then
#' takes the median of all the frames in the stack image, uses this to create a
#' mask with the found threshold and then applies this mask to every frame in
#' the stack (so for a given pillar in the image stack, either all the pixels
#' therein are thresholded away or all are untouched, where pillar x,y of
#' channel ch is img[y, x, ch, ]).
#'
#' \itemize{\item Values greater than or equal to the found threshold
#' \emph{pass} the thresholding and values less than the threshold \emph{fail}
#' the thresholding.
#'
#' \item{For ignore_white = TRUE, if the maximum value in the array is one of
#' 2^8-1, 2^16-1 or 2^32-1, then those max values are ignored.
#' That's because they're the white values in 8, 16 and 32-bit images
#' respectively (and these are the common image bit sizes to work with). This
#' guesswork has to be done because R does not know how many bits the image
#' was on disk. This guess is very unlikely to be wrong, and if it is, the
#' consequences are negligible anyway. If you're very concerned, then just
#' specify the white value as an integer in this ignore_white argument.}
#'
#' \item{If you have set ignore_black = TRUE and/or ignore_white = TRUE but
#' you are still getting error/warning messages telling you to try them, then
#' your chosen method is not working for the given array, so you should try a
#' different method.}
#'
#' \item For a given array, if all values are less than 2^8, saturated value
#' is 2^8 - 1, otherwise, saturated value is 2^16 - 1. }
#'
#' @param img A 3-dimensional array (the image stack, possibly a time-series of
#'   images) where the \eqn{n}th slice is the \eqn{n}th image in the stack.
#' @param method The name of the thresholding method you wish to use. The
#'   available methods are "IJDefault", "Huang", "Huang2", "Intermodes",
#'   "IsoData", "Li", "MaxEntropy", "Mean", "MinErrorI", "Minimum",
#'   "Moments", "Otsu", "Percentile", "RenyiEntropy", "Shanbhag",
#'   "Triangle" and "Yen". Partial matching is performed i.e. method = "h"
#'   is enough to get you "Huang" and method = "in" is enough to get you
#'   "Intermodes". To perform \emph{manual} thresholding (where you set the
#'   threshold yourself), supply the threshold here as a number e.g. method =
#'   3 (so note that this would \emph{not} select the third method in the above
#'   list of methods). This manual threshold will then be used to threshold the
#'   median stack to create a 2D mask and then this mask will be applied to all
#'   frames in the stack. If you want a different method for each channel,
#'   specify this parameter as a vector or list, one element per channel.
#' @param ignore_black Ignore black pixels/elements (zeros) when performing the
#'   thresholding?
#' @param ignore_white Ignore white pixels when performing the thresholding? If
#'   set to TRUE, the function makes a good guess as to what the white
#'   (saturated) value would be (see 'Details').
#' @param fail When using auto_thresh_apply_mask(), to what value do you wish
#'   to set the pixels which fail to exceed the threshold? fail = 'saturate'
#'   sets them to saturated value (see 'Details'). fail = 'zero' sets them to
#'   zero. You can also specify directly here a natural number (must be between
#'   0 and 2^32 - 1) to use.
#' @param ignore_na This should be TRUE if NAs in int_arr should be
#'   ignored or FALSE if you want the presence of NAs in int_arr to throw
#'   an error.
#'
#' @return An object of class [stack_threshed_img] which is the thresholded
#'   image (an array in the style of an [ijtiff_img][ijtiff::ijtiff_img]).
#'   Pillars not exceeding the threshold are set to the fail value (default
#'   NA).
#'
#' @examples
#' \donttest{
#'   package = "autothresholdr"
#' ))
#' ijtiff::display(img[, , 1, 1])
#' ijtiff::display(img[, , 1, 1])
#' }
#' @export
med_stack_thresh <- function(img, method, fail = NA,
ignore_black = FALSE, ignore_white = FALSE,
ignore_na = FALSE) {
checkmate::assert_array(img, min.d = 3, max.d = 4)
checkmate::assert_numeric(img, lower = 0)
if (all(is.na(img))) {
custom_stop(
"img cannot be all NAs.",
"Every element of your img is NA."
)
}
if (max(img, na.rm = TRUE) > .Machine$integer.max) { custom_stop( "All elements of img must be in the integer range.", " Your maximum element is {format(max(img, na.rm = TRUE), scientific = FALSE)}, which is greater than the maximum allowed value for an int, {format(.Machine$integer.max, scientific = FALSE)}.
"
)
}
if (!isTRUE(all.equal(floor(img), img, check.attributes = FALSE))) {
img[is.na(img)] <- 0
FALSE,
as.vector(img) %>% {
. != floor(.)
}
)
custom_stop(
"img must be an array of integers.",
"
Element {format(bad_index, scientific = FALSE)} of img is
{format(img[bad_index], scientific = FALSE)}, which is not an integer.
"
)
}
if (length(dim(img)) == 3) {
dim(img) %<>% {
c(.[1:2], 1, .[3])
}
}
d <- dim(img)
n_ch <- dim(img)[3]
out <- array(as.vector(img), dim = d)
if (length(method) == 1) method %<>% rep(n_ch)
thresh <- as.list(seq_len(n_ch))
for (i in seq_len(n_ch)) {
if (length(unique(as.vector(img[, , i, ]))) == 1) {
custom_stop(
"Constant arrays cannot bbe thresholded.",
"The array given for thresholding is constant. ",
"All values are equal to {format(img[[1]], scientific = FALSE)}."
)
}
fail <- translate_fail(img, fail)
thresh[[i]] <- auto_thresh(img[, , i, ], method[[i]],
ignore_black = ignore_black,
ignore_white = ignore_white,
ignore_na = ignore_na
)
med_stack <- median_pillars(structure(img[, , i, ], dim = dim(img)[-3]))
out[, , i, ][set_indices] <- fail
}
stack_threshed_img(
img = out, thresh = thresh, fail_value = fail,
stack_thresh_method = "median"
)
}


## Try the autothresholdr package in your browser

Any scripts or data that you put into this service are public.

autothresholdr documentation built on Jan. 4, 2022, 5:11 p.m.