Efficient implementation of summed area table/integral image in R

Question

I am trying to construct a summed area table or integral image given an image matrix. For those of you who dont know what it is, from wikipedia:

A summed area table (also known as an integral image) is a data structure and algorithm for quickly and efficiently generating the sum of values in a rectangular subset of a grid

In other words, its used to sum up values of any rectangular region in the image/matrix in constant time.

I am trying to implement this in R. However, my code seems to take too long to run.

Here is the pseudo code from this link. in is the input matrix or image and intImg is whats returned

for i=0 to w do
   sum←0

   for j=0 to h do
      sum ← sum + in[i, j]

      if i = 0 then
         intImg[i, j] ← sum
      else
         intImg[i, j] ← intImg[i − 1, j] + sum
      end if
   end for 
end for

And here is my implementation

w = ncol(im)
h = nrow(im)
intImg = c(NA)
length(intImg) = w*h

for(i in 1:w){ #x
  sum = 0;
  for(j in 1:h){ #y
    ind = ((j-1)*w)+ (i-1) + 1 #index
    sum = sum + im[ind]
    if(i == 1){
      intImg[ind] = sum
    }else{
      intImg[ind] = intImg[ind-1]+sum
    }
  }
}
intImg = matrix(intImg, h, w, byrow=T)

Example of input and output matrix:

However, on a 480x640 matrix, this takes ~ 4 seconds. In the paper they describe it to take on the order of milliseconds for those dimensions.

Am I doing something inefficient in my loops or indexing?

I considered writing it in C++ and wrapping it in R, but I am not very familiar with C++.

Thank you

Answer 1

You could try to use apply (isn't faster than your for-loops if you pre-allocating the memory):

areaTable <- function(x) {
  return(apply(apply(x, 1, cumsum), 1, cumsum))
}

areaTable(m)
#      [,1] [,2] [,3] [,4]
# [1,]    4    5    7    9
# [2,]    4    9   12   17
# [3,]    7   13   16   25
# [4,]    9   16   22   33