README.md

Thermimage: Thermal Image Analysis

cran version downloads total downloads Research software impact

This is a collection of functions for assisting in converting extracted raw data from infrared thermal images and converting them to estimate temperatures using standard equations in thermography. Provides an open source proxy tool for assisting with infrared thermographic analysis.

Recent/release notes

Features

Installation

On current R (>= 3.0.0)

install.packages("Thermimage")
library("devtools"); install_github("gtatters/Thermimage",dependencies=TRUE)

Package Imports

OS Requirements

Import, Export, Image Processing

A typical thermal image

Galapagos Night Heron

Normally, these thermal images require access to software that only runs on Windows operating system. This package will allow you to import certain FLIR jpgs and videos and process the images.

Import FLIR JPG

To load a FLIR JPG, you first must install Exiftool as per instructions above. Open sample flir jpg included with Thermimage package:

library(Thermimage)
f<-paste0(system.file("extdata/IR_2412.jpg", package="Thermimage"))
img<-readflirJPG(f, exiftoolpath="installed")
dim(img)

> [1] 480 640

The readflirJPG function has used Exiftool to figure out the resolution and properties of the image file. Above you can see the dimensions are listed as 480 x 640. Before plotting or doing any temperature assessments, let's extract the meta-tages from the thermal image file.

Extract meta-tags from thermal image file

cams<-flirsettings(f, exiftoolpath="installed", camvals="")

This produes a rather long list of meta-tags. If you only want to see your camera calibration constants, type:

plancks<-flirsettings(f, exiftoolpath="installed", camvals="-*Planck*")
unlist(plancks$Info)

> PlanckR1       PlanckB       PlanckF       PlanckO      PlanckR2 
> 2.110677e+04  1.501000e+03  1.000000e+00 -7.340000e+03  1.254526e-02 

If you want to check the file data information, type:

cbind(unlist(cams$Dates))

> FileModificationDateTime "2017-03-22 22:15:09"
> FileAccessDateTime       "2017-03-22 23:27:31"
> FileInodeChangeDateTime  "2017-03-22 22:15:10"
> ModifyDate               "2013-05-09 16:22:23"
> CreateDate               "2013-05-09 16:22:23"
> DateTimeOriginal         "2013-05-09 22:22:23"

or just:

cams$Dates$DateTimeOriginal

> [1] "2013-05-09 22:22:23"

The most relevant variables to extract for calculation of temperature values from raw A/D sensor data are listed here. These can all be extracted from the cams output as above. I have simplified the output below, since dealing with lists can be awkward.

Emissivity<-  cams$Info$Emissivity                    # Image Saved Emissivity - should be ~0.95 or 0.96
dateOriginal<-cams$Dates$DateTimeOriginal             # Original date/time extracted from file
dateModif<-   cams$Dates$FileModificationDateTime     # Modification date/time extracted from file
PlanckR1<-    cams$Info$PlanckR1                      # Planck R1 constant for camera  
PlanckB<-     cams$Info$PlanckB                       # Planck B constant for camera  
PlanckF<-     cams$Info$PlanckF                       # Planck F constant for camera
PlanckO<-     cams$Info$PlanckO                       # Planck O constant for camera
PlanckR2<-    cams$Info$PlanckR2                      # Planck R2 constant for camera
OD<-          cams$Info$ObjectDistance                # object distance in metres
FD<-          cams$Info$FocusDistance                 # focus distance in metres
ReflT<-       cams$Info$ReflectedApparentTemperature  # Reflected apparent temperature
AtmosT<-      cams$Info$AtmosphericTemperature        # Atmospheric temperature
IRWinT<-      cams$Info$IRWindowTemperature           # IR Window Temperature
IRWinTran<-   cams$Info$IRWindowTransmission          # IR Window transparency
RH<-          cams$Info$RelativeHumidity              # Relative Humidity
h<-           cams$Info$RawThermalImageHeight         # sensor height (i.e. image height)
w<-           cams$Info$RawThermalImageWidth          # sensor width (i.e. image width)

