Skip to content
Snippets Groups Projects
Commit 6858c2d7 authored by dschnei5's avatar dschnei5
Browse files

Upload New File

parent addf3950
Branches patch-1
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
ThermalImageScript_1.1.3.R 0 → 100644
+ 480
− 0
View file @ 6858c2d7
################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());
#########################################################################################################################
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment