From 6858c2d744181ab8c421ebd22016154632fd62e7 Mon Sep 17 00:00:00 2001 From: dschnei5 <dschnei5@une.edu.au> Date: Fri, 22 Feb 2019 12:32:33 +1100 Subject: [PATCH] Upload New File --- ThermalImageScript_1.1.3.R | 480 +++++++++++++++++++++++++++++++++++++ 1 file changed, 480 insertions(+) create mode 100644 ThermalImageScript_1.1.3.R diff --git a/ThermalImageScript_1.1.3.R b/ThermalImageScript_1.1.3.R new file mode 100644 index 0000000..59a3631 --- /dev/null +++ b/ThermalImageScript_1.1.3.R @@ -0,0 +1,480 @@ +################Thermal Image Analysis###################### + +#1.0 - Does it all with the user measured temps and humidity values and crop coordinates file +#1.1 - For images that have correct endian - Resizes small images as well +# - Resaves .png files from flir export when binary image is .tiff +#1.1.1 - Attempts to do the same as imageJ - auto thresholding for ROI - along with full ROI +#1.1.2 - Modified for Surantha's tree canopies - FLIR ONE Camera +# - Corrects incorrect makers notes +#1.1.3 - Writes the images to ascii file converted to temp + + +####Load libraries#### + +library(exifr) +library(EBImage) +library(plyr) +library(Thermimage) +library(autothresholdr) +library(raster) + +#####Add function for later#### + +as.numeric.factor <- function(x) {as.numeric(levels(x))[x]} + + +####Set wd and get a list of images#### + +setwd(choose.dir(caption = "Select a folder full of raw and original FLIR images", default = "C:/R_2/2017-07-12_SuranthaThermal/Images_01")); +wd <- getwd(); +wd; + +files.jpg <- list.files(getwd(), pattern = ".jpg", ignore.case = TRUE, recursive = FALSE); + +####Use exifTool to extract the raw thermal image#### + +#exiftool.call(args = c("-rawthermalimage", "-b", "-w png"), fnames = files.jpg, wait = TRUE); + +####Get a list of raw thermal images#### + +files.png <- list.files(getwd(), pattern = ".png", ignore.case = TRUE, recursive = FALSE); +#file.remove(files.png); + +####Get a list of the image masks#### + +#files.tiff <- list.files(getwd(), pattern = ".tiff", ignore.case = TRUE, recursive = FALSE); + +##Read in user temps and humidity + +#files.csv <- list.files(getwd(), pattern = ".csv", ignore.case = TRUE, recursive = FALSE); + +#Uhumtemps <- read.csv(files.csv[1], header = TRUE); + + +####Get the exif info from all original jpegs#### + +Exif.info.db <- data.frame(); + + +Exif.info.db <- exifr(paste0(wd,"/",files.jpg)); + + + +####Place user relative humidity and temperature into Exif info#### + +#Exif.info.db$RelativeHumidity <- Uhumtemps$Humidity +#Exif.info.db$AtmosphericTemperature <- Uhumtemps$Temp + +#Exif.info.db$FocusDistance[8:55] <- 1.8 #Fix focus distance error in Hadriana's dataset + +####Set Object Distance into Exif Info (metres)#### Not used in this script (using focal distance data) + +#Exif.info.db$ObjectDistance <- as.numeric(readline("What is the distance between the camera and target object in metres? (1.5) ")); + +#Exif.info.db$ObjectDistance <- as.numeric(1.0) #Not using because found that focal distance is better + + + +####Create the output data frame#### + +#output.df1 <- data.frame(); + +#output.df2 <- data.frame(); + +#output.df3 <- data.frame(); + +#output.df4 <- data.frame(); + +#output.df5 <- data.frame(); + +#output.df6 <- data.frame(); + +#output.df7 <- data.frame(); + +#output.df8 <- data.frame(); + +#output.df9 <- data.frame(); + + + + +####Do it all with a loop#### + +for (j in seq_along(files.png)){ + + #j = 47; + print(paste("Image",j)); + + ####Use imagemagick to resave .png file (used if binary image is a .tiff) + + cd <- paste0("cd ",wd); + shell(cd,translate = TRUE, wait = TRUE); + cmd <- paste0("magick convert ", files.png[j]," ", files.png[j]) + shell(cmd,translate = FALSE, wait = TRUE); + + + ####Get the correct endian raw image#### + #cd <- paste0("cd ",wd); + #shell(cd,translate = TRUE, wait = TRUE); + #cmd <- paste0("exiftool -b -RawThermalImage -F ", files.jpg[j], " | magick convert - gray:- |magick convert -depth 16 -endian msb -size ", Exif.info.db$RawThermalImageWidth[j], "x", Exif.info.db$RawThermalImageHeight[j], " gray:- ",files.png[j]) + #cmd + #shell(cmd,translate = FALSE, wait = TRUE); + + + ##Use imageMagick to do the wavelet denoise - if necessary## + + #cd <- paste0("cd ",wd); + #cmd <- paste0("magick convert ",files.png[j]," -wavelet-denoise 50% ", files.png[j]); + #shell(cd,translate = TRUE, wait = TRUE); + #shell(cmd,translate = FALSE, wait = TRUE); + + ##read in image## + + #im <- readImage(paste0(wd,"/",files.png[j])); + + #print(dim(im)); + + #im <- EBImage::resize(im,640,480, filter = "bilinear"); + + #print(dim(im)); + + #EBImage::display(im); + #im.thresh <- thresh(im,w=5,h=5,offset = 0); + #EBImage::display(im.thresh) + #im.thresh.i <- ifelse(im.thresh == 1,0,1) + #display(im.thresh.i) + #im.thresh.i <- opening(im.thresh.i,makeBrush(5,shape = "disc")) + #im.thresh <- opening(im.thresh,makeBrush(5,shape = "disc")) + #oi.i <- bwlabel(im.thresh.i) + #display(normalize(oi)) + #cols.i <- c('black', sample(rainbow(max(oi.i)))) + #oirainbow.i <- Image(cols[1+oi.i], dim=dim(oi.i)) + #display(oirainbow.i) + #oi <- bwlabel(im.thresh) + #display(normalize(oi)) + #cols <- c('black', sample(rainbow(max(oi)))) + #oirainbow <- Image(cols[1+oi], dim=dim(oi)) + #display(oirainbow) + + + ##scale image## + + #max(im) + #im.ran <- Exif.info.db[j,73]; + #im.med <- Exif.info.db[j,72]; + #im.min <- Exif.info.db[j,72]-0.5*Exif.info.db[j,73]; + #im.max <- Exif.info.db[j,72]+0.5*Exif.info.db[j,73]; + #im <- im * im.ran + im.min; + #im <- im * 2^16; + + ##Threshold Image## + + #im.thresh <- as.integer(im) + + #im.thresh <- auto_thresh(im.thresh,"Otsu"); + + #im.masked <- readImage(paste0(wd,"/",files.tiff[j])) + + #mask <- ifelse(im >= im.thresh, 1,0); + + #im.masked <- as.array(im.masked); + + #im.masked.tree <- ifelse(im.masked == 1,im,NA); + + #im.masked.sky <- ifelse(im.masked == 0, im,NA); + + #hist(im.masked.tree) + #hist(im.masked.sky) + + ##Determine the number of ROIs for that image; + + #ROI.n <- (dim(Uhumtemps[j,-c(1,2,3,4)])[2] - sum(is.na(Uhumtemps[j,-c(1,2,3,4)]))) / 4; + #ROI.n <- 1 + + + ##Extract 1 Regions of Interest - This is where you need to customise for your images## + + #crop.ims <- c(); + #averages <- c(); + #maximums <- c(); + #minimums <- c(); + + #for (l in 1:ROI.n){ + + # print(paste("Entire Image ",l)); + + # nam <- ifelse(l <= 9, paste0("im0",l),paste0("im",l)); + # #assign(nam, im[Uhumtemps[j,(l-1)*4+5]:Uhumtemps[j,(l-1)*4+7],Uhumtemps[j,(l-1)*4+6]:Uhumtemps[j,(l-1)*4+8]]); + # assign(nam, im); + #a <- ifelse(dim(get(nam))[1] > dim(get(nam))[2],dim(get(nam))[1],as.numeric(dim(get(nam))[2] - dim(get(nam))[1])); + #b <- ifelse(dim(get(nam))[2] < dim(get(nam))[1],as.numeric(dim(get(nam))[1] - dim(get(nam))[2]),dim(get(nam))[2]); + #a <- ifelse(a == 0, dim(get(nam))[1], a); + #b <- ifelse(b == 0, dim(get(nam))[2], b); + #assign(nam,(get(nam)[-c(b),-c(a)])); + # d <- round(dim(get(nam))[2] * dim(get(nam))[1] * 0.001, digits = 0); + # crop.ims <- c(crop.ims,nam); + # averages <- c(averages, mean(get(nam))); + # max0.001 <- data.frame(Breaks = hist(get(nam))[[1]], Counts = c(0,hist(get(nam))[[2]])); + # max0.001 <- max0.001[with(max0.001, order(-max0.001[1])),]; + # max0.001[,2] <- ifelse(max0.001[,2] <= d, NA, max0.001[,2]); + # max0.001 <- max0.001[complete.cases(max0.001),]; + # maximums <- c(maximums,max(max0.001[,1])); + + # min0.001 <- data.frame(Breaks = hist(get(nam))[[1]], Counts = c(0,hist(get(nam))[[2]])); + # min0.001 <- min0.001[with(min0.001, order(min0.001[1])),]; + # min0.001[,2] <- ifelse(min0.001[,2] <= d, NA, min0.001[,2]); + # min0.001 <- min0.001[complete.cases(min0.001),]; + # minimums <- c(minimums,min(min0.001[,1])); + + #}; #END l loop + + #rm(l, a, b, d); + + + ########################################################################################################## + ###Tree Only### + + #crop.ims.tree <- c(); + #averages.tree <- c(); + #maximums.tree <- c(); + #minimums.tree <- c(); + +# for (l in 1:ROI.n){ + +# print(paste("Tree Crop ",l)); + +# nam.tree <- ifelse(l <= 9, paste0("im0",l,"_tree"),paste0("im",l,"_tree")); +# assign(nam.tree, im.masked.tree); +# d.tree <- round(dim(get(nam.tree))[2] * dim(get(nam.tree))[1] * 0.001, digits = 0); +# crop.ims.tree <- c(crop.ims.tree,nam.tree); +# averages.tree <- c(averages.tree, mean(get(nam.tree),na.rm = TRUE)); +# max0.001.tree <- data.frame(Breaks = hist(get(nam.tree))[[1]], Counts = c(0,hist(get(nam.tree))[[2]])); +# max0.001.tree <- max0.001.tree[with(max0.001.tree, order(-max0.001.tree[1])),]; +# max0.001.tree[,2] <- ifelse(max0.001.tree[,2] <= d.tree, NA, max0.001.tree[,2]); +# max0.001.tree <- max0.001.tree[complete.cases(max0.001.tree),]; +# maximums.tree <- c(maximums.tree,max(max0.001.tree[,1])); + +# min0.001.tree <- data.frame(Breaks = hist(get(nam.tree))[[1]], Counts = c(0,hist(get(nam.tree))[[2]])); + # min0.001.tree <- min0.001.tree[with(min0.001.tree, order(min0.001.tree[1])),]; + # min0.001.tree[,2] <- ifelse(min0.001.tree[,2] <= d.tree, NA, min0.001.tree[,2]); + # min0.001.tree <- min0.001.tree[complete.cases(min0.001.tree),]; + # minimums.tree <- c(minimums.tree,min(min0.001.tree[,1])); + # + # }; #END l loop + # + # rm(l, a.f, b.f, d.f); + # + # ########################################################################################################## + # ###Dag### + # + # crop.ims.sky <- c(); + # averages.sky <- c(); + # maximums.sky <- c(); + # minimums.sky <- c(); + # + # for (l in 1:ROI.n){ + # + # print(paste("Sky Crop ",l)); + # + # nam.sky <- ifelse(l <= 9, paste0("im0",l,"_tree"),paste0("im",l,"_tree")); + # assign(nam.sky, im.masked.sky); + # d.sky <- round(dim(get(nam.sky))[2] * dim(get(nam.sky))[1] * 0.001, digits = 0); + # crop.ims.sky <- c(crop.ims.sky,nam.sky); + # averages.sky <- c(averages.sky, mean(get(nam.sky),na.rm = TRUE)); + # max0.001.sky <- data.frame(Breaks = hist(get(nam.sky))[[1]], Counts = c(0,hist(get(nam.sky))[[2]])); + # max0.001.sky <- max0.001.sky[with(max0.001.sky, order(-max0.001.sky[1])),]; + # max0.001.sky[,2] <- ifelse(max0.001.sky[,2] <= d.sky, NA, max0.001.sky[,2]); + # max0.001.sky <- max0.001.sky[complete.cases(max0.001.sky),]; + # maximums.sky <- c(maximums.sky,max(max0.001.sky[,1])); + # + # min0.001.sky <- data.frame(Breaks = hist(get(nam.sky))[[1]], Counts = c(0,hist(get(nam.sky))[[2]])); + # min0.001.sky <- min0.001.sky[with(min0.001.sky, order(min0.001.sky[1])),]; + # min0.001.sky[,2] <- ifelse(min0.001.sky[,2] <= d.sky, NA, min0.001.sky[,2]); + # min0.001.sky <- min0.001.sky[complete.cases(min0.001.sky),]; + # minimums.sky <- c(minimums.sky,min(min0.001.sky[,1])); + # + # }; #END l loop + # + # rm(l, a.d, b.d, d.d); + # + # + # + # ##Get Temperatures - check your exif info to ensure column indexes are correct## + # + # + # output.df1 <- rbind.fill(output.df1,as.data.frame(t(c(Exif.info.db[j,3],raw2temp(averages, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j]))))); + # + # + # + # output.df2 <- rbind.fill(output.df2,as.data.frame(t(raw2temp(maximums, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # output.df3 <- rbind.fill(output.df3,as.data.frame(t(raw2temp(minimums, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # + # output.df4 <- rbind.fill(output.df4,as.data.frame(t(raw2temp(averages.tree, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # + # output.df5 <- rbind.fill(output.df5,as.data.frame(t(raw2temp(maximums.tree, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # output.df6 <- rbind.fill(output.df6,as.data.frame(t(raw2temp(minimums.tree, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # output.df7 <- rbind.fill(output.df7,as.data.frame(t(raw2temp(averages.sky, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # + # output.df8 <- rbind.fill(output.df8,as.data.frame(t(raw2temp(maximums.sky, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + # + # output.df9 <- rbind.fill(output.df9,as.data.frame(t(raw2temp(minimums.sky, + # E = Exif.info.db$Emissivity[j], + # OD = Exif.info.db$FocusDistance[j], + # RTemp = Exif.info.db$ReflectedApparentTemperature[j], + # ATemp = Exif.info.db$AtmosphericTemperature[j], + # IRWTemp = Exif.info.db$IRWindowTemperature[j], + # IRT = Exif.info.db$IRWindowTransmission[j], + # RH = Exif.info.db$RelativeHumidity[j], + # PR1 = Exif.info.db$PlanckR1[j], + # PB = Exif.info.db$PlanckB[j], + # PF = Exif.info.db$PlanckF[j], + # PO = Exif.info.db$PlanckO[j], + # PR2 = Exif.info.db$PlanckR2[j])))); + # + + im.tmp <- raster(paste0(wd,"/",files.png[j])) + + im.tmp <- raw2temp(im.tmp,E = Exif.info.db$Emissivity[j], + OD = Exif.info.db$FocusDistance[j], + RTemp = Exif.info.db$ReflectedApparentTemperature[j], + ATemp = Exif.info.db$AtmosphericTemperature[j], + IRWTemp = Exif.info.db$IRWindowTemperature[j], + IRT = Exif.info.db$IRWindowTransmission[j], + RH = Exif.info.db$RelativeHumidity[j], + PR1 = Exif.info.db$PlanckR1[j], + PB = Exif.info.db$PlanckB[j], + PF = Exif.info.db$PlanckF[j], + PO = Exif.info.db$PlanckO[j], + PR2 = Exif.info.db$PlanckR2[j]) + + + writeRaster(im.tmp,paste0("out_",files.jpg[j]),format = "ascii") + rm(im.tmp) + +}; # END j Loop + +#rm(j); + + +##Add some column names, write a .csv and cleanup + +# output.df10 <- cbind(output.df1,output.df2,output.df3,output.df4,output.df5,output.df6,output.df7,output.df8,output.df9) +# col.ns <- 1:((dim(output.df10)[2] - 1)/9) +# colnames(output.df10) <- c("ImageName",ifelse(col.ns <= 9,paste0("All_ROI_Average_0",col.ns), paste0("All_ROI_Average_",col.ns)),ifelse(col.ns <= 9,paste0("All_ROI_Maximum_0",col.ns), paste0("All_ROI_Maximums_",col.ns)),ifelse(col.ns <= 9,paste0("All_ROI_Minimum_0",col.ns), paste0("All_ROI_Minimums_",col.ns)),ifelse(col.ns <= 9,paste0("Tree_ROI_Average_0",col.ns), paste0("Tree_ROI_Average_",col.ns)),ifelse(col.ns <= 9,paste0("Tree_ROI_Maximum_0",col.ns), paste0("Tree_ROI_Maximums_",col.ns)),ifelse(col.ns <= 9,paste0("Tree_ROI_Minimum_0",col.ns), paste0("Tree_ROI_Minimums_",col.ns)), ifelse(col.ns <= 9,paste0("Sky_ROI_Average_0",col.ns), paste0("Sky_ROI_Average_",col.ns)),ifelse(col.ns <= 9,paste0("Sky_ROI_Maximum_0",col.ns), paste0("Sky_ROI_Maximums_",col.ns)),ifelse(col.ns <= 9,paste0("Sky_ROI_Minimum_0",col.ns), paste0("Sky_ROI_Minimums_",col.ns))) +# +# #####Stick them all together and export############# +# +# write.table(t(output.df10), file.choose(new = TRUE),sep = ",", col.names = FALSE); + +rm(list = ls()); +######################################################################################################################### + -- GitLab