R test script for calibrating JPEG images - this is a work in progress
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment