############################################################ # cleans and thresholds image clean <- function(x, thresh_win = 30, thresh_sense = 0.01, schmutz = 81, element_diam = 9, edgefill = TRUE, verbose = FALSE){ require(EBImage) if(!is.Image(x)) stop('x is not an image') col <- x # convert to greyscale and normalize values gs <- normalize(channel(col, mode = 'gray')) # adaptive thresholding; # offset is 'sensitivity'; w and h give window size bin <- thresh(gs, w = thresh_win, h = thresh_win, offset = thresh_sense) # Note: # white = 1 (background, areoles) # black = 0 (foreground, veins) # opening and closing # kernel <- makeBrush(element_diam, shap = 'disc') # bin <- opening(bin, kern = kernel) # bin <- closing(bin, kern = kernel) # remove unconnected black in white invert <- !bin # invert image so vein network is white/foreground invert.lab <- bwlabel(invert) # lable inverted image tab <- tabulate(invert.lab) biggest.index <- (1:length(tab))[tab == max(tab)] biggest <- invert.lab == biggest.index clean <- as.Image(!biggest) # invert image so veins are again black/background storage.mode(clean) <- 'numeric' if(verbose) cat('removing all background pixels (white) unconnected to vein network in foreground (black)\n') if(edgefill){ # remove edge effects SLOW! for(i in seq(from = 1, to = nrow(clean), by = sqrt(schmutz))){ clean <- floodFill(clean, pt = c(i,1), col = 0) clean <- floodFill(clean, pt = c(i,ncol(clean)), col = 0) } for(j in seq(from = 1, to = ncol(clean), by = sqrt(schmutz))){ clean <- floodFill(clean, pt = c(1,j), col = 0) clean <- floodFill(clean, pt = c(nrow(clean),j), col = 0) } } #remove small white in black clean.lab <- bwlabel(clean) tab <- tabulate(clean.lab) fts <- computeFeatures.shape(clean.lab) small_ids <- (1:max(clean.lab))[fts[,1] < schmutz] clean[clean.lab %in% small_ids] <- 0 return (clean) } # End of function ############################################################ # sizing transform st <- function(x, max_mask = max(dim(x)), draw = FALSE, negate = FALSE){ # x is an image object # max_mask is an integer giving an upper bound for mask size # draw is FALSE for no image output; TRUE for screen output; # a string to produce an output file # negate = TRUE reverses foreground and background require(EBImage) if(!is.Image(x)){ stop('x is not an image') } if(length(table(x)) > 2){ warning('input image does not seem to be binary') } if(negate){ x <- !x } # initialize variables to hold the sequentially larger masks masks <- list() masks.im <- Image(x) # initialize index i <- 1 # initialize working mask thismask <- x # initialize size transform matrix sizes <- matrix(0, ncol = ncol(x), nrow = nrow(x)) # loop through sequentially larger odd circular kernals # until all foreground is eliminated while(sum(thismask > 0) > 0){ # while there are foreground pixels left.... thisdisc <- 2 * i + 1 #cat(thisdisc) #cat(max_mask) if(thisdisc > max_mask) break() cat(paste(thisdisc, '\n')) thismask <- opening(x, kern = makeBrush(thisdisc, shape = 'disc')) # perform an opening if(draw == TRUE) image(thismask) masks.im <- combine(masks.im, thismask) # save this mask as an image masks[[i]] <- apply(thismask, 2, function(x){as.numeric(as.logical(x))}) # and as a matrix names(masks)[i] <- paste('d', thisdisc, sep = '') i <- i + 1 sizes[apply(thismask > 0, 2, as.logical)] <- thisdisc # where thismask still has foreground pixels, put in the size of the current disk } if(draw == TRUE){ par(mfrow = c(1,2)) image(sizes) plot(as.table(table(sizes)[-1]), type = 'h', ylab = 'Frequency', xlab = 'Diameter in Pixels', main = 'Table of Image Sizing Transform') }else if(is.character(draw)){ nio <- dim(masks.im)[3] writeImage(masks.im, files = paste(draw, 'mask', 1:nio, '.jpg', sep = '')) pdf(paste(draw, 'Rplot.pdf', sep = '')) plot(as.table(table(sizes)[-1]), type = 'h', ylab = 'Frequency', xlab = 'Diameter in Pixels', main = 'Table of Image Sizing Transform') dev.off() } return(list(kernels = seq(3, i, by = 2), sizes = sizes, masks = masks, masks.im = masks.im)) } # End of function ############################################################ # a function to perform a morphological thinning on an image # algorithm from Glasbey and Horgan, Image Analysis for the # Biological Sciences; Chapt. 5, p.15f thin <- function(x, verbose = FALSE){ dims <- dim(x) #if(is.Image(x)) # x <- x@.Data[,,1] # edge effects are dealt with by cloning the first and last # rows and columns to produce a matrix two pixels wider than # the input matrix and then truncating the output matrix (by # one pixel all around) to return to the initial size. x <- rbind(x[1,], x, x[nrow(x),]) x <- cbind(x[,1], x, x[,ncol(x)]) for(i in 2:dims[1]){ # loop through each pixel in image for(j in 2:dims[2]){ # comparing with shape elements 1--8 if(verbose) cat(i, j, ':\t', sep = '\t') # se1 if(x[i,j] == 0 && (x[i-1,j-1] == 0 && length(x[i-1,j-1])) && (x[i,j-1] == 0 && length(x[i,j-1])) && (x[i+1,j-1] == 0 && length(x[i,j])) && (x[i-1,j+1] == 1 && length(x[i-1,j+1])) && (x[i,j+1] == 1 && length(x[i,j+1])) && (x[i+1,j+1] == 1 && length(x[i+1,j+1]))){ x[i,j] <- 1 if(verbose) cat('1') }else{ if(verbose) cat('.') } # se2 if(x[i,j] == 0 && (x[i,j-1] == 0 && length(x[i,j-1])) && (x[i+1,j-1] == 0 && length(x[i+1,j-1])) && (x[i+1,j] == 0 && length(x[i+1,j])) && (x[i-1,j+1] == 1 && length(x[i-1,j+1])) && (x[i,j+1] == 1 && length(x[i,j+1])) && (x[i-1,j] == 1 && length(x[i-1,j]))){ x[i,j] <- 1 if(verbose) cat('2') }else{ if(verbose) cat('.') } # se3 if(x[i,j] == 0 && (x[i+1,j-1] == 0 && length(x[i+1,j-1])) && (x[i+1,j] == 0 && length(x[i+1,j])) && (x[i+1,j+1] == 0 && length(x[i+1,j+1])) && (x[i-1,j-1] == 1 && length(x[i-1,j-1])) && (x[i-1,j] == 1 && length(x[i-1,j])) && (x[i-1,j+1] == 1 && length(x[i-1,j+1]))){ x[i,j] <- 1 if(verbose) cat('3') }else{ if(verbose) cat('.') } # se4 if(x[i,j] == 0 && (x[i+1,j] == 0 && length(x[i+1,j])) && (x[i+1,j+1] == 0 && length(x[i+1,j+1])) && (x[i,j+1] == 0 && length(x[i,j+1])) && (x[i-1,j] == 1 && length(x[i-1,j])) && (x[i-1,j-1] == 1 && length(x[i-1,j-1])) && (x[i,j-1] == 1 && length(x[i,j-1]))){ x[i,j] <- 1 if(verbose) cat('4') }else{ if(verbose) cat('.') } # se5 if(x[i,j] == 0 && (x[i-1,j+1] == 0 && length(x[i-1,j+1])) && (x[i,j+1] == 0 && length(x[i,j+1])) && (x[i+1,j+1] == 0 && length(x[i+1,j+1])) && (x[i-1,j-1] == 1 && length(x[i-1,j-1])) && (x[i,j-1] == 1 && length(x[i,j-1])) && (x[i+1,j-1] == 1 && length(x[i+1,j-1]))){ x[i,j] <- 1 if(verbose) cat('5') }else{ if(verbose) cat('.') } # se6 if(x[i,j] == 0 && (x[i-1,j] == 0 && length(x[i-1,j])) && (x[i-1,j+1] == 0 && length(x[i-1,j+1])) && (x[i,j+1] == 0 && length(x[i,j+1])) && (x[i+1,j] == 1 && length(x[i+1,j])) && (x[i+1,j-1] == 1 && length(x[i+1,j-1])) && (x[i,j-1] == 1 && length(x[i,j-1]))){ x[i,j] <- 1 if(verbose) cat('6') }else{ if(verbose) cat('.') } # se7 if(x[i,j] == 0 && (x[i-1,j-1] == 0 && length(x[i-1,j-1])) && (x[i-1,j] == 0 && length(x[i-1,j])) && (x[i-1,j+1] == 0 && length(x[i-1,j+1])) && (x[i+1,j-1] == 1 && length(x[i+1,j-1])) && (x[i+1,j] == 1 && length(x[i+1,j])) && (x[i+1,j+1] == 1 && length(x[i+1,j+1]))){ x[i,j] <- 1 if(verbose) cat('7') }else{ if(verbose) cat('.') } # se8 if(x[i,j] == 0 && (x[i-1,j] == 0 && length(x[i-1,j])) && (x[i-1,j-1] == 0 && length(x[i-1,j-1])) && (x[i,j-1] == 0 && length(x[i,j-1])) && (x[i+1,j] == 1 && length(x[i+1,j])) && (x[i+1,j+1] == 1 && length(x[i+1,j+1])) && (x[i,j+1] == 1 && length(x[i,j+1]))){ x[i,j] <- 1 if(verbose) cat('8\n') }else{ if(verbose) cat('.\n') } } } x <- x[-nrow(x),-ncol(x)] x <- x[-1,-1] thinned <- Image(x) #thinned@compression <- "NONE" return(thinned) } # End of function ############################################################ # thinning to a 1 pixel thick skeleton skeleton <- function(x, draw = TRUE, verbose = FALSE){ prev.iter <- x skel <- thin(prev.iter) if(draw) image(skel) ii <- 1 while(any(skel != prev.iter)){ cat('Iteration', ii, '\n') prev.iter <- skel skel <- thin(skel) ii <- ii + 1 if(draw) image(skel) } return(skel) } # End of function ############################################################ # autocrossproduct function for 4 angles (0, pi/4, pi/2, 3*pi/4) axp <- function(x, maxkern = min(nrow(x), ncol(x))/2 - 1, draw = TRUE){ require(EBImage) kerns <- seq(3, maxkern, by = 2) xpixct <- sum(x > 0) vert <- rep(0, length(kerns)) horiz <- rep(0, length(kerns)) nesw <- rep(0, length(kerns)) # NE-SW diagonal nwse <- rep(0, length(kerns)) # NW-SE diagonal for(i in kerns){ if(draw) par(mfrow = c(2,2)) ekern <- matrix(as.integer(0), nrow = i, ncol = i) ekern[1,1] <- as.integer(1) ekern[i,i] <- as.integer(1) nwse[kerns == i] <- sum(erode(x, kern = ekern) > 0) if(draw){ image(erode(x, kern = ekern), main = paste( 'NW-SE diagonal, kernel size = ', i, 'pixels')) } ekern <- matrix(as.integer(0), nrow = i, ncol = i) ekern[1,i] <- as.integer(1) ekern[i,1] <- as.integer(1) nesw[kerns == i] <- sum(erode(x, kern = ekern) > 0) if(draw){ image(erode(x, kern = ekern), main = paste( 'NE-SW diagonal, kernel size = ', i, 'pixels')) } ekern <- matrix(as.integer(0), nrow = i, ncol = i) ekern[1,(i+1)/2] <- as.integer(1) ekern[i,(i+1)/2] <- as.integer(1) vert[kerns == i] <- sum(erode(x, kern = ekern) > 0) if(draw){ image(erode(x, kern = ekern), main = paste( 'Vertical, kernel size = ', i, 'pixels')) } ekern <- matrix(as.integer(0), nrow = i, ncol = i) ekern[(i+1)/2,1] <- as.integer(1) ekern[(i+1)/2,i] <- as.integer(1) horiz[kerns == i] <- sum(erode(x, kern = ekern) > 0) if(draw){ image(erode(x, kern = ekern), main = paste( 'Horizontal, kernel size = ', i, 'pixels')) } } axps <- list(kerns = kerns, vert = vert, horiz = horiz, nesw = nesw, nwse = nwse) if (draw){ plot.axp(axps) } return(axps) } # End of function plot.axp <- function(x){ kerns <- x$kerns axps <- x par(mfrow = c(1,2)) yrng <- range(unlist(axps[-1])) plot(kerns, axps$vert, xlab = 'kernel length in pixels', ylab = 'autocrossproduct (pixel count)', ylim = yrng) text(kerns[1], axps$vert[1], 'vert', adj = c(-0.2, 0.5)) points(kerns, axps$horiz, col = 'red') text(kerns[1], axps$horiz[1], 'horiz', adj = c(-0.2, 0.5), col = 'red') points(kerns, axps$nesw, col = 'green') text(kerns[1], axps$nesw[1], 'nesw', adj = c(-0.2, 0.5), col = 'green') points(kerns, axps$nwse, col = 'blue') text(kerns[1], axps$nwse[1], 'nwse', adj = c(-0.2, 0.5), col = 'blue') lengths = vector(mode = 'list', length = 4) names(lengths) <- names(axps)[-1] for(i in 2:5){ lengths[i-1] <- axps$kerns[axps[[i]]==min(axps[[i]])] } rg <- max(unlist(lengths)) plot(0, xlim = c(-rg,rg), ylim = c(-rg,rg), type = 'n') segments(0, 0, 0, lengths$vert) segments(0, 0, 0, -lengths$vert) segments(0, 0, lengths$horiz, 0, col = 'red') segments(0, 0, -lengths$horiz, 0, col = 'red') segments(0, 0, lengths$nesw * cos(pi/4), lengths$nesw * cos(pi/4), col = 'green') segments(0, 0, -lengths$nesw * cos(pi/4), -lengths$nesw * cos(pi/4), col = 'green') segments(0, 0, lengths$nwse * cos(pi/4), -lengths$nwse * cos(pi/4), col = 'blue') segments(0, 0, -lengths$nwse * cos(pi/4), lengths$nwse * cos(pi/4), col = 'blue') } # End of function # 08-07-08 working on the axp() (autocrossproduct) function #im.syn <- readImage(paste(PATH, 'synthetic.binary.image.tif', # sep = '')) # #im.h <- readImage(paste(PATH, 'h.tif', sep = '')) #im.v <- readImage(paste(PATH, 'v.tif', sep = '')) #im.b <- readImage(paste(PATH, 'b.tif', sep = '')) #im.sp <- readImage(paste(PATH, 'sp.tif', sep = '')) #im.h.axp <- axp(im.h, draw = FALSE, maxkern = 30) #im.v.axp <- axp(im.v, draw = FALSE, maxkern = 30) #im.b.axp <- axp(im.b, draw = FALSE, maxkern = 30) #im.sp.axp <- axp(im.sp, draw = FALSE, maxkern = 30) ############################################################ # 08-07-22 # working on ferets() ferets <- function(x){ } #im.mask.ws <- watershed(distmap(im.mask), tolerance = 3, ext = 1) #profiles <- edgeProfile(im.mask.ws, ref = im.real, n = 32, # fft = FALSE, # scale = FALSE, rotate = FALSE) # #hF <- hullFeatures(im.mask.ws) #centers <- hF[,c(1,2)] # #image(im.mask.ws@.Data[,,1]) #image(im.mask.ws) # #points(centers[,1], pretty(max(centers[,2]))[2] - # centers[,2], col = 'red') #for(i in 1:nrow(profiles)){ # points(centers[i,1] + profiles[i,] * # cos(seq(-pi, pi, length.out = 32)), # pretty(max(centers[,2]))[2] - # centers[i,2] - profiles[i,] * # sin(seq(-pi, pi, length.out = 32)), # type = 'l', col = 'green') #} # trying to clean up the image.... #mk3 <- morphKern(3) #mk5 <- morphKern(5) #im.mask.clean <- erode(im.mask, mk3) #dilate(erode(closing(mask, mk5), mk3), mk5) #mask.bw <- im.mask #mask.bw[im.mask < 1] <- 0 #par(mfrow = c(2,2)) #image(im.mask) #image(im.mask.clean) #image(mask.bw) #mask.bw <- watershed(distmap(mask.bw), tolerance = 3, ext = 1) #profiles <- edgeProfile(mask.bw, ref = im.real, n = 32, # fft = FALSE, # scale = FALSE, rotate = FALSE) #hF <- hullFeatures(mask.bw) #centers <- hF[,c(1,2)] #image(mask.bw) #points(centers[,1], pretty(max(centers[,2]))[2] - # centers[,2], col = 'red') #for(i in 1:nrow(profiles)){ # points(centers[i,1] + profiles[i,] * # cos(seq(-pi, pi, length.out = 32)), # pretty(max(centers[,2]))[2] - # centers[i,2] - profiles[i,] * # sin(seq(-pi, pi, length.out = 32)), # type = 'l', col = 'green') #}