################################################################ ### FUNCTION TO CALCULATE THE MEAN LEAF AREA FROM A JPG IMAGE ################################################################ # Author: Gabriel Arellano # Warranty: none # The leaf area is usually calculated with ImageJ "by hand". # ImageJ is a extremely powerful and free software that allows lots of calculations. Nothing personal against it. # However, it requires some interaction with the user. This takes time and make results non-replicable. # The algorithm for the leaf area calculation is, in principle, very simple. Why no-one uses R to do that? # The problem comes if one wants to establish the scale at which the photograph (or scanned image) has been taken # ImageJ has some tools to do that "by hand". # This R function does that automatically. # To do so, it requires the picture to be taken over a particular white-and-black frame # That frame is cheap and super-fast to build, I invested 4 minutes and 1$ in doing mine. # See some example of this frame at my webpage (www.gabrielareto.weebly.com) # The most important two things: # (1) do not place anything except leaves in the white "cuore" (the center of the frame); do not write there. # (2) do not place anything (at least anything clear) in the black strips that define the frame. leafArea <- function(file, bw.threshold=0.5, big.peak=0.95, area.cuore=30*45, number.of.leaves=3, trace=T) { # ARGUMENTS: # file: the path to the file with the picture # bw.threshold: the threshold to classify a pixel as black or white; 0.5 by default # big.peak is the quantile above which a peak is considered a peak; 0.95 by default # area.cuore is the area of the white cuore, in cm^2 # number.of.leaves: defined by the user; by default 3 # trace: if TRUE, shows how the image is processed. Note that this can be very slow!!!! # RETURNS: # the mean are per leaf, in cm^2 # packages required require(jpeg) require(raster) # Some functions to wrap: matrixToRaster <- function(x) raster(list(x=1:nrow(x), y=1:ncol(x), z=x)) ReduceResolution <- function(x, p=200000) # p is the desired sixe in pixels, approximately { require(raster) size.now = dim(x)[1]*dim(x)[2] r <- raster(list(x=1:nrow(x), y=1:ncol(x), z=x)) # creates a raster return(t(as.matrix(aggregate(r, fact=round(sqrt(size.now/p)))))) } # Function to define approximate peaks of a vector, one per half: findTwoPeaks <- function(x, q = big.peak) { # ARGUMENTS: # x is a vector # q is the quantile above which we consider "high peak" or not jumps <- abs(c(x, NA) - c(NA, x)) # the contrast, the "jump" from clear to dark or from dark to clear interesting.indexes <- which(jumps > quantile(jumps, probs=q, na.rm=T)) # places where the jumps are very large distances.to.center <- interesting.indexes - length(x)/2 interesting.indexes.first.half <- interesting.indexes[distances.to.center < 0] interesting.indexes.second.half <- interesting.indexes[distances.to.center > 0] closest.to.center.first.half <- interesting.indexes.first.half[length(interesting.indexes.first.half)] closest.to.center.second.half <- interesting.indexes.second.half[1] #plot(x); abline(v= closest.to.center.first.half); abline(v= closest.to.center.second.half) # check it! it works nicely! return(c(closest.to.center.first.half=closest.to.center.first.half, closest.to.center.second.half=closest.to.center.second.half)) } # Function to explore the surroundings of a given position of a vector # it returns the closest point to the left of the peak and to the right of a peak closestNonPeaks <- function(x, index) { # x is a vector # index is the interesting area (a proxy for a peak) # plot: logical k = round(length(x)/10) # how far we will go from the given index alrededores <- (index-k):(index +k) sub.x <- x[alrededores] black <- sub.x > (max(sub.x)-min(sub.x)) first.white.before <- alrededores [which(black)[1]-1] first.white.after <- alrededores [which(black)[sum(black)]+1] #plot(x); abline(v=first.white.before); abline(v=first.white.after) # check it! it works nicely! return(c(before= first.white.before, after= first.white.after)) } # Reads image: image <- readJPEG(file) # Promediates the "darkness" in all channels image <- 1-(image[,,1] + image[,,2] + image[,,3])/3 if(trace==TRUE) { par(mfrow=c(2, 3)) image(image, main="Gray scale") } # Reduces resolution: image <- ReduceResolution(x = image) if(trace==TRUE) image(image, main="Reduced resolution") # Binary classification: black or white image[image >= bw.threshold] <- 1 image[image < bw.threshold] <- 0 if(trace==TRUE) image(image, main="Black and white") # Makes the first subset of the image: x1 <- rowSums(image) # two vectors, darkness along X and darkness along Y x2 <- colSums(image) # peaks (approximately) a <- findTwoPeaks(x=x1, q= big.peak) b <- findTwoPeaks(x=x2, q= big.peak) # exact inner borders of those peaks crop1 <- closestNonPeaks(x1, a["closest.to.center.first.half"])["after"] crop2 <- closestNonPeaks(x1, a["closest.to.center.second.half"])["before"] crop3 <- closestNonPeaks(x2, b["closest.to.center.first.half"])["after"] crop4 <- closestNonPeaks(x2, b["closest.to.center.second.half"])["before"] # Subsets the matrix: image <- image[crop1:crop2, crop3:crop4] if(trace==TRUE) image(image, main="Selects the center (of known area)") # Gives the scale: area.of.one.pixel <- area.cuore/length(image) # Second round of subsetting the cuore a <- rowSums(image) # darkness b <- colSums(image) a2 <- round(a/(length(a)/10))*(length(a)/10) # "takes to the ground" b2 <- round(b/(length(b)/10))*(length(b)/10) crop1 <- which(a2==min(a2))[1] crop2 <- which(a2==min(a2))[length(which(a2==min(a2)))] crop3 <- which(b2==min(b2))[1] crop4 <- which(b2==min(b2))[length(which(b2==min(b2)))] image <- image[crop1:crop2, crop3:crop4] if(trace==TRUE) image(image, main="Remove borders: first attemp") # And a third round, and everything should be perfect... a <- rowSums(image) # darkness b <- colSums(image) crop1 <- which(a==0)[1] crop2 <- which(a==0)[length(which(a==0))] crop3 <- which(b==0)[1] crop4 <- which(b==0)[length(which(b==0))] image <- image[crop1:crop2, crop3:crop4] if(trace==TRUE) image(image, main="THIS MUST BE ONLY LEAVES") # Makes the final calculations: leaf.area.in.pixels <- sum(image)/number.of.leaves return(leaf.area.in.pixels * area.of.one.pixel) }