Created
July 10, 2015 19:42
-
-
Save nedhorning/f384d9afcc043741def3 to your computer and use it in GitHub Desktop.
R test script for calibrating JPEG images - this is a work in progress
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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