Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
You can’t perform that action at this time.