R test script for calibrating JPEG images - this is a work in progress

autoCalibrateJPEG.R

setwd("/media/nedhorning/684EE5FF4EE5C642/AMNH/PhotoMonitoring/FilterTests/ChrisPhotos")
library(raster)
library(maptools)
library(stringr)

#############################   SET VARIABLES HERE  ###################################
# Enter the value for gamma *** Vout = Vin ^ gamma : Vin = Vout ^ (1/gamma)  : gamma = log(Vout, Vin)  : Vin = RAW DN and Vout = JPEG DN
blueGamma <- 0.8
redGamma <- 0.8

# Flag to subtract the blue band from the red band
subtractBlue <- TRUE

# Value to specify percent of the blue band to use (ignored if subtractBlue is FALSE)
percentBlue <- 0.8

# Flag to pause to script after printing first graph
useReturn <- FALSE

# Name of the attribute that holds the integer target type identifyer
# attName <- 'targetName'

# Input image name
imageName <- "DB660-850Rori.JPG"

# Set wavelength for camera visible and NIR
refVisWavelength <- 660
refNIR_Wavelength <- 850

createPDF <- TRUE
########################################################################################3
# Create output base file name
baseName <- str_sub(imageName, 1, str_locate(imageName, ignore.case(".jpg"))[1]-1)
if (subtractBlue) {
  outBaseName <- paste(baseName, refVisWavelength, "_", refNIR_Wavelength, "SubtractBlue", percentBlue*100, sep="")
} else {
  outBaseName <- paste(baseName, refVisWavelength, "_", refNIR_Wavelength, sep="")
}

# Read the image bands 
cat("Reading image\n")
inImage <- brick(imageName)
names(inImage) <- c("red", "green", "blue")

redBand <- raster(inImage, layer=1)
greenBand <- raster(inImage, layer=2)
blueBand <- raster(inImage, layer=3)

blueBand <- blueBand ^ (1/blueGamma)
redBand <- redBand ^ (1/redGamma)

rgbImage <- stack(redBand, greenBand, blueBand)

# PLot single band and define white and black target locatinos
plot(greenBand, col=gray(seq(0, 1, length.out=256)))
cat("Click two points to define white target rectangle\n")
whiteTargetExtent <- drawExtent()

cat("Click two points to define black target rectangle\n")
blackTargetExtent <- drawExtent()

# Create data structures for storing values
rgbValues <- data.frame()
targetNumbers <- character()

# Calculate mean pixel values under each target polygon for each image band
cat("Extracting target pixels\n")
whiteTargetValues <- extract(rgbImage, whiteTargetExtent, df=TRUE)
whiteTargetMeans <- apply(whiteTargetValues, 2, mean)
rgbValues <- rbind(rgbValues, whiteTargetMeans)
blackTargetValues <- extract(rgbImage, blackTargetExtent, df=TRUE)
blackTargetMeans <- apply(blackTargetValues, 2, mean)
rgbValues <- rbind(rgbValues, blackTargetMeans)

photoRGB <- cbind(c("printerPaper", "tarPaper"), rgbValues) 
names(photoRGB) <- c("targetName", "red", "green", "blue")

# Read CSV file with target reference spectra
specDataRaw <- read.csv("/media/nedhorning/684EE5FF4EE5C642/AMNH/PhotoMonitoring/ColorTests/CalibrationTargets/InfragramCalSamples.csv")
tarPaper <- (specDataRaw[,8] + specDataRaw[,9]) / 2
printerPaper <- (specDataRaw[,10] + specDataRaw[,11]) / 2
avgAll <- cbind(specDataRaw[,1], tarPaper, printerPaper)

# Pull out reference data for the specify visible and NIR wavelengths
refVis <- avgAll[which(avgAll[,1]==refVisWavelength),-1] 
refNIR <- avgAll[which(avgAll[,1]==refNIR_Wavelength),-1]

# Put sample data into a data frame and order it so photo and reference data match
targetMapping <- match(names(refVis), photoRGB[,1])
photoRGB <- photoRGB[targetMapping,]
dfVis <- data.frame(refVis, photoRGB)
dfNIR <- data.frame(refNIR, photoRGB)

if (subtractBlue) {
  fitRedBlue <- lm(blue ~ red , data=dfVis)
  dfVis$red <- c(dfVis$red - (dfVis$blue*percentBlue))
  predBlue <- predict(redBand, fitRedBlue)
  redBand <- redBand - (predBlue * percentBlue)
  names(redBand) <- "red"
  rgbImage <- stack(redBand, greenBand, blueBand)
}

blue <- photoRGB$blue
green <- photoRGB$green
red = photoRGB$red

# Calculate visible linear regression
cat("Predicting vis band regression\n")
if (refVisWavelength >550) {
  fitVis <- lm(refVis ~ red , data=dfVis)
 
  visBandCal <- predict(redBand, fitVis)
  print(summary(fitVis))
  cat("Predicting vis band regression\n")

  # Calculate NIR linear regression using single and two bands
  fitNIR <- lm(refNIR ~ blue, data=dfNIR)
  print(summary(fitNIR))
  cat("Predicting NIR band single regression\n")
  NIRBandCal <- predict(blueBand, fitNIR)
} else {
  fitVis <- lm(refVis ~ blue , data=dfVis)
  print(summary(fitVis))
  cat("Predicting vis band single regression\n")
  visBandCal <- predict(blueBand, fitVis)

  # Calculate NIR linear regression using single and two bands
  fitNIR <- lm(refNIR ~ red, data=dfNIR)
  print(summary(fitNIR))
  cat("Predicting NIR band single regression\n")
  NIRBandCal <- predict(redBand, fitNIR)
}

cat("Calculating NDVI\n")
ndviCal <- (NIRBandCal - visBandCal) / (NIRBandCal + visBandCal)

# Clip NDVI values under 0 and NA to 0 and over 1 to 1
ndviCal[(ndviCal < -1) | is.na(ndviCal)] <- -1
ndviCal[ndviCal > 1] <- 1

# Determine output image names
outName <- paste(outBaseName, "_NDVI.tif", sep="")

cat("Writing image\n")
writeRaster(ndviCal, filename=outName, format='GTiff', datatype='FLT4S', overwrite=TRUE)

plot(ndviCal)