#' Extracts color pattern from landmark transformed image using watershed segmentation. This function
#' works interactively by allowing to pick a starting pixel within each pattern element from which the
#' watershed will extract the pattern. This function works best for patterns with sharp boundaries.
#'
#' @param sampleList List of RasterStack objects.
#' @param landList Landmark list as returned by \code{\link[patternize]{makeList}}.
#' @param IDlist List of sample IDs should be specified when masking outline and transformRef
#' is 'meanshape'.
#' @param adjustCoords Adjust landmark coordinates in case they are reversed compared to pixel
#' coordinates (default = FALSE).
#' @param transformRef ID or landmark matrix of reference sample for shape to which color patterns
#' will be transformed to. Can be 'meanshape' for transforming to mean shape of Procrustes
#' analysis.
#' @param resampleFactor Integer for downsampling image used by \code{\link{redRes}}.
#' @param transformType Transformation type as used by \code{\link[Morpho]{computeTransform}}
#' (default ='tps').
#' @param maskOutline When outline is specified, everything outside of the outline will be masked for
#' the color extraction (default = NULL).
#' @param cartoonID ID of the sample for which the cartoon was drawn and will be used for masking
#' (should be set when transformRef = 'meanShape').
#' @param correct Correct image illumination using a linear model (default = FALSE).
#' @param blur Blur image for priority map extraction (default = TRUE).
#' @param sigma Size of sigma for Gaussian blurring (default = 5).
#' @param bucketfill Use a bucket fill on the background to fill holes (default = TRUE).
#' @param cleanP Integer to remove spurious areas with width smaller than cleanP (default = NULL).
#' @param splitC Integer to split selected patterns into connected components and remove ones with
#' areas smaller than splitC (default = NULL).
#' @param plotTransformed Plot transformed image (default = FALSE).
#' @param plotCorrect Plot corrected image, corrected for illumination using a linear model
#' (default = FALSE).
#' @param plotEdges Plot image gradient (default = FALSE).
#' @param plotPriority Plot priority map (default = FALSE).
#' @param plotWS Plot watershed result (default = FALSE).
#' @param plotBF Plot bucketfill (default = FALSE).
#' @param plotFinal Plot extracted patterns (default = FALSE).
#'
#' @return List of raster objects.
#'
#' @examples
#'
#' \dontrun{
#' IDlist <- c('BC0077','BC0071','BC0050','BC0049','BC0004')
#' prepath <- system.file("extdata", package = 'patternize')
#' extension <- '_landmarks_LFW.txt'
#'
#' landmarkList <- makeList(IDlist, 'landmark', prepath, extension)
#'
#' extension <- '.jpg'
#' imageList <- makeList(IDlist, 'image', prepath, extension)
#'
#' outline_BC0077 <- read.table(paste(system.file("extdata", package = 'patternize'),
#' '/BC0077_outline.txt', sep=''), header = FALSE)
#'
#' rasterList_W <- patLanW(imageList, landmarkList, IDlist, transformRef = 'meanshape',
#' adjustCoords = TRUE, plotTransformed = FALSE, correct = TRUE, plotCorrect = FALSE, blur = FALSE,
#' sigma = 2, bucketfill = FALSE, cleanP = 0, splitC = 10, plotPriority = TRUE, plotWS = TRUE,
#' plotBF = TRUE, plotFinal = TRUE, maskOutline = outline_BC0077, cartoonID = 'BC0077')
#' }
#'
#' @export
#' @import raster
#' @importFrom imager as.cimg imfill isoblur imgradient imsplit watershed bucketfill clean split_connected highlight enorm parany add
#' @importFrom purrr discard
#' @importFrom dplyr sample_n
#' @importFrom magrittr %>%
#' @importFrom graphics locator
#' @importFrom stats lm
#' @importFrom Morpho procSym computeTransform applyTransform
patLanW <- function(sampleList,
landList,
IDlist = NULL,
adjustCoords = FALSE,
transformRef = 'meanshape',
resampleFactor = NULL,
transformType = 'tps',
maskOutline = NULL,
cartoonID = NULL,
correct = FALSE,
blur = TRUE,
sigma = 3,
bucketfill = TRUE,
cleanP = NULL,
splitC = NULL,
plotTransformed = FALSE,
plotCorrect = FALSE,
plotEdges = FALSE,
plotPriority = FALSE,
plotWS = FALSE,
plotBF = FALSE,
plotFinal = FALSE){
declare( . <- as_nonstandard())
rasterList <- list()
# Check whether sampleList and landList have the same length
if(length(sampleList) != length(landList)){
stop("sampleList is not of the same length as lanArray")
}
for(n in 1:length(sampleList)){
if(names(sampleList)[n] != names(landList)[n]){
stop("samples are not in the same order in sampleList and lanArray")
}
}
# Make landmark array
lanArray <- lanArray(landList, adjustCoords, sampleList)
# Set the reference shape
if(is.matrix(transformRef)){
refShape <- transformRef
}
if(!is.matrix(transformRef)){
if(transformRef == 'meanshape'){
invisible(capture.output(transformed <- Morpho::procSym(lanArray)))
refShape <- transformed$mshape
}
if(transformRef %in% names(landList)){
e <- which(names(landList) == transformRef)
refShape <- lanArray[,,e]
}
}
# Transform the outline for masking if 'meanShape'
if(!is.null(cartoonID)){
if(any(c(!is.null(maskOutline), transformRef == 'meanshape'))){
indx <- which(names(sampleList) == cartoonID)
maskOutlineNew <- maskOutline
extPicture <- raster::extent(sampleList[[indx]])
maskOutlineNew[,2] <- extPicture[4]-maskOutlineNew[,2]
}
}
if(is.null(cartoonID)){
maskOutlineNew <- maskOutline
}
if(all(c(!is.null(maskOutline), transformRef == 'meanshape'))){
invisible(capture.output(cartoonLandTrans <- Morpho::computeTransform(refShape,
as.matrix(lanArray[,,indx]),
type="tps")))
maskOutlineMean <- Morpho::applyTransform(as.matrix(maskOutlineNew), cartoonLandTrans)
}
n <- 1
todo <- 0
# Run the loop for each sample
while(n <= length(sampleList)){
if(todo != 'r' && todo != 'a'){
image <- sampleList[[n]]
extRaster <- raster::extent(image)
# Reduce resolution
if(!is.null(resampleFactor)){
image <- redRes(image, resampleFactor)
}
# Transform image using landmarks
invisible(capture.output(transMatrix <- Morpho::computeTransform(refShape,
as.matrix(lanArray[,,n]),
type = transformType)))
imageDF1 <- raster::as.data.frame(image[[1]], xy = TRUE)
imageDF2 <- raster::as.data.frame(image[[2]], xy = TRUE)
imageDF3 <- raster::as.data.frame(image[[3]], xy = TRUE)
invisible(capture.output(imageT <- Morpho::applyTransform(as.matrix(imageDF1)[,1:2], transMatrix)))
r <- raster::raster(nrow = dim(image)[1], ncol = dim(image)[2])
raster::extent(r) <- c(min(imageT[,1]),max(imageT[,1]),min(imageT[,2]),max(imageT[,2]))
# Rasterize the transformed image and fill in NA values using
imageT1r <- raster::rasterize(imageT, field = imageDF1[,3], r, fun = mean)
imageT1rf <- focal(imageT1r, w=matrix(1,nrow=3,ncol=3), fun=fill.na, pad = TRUE, na.rm = FALSE)
imageT2r <- raster::rasterize(imageT, field = imageDF2[,3], r, fun = mean)
imageT2rf <- focal(imageT2r, w=matrix(1,nrow=3,ncol=3), fun=fill.na, pad = TRUE, na.rm = FALSE)
imageT3r <- raster::rasterize(imageT, field = imageDF3[,3], r, fun = mean)
imageT3rf <- focal(imageT3r, w=matrix(1,nrow=3,ncol=3), fun=fill.na, pad = TRUE, na.rm = FALSE)
imageTr <- raster::stack(imageT1rf, imageT2rf, imageT3rf)
imageTr[is.na(imageTr)] <- 255
# Plot transformed raster
if(plotTransformed){
# imageTr <- raster::flip(imageTr, 'y')
x <- as.array(imageTr)/255
cols <- rgb(x[,,1], x[,,2], x[,,3], maxColorValue=1)
uniqueCols <- unique(cols)
x2 <- match(cols, uniqueCols)
dim(x2) <- dim(x)[1:2]
raster::image(t(x2), col=uniqueCols, yaxt='n', xaxt='n', main = paste(names(landList)[n],'transformed', sep=' '))
# imageTr <- raster::flip(imageTr, 'y')
}
# MaskOutline
if(!is.null(maskOutline)){
if(transformRef[1] != 'meanshape'){
imageTr <- maskOutline(imageTr, maskOutlineNew, refShape = 'target', crop = c(0,0,0,0),
maskColor = 255, imageList = sampleList)
cropEx <- c(min(maskOutlineNew[,1]), max(maskOutlineNew[,1]), min(maskOutlineNew[,2]), max(maskOutlineNew[,2]))
}
if(transformRef[1] == 'meanshape'){
imageTr <- maskOutline(imageTr, outline = maskOutline, refShape = 'mean',
IDlist = IDlist, landList = landList, imageList = sampleList,
adjustCoords = TRUE, cartoonID = cartoonID)
cropEx <- c(min(maskOutlineMean[,1]), max(maskOutlineMean[,1]), min(maskOutlineMean[,2]), max(maskOutlineMean[,2]))
}
imageTr <- raster::crop(imageTr, cropEx)
raster::extent(imageTr) <- cropEx
}
# Transform to imager format
imA <- raster::as.array(imageTr)
imR <- as.raster(imA, nrow = dim(imageTr)[1], ncol = dim(imageTr)[2], max = 255)
im <- imager::as.cimg(imR)
if(n == 1){
im1 <- imager::as.cimg(imR)
}
# Correct illumination using linear model
if(correct){
if(plotCorrect){
layout(t(1:2))
par(mar=c(2, 2, 2, 2))
plot(im, main = paste(names(landList)[n], 'original', sep=' '))
}
d <- as.data.frame(im)
# Subsamble, fit a linear model
m <- sample_n(d, 1e4) %>% lm(value ~ x*y, data=.)
# Correct by removing the trend
im <- im-predict(m,d)
if(plotCorrect) plot(im, main = paste(names(landList)[n], 'corrected', sep=' '))
}
# # Blur image for color extraction
# if(focal){
#
# gf <- focalWeight(image, sigma, "Gauss")
#
# rrr1 <- raster::focal(image[[1]], gf)
# rrr2 <- raster::focal(image[[2]], gf)
# rrr3 <- raster::focal(image[[3]], gf)
#
# image <- raster::stack(rrr1, rrr2, rrr3)
#
# image[is.na(image)] <- 255
# }
# Blur image for priority map extraction
if(blur){
edges <- isoblur(im, sigma) %>% imgradient("xy") %>% enorm %>% imsplit("c") %>% add %>% sqrt
}
else{
edges <- imgradient(im,"xy") %>% enorm
}
layout(t(1:1))
if(plotEdges) plot(edges)
priority <- 1/(255+edges)
if(plotPriority) plot(priority)
# Extract seed points from image
seed <- imager::imfill(dim=dim(im)) #Empty image
}
layout(t(1:1))
par(mar=c(2, 2, 2, 2))
plot(im)
print("Choose points to identify patterns. Click outside image area to stop.")
while(1){
xy <- locator(n=1)
if(any(c(as.numeric(xy)[1] > dim(im)[1], as.numeric(xy)[1] < 0, as.numeric(xy)[2] > dim(im)[2], as.numeric(xy)[2] < 0))){
break
}
print(paste('x: ', as.character(xy)[1], 'y: ', as.character(xy)[2]))
seed[as.numeric(xy[1]),as.numeric(xy[2]),1,1] <- 2
}
print("Choose at least one point to identify background. Click outside image area to stop.")
while(1){
xy <- locator(n=1)
if(any(c(as.numeric(xy)[1] > dim(im)[1], as.numeric(xy)[1] < 0, as.numeric(xy)[2] > dim(im)[2], as.numeric(xy)[2] < 0))){
break
}
print(paste('x: ', as.character(xy)[1], 'y: ', as.character(xy)[2]))
xyBF <- xy
seed[as.numeric(xy[1]),as.numeric(xy[2]),1,1] <- 1
}
# Run watershed transform
ws <- (imager::watershed(seed, priority) == 2) # multiply by 255 if RGB extracted
if(plotWS) plot(ws)
# Use a bucket fill on the background to fill holes (starts from last selected bg pixel)
if(bucketfill){
ws <- imager::bucketfill(ws, as.numeric(xyBF[1]), as.numeric(xyBF[2]), color=2) %>% {!( . == 2)}
}
if(plotBF) plot(ws)
# Remove spurious areas with width < cleanP
if(!is.null(cleanP)){
ws <- imager::clean(ws,cleanP)
}
# Split into connected components and remove ones with areas < splitC
if(!is.null(splitC)){
ws <- imager::split_connected(ws) %>% purrr::discard(~ sum(.) < splitC) %>% parany
}
# Transform to raster
wsR <- as.raster(ws)
wsM <- matrix(t(col2rgb(wsR))[,1]/255, nrow = dim(im)[1], ncol = dim(im)[2])
patternRaster <- raster::raster(t(wsM))
raster::extent(patternRaster) <- raster::extent(0, dim(im)[1], 0, dim(im)[2])
if(plotFinal){
layout(t(1:2))
plot(im, main="Watershed")
highlight(ws, col="blue")
plot(patternRaster, legend = FALSE, bty = 'n', main = 'Extracted pattern')
}
todo <- readline(prompt="Press [enter] to continue, 'a' to add points 'r' to redo or 'q' to quit >>> ")
if(todo == 'q') break
if(todo == 'a') next
if(todo == 'r'){
seed <- imager::imfill(dim=dim(im))
next
}
else{
y <- raster::raster(ncol = dim(im1)[2], nrow = dim(im1)[1])
extent(y) <- raster::extent(0, dim(im)[1], 0, dim(im)[2])
patternRaster <- resample(patternRaster, y, method = 'ngb')
raster::extent(patternRaster) <- c(min(refShape[,1]), max(refShape[,1]), min(refShape[,2]), max((refShape[,2])))
rasterList[[names(landList)[n]]] <- patternRaster
print(paste('sample', names(landList)[n], 'done and added to rasterList', sep=' '))
n = n + 1
}
}
return(rasterList)
}
# Function used in focal to calculate average for NA values in tranformed matrix
fill.na <- function(x, i=5) {
if( is.na(x)[i] ) {
return( round(mean(x, na.rm=TRUE),0) )
} else {
return( round(x[i],0) )
}
}
###
# The following code is added to suppress raising a R CMD check note resulting from magrittr '.'
###
# Placeholder function.
as_nonstandard <- function()
{
stop("This is a placeholder function, and should not be called directly.")
}
# Declare symbols for nonstandard evalution.
declare <- function(declaration)
{
d <- substitute(declaration)
p <- parent.frame()
if (!identical(d[[1]], quote(`<-`)))
stop("Declarations follow the <- syntax.", call. = FALSE)
if (!is.symbol(d[[2]]))
stop("LHS needs to be a symbol", call. = FALSE)
sym <- as.character(d[[2]])
msg <- sprintf("Symbol '%s' is declared for non-standard use.", sym)
eval(
substitute(
delayedAssign(x, stop(m, call. = FALSE), environment(), environment()),
list(m = msg, x = sym)), p, p)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.