We assume the user has R installed on their computer and has a basic familiarity with R syntax, base functions, and coding in general. To install the imagefx package from the R console use install.packages('imagefx') and then load with:
Read in image data
Data files can have various formats including jpg, png, tiff, etc… and it is important that appropriate packages are installed to read in the desired data set. For example, if you are trying to read in images with the extension .jpeg, you will need to install the jpeg package using
Image data should be saved within a well organized directory. This will ease reading in and analyzing sequences of images (i.e. video), which is discussed later. In the case of webcam data, raw image frames may be located in /name_of_camera/year/month/day/hour/. Please keep in mind the difference in directory nomenclature between Apple, Windows, and Linux operating systems.
The below code is an example of how to read in an image file (in this case a .jpeg) once the jpeg package is installed.
Once the data are loaded in the R workspace, you can get a sense for the structure and the image itself with:
In this case the dimensions of the image correspond to the number of rows, number of columns, and number of color channels, respectively. In other words, this image has 300 rows, 400 columns and a red, green, and blue color channel. If the image your read in is gray scaled, it is likely the image will be read in as a matrix with only rows and columns.
Note the origin in the image is the bottom left corner and that the 400 rows correspond to the width of the image and the 300 columns correspond to the height of the image. This is an important, and sometimes confusing, aspect of how R reads in image data and should be considered when cropping and further processing the image.
Preprocess image
Before analyzing volcano image data, preprocessing is helpful, and sometimes necessary, to highlight the plume region. To this end, we first crop the image and then remove an color trends not associated with volcanic activity.
Crop region above the vent
We are interested in monitoring the plume activity above the crater rim and therefore must crop the image in one of two ways. If you do not already know the bottom left and top right corners of the crop region you can pick them manually by supplying the name of the image as the single input into the crop.image function:
The image will be plotted and the user must first select the bottom left corner then, the top right corner. The crop area will be indicated by a dashed rectangle and if deemed sufficient, the user may select crop in the top right margin of the plot. If the crop region needs to be modified, the user may click repick in the top left margin and repeat the picking process.
If you know the bottom left and top right corners of the crop region, you may crop the image automatically as such:
The output from crop.image is a list of length two whose first component is the cropped image and second is a vector indicating the locations of the crop region’s bottom left and top right corners in the original image. We can disregard the crop locations and keep the crop image and plot the results with:
Remove trends
Light sources outside the frame (e.g. the Sun or moon) can introduce color gradients in the image. These trends may detract from any plume activity and should be removed. We accomplish this by subtracting the mean and best fit plane from each color channel in the image.
Unless the image is gray scaled (i.e. one color channel) each RGB channel must be processed separately as such:
## separate the image into its RGB color channels
img.r <- img.crop[,,1]
img.g <- img.crop[,,2]
img.b <- img.crop[,,3]
## subtract the mean from each color channel matrix
img.r.dmean <- img.r - mean(img.r)
img.g.dmean <- img.g - mean(img.g)
img.b.dmean <- img.b - mean(img.b)
## Find the best fit plane (trend) in each color channel matrix
img.r.trend <- fit3d(img.r.dmean)
img.g.trend <- fit3d(img.g.dmean)
img.b.trend <- fit3d(img.b.dmean)
## subtract the fitted plane from each color channel (i.e. detrend)
img.r.dtrend <- img.r.dmean - img.r.trend
img.g.dtrend <- img.g.dmean - img.g.trend
img.b.dtrend <- img.b.dmean - img.b.trend
## plot the red channel detrend and original
plot.new()
split.screen(c(1,2))
#> [1] 3 4
screen(1)
image2(img.r,asp=1,xlab="",ylab="",main="Red Channel")
screen(2)
image2(img.r.dtrend,asp=1,xlab="",ylab="",main="Red Channel Detrended")
Notice in the above detrended image the colors in the top left and bottom right regions are both light gray. However in the original image, the top left hand region is significantly darker than the bottom right hand region. Although the difference may seem insignificant, the following blob detection algorithm often fails without this preprocessing.
Run algorithm on series of images
As mentioned earlier, a well organized directory is critical for proper time series analysis of images. Directories in which images are stored should be named in a way so that files are sorted with the oldest listed first and the newest last. For example, file names may be specific and include the hour, min, and second in the file name (e.g. 010101.jpeg) or, if the sampling rate is known, as a general numeric sequence with zero padding (e.g. 0001.jpeg)
## name the directory which holds the image files
img.dir <- 'raw_images/year/month/day/hour/'
## list the files in this directory. Note they should be sorted from oldest to newest
img.files <- list.files(img.dir)
## how many statistics will you save.
## this will change if you modify the calc.blob.stats function
num.stats=14
## set up matrices to hold all the blob statistics for each RGB channel
all.blob.stats.r = matrix(NA,nrow=length(img.files),ncol=num.stats)
all.blob.stats.g = matrix(NA,nrow=length(img.files),ncol=num.stats)
all.blob.stats.b = matrix(NA,nrow=length(img.files),ncol=num.stats)
## loop over each image and perform the steps outlined above
i=1
while(i<=length(img.files)) {
## what is the current image file
cur.img.file <- paste(img.dir,img.files[i],sep="")
## read in the current image
cur.img <-readJPEG(cur.img.file)
## Preprocess the image ...
## Perform blob detection ...
## Calculate blob statistics ...
## save the statistics from the current frame to the all blob statistics matrix
all.blob.stats.r[i,] <- blob.stats.r
all.blob.stats.g[i,] <- blob.stats.g
all.blob.stats.b[i,] <- blob.stats.b
## move onto the next file
i=i+1
}
## combine the blob statistis associated with each RGB channel into a list
all.stats <- list(r=all.blob.stats.r, g=all.blob.stats.g, b=all.blob.stats.b)
## save all the statistics in an appropriate directory
save(all.stats,file='blob_stats/year/month/day/hour/blob_stats.RData')
We perform the outlined routine above on a series of image frames collected at Erebus Volcano, Antarctica in order to track eruptive activity through time. These data come preloaded with imagefx and show a single bubble burst from the lava lake.
## pick a single statistic from all the blob statistics
## in this case, the sum of the red channel blob region
blob.sum <- blob.stats$r[,1]
## PLOTTING ##
## set up the plot
close.screen(all.screens=TRUE)
split.screen(c(2,1))
#> [1] 1 2
split.screen(c(1,3),screen=1)
#> [1] 3 4 5
## plot some example images
screen(3)
par(mar=c(0,0,0,0))
image2(erebus.40,axes=FALSE,xlab="",ylab="",asp=1)
text(40,nrow(erebus.40),labels='t1',pos=1,font=2,col='white')
screen(4)
par(mar=c(0,0,0,0))
image2(erebus.70,axes=FALSE,xlab="",ylab="",asp=1)
text(40,nrow(erebus.70),labels='t2',pos=1,font=2,col='white')
screen(5)
par(mar=c(0,0,0,0))
image2(erebus.90,axes=FALSE,xlab="",ylab="",asp=1)
text(40,nrow(erebus.90),labels='t3',pos=1,font=2,col='white')
## plot raw blob stats
screen(2)
plot(blob.sum,type='o',bg='gray80',xlab="",ylab='',
pch=21, cex=0.5, ylim=c(min(blob.sum),max(blob.sum)+500),axes=FALSE)
## add a title
title(main='Blob Sum in Red Channel')
## add some axis information
mtext('normalized amplitude',side=2,line=1)
mtext('time (frames)',side=1,line=3)
axis(1)
## add text and lines indicating where the images occur
abline(v=c(40,70,90),lty=2,col='gray80')
text(x=c(40,70,90),y=rep(max(blob.sum)+400,3),labels=c('t1','t2','t3'),pos=2)
Smooth blob statistics
Raw images are often noisy from sudden changes in acquisition settings (e.g. exposure time), camera shaking, instrument noise, etc… Noise may be smoothed by applying a filter in the time domain.
## pick a single statistic from (i.e. the sum of the red channel blob region)
blob.sum <- blob.stats$r[,1]
## choose the window length for the running average operator
win.length = 15
## smooth the blob sum stat according to window length
blob.sum.filt <- run.avg(blob.sum,win.length)
## PLOTTING ##
plot(blob.sum,type='o',bg='gray80',xlab="time (frame)",
ylab='blob sum in red channel',pch=21,cex=0.5)
lines(blob.sum.filt,col='gray30',lwd=4)
lines(blob.sum.filt,col='red',lwd=2)
## add a legend
legend('topright',legend=c('raw','filtered'),col=c('gray80','red'),
pch=c(21,NA),lty=c(1,1),pt.bg=c('gray30',NA))
The smoothed bob statistics (and to some degree the raw blob statistics) serve as additional time series information to supplement other datasets including seismic and infrasound. Additionally these statistics may serve as feature inputs in machine learning algorithms to classify volcano video data.
