From b56c370dfaf9a3203259ca4ae26959bea9ca9674 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sat, 17 Sep 2011 09:53:08 -0700
Subject: [PATCH] Add summed area table, using Nicolas Pinto's suggestion for a
 pure Python implementation.

---
 scikits/image/transform/__init__.py       |  8 +--
 scikits/image/transform/sum_area_table.py | 61 +++++++++++++++++++++++
 scikits/image/transform/tests/test_sat.py | 28 +++++++++++
 3 files changed, 93 insertions(+), 4 deletions(-)
 create mode 100644 scikits/image/transform/sum_area_table.py
 create mode 100644 scikits/image/transform/tests/test_sat.py

diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py
index 91711c6e..407ab0d4 100644
--- a/scikits/image/transform/__init__.py
+++ b/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 *
diff --git a/scikits/image/transform/sum_area_table.py b/scikits/image/transform/sum_area_table.py
new file mode 100644
index 00000000..6811cdc5
--- /dev/null
+++ b/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
diff --git a/scikits/image/transform/tests/test_sat.py b/scikits/image/transform/tests/test_sat.py
new file mode 100644
index 00000000..481d8e28
--- /dev/null
+++ b/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()