Band pass filtering of n-dimensional arrays

Ana Todorovic had a matlab script to band pass filter some images. A very rudimentary script, and it didn't seem to work. It used some functions that seem missing in R. Also the available image libraries for R aren't that great yet, and the package rimage doesn't even compile out of the box on Mac OS X... That's a bit of a pity because many of the functions in there are compiled C code and thereby aren't available when the package won't install. Here's some alternative R code.

# shift wraps the content of an n-D array half way arround the borders,
# compare with 'fftshift' in matlab
shift <- fftshift <- function (x)
{
d = dim(as.array(x))
idx = sapply(d, function(n) list((1:n + n%/%2)%%n + 1))
do.call("[", c(list(x), idx))
}

slab <- function (v, i) # generalization of 'row' and 'col' functions to n-D arrays
{
v <- as.array(v)
if (i > length(dim(v)))
stop("wrong number of dimensions")
else array(rep(1:dim(v)[i], each = prod(c(1, dim(v))[1:i])),
dim(v))
}

lowpass2d <- function(x, radius){ # low-pass filter 2d images at cutoff frequency radius
x = as.array(x)
m = matrix(0,nrow(x),ncol(x))
m = sqrt((row(m)-nrow(m)/2)^2+(col(m)-ncol(m)/2)^2)
Mod(fft(fft(x) * shift(m < radius) / length(x), inverse=TRUE));
}

lowpass <- # generalizes lowpass2d from 2D to nD
function(x, radius, taper=function(x) (sign(x)+1)/2){
x = as.array(x)
m = sapply(1:length(dim(x)), function(i) (slab(x,i)-dim(x)[i]/2)^2)
m = array(sqrt(rowSums(m)), dim(x))
Mod(fft(fft(x) * shift(taper(radius-m)) / length(x), inverse=TRUE));
}

highpass2d <- function(x, radius){
x = as.array(x)
m = matrix(0,nrow(x),ncol(x))
m = sqrt((row(m)-nrow(m)/2)^2+(col(m)-ncol(m)/2)^2)
Mod(fft(fft(x) * shift(m >= radius) / length(x), inverse=TRUE));
}

highpass <- # generalizes highpass2d from 2D to nD
function(x, radius, taper=function(x) (sign(x)+1)/2){
x = as.array(x)
m = sapply(1:length(dim(x)), function(i) (slab(x,i)-dim(x)[i]/2)^2)
m = array(sqrt(rowSums(m)), dim(x))
Mod(fft(fft(x) * shift(taper(m-radius)) / length(x), inverse=TRUE));
}

bandpass <- # bandpass between lower and upper frequencies
function(x, lower, upper, taper=function(x) (sign(x)+1)/2)
lowpass(highpass(x, lower, taper=taper), upper, taper=taper)

To test these rudimentary filters it useful to generate a grating image.

grating <- # this function is not very satisfactory as you will soon notice
function (nrow, ncol, deg = 0, ncycles = nrow/2, shift = 0)
{
if (!deg %in% c(0, 90, 45, 135))
stop("'deg' should be one of ", c(0, 90, 45, 135))
im = matrix(, nrow, ncol)
im = if (deg == 0)
row(im)
else if (deg == 90)
col(im)
else if (deg == 45)
row(im) - col(im)
else row(im) + col(im)
im = im - min(im)
len = round(2 * ncycles)
s = matrix(seq(min(im), max(im), len = len)[1:ifelse(len%%2 ==
0, len, len - 1)] - shift, 2)
apply(outer(im, s[1, ], ">") & outer(im, s[2, ], "<="), 1:2,
any)
}

Here's an example of filtering an image Is the band-passed image improved? No, of course not, but that was not the point! Here's an example where some noise was added to the original and then partly removed again with a low-pass filter. I suppose this type of filtering is not very useful to clean up degraded images....

cbind, rbind, zbind

I needed to bind a couple of matrices in a three dimensional array. Surprisingly, there was no matrix equivalent of cbind/rbind, which neatly bind vectors and matrices into matrices. So I defined my own: zbind. It is slightly more general, as it binds any number of similar shaped arrays along a new dimension. Here's the code:

function (X, Y, ..., higher.dim = TRUE) 
# Binds arrays X, Y, ... of the same shape in an array
# of one higher dimension: zbind(1:10,15:25) == cbind(1:10,15:25).
# If higher.dim=FALSE the arguments are bind in the last dimension
{
    if (missing(Y)) 
        return(X)
    X = as.array(X)
    Y = as.array(Y)
    if (higher.dim) {
        dim(X) = c(dim(X), 1)
        dim(Y) = c(dim(Y), 1)
    }
    dx = dim(X)
    nx = length(dx)
    dy = dim(Y)
    ny = length(dy)
    idx = 1:(nx - 1)
    if (any(dx[idx] != dy[idx])) 
        stop("Arguments should have the same shape")
    z = c(X, Y)
    zz = array(z, c(dx[-nx], length(z)/prod(dx[-nx])))
    Recall(zz, ..., higher.dim = FALSE)
}