## New features in imager 0.20

imager, an R package for image processing, has been updated to v0.20 on CRAN. It’s a major upgrade with a lot of new features, better documentation and a more consistent API.

imager now has 130 functions, and I myself keep forgetting all that’s in there. I’ve added a tutorial vignette that should help you get started. It goes through a few basic tasks like plotting and histogram equalisation and builds up to a multi-scale blob detector. It also covers plotting with ggplot2 and has a thematic list of functions.

New features added in the last months include new assignment functions, a utility for getting information on image files (iminfo), auto-thresholding based on k-means, much better array subset operators, updated docs and a reorganised codebase. Windows support should also have improved. Last but not least, you can now interrupt lengthy computations by hitting Ctrl+c or the stop button in RStudio.

imager now has some easy-to-use replacement functions, meaning you can now do set image channels or change frames using a convenient R-like syntax:

library(imager)
boats.cp = boats #Make a copy of the boats image
R(boats.cp) = 0 #Set red channel to 0
G(boats.cp) =  0 #Set blue channel to 0
plot(boats.cp,main="Just the blue channel")
R(boats.cp) = G(boats)
G(boats.cp) = R(boats)
plot(boats.cp,main="Swapping red and green channels")


see ?imager.replace for more.

Auto-thresholding finds an optimal threshold for converting an image to binary values, based on k-means (it’s essentially a variant of Otsu’s method).
Here’s an illustration on a sketch by Thomas Gainsborough:

url = "https://upload.wikimedia.org/wikipedia/commons/thumb/3/30/Study_of_willows_by_Thomas_Gainsborough.jpg/375px-Study_of_willows_by_Thomas_Gainsborough.jpg"
layout(t(1:2))
plot(im)
grayscale(im) %>% threshold %>% plot


Point-wise reductions are useful for combining a list of images into a single output image. For example, enorm(list(A,B,C)) computes $\sqrt{A^2+B^2+C^2}$, ie. the Euclidean norm. Here’s how you can use it to compute gradient magnitude:

imgradient(im,"xy") %>% enorm %>% plot("Gradient magnitude")


A note on compiling imager: if for some reason R tries to install imager from source (Linux or Mac), you will need the fftw library. On a Mac the easiest way is to grab it via Homebrew (“brew install fftw”), in Ubuntu “sudo apt-get install libfftw3-dev” should do it.

## New R package for Eyelink eye-trackers

Eyelink eye-trackers output an avalanche of disorganised crap. I’ve written an R package that will hopefully filter that crap for you. It’s called eyelinker and it’s on Github.

It outputs a set of dataframes containing raw traces, saccades, fixations and blinks, meaning it’s easy to produce plots like this one:

There’s a vignette explaining everything, just hit vignette(“basics”,package=”eyelinker”).

I’ve tested it on some of our local datasets but given the relatively free-form nature of Eyelink asc files, there’s no guarantee it will work everywhere.

Bug reports are welcome on the github issues page.

## imager now on CRAN, and a non-linear filtering example

imager is an R package for image processing that’s fairly fast and now quite powerful (if I may say so myself). It wraps a neat C++ library called CImg, by David Tschumperlé (CNRS). It took quite a bit of work, but imager is now on CRAN, so that installing it is as easy as:

install.packages("imager")


Here’s an example of using imager for max-filtering. A max-filter replaces each pixel value with the maximum value of its neighbours. Usually you’d write a loop, but we want to do things the R way (warning: this only works for small neighbourhoods, for reasons that will become obvious).

library(imager)
nhood <- expand.grid(dx=-2:2,dy=-2:2) #We want to include all pixels in a square 5x5 neighbourhood
im.s <- alply(nhood,1,function(d) imshift(im,d$dx,d$dy))


The result is a list of shifted versions of the image im (there are 5×5 = 25 different shifts, so you can imagine that you wouldn’t want to do this with a 20×20 neighbourhood and a large image!)
Now running a max (or min) filter is just a matter of calling pmax (or pmin):

max.filt <- do.call(pmax,im.s)
min.filt <- do.call(pmin,im.s)


Here’s our max-filtered image:

and here’s the min-filtered one:

## New package for image processing in R

I’ve written a package for image processing in R, with the goal of providing a fast API in R that lets you do things in C++ if you need to. The package is called imager, and it’ on Github.
The whole thing is based on CImg, a very nice C++ library for image processing by David Tschumperlé.

Features:

• Handles images in up to 4 dimensions, meaning you can use it for volumetric/hyperspectral/data or short videos
• Facilities for taking subsets of images, pixel neighbourhoods, etc.
• All the usual image processing stuff (filters, morphology, transformations, interpolation, etc.)

The package is still in an early phase but it can already do a lot of useful things as you’ll see from the documentation.

Example code:

library(imager)
layout(t(1:3))
plot(im,main="Original image")


## Neurostats 2014 Highlights

Last week the Neurostats 2014 workshop took place at the University of Warwick (co-organised by Adam Johansen, Nicolas Chopin, and myself). The goal was to put some neuroscientists and statisticians together to talk about neural data and what to do with it. General impressions:

• The type of Bayesian hierarchical modelling that Andrew Gelman has been advocating for years is starting to see some use in neuroimaging. On the one hand it makes plenty of sense since the data at the level of individual subjects can be cr*p and so one could really use a bit of clever pooling. On the other, imaging data is very high-dimensional, running a Gibbs sampler can take days, and it’s not easy making the data comparable across subjects.
• You have to know your signals. Neural data can be unbelievably complicated and details matter a lot, as Jonathan Victor showed in his talk. A consequence if that if you as a neuroscientist have a data analysis problem, it’s not enough to go see a statistician and ask for advice. If you have EEG data you need to find someone who knows *specifically* about all the traps and pitfalls of EEG, or else someone who’s willing to learn about these things. A consequence is that we should think about training neurostatisticians, the way we already have biostatisticians, econometricians and psychometricians.

There were plenty of interesting talks, but below are some of my personal highlights.

## Poisson transform – update

Michael Gutmann (University of Helsinki) recently wrote me with some comments on the Poisson transform paper (here). It turns out that the Poisson likelihood we define in the paper is a special case of more general frameworks he has worked on, the most recent being:
M.U. Gutmann and J.Hirayama (2011). Bregman Divergence as General Framework to Estimate Unnormalized Statistical Models,UAI.
available at arxiv.org/abs/1202.3727.

The paper gives a very general (and interesting) framework for estimation using divergences between the empirical distribution of the data and a theoretical model that is not necessarily normalised.
What we call the Poisson transform appears when taking $\Psi(x) = x\log x$ as the generating function for the Bregman divergence. The same choice of Bregman divergence also corresponds to the generalised Kullback-Leibler divergence used in Minka (2005) Divergence measures and message passing. Presumably there are other connections we hadn’t seen either.

Michael also points out the following paper by Mnih & Teh (ICML 2012), who use noise-contrastive learning in a sequential unnormalised model: http://arxiv.org/abs/1206.6426. They ignore normalisation constants, which I wouldn’t recommend as a general strategy (it generally leads to biased estimates). See our paper for a solution that uses semiparametric inference.

## Statistical Challenges in Neuroscience

A workshop on statistics and neuroscience, to take place at the University of Warwick, UK, Sept. 3-5 2014. We’ll talk spikes, voxels, pixels, MCMC, and so on.Official call for posters below the fold.