Convert raw binary to temperature

Now you have the img loaded, look at the values:

str(img)

> int [1:480, 1:640] 18090 18074 18064 18061 18081 18057 18092 18079 18071 18071 ...

If stored with a TIFF header, the data load in as a pre-allocated matrix of the same dimensions of the thermal image, but the values are integers values, in this case ~18000. The data are stored as in binary/raw format at 2^16 bits of resolution = 65535 possible values, starting at 1. These are not temperature values. They are, in fact, radiance values or absorbed infrared energy values in arbitrary units. That is what the calibration constants are for. The conversion to temperature is a complicated algorithm, incorporating Plank's law and the Stephan Boltzmann relationship, as well as atmospheric absorption, camera IR absorption, emissivity and distance to namea few. Each of these raw/binary values can be converted to temperature, using the raw2temp function:

temperature<-raw2temp(img, ObjectEmissivity, OD, ReflT, AtmosT, IRWinT, IRWinTran, RH,
                      PlanckR1, PlanckB, PlanckF, PlanckO, PlanckR2)
str(temperature)      

> num [1:480, 1:640] 23.7 23.6 23.6 23.6 23.7 ...

The raw binary values are now expressed as temperature in degrees Celsius (apologies to Lord Kelvin). Let's plot the temperature data:

library(fields) # should be loaded imported when installing Thermimage
plotTherm(t(temperature), h, w)

FLIR JPG on import

The FLIR jpg imports as a matrix, but default plotting parameters leads to it being rotated 270 degrees (counter clockwise) from normal perspective, so you should either rotate the matrix data before plotting, or include the rotate270.matrix transformation in the call to the plotTherm function:

plotTherm(temperature, w=w, h=h, minrangeset = 21, maxrangeset = 32, trans="rotate270.matrix")

FLIR JPG rotate 270

If you prefer a different palette:

plotTherm(temperature, w=w, h=h, minrangeset = 21, maxrangeset = 32, trans="rotate270.matrix", 
          thermal.palette=rainbowpal)
plotTherm(temperature, w=w, h=h, minrangeset = 21, maxrangeset = 32, trans="rotate270.matrix", 
          thermal.palette=glowbowpal)
plotTherm(temperature, w=w, h=h, minrangeset = 21, maxrangeset = 32, trans="rotate270.matrix", 
          thermal.palette=midgreypal)
plotTherm(temperature, w=w, h=h, minrangeset = 21, maxrangeset = 32, trans="rotate270.matrix", 
          thermal.palette=midgreenpal)

FLIR JPG rotate 270 rainbow palette

FLIR JPG rotate 270 glowbow palette

FLIR JPG rotate 270 midgrey palette

Export Image or Video

Finding a way to quantitatively analyse thermal images in R is a challenge due to limited interactions with the graphics environment. Thermimage has a function that allows you to write the image data to a file format that can be imported into ImageJ.

First, the image matrix needs to be transposed (t) to swap the row vs. column order in which the data are stored, then the temperatures need to be transformed to a vector, a requirement of the writeBin function. The function writeFlirBin is a wrapper for writeBin, and uses information on image width, height, frame number and image interval (the latter two are included for thermal video saves) but are kept for simplicity to contruct a filename that incorporates image information required when importing to ImageJ:

writeFlirBin(as.vector(t(temperature)), templookup=NULL, w=w, h=h, I="", rootname="FLIRjpg")

The raw file can be found here: https://github.com/gtatters/Thermimage/blob/master/READMEimages/FLIRjpg_W640_H480_F1_I.raw?raw=true

Import Raw File into ImageJ

