#' @title Neighborhood Correlation Image
#' @description Calculates the neighborhood of 2 registered images, and then
#' does the correlation of each neighbor with itself in the other imaging modality
#'
#' @param img1 image of class \code{nifti}, \code{antsImage}, \code{character}
#' @param img2 image of class \code{nifti}, \code{antsImage}, \code{character}
#' @param mask Binary image of class \code{nifti},
#' \code{antsImage}, \code{character}. Only neighborhoods inside the mask
#' will be taken
#' @param radius vector of length 3 for number of voxels to go in each direction.
#' Default is 27 neighbors (including voxel at center).
#' Passed to \code{\link{getNeighborhoodInMask}}.
#' @param method Type of correlation.
#' @param verbose print diagnostic messages
#'
#' @return Object of class \code{nifti}
#' @export
#'
#' @examples \dontrun{
#' corr_img("T1_Image.nii.gz", "FLAIR_Image.nii.gz", mask = "Mask.nii.gz")
#' }
#' @importFrom matrixStats colCounts
corr_img = function(
img1,
img2,
mask = NULL,
radius = rep(1,3),
method = c("pearson", "spearman"),
verbose = TRUE){
method = match.arg(method)
img1 = check_ants(img1)
img2 = check_ants(img2)
if (is.null(mask)) {
mask = array(1, dim = dim(img1))
mask = as.antsImage(mask)
mask = antsCopyImageInfo(
target = mask,
reference = img1)
}
xmask = mask
mask = check_ants(mask)
if (verbose) {
message("Checking dimensions of images")
}
stopifnot(
same_dims(img1, img2, mask)
)
# col_scale = function(x){
# cm = colMeans(x, na.rm = TRUE)
# csd = colSds(x, center = cm,
# na.rm = TRUE)
# x = t( (t(x) - cm) / csd )
# }
# neigh = function(img, mask, radius, method){
# grads = getNeighborhoodInMask(
# image = img,
# mask = mask,
# radius = radius,
# spatial.info = TRUE)
# # not_zero = grads1$offsets != 0
# # neighbor = rowSums(not_zero) > 0
# if (method == "spearman") {
# grads$values = colRanks(grads$values,
# ties.method = "average",
# preserveShape = TRUE)
# }
# grads$values = col_scale(
# grads$values)
# return(grads)
# }
if (verbose) {
message("Getting Neighborhood for Image 1")
}
neigh1 = scaled_neighborhood(img = img1,
mask = mask,
radius = radius,
method = method,
center = TRUE,
scale = TRUE,
verbose = verbose)
# neigh1 = neigh(img = img1,
# mask = mask,
# radius = radius,
# method = method)
if (verbose) {
message("Getting Neighborhood for Image 2")
}
neigh2 = scaled_neighborhood(img = img2,
mask = mask,
radius = radius,
method = method,
center = TRUE,
scale = TRUE,
verbose = verbose)
# neigh2 = neigh(img = img2,
# mask = mask,
# radius = radius,
# method = method)
stopifnot(
identical(neigh1$indices,
neigh2$indices)
)
if (verbose) {
message("Calculating Correlations")
}
corrs = neigh1$values * neigh2$values
nNA <- colCounts(corrs,
value = NA_real_,
na.rm = FALSE)
n = nrow(corrs) - nNA
corrs = colSums(corrs, na.rm = TRUE)
corrs = corrs / (n - 1)
# divide by n - 1
# corrs = colMeans(corrs, na.rm = TRUE)
# #####################
# # Must ADD 1 because they are 0-indexed!
# #####################
# if (verbose) {
# message("Reordering indices from ANTsR output")
# }
# out = reorder_neighborhood(neigh1,
# img.dim = dim(img1))
# ord = out$order_ind
# corrs = corrs[ ord ]
#
if (verbose) {
message("Creating Output Image")
}
xmask = check_nifti(xmask)
ximg = xmask
ximg[xmask == 1] = corrs
ximg = cal_img(ximg)
ximg = datatyper(ximg,
datatype =
convert.datatype()$FLOAT32,
bitpix =
convert.bitpix()$FLOAT32)
rm(list = c("img1", "img2", "mask", "xmask"))
rm(list = c("neigh1", "neigh2"))
for (i in 1:10) {
gc()
}
return(ximg)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.