Simon Barthelmé (GIPSA-lab, CNRS)

The following is a implemention in imager of Quads by Michael Fogleman. It illustrates the use of imager, purrr and recursion to work with recursive data structures. A quadtree is a recursive data-structure that subdivides an image into four quadrants, which can themselves be subdivided into sub-quadrants, etc. For compression, we build a quadtree such that image regions that are unhomogeneous are divided more finely. Homogenous image regions are simply assigned their mean colour.

Subdividing can be done using “imsplit”

library(imager)
library(purrr)

#Divide along x, then y
qsplit <- function(im)
{
imsplit(im,"x",2) %>% map(~ imsplit(.,"y",2)) %>%
flatten
}

qsplit(im) %>% as.imlist %>% plot The inverse operation uses “imappend”:

qunsplit <- function(l)
{
list(l[1:2],l[3:4]) %>% map(~ imappend(.,"y")) %>% imappend("x")
}

qsplit(im) %>% qunsplit %>% plot The algorithm works by iterative refinement: at each iteration we find the least homogeneous region, and split it into four quadrants. Homogeneity is measured by the standard deviation of pixel values:

#Max std. dev across channels
imsd <- function(im)
{
imsplit(im,"c") %>% map_dbl(sd) %>% max
}

The refinement function is recursive. It traverses the tree depth-first, looking for the least-homogeneous leaf (leaves are image regions).

refine <- function(l)
{
if (is.cimg(l)) #We have a leaf
{
qs <- qsplit(l) #Split
if (any(dim(l)[1:2] <= 4)) #Quadrants are very small
{
qs$sds <- rep(0,4) #Prevent further refinement } else { qs$sds <- map_dbl(qs,imsd) #Store std.dev of children
}
qs
}
else #Not a leaf, explore further
{
indm <- which.max(l$sds) #Find child with max. std. dev l[[indm]] <- refine(l[[indm]]) #Refine l$sds[indm] <- max(l[[indm]]\$sds) #Update std. dev
l
}
}

Refinement will produce ever-deeper trees. To visualise them, we need to reconstruct an image from the tree. That’s the job of rebuild, which again is recursive:

rebuild <- function(l,borders=FALSE)
{
map(l[-5],~ if (is.cimg(.)) meanim(.,borders=borders) else rebuild(.,borders=borders)) %>% qunsplit
}

#Produce an image that's just the average of image im
meanim <- function(im,borders=FALSE)
{
im <- imsplit(im,"c") %>% map(~ 0*. + mean(.)) %>% imappend("c")
if (borders)
{
im[px.borders(im)] <- 0
}
im
}

The way rebuild works is as follows: look at the current node. If the node is a leaf, approximate it. If the node has children, call rebuild on the children and recombine them using qunsplit.

iter.refine <- function(im,nIter)
{
for (i in seq_len(nIter)) { im <- refine(im) };
im
}

#The first four iterations of the process
map_il(1:4,~ iter.refine(im,.) %>% rebuild) %>% plot #After 200 iterations
iter.refine(im,200) %>% rebuild(borders=T) %>% plot The next image is an animation showing the first 1,000 steps of the process: 