The .raw file is simply the pixel data saved in raw format but with real 32-bit precision. This means that the temperature data (negative or positive values) are encoded in 4 byte chunks. ImageJ has a plethora of import functions, and the File-->Import-->Raw option provides great flexibility. Once opening the .raw file in ImageJ, set the width, height, number of images (i.e. frames or stacks), byte storage order (little endian), and hyperstack (if desired):

ImageJ Import Settings

The image imports clearly just as it would in a thermal image program. Each pixel stores the calculated temperatures as provided from the raw2temp function above.

Image Imported into ImageJ

Importing Thermal Videos

Importing thermal videos (March 2017: still in development) is a little more involved and less automated, but below are steps that have worked for seq and fcf files tested.

Set file info and extract meta-tags as done above:

# set filename as v
v<-paste0(system.file("extdata/SampleSEQ.seq", package="Thermimage"))

# Extract camera values using Exiftool (needs to be installed)
camvals<-flirsettings(v)
w<-camvals$Info$RawThermalImageWidth
h<-camvals$Info$RawThermalImageHeight

Create a lookup variable to convert the raw binary to actual temperature estimates, use parameters relevant to the experiment. You could use the values stored in the FLIR meta-tags, but these are not necessarily correct for the conditions of interest. suppressWarnings() is used because of NaN values returned for binary values that fall outside the range.

suppressWarnings(
templookup<-raw2temp(raw=1:65535, E=camvals$Info$Emissivity, OD=camvals$Info$ObjectDistance, RTemp=camvals$Info$ReflectedApparentTemperature, ATemp=camvals$Info$AtmosphericTemperature, IRWTemp=camvals$Info$IRWindowTemperature, IRT=camvals$Info$IRWindowTransmission, RH=camvals$Info$RelativeHumidity, PR1=camvals$Info$PlanckR1,PB=camvals$Info$PlanckB,PF=camvals$Info$PlanckF,PO=camvals$Info$PlanckO,PR2=camvals$Info$PlanckR2)
)
plot(templookup, type="l", xlab="Raw Binary 16 bit Integer Value", ylab="Estimated Temperature (C)")

Binary to Temperature Conversion

The advantage of using the templookup variable is in its index capacity. For computations involving large files, this is most efficient way to convert the raw binary values rapidly without having to call the raw2temp function repeatedly. Thus, for a raw binary value of 17172, 18273, and 24932:

templookup[c(17172, 18273, 24932)]

> [1] 18.30964 24.77935 57.07821

We will use the templookup later on, but first to detect where the image frames can be found in the video file. Using the width and height information, we use this to find where in the video file these are stored. This corresponds to reproducible locations in the frame header:

fl<-frameLocates(v, w, h)
n.frames<-length(fl$f.start)
n.frames; fl

> [1] 2
> $h.start
> [1]    162 308688
> 
> $f.start
> [1]   1391 309917

The relative positions of the header start (h.start) are 162 and 308688, and the frame start (f.start) positions are 1391 and 309917. The video file is a short, two frame (n.frames) sequence from a thermal video.

Then pass the fl data to two different functions, one to extract the time information from the header, and the other to extract the actual pixel data from the image frame itself. The lapply function will have to be used (for efficiency), but to wrap the function across all possible detected image frames. Note: For large files, the parallel function, mclapply, is advised (?getFrames for an example):

extract.times<-do.call("c", lapply(fl$h.start, getTimes, vidfile=v))
data.frame(extract.times)

> 1 2012-06-13 15:52:08.698
> 2 2012-06-13 15:52:12.665

Interval<-signif(mean(as.numeric(diff(extract.times))),3)
Interval

> [1] 3.97

This particluar sequence was actually captured at 0.03 sec intervals, but the sample file in the package was truncated to only two frames to minimise online size requirements for CRAN. At present, the getTimes function cannot accurately render the time on the first frame. On the original 100 frame file, it accurately captures the real time stamps, so the error is appears to be how FLIR saves time stamps (save time vs. modification time vs. original time appear highly variable in .seq and .fcf files). Precise time capture is not crucial but is helpful for verifying data conversion.

