Calculate 3D variant for summed area table using numpy cumsum

Question

In case of a 2D array array.cumsum(0).cumsum(1) gives the Integral image of the array.

What happens if I compute array.cumsum(0).cumsum(1).cumsum(2) over a 3D array?

Do I get a 3D extension of Integral Image i.e, Integral volume over the array?

Its hard to visualize what happens in case of 3D.

I have gone through this discussion. 3D variant for summed area table (SAT)

This gives a recursive way on how to compute the Integral volume. What if I use the cumsum along the 3 axes. Will it give me the same thing?

Will it be more efficient than the recursive method?

Answer 1

Yes, the formula you give, array.cumsum(0).cumsum(1).cumsum(2), will work.

What the formula does is compute a few partial sums so that the sum of these sums is the volume sum. That is, every element needs to be summed exactly once, or, in other words, no element can be skipped and no element counted twice. I think going through each of these questions (is any element skipped or counted twice) is a good way to verify to yourself that this will work. And also run a small test:

x = np.ones((20,20,20)).cumsum(0).cumsum(1).cumsum(2)

print x[2,6,10]   # 231.0
print 3*7*11      # 231

Of course, with all ones there could two errors that cancel each other out, but this wouldn't happen everywhere so it's a reasonable test.

As for efficiency, I would guess that the single pass approach is probably faster, but not by a lot. Also, the above could be sped up using an output array, eg, cumsum(n, out=temp) as otherwise three arrays will be created for this calculation. The best way to know is to test (but only if you need to).