imager contains a large array of functions for working with image data, with most of these functions coming from the CImg library by David Tschumperlé. This vignette is just a short tutorial, you’ll find more information and examples on the website. Each function in the package is documented and comes with examples, so have a look at package documentation as well.

1 Plotting and loading images

imager comes with an example picture of boats. Let’s have a look:

library(imager)
plot(boats)

Note the y axis running downwards: the origin is at the top-left corner, which is the traditional coordinate system for images. imager uses this coordinate system consistently. Image data has class “cimg”:

class(boats)
[1] "cimg"         "imager_array" "numeric"     

and we can get some basic info by typing:

boats
Image. Width: 256 pix Height: 384 pix Depth: 1 Colour channels: 3 

Width and height should be self-explanatory. Depth is how many frames the image has: if depth > 1 then the image is actually a video. Boats has three colour channels, the usual RGB. A grayscale version of boats would have only one:

grayscale(boats)
Image. Width: 256 pix Height: 384 pix Depth: 1 Colour channels: 1 

An object of class cimg is actually just a thin interface over a regular 4D array:

dim(boats)
[1] 256 384   1   3

We’ll see below how images are stored exactly. For most intents and purposes, they behave like regular arrays, meaning the usual arithmetic operations work:

log(boats)+3*sqrt(boats)
Image. Width: 256 pix Height: 384 pix Depth: 1 Colour channels: 3 
mean(boats)
[1] 0.5089061
sd(boats)
[1] 0.144797

Now you might wonder why the following two images look exactly the same:

layout(t(1:2))
plot(boats)
plot(boats/2)

That’s because the plot function automatically rescales the image data so that the whole range of colour values is used. There are two reasons why that’s the default behaviour:

  1. There’s no agreed-upon standard for how RGB values should be scaled. Some software, like CImg, uses a range of 0-255 (dark to light), other, like R’s rgb function, uses a 0-1 range.
  2. Often it’s just more convenient to work with a zero-mean image, which means having negative values.

If you don’t want imager to rescale the colours automatically, set rescale to FALSE, but now imager will want values that are in the \([0,1]\) range.

layout(t(1:2))
plot(boats,rescale=FALSE)
plot(boats/2,rescale=FALSE)

If you’d like tighter control over how imager converts pixel values into colours, you can specify a colour scale. R likes its colours defined as hex codes, like so:

rgb(0,1,0)
[1] "#00FF00"

The function rgb is a colour scale, i.e., it takes pixel values and returns colours. We can define an alternative colour scale that swaps the red and green values:

cscale <- function(r,g,b) rgb(g,r,b)
plot(boats,colourscale=cscale,rescale=FALSE)

In grayscale images pixels have only one value, so that the colour map is simpler: it takes a single value and returns a colour. In the next example we convert the image to grayscale

#Map grayscale values to blue
cscale <- function(v) rgb(0,0,v)
grayscale(boats) %>% plot(colourscale=cscale,rescale=FALSE)

The scales package has a few handy functions for creating colour scales, for example by interpolating a gradient:

cscale <- scales::gradient_n_pal(c("red","purple","lightblue"),c(0,.5,1))
#cscale is now a function returning colour codes
cscale(0)
[1] "#FF0000"
grayscale(boats) %>% plot(colourscale=cscale,rescale=FALSE)

See the documentation for plot.cimg and as.raster.cimg for more information and examples.

The next thing you’ll probably want to be doing is to load an image, which can be done using load.image. imager ships with another example image, which is stored somewhere in your R library. We find out where using system.file

fpath <- system.file('extdata/parrots.png',package='imager')

We’re now ready to load the image:

parrots <- load.image(fpath)
plot(parrots)

imager supports JPEG, PNG, TIFF and BMP natively - for other formats you’ll need to install ImageMagick.

2 Example 1: Histogram equalisation

Histogram equalisation is a textbook example of a contrast-enhancing filter. It’s also a good topic for an introduction to what you can do with imager.

Image histograms are just histogram of pixel values, which are of course pretty easy to obtain in R:

grayscale(boats) %>% hist(main="Luminance values in boats picture")

Since images are stored essentially as arrays, here we’re just using R’s regular hist function, which treats our array as a vector of values. If we wanted to look only at the red channel, we could use:

R(boats) %>% hist(main="Red channel values in boats picture")

#Equivalently:
#channel(boats,1) %>% hist(main="Red channel values in boats picture")

Another approach is to turn the image into a data.frame, and use ggplot to view all channels at once:

library(ggplot2)
bdf <- as.data.frame(boats)
head(bdf,3)
  x y cc     value
1 1 1  1 0.3882353
2 2 1  1 0.3858633
3 3 1  1 0.3849406
bdf <- plyr::mutate(bdf,channel=factor(cc,labels=c('R','G','B')))
ggplot(bdf,aes(value,col=channel))+geom_histogram(bins=30)+facet_wrap(~ channel)

What we immediately see from these histograms is that the middle values are in a sense over-used: there’s very few pixels with high or low values. Histogram equalisation solves the problem by making histograms flat: each pixel’s value is replaced by its rank, which is equivalent to running the data through their empirical cdf.

As an illustration of what this does, see the following example:

x <- rnorm(100)
layout(t(1:2))
hist(x,main="Histogram of x")
f <- ecdf(x)
hist(f(x),main="Histogram of ecdf(x)")

We can apply it directly to images as follows:

boats.g <- grayscale(boats)
f <- ecdf(boats.g)
plot(f,main="Empirical CDF of luminance values")

Again we’re using a standard R function (ecdf), which returns another function corresponding to the ECDF of luminance values in boats.g.

If we run the pixel data back through f we get a flat histogram:

f(boats.g) %>% hist(main="Transformed luminance values")

Now the only problem is that ecdf is base R, and unaware of our cimg objects. The function f took an image and returned a vector:

f(boats.g) %>% str
 num [1:98304] 0.171 0.165 0.163 0.164 0.165 ...

If we wish to get an image back we can just use as.cimg:

f(boats.g) %>% as.cimg(dim=dim(boats.g)) %>% plot(main="With histogram equalisation")

So far we’ve run this on a grayscale image. If we want to do this on RGB data, we need to run the equalisation separately in each channel. imager enables this using its split-apply-combine tricks:

#Hist. equalisation for grayscale
hist.eq <- function(im) as.cimg(ecdf(im)(im),dim=dim(im))

#Split across colour channels, 
cn <- imsplit(boats,"c")
cn #we now have a list of images
Image list of size 3 
cn.eq <- llply(cn,hist.eq) #run hist.eq on each
imappend(cn.eq,"c") %>% plot(main="All channels equalised") #recombine and plot

There’s even a one-liner to do this:

iiply(boats,"c",hist.eq) 
Image. Width: 256 pix Height: 384 pix Depth: 1 Colour channels: 3 

We can use it to check that all channels have been properly normalised:

iiply(boats,"c",hist.eq) %>% as.data.frame %>% ggplot(aes(value))+geom_histogram(bins=30)+facet_wrap(~ cc)

Our trick worked.

3 Example 2: Edge detection

Edge detection relies on image gradients, which imager returns via:

gr <- imgradient(boats.g,"xy")
gr
Image list of size 2 
plot(gr,layout="row")

The object “gr” is an image list, with two components, one for the gradient along \(x\), the other for the gradient along \(y\). “gr” is an object with class “imlist”, which is just a list of images but comes with a few convenience functions (for example, a plotting function as used above).

To be more specific, noting \(I(x,y)\) the image intensity at location \(x,y\), what imager returns is an approximation of: \[ \frac{\partial}{\partial x}I \] in the first panel and: \[ \frac{\partial}{\partial y}I \] in the second.

The magnitude of the gradients thus tell us how fast the image changes around a certain point. Image edges correspond to abrubt changes in the image, and so it’s reasonable to estimate their location based on the norm of the gradient \[ \sqrt{\left(\frac{\partial}{\partial x}I\right)^{2}+\left(\frac{\partial}{\partial y}I\right)^{2}} \]

In imager:

dx <- imgradient(boats.g,"x")
dy <- imgradient(boats.g,"y")
grad.mag <- sqrt(dx^2+dy^2)
plot(grad.mag,main="Gradient magnitude")

Here’s a handy shortcut:

imgradient(boats.g,"xy") %>% enorm %>% plot(main="Gradient magnitude (again)")