After extracting times, then extract the frame data, with the getFrames function:

alldata<-unlist(lapply(fl$f.start, getFrames, vidfile=v, w=w, h=h))
class(alldata); length(alldata)/(w*h)

> [1] "integer"
> [1] 2

The raw binary data are stored as an integer vector. length(alldata)/(w*h) verifies the total # of frames in the video file is 2.

It is best to convert the temperature data in the following manner, although depending on file size and system limits, you may wish to delay converting to temperature until writing the file.

alltemperature<-templookup[alldata]

I recommend converting the binary and/or temperature variables to a matrix class, where each column represents a separate image frame, while the individual rows correspond to unique pixel positions. Pixels are filled into the row values the same way across all frames. Dataframes and arrays are much slower for processing large files.

alldata<-unname(matrix(alldata, nrow=w*h, byrow=FALSE))
alltemperature<-unname(matrix(alltemperature, nrow=w*h, byrow=FALSE))
dim(alltemperature)

> [1] 307200      2

Frames extracted from thermal vids are upside down, so use the mirror.matrix function inside the plotTherm function.

plotTherm(alltemperature[,1], w=w, h=h, trans="mirror.matrix")
plotTherm(alltemperature[,2], w=w, h=h, trans="mirror.matrix")

Sample Sequence Frame 1

Sample Sequence Frame 2

Now, export entire sequence to a raw bin for opening in ImageJ - smallish file size

writeFlirBin(bindata=alldata, templookup, w, h, Interval, rootname="SampleSEQ")

