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)
im <- load.example('parrots') %>% imresize(.5)
#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
#Optionally, add borders
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: