From 0b7b225660b952d9c39ce1c88959f130c38d0ca6 Mon Sep 17 00:00:00 2001
From: Tabish <syedtabish.a@gmail.com>
Date: Sat, 22 Mar 2014 22:02:01 +0530
Subject: [PATCH] Modify 'integral_image' to support nD

---
 skimage/transform/integral.py            | 72 ++++++++++++++++--------
 skimage/transform/tests/test_integral.py | 27 ++-------
 2 files changed, 53 insertions(+), 46 deletions(-)

diff --git a/skimage/transform/integral.py b/skimage/transform/integral.py
index 3503a652..c58ca9ea 100644
--- a/skimage/transform/integral.py
+++ b/skimage/transform/integral.py
@@ -18,8 +18,8 @@ def integral_image(x):
 
     Returns
     -------
-    S : ndarray
-        Integral image / summed area table.
+    S : scalar value
+        summed area table.
 
     References
     ----------
@@ -27,44 +27,66 @@ def integral_image(x):
            ACM SIGGRAPH Computer Graphics, vol. 18, 1984, pp. 207-212.
 
     """
-    return x.cumsum(1).cumsum(0)
+    dim = len(x.shape)
+    S = x
+    for i in range(dim):
+        S = S.cumsum(i)
+    return S
 
 
-def integrate(ii, r0, c0, r1, c1):
+def integrate(ii, start, end):
     """Use an integral image to integrate over a given window.
 
     Parameters
     ----------
     ii : ndarray
         Integral image.
-    r0, c0 : int or ndarray
-        Top-left corner(s) of block to be summed.
-    r1, c1 : int or ndarray
-        Bottom-right corner(s) of block to be summed.
+    start : int or ndarray or list
+        Top-left corner of block to be summed.
+    end  : int or ndarray or list
+        Bottom-right corner of block to be summed.
 
     Returns
     -------
-    S : scalar or ndarray
-        Integral (sum) over the given window(s).
-
+    S : scalar 
+        Integral (sum) over the given window.
+    Notes
+    -----
+    Explination:
+        For a 2D array say(10 x 10) intergral from start=(2,3) to end=(5,6) is
+        #replace 'zero' elements from end -> permutation('00')
+        +Intgral_array[5,6]
+        #replace 'one' elements from end by 'start coorinate - 1' -> permutation('10','01')
+        -(Integral_array[5,(3 - 1)] + integral_array[(2 - 1), 6])
+        #replace 'two' elements from end by 'start coordinate - 1' -> permutation('11')
+        +(Integral_array[(2-1),(3-1)])
     """
-    if np.isscalar(r0):
-        r0, c0, r1, c1 = [np.asarray([x]) for x in (r0, c0, r1, c1)]
+    #make sure start and end both are arrays	
+    start = np.asarray(start)
+    end = np.asarray(end)
 
-    S = np.zeros(r0.shape, ii.dtype)
+    if(np.any(start < 0) or np.any(end < 0)):
+        raise IndexError('cordinates must be non negative')
 
-    S += ii[r1, c1]
+    if(np.any((end - start) < 0)):
+        raise IndexError('end coordinates must be greater or equal to start')
 
-    good = (r0 >= 1) & (c0 >= 1)
-    S[good] += ii[r0[good] - 1, c0[good] - 1]
+    dim = len(ii.shape)  #No. of dimensions of input nd-array 
+    S = 0
+    bit_perm = 2**dim  #bit_perm is the total number of elements in expression of S
+    width = len(bin(bit_perm-1)[2:])
 
-    good = r0 >= 1
-    S[good] -= ii[r0[good] - 1, c1[good]]
-
-    good = c0 >= 1
-    S[good] -= ii[r1[good], c0[good] - 1]
-
-    if S.size == 1:
-        return np.asscalar(S)
+    for i in range(bit_perm):  #for all permutations
+        #generate boolean array corresponding to permutation eg [True, False] for '10'		      
+        binary = bin(i)[2:].zfill(width)
+        bool_mask = [bit == '1' for bit in binary]
+        
+        sign = (-1)**sum(bool_mask)  #determine sign of permutation
+        bad = np.any(((start - 1)*bool_mask) < 0)
+        if(bad):
+            continue
 
+        corner_point = (end * (np.invert(bool_mask))) + ((start - 1) * bool_mask)
+        
+        S += sign*ii[tuple(corner_point)]
     return S
diff --git a/skimage/transform/tests/test_integral.py b/skimage/transform/tests/test_integral.py
index 4a641e71..307728f4 100644
--- a/skimage/transform/tests/test_integral.py
+++ b/skimage/transform/tests/test_integral.py
@@ -16,30 +16,15 @@ def test_validity():
 
 
 def test_basic():
-    assert_equal(x[12:24, 10:20].sum(), integrate(s, 12, 10, 23, 19))
-    assert_equal(x[:20, :20].sum(), integrate(s, 0, 0, 19, 19))
-    assert_equal(x[:20, 10:20].sum(), integrate(s, 0, 10, 19, 19))
-    assert_equal(x[10:20, :20].sum(), integrate(s, 10, 0, 19, 19))
+    assert_equal(x[12:24, 10:20].sum(), integrate(s, [12, 10], [23, 19]))
+    assert_equal(x[:20, :20].sum(), integrate(s, [0, 0], [19, 19]))
+    assert_equal(x[:20, 10:20].sum(), integrate(s, [0, 10], [19, 19]))
+    assert_equal(x[10:20, :20].sum(), integrate(s, [10, 0], [19, 19]))
 
 
 def test_single():
-    assert_equal(x[0, 0], integrate(s, 0, 0, 0, 0))
-    assert_equal(x[10, 10], integrate(s, 10, 10, 10, 10))
-
-def test_vectorized_integrate():
-    r0 = np.array([12, 0, 0, 10, 0, 10, 30])
-    c0 = np.array([10, 0, 10, 0, 0, 10, 31])
-    r1 = np.array([23, 19, 19, 19, 0, 10, 49])
-    c1 = np.array([19, 19, 19, 19, 0, 10, 49])
-
-    expected = np.array([x[12:24, 10:20].sum(),
-                         x[:20, :20].sum(),
-                         x[:20, 10:20].sum(),
-                         x[10:20, :20].sum(),
-                         x[0,0],
-                         x[10, 10],
-                         x[30:, 31:].sum()])
-    assert_equal(expected, integrate(s, r0, c0, r1, c1))
+    assert_equal(x[0, 0], integrate(s, [0, 0], [0, 0]))
+    assert_equal(x[10, 10], integrate(s, [10, 10], [10, 10]))
 
 
 if __name__ == '__main__':