## Art, science, and fun of image analysis in R

A long flight and I’ve already finished reading my Star Trek novel, so I turn to **R**. I found myself in need of a graphic for use as a background image in a web application I was coding in **R** that is used in the context of wildfire modeling. Naturally, I thought to use a picture of flames. They are all over the internet (statisticians like populations that are easy to sample).

But okay, wow, really? Who hasn’t littered a webpage with flames before? This is reminding me of Boyd’s Toast Store. This could turn out one of two ways, nerdy and hokey or nerdy and awesome. It was important to me that the image convey analytic prowess, methodological innovation, computational technique, and scientific investigation. Because that’s how we roll at SNAP. And because for me, sitting in front of an **R** terminal leaves me subject to the ‘one more turn’ syndrome. So I had to up the nerd a notch, in a way that would require as much imagination as skill, which is important for the ‘fun’ factor. Here I will discuss how I go from the first picture below to the second.

## Goals and ideas

But first a quick note on some limitations. This is just a fun example. I did no survey whatsoever of existing packages on image analysis, pattern recognition, or anything similar available in **R**. If you are looking for such, do keep looking, because I did not look, and this post is not intended in any way to represent image analysis in general using **R**. I do not expect my code to even generalize well to other images. This is really a case of making a nice graphic for a particular purpose, which I’ve written about here to provide a sense of how I sometimes think about specific data visualization tasks.

Clearly from the outset this is about achieving an aesthetically pleasing graphic, not the goal of most projects in their own right. Usually I tell people that a good plot is one that draws attention to the underlying data, not to itself. But here we have an exception. I just want this graphic to look good! So rather than merely adding an image of flames to the background of my web application, I decided that I would load the image into **R** and first do some image processing on it.

My goal was to accurately map semi-transparent overlays of common statistical plots and graphs to the flames in the image in a way suggestive of pattern recognition. I also felt like doing something somewhat from scratch, by which I mean, without the use of **R** packages made specifically for image analysis, e.g., `ripa`

. It’s puzzle-solving time!

The easiest thing for me to imagine at first was a semi-transparent overlay of histogram bars that tracked the heights of the flames in the image. This was followed by the idea of a gridded heat map, on a somewhat coarser scale (an aggregation by some factor of the rows and columns of the image matrix). I decided the heat map would be the first layer on top of the flames image, and then the histogram bars, and both would be semi-transparent. Later I decided to add opaque smoothing splines and other statistically derived plotting features that tracked the pattern of flame peaks in the image.

This was all quite fitting (forgive the pun), as one statistic we routinely compute at SNAP with our spatial modeling of wildfire is the area burned on the landscape, and although this image matrix of flames has a two-dimensional space with a vertical orientation orthogonal to such a landscape, I quite literally found myself estimating the boundaries of the area burned.

The heat map easy. It also servers as a backdrop to other more complex plotting elements. So let’s begin there. For this all I have to do is choose a level of aggregation. I convert the matrix to a raster layer using the `raster`

package, do my aggregation by mean, and obtain a coarser, more blocky image that of course corresponds to the flames image from which it is derived. I simply aggregate it to comical levels so it doesn’t look like fire anymore and plot it over top of the original with semi-transparency using a similar color palette (but not quite `heat.colors`

).

## Take time to think about the problem

So how did I implement everything I saw in my head? As usual, the gap between imagination and reality was bridged by **R** code. And the gap between my imagination and my **R** code was bridged by statistics and mathematics. Look at the original image. What do you see? I see a curve that rides along the tops of the flames.

Actually, I don’t see it. I infer it. It’s an abstraction. The human brain is very good at pattern recognition. In fact we easily see patterns even when they do not exist. And thank goodness because combating this failing, when it’s a failing, is partly why statisticians exist. Humans have a very high rate of false positives; high sensitivity and low specificity if you will. We’ll find a pattern in anything. I’m not an expert on computers by any means, but I know computers are not nearly as good at seeing patterns as we are. But regarding using software to analyze this image, we’re not quite there yet. For the moment we’re still on our ability to see patterns. And this inferred curve is very easy and completely natural for us to see abstractly.

But the neat thing is even though I can’t see any explicit curve on the graphic, this abstract idea allows me to alter how I see what is actually in the image, the flames themselves. The curve is a mathematical function in the most basic sense. Imagine this image has standard `x`

and `y`

axes and `y = f(x)`

. When I say it is a function, the idea I am thinking of is that each unique value of `x`

is mapped to a unique value of `y`

whereas the reverse need not be true. If I turned this image of flames on its side, the inferred curve would no longer be a valid mathematical function in this sense, at least not unless I then approached it from the side in my code.

Next let’s think about calculus. In the image, the flames represent the area under the curve. Although the curve is not directly observable, I can observe the area under it, albeit not perfectly. I don’t have the function, but I have observations which describe the function’s integral. The calculus speak is just for conceptualization of course. I do not possess any kind of deterministic mathematical function, equation, or otherwise closed-form expression of which I could take a derivative, plot it, and be done with it; this is not math class.

## Framing the problem frames the solutions

So now let’s think about this problem in a more statistical framework. I have a sample of observations of the integral, or area under the curve, of the function which describes this curve I seek to fit over top of the flames. And it is a very noisy sample at that. If I pretend the general curve that fits over the flames well (in the ‘Goldilocks zone’, neither underfitting nor overfitting) is the ‘truth’, then the observed sample clearly is noisy. In fact I can treat the observations as having substantial measurement error, in the sense that to overfit a very flexible curve perfectly to the flame peaks would require using a curve that is not a valid function. When a lick of flame curls around or doubles back I have to ignore that. The trick of implementation will be in using this noisy sample of observations of the area under the curve to estimate the “average” curve.

