Add summed area table, using Nicolas Pinto's suggestion for a pure

Python implementation.

scikits/image/transform/__init__.py

@@ -1,4 +1,4 @@
-from hough_transform import *
-from finite_radon_transform import *
-from project import *
-
+from .hough_transform import *
+from .finite_radon_transform import *
+from .project import *
+from .sum_area_table import *

scikits/image/transform/sum_area_table.py

@@ -0,0 +1,61 @@
+def sat(x):
+    """Summed area table / integral image.
+
+    The integral image contains the sum of all elements above and to the
+    left of it, i.e.:
+
+    .. math::
+
+       S[m, n] = \sum_{i \leq m} \sum_{j \leq n} X[i, j]
+
+    Parameters
+    ----------
+    X : ndarray of uint
+        Input image.
+
+    Returns
+    -------
+    S : ndarray
+        Summed area table.
+
+    References
+    ----------
+    .. [1] F.C. Crow, "Summed-area tables for texture mapping,"
+           ACM SIGGRAPH Computer Graphics, vol. 18, 1984, pp. 207-212.
+
+    """
+    return x.cumsum(1).cumsum(0)
+
+def sat_sum(sat, r0, c0, r1, c1):
+    """Using a summed area table / integral image, calculate the sum
+    over a given window.
+
+    Parameters
+    ----------
+    sat : ndarray of uint64
+        Summed area table / integral image.
+    r0, c0 : int
+        Top-left corner of block to be summed.
+    r1, c1 : int
+        Bottom-right corner of block to be summed.
+
+    Returns
+    -------
+    S : int
+        Sum over the given window.
+
+    """
+    S = 0
+
+    S += sat[r1, c1]
+
+    if (r0 - 1 >= 0) and (c0 - 1 >= 0):
+        S += sat[r0 - 1, c0 - 1]
+
+    if (r0 - 1 >= 0):
+        S -= sat[r0 - 1, c1]
+
+    if (c0 - 1 >= 0):
+        S -= sat[r1, c0 - 1]
+
+    return S

scikits/image/transform/tests/test_sat.py

@@ -0,0 +1,28 @@
+import numpy as np
+from numpy.testing import *
+
+from scikits.image.transform import sat, sat_sum
+
+x = (np.random.random((50, 50)) * 255).astype(np.uint8)
+s = sat(x)
+
+def test_validity():
+    y = np.arange(12).reshape((4,3))
+
+    y = (np.random.random((50, 50)) * 255).astype(np.uint8)
+    assert_equal(sat(y)[-1, -1],
+                 y.sum())
+
+def test_basic():
+    assert_equal(x[12:24, 10:20].sum(), sat_sum(s, 12, 10, 23, 19))
+    assert_equal(x[:20, :20].sum(), sat_sum(s, 0, 0, 19, 19))
+    assert_equal(x[:20, 10:20].sum(), sat_sum(s, 0, 10, 19, 19))
+    assert_equal(x[10:20, :20].sum(), sat_sum(s, 10, 0, 19, 19))
+
+def test_single():
+    assert_equal(x[0, 0], sat_sum(s, 0, 0, 0, 0))
+    assert_equal(x[10, 10], sat_sum(s, 10, 10, 10, 10))
+
+
+if __name__ == '__main__':
+    run_module_suite()