The newly written 32=bit video file (https://github.com/gtatters/Thermimage/blob/master/READMEimages/SampleSEQ_W640_H480_F2_I3.97.raw?raw=true) can now be imported into ImageJ, as desribed above for the single image. Each frame is converted into a stack in ImageJ.

Heat Transfer Calculations

The other functions in Thermimage relate to steady state modelling of heat exchange using surface temperatures extracted from thermal images. Provided below are basic steps and other information that is required to make use of the heat transfer calculations. References for these calculations are found within the Thermimage package.

Minimum required information (units)

Required if working outdoors with solar radiation

Can be estimated or provided

Ground Temperature Estimation

Assemble paramters

# Create a randomly generated sample:

Ta<-rnorm(20, 25, sd=10)
Ts<-Ta+rnorm(20, 5, sd=1)
RH<-rep(0.5, length(Ta))
SE<-rnorm(20, 400, sd=50)
Tg<-Tground(Ta,SE)
A<-rep(0.4,length(Ta))
L<-rep(0.1, length(Ta))
V<-rep(1, length(Ta))
shape<-rep("hcylinder", length(Ta))
c<-forcedparameters(V=V, L=L, Ta=Ta, shape=shape)$c
n<-forcedparameters(V=V, L=L, Ta=Ta, shape=shape)$n
a<-freeparameters(L=L, Ts=Ts, Ta=Ta, shape=shape)$a
b<-freeparameters(L=L, Ts=Ts, Ta=Ta, shape=shape)$b
m<-freeparameters(L=L, Ts=Ts, Ta=Ta, shape=shape)$m
type<-rep("forced", length(Ta))
rho<-rep(0.1, length(Ta))
cloud<-rep(0, length(Ta))

d<-data.frame(Ta, Ts, Tg, SE, RH, rho, cloud, A, V, L, c, n, a, b, m, type, shape)
head(d)

>          Ta        Ts       Tg       SE  RH rho cloud   A V   L     c     n a    b    m   type     shape
> 1 28.322047 34.426894 32.67210 359.5081 0.5 0.1     0 0.4 1 0.1 0.174 0.618 1 0.58 0.25 forced hcylinder
> 2 19.295451 23.458105 23.72816 366.3394 0.5 0.1     0 0.4 1 0.1 0.174 0.618 1 0.58 0.25 forced hcylinder
> 3 23.640834 26.932211 28.29766 384.8615 0.5 0.1     0 0.4 1 0.1 0.174 0.618 1 0.58 0.25 forced hcylinder
> 4  6.971665  8.822035 12.14272 427.3600 0.5 0.1     0 0.4 1 0.1 0.174 0.618 1 0.58 0.25 forced hcylinder
> 5 32.594745 39.277282 38.40423 480.1226 0.5 0.1     0 0.4 1 0.1 0.174 0.618 1 0.58 0.25 forced hcylinder
> 6 22.613530 28.058783 27.91851 438.4282 0.5 0.1     0 0.4 1 0.1 0.174 0.618 1 0.58 0.25 forced hcylinder

Applying the heat transfer calculations

(qtotal<-qrad() + qconv())

> [1] -186.8849

Notice how the above example yielded an estimate without any inputs. This is because there are default values in all the functions. In this case, the estimate is negative, meaning a net loss of heat to the environment. It's units are in W/m2. To convert the above measures into total heat flux (Q), the Area of each part is required.

Area1<-0.2 # units are in m2
Area2<-0.3 # units are in m2
Qtotal1<-qrad()*Area1 + qconv()*Area1
Qtotal2<-qrad()*Area2 + qconv()*Area2
(QtotalAll<-Qtotal1 + Qtotal2)

> [1] -93.44246

This approach is used in animal thermography, such that all component shapes sum to estimate entire body heat exchange: WholeBody Heat Exchange = Qtotal1 + Qtotal2 + Qtotal3 ... Qtotaln, where 1-n represent the different regions of the body modelled as geomtric shapes.

Qtotal is made up of two components: qrad + qconv

What is qabs

qabs(Ta = 20, Tg = NULL, RH = 0.5, E = 0.96, rho = 0.1, cloud = 0, SE = 400)

> [1] 720.2545

compare to a shaded environment with lower SE:

qabs(Ta = 20, Tg = NULL, RH = 0.5, E = 0.96, rho = 0.1, cloud = 0, SE = 100)

> [1] 440.2954

What is qrad

qrad(Ts = 27, Ta = 20, Tg = NULL, RH = 0.5, E = 0.96, rho = 0.1, cloud = 0, SE = 100)

> [1] -1.486309

Ground temperature

qrad(Ts = 30, Ta = 25, Tg = 28, RH = 0.5, E = 0.96, rho = 0.1, cloud = 0, SE = 100)

> [1] 14.29534

What is hconv

What is qconv

qconv(Ts = 30, Ta = 20, V = 1, L = 0.1, type = "forced", shape="hcylinder")
> [1] -102.6565

qconv(Ts = 30, Ta = 20, V = 1, L = 0.1, type = "forced", shape="hplate")
> [1] -124.323

qconv(Ts = 30, Ta = 20, V = 1, L = 0.1, type = "forced", shape="sphere")
> [1] -186.3256

Which is higher: convection or radiation?

qconv(Ts = 30, Ta = 20, V = 0.1, L = 0.1, type = "forced", shape="hcylinder")

> [1] -32.58495
qrad(Ts = 30, Ta = 20, Tg = NULL, RH = 0.5, E = 0.96, rho = 0.1, cloud = 0, SE = 0)

> [1] -112.6542
qconv(Ts = 30, Ta = 20, V = 0.1, L = 0.01, type = "forced", shape="hcylinder")

> [1] -111.4338

qrad(Ts = 30, Ta = 20, Tg = NULL, RH = 0.5, E = 0.96, rho = 0.1, cloud = 0, SE = 0)

> [1] -112.6542

Sample Calculations

qrad.A<-with(d, qrad(Ts, Ta, Tg, RH, E=0.96, rho, cloud, SE))
qconv.free.A<-with(d, qconv(Ts, Ta, V, L, c, n, a, b, m, type="free", shape))
qconv.forced.A<-with(d, qconv(Ts, Ta, V, L,  c, n, a, b, m, type, shape))

qtotal<-A*(qrad.A + qconv.forced.A)

d<-data.frame(d, qrad=qrad.A*A, qconv=qconv.forced.A*A, qtotal=qtotal)
head(d[,-c(3:17)]) # remove some of the extraneous constants for display purposes:

>          Ta        Ts      qrad      qconv    qtotal
> 1 28.322047 34.426894  97.84887 -24.827576  73.02129
> 2 19.295451 23.458105 105.43127 -17.107846  88.32342
> 3 23.640834 26.932211 114.29364 -13.456416 100.83722
> 4  6.971665  8.822035 133.51048  -7.733644 125.77684
> 5 32.594745 39.277282 141.46590 -27.054586 114.41132
> 6 22.613530 28.058783 129.24793 -22.289189 106.95875

A<-0.0097169
L<-0.0587
Ta<-10
Tg<-Ta
Ts<-15.59
SE<-0
rho<-0.1
E<-0.96
RH<-0.5
cloud<-1
V<-5
type="forced"
shape="hcylinder"
(qrad.t<-qrad(Ts=Ts, Ta=Ta, Tg=Tg, RH=RH, E=E, rho=rho, cloud=cloud, SE=SE))

> [1] -37.90549
(qrad.t<-qrad(Ts=Ts, Ta=Ta, Tg=Tg, RH=RH, E=E, rho=rho, cloud=0, SE=SE))

> [1] -83.03191
(qconv.forced.t<-qconv(Ts, Ta, V, L, type=type, shape=shape))

> [1] -192.7086
(qtotal.t<-(qrad.t + qconv.forced.t))

> [1] -230.6141
qtotal.t*A

> [1] -2.240854
qrad(Ts=Ts, Ta=Ta, Tg=Tg, RH=0.5, E=0.96, rho=rho, cloud=1, SE=0) + qconv(Ts, Ta, V=0, L, type="free", shape=shape)

> [1] -64.83292

Estimating Operative Temperature from Environmental Data

Operative temperature with varying reflectances:

Ts<-40
Ta<-30
SE<-seq(0,1100,100)
Toperative<-NULL
for(rho in seq(0, 1, 0.1)){
  temp<-Te(Ts=Ts, Ta=Ta, Tg=NULL, RH=0.5, E=0.96, rho=rho, cloud=1, SE=SE, V=1, 
           L=0.1, type="forced", shape="hcylinder")
  Toperative<-cbind(Toperative, temp)
}
rho<-seq(0, 1, 0.1)
Toperative<-data.frame(SE=seq(0,1100,100), Toperative)
colnames(Toperative)<-c("SE", seq(0,1,0.1))
matplot(Toperative$SE, Toperative[,-1], ylim=c(25, 50), type="l", xlim=c(0,1000),
        main="Effects of Altering Reflectance from 0 to 1",
        ylab="Operative Temperature (°C)", xlab="Solar Radiation (W/m2)", lty=1,
        col=flirpal[rev(seq(1,380,35))])
for(i in 2:12){
    ymax<-par()$yaxp[2]
    xmax<-par()$xaxp[2]  
    x<-Toperative[,1]; y<-Toperative[,i]
    lm1<-lm(y~x)
    b<-coefficients(lm1)[1]; m<-coefficients(lm1)[2]
    if(max(y)>ymax) {xpos<-(ymax-b)/m; ypos<-ymax}
    if(max(y)<ymax) {xpos<-xmax; ypos<-y[which(x==1000)]}
    text(xpos, ypos, labels=rho[(i-1)])
}

Modelling Operative Temperature as a function of Reflectance



Try the Thermimage package in your browser

Any scripts or data that you put into this service are public.

Thermimage documentation built on May 29, 2017, 6:47 p.m.