Keep in mind I don’t quite want to plot the curve itself yet. My initial goals are the heat map and histogram bars. But if I estimate a smooth function that describes this curve, then I will have an estimate of `y`

for each value of `x`

. I will effectively have estimated the heights of the flames across the columns of the image matrix and get the histogram for free. I could do something like draw bars from zero up to those heights, but that’s not quite what I do. I have a sample. I intend to do some kind of estimation. Let’s continue thinking of this as a statistical problem, specifically a sampling problem.

I decided to sample column indices `1:n`

. Look at the bottom of the image. The last row has fire in every column. It varies in intensity, but is uniformly present enough and intense enough across all cells in the bottom row that the values are not very different from cell to cell. If I were to sample with replacement a large number of these column indices with probability proportional to the intensity of the cell values, I would expect to draw similar frequencies of all column indices. However, doing the same sampling of column indices in the top row with probability proportional to cell values would result in sampling largely from a small number of column indices corresponding to the high flame peak near the right side of the image and relatively no samples of other column indices being drawn.

I do this sampling for each row and then pool my samples of column indices, resulting in a frequency table of how many times each column number, `1:n`

, was drawn. This allows me to make a semi-transparent histogram overlay for the flames in the spirit of pattern recognition. My actual technique was a bit more complex, involving weighting the probabilities of sampling each cell in a row by more than strictly their cell values, as well as altering this weighting scheme monotonically across the rows. But you have the general idea for now. You can get specifics from the **R** code.

Note that I do not actually use `hist`

to make a histogram. Instead I use `barplot`

, but with no spacing between bars. They are not equivalent, and a straight up histogram plot would not give the desired result. Also, for the bar plot, I have to rescale the bars to `[0,1]`

, the `ylim`

of the `image`

plot of the flames matrix, rather than plot the values from the frequency table obtained by the sampling procedure because any count beyond one would be completely off the scale of the plot.

Next comes `lines`

. But I don’t like the lines, too jagged. I try `loess`

but it doesn’t work well with so few `x`

values either. But then I go with smoothing splines and am very satisfied.

Figuring I’m on a roll, I decide to plot more smoothing splines. I want one to suggest a sort of upper bound for the hypothetical curve I am estimating. An interval estimate trumps a point estimate, so I would rather have something akin to confidence bands than a single line that is misleadingly suggestive to the casual observer as implying high precision.

## Turn the problem inside out

Specifically in the context of this image of fire, I think about how rapidly flames flicker and move, shrink and grow, from the human perception of time. It just seems natural to have a bounded interval estimate of the curve using bands rather than a single line. I make the second line by the same methods as the first, with a couple tweaks. As with the first sampling technique, for each row I sample with probability somewhat proportional to cell value, but this time inversely proportional, because I want to estimate the boundary of the darkness, not the fire. I also tweak the weighting schemes so that the two estimated curves do not completely coincide. One needs to be more tuned to the darkness and one to the light so that they will estimate different curves and leave a gap of uncertainty or margin of error between them. I make a third line to describe the mean by averaging the two outer curves.

## Poised to plot

The hard work has been done and I am in a position to draw just about anything I want on the image with respect to the flames. At this point I’m just getting fancy. After all, the goal is largely aesthetics.

It’s also important to know when to stop. I thought of placing an upside down histogram that begins at the top of the image, mirroring the bottom histogram somewhat and terminating at the upper smoothing spline curve. But this would be too much visual clutter, especially with consideration of using my graphic as a background image for my web application. I know there is more room for visual noise at the bottom of the app screen and less at the top where there will tend to be text present. In fact this is another reason I chose this particular image of fire. It occurs predominantly at the bottom of the image and fades to black as it nears the top, leaving room for text.

I go with line segments where the upside down histogram bars would have terminated.

This works well, but something in between the two relatively extreme options of histogram bars and flat horizontal lines may look better. So I add points (not vertical line segments) above these horizontal lines. They look like line segments of varying length but actually they are points using `pch="|"`

with randomly varying `cex`

values using `runif`

. I also use `runif`

to sample the `x`

values of `(x,y)`

pairs and `rnorm`

to sample the `y`

values. I sample the latter with mean zero, then take the absolute value, and then add the `y`

value of the estimated curve for the given `x`

. And there it is; the final plot.

I could have sampled from a distribution of a random variable that only takes non-negative values, but I was too lazy to use a log-normal and didn’t feel like thinking about parameters for gamma or exponential, etc. Sometimes it’s just nice to use `rnorm`

and not get needlessly complicated. As if this isn’t complicated enough already. For one thing, thinking statistically is helpful, but at this stage I am essentially doing art, not math, so the details are not critical.

I also scale the opacity of these points based on their position in the vertical range from the top of the image to the lowest reach among the points making the upper bound curve. This ensures they disappear fast enough to be out of the way of any text I add in the app. They also add a nice touch without adding too much. It wasn’t my goal to make it look like a battle between fire and rain, that’s just inadvertently how it turned out. And it seems apropos.

Below are comparisons of the original, the final with and without the original in the background, and a screenshot from the web application which uses this plot as a background image.