# Leaf Venation Analysis # W. A. Green # # last modified: 2016-02-03 cat(date()) # start time # load libraries library(EBImage) source('http://bricol.net/research/imageanalysis/RFunctions/image_analysis_functions.R') # set global variables PATH <- '~/Desktop/vein_density/forSLW/v2/' THRESH_WIN = 30 # #(threshold window size) THRESH_SENSE = 0.01 # #(threshold sensitivity) SCHMUTZ_MIN <- 81 #(large-pass filter size; 9 pixels sq.) # Read in image(s) ims <- dir(paste(PATH, 'images', sep = '')) ims_clean <- gsub('\\.\\w*', '', ims, perl = TRUE) #image names without file extensions # Make a frame to collect data dat <- as.data.frame(cbind(rep(0, length(ims)), rep(0, length(ims)))) names(dat) <- c('mode_area', 'mode_log_dm') rownames(dat) <- ims_clean for(i in seq_along(ims)){ #loop through all images #i <- 3 #just one image raw <- readImage(paste(PATH, 'images/', ims[i], sep = '')) fnm <- paste(ims_clean[i], '_raw.tiff', sep = '') #fnm for 'filename' writeImage(raw, fnm, bits.per.sample = 16) bin.cl <- clean(raw, thresh_win = THRESH_WIN, thresh_sense = THRESH_SENSE, schmutz = SCHMUTZ_MIN) fnm <- paste(ims_clean[i], '_cl.tiff', sep = '') writeImage(bin.cl, fnm, bits.per.sample = 16) # Modal areole size bin.lab <- bwlabel(bin.cl) bin.fts <- computeFeatures.shape(bin.lab) ar_sizes <- bin.fts[,'s.area'] dens <- density(ar_sizes) mod <- dens$x[dens$y == max(dens$y)] dat[i, 1] <- mod fnm <- paste(ims_clean[i], '_area.pdf', sep = '') pdf(file = fnm, width = 7, height = 5) hist(ar_sizes, main = 'Areole size distribution', xlab = 'Areole area in pixels', freq = FALSE) lines(dens) dev.off() # Size transform for areoles: # # bin.st <- st(bin.cl, max_mask = floor(sqrt(mean(ar_sizes)))) # # save(bin.st, file = 'areole.st.R') # display(bin.st$sizes/max(bin.st$sizes)) # fnm <- paste(ims_clean[i], '_st.tiff', sep = '') # writeImage(bin.st$sizes/max(bin.st$sizes), fnm, bits.per.sample = 16) # fnm <- paste(ims_clean[i], '_st.pdf', sep = '') # pdf(file = fnm, width = 7, height = 5) # plot(as.table(table(bin.st$sizes)[-1]), type = 'h', # ylab = 'Frequency', xlab = 'Diameter in Pixels', # main = 'Sizing transform') # dev.off() # Distance map for areoles: bin.dm <- distmap(bin.cl) # save(bin.dm, file = 'dm.R') dm.bar <- mean(bin.dm) #dm.log.bar <- log(bin.dm) #dm.log.bar[is.infinite(dm.log.bar)] <- 0 #dm.log.bar <- mean(dm.log.bar) dens <- density(log(as.numeric(bin.dm)), adjust = 10) mod <- dens$x[dens$y == max(dens$y)] dat[i, 2] <- mod # display(bin.dm/max(bin.dm)) fnm <- paste(ims_clean[i], '_dm.tiff', sep = '') writeImage(bin.dm/max(bin.dm), fnm, bits.per.sample = 16) fnm <- paste(ims_clean[i], '_dm.pdf', sep = '') pdf(file = fnm, height = 5, width = 12) par(mfrow = c(1,2)) # not log tranformed hist(as.numeric(bin.dm), main = 'Areole distance map', freq = FALSE, xlab = 'Euclidean distance from nearest vein in pixels', ylab = 'Frequency') lines(density(as.numeric(bin.dm), adjust = 10), col = 'red') abline(v = dm.bar, col = 'grey', lwd = 3) text(x = dm.bar, y = 0.55, 'Mean', col = 'grey', xpd = NA) # log tranformed hist(log(as.numeric(bin.dm)), main = 'Areole distance map', freq = FALSE, border = 'grey', xlab = 'Euclidean distance from nearest vein in log pixels') lines(dens, col = 'red') #abline(v = dm.log.bar, col = 'green', lwd = 3) abline(v = mod, col = 'red', lwd = 3) text(x = mod, y = 0.73, 'Mode', col = 'red', xpd = NA) abline(v = log(dm.bar), col = 'grey', lwd = 3) text(x = log(dm.bar), y = 0.73, 'Mean', col = 'grey', xpd = NA) dev.off() } #end loop through images cat(date()) # end time dat <- cbind(dat, scale = gsub('\\w*\\_', '', rownames(dat))) write.table(dat, file = 'data_out.csv', sep = '\t') pdf(file = 'area_vs_dm.pdf', width = 7, height = 5) plot(dat[,1], exp(dat[,2]), log = 'xy', pch = substr(dat[,3], 1, 1), main = 'Area versus distance map', xlab = 'Modal areole area', ylab = 'Modal distance map') dev.off() pdf(file = 'dm_pairs.pdf', width = 7, height = 7) pairs(cbind(macro = dat[,2][dat$scale == 'macro'], crop = dat[,2][dat$scale == 'crop'], crop2 = dat[,2][dat$scale == '2']), main = 'Mode of log-transformed distance map density') dev.off() pdf(file = 'area_pairs.pdf', width = 7, height = 7) pairs(cbind(macro = dat[,1][dat$scale == 'macro'], crop = dat[,1][dat$scale == 'crop'], crop2 = dat[,1][dat$scale == '2']), main = 'Mode of areole area density', log = 'xy') dev.off() pdf(file = 'macro_vs_2.pdf', width = 7, height = 5) plot(dat[,2][dat$scale == 'macro'], dat[,2][dat$scale == '2'], main = 'Macro versus Crop 2 modal log distance map', xlab = 'Macro', ylab = 'Crop 2', pch = '') labs <- gsub('\\_\\w*', '', rownames(dat)[dat$scale == 'macro']) text(dat[,2][dat$scale == 'macro'], dat[,2][dat$scale == '2'], labs, cex = 0.5) dev.off()