1import numpy as np
2
3
4def integral_image(image, *, dtype=None):
5    r"""Integral image / summed area table.
6
7    The integral image contains the sum of all elements above and to the
8    left of it, i.e.:
9
10    .. math::
11
12       S[m, n] = \sum_{i \leq m} \sum_{j \leq n} X[i, j]
13
14    Parameters
15    ----------
16    image : ndarray
17        Input image.
18
19    Returns
20    -------
21    S : ndarray
22        Integral image/summed area table of same shape as input image.
23
24    Notes
25    -----
26    For better accuracy and to avoid potential overflow, the data type of the
27    output may differ from the input's when the default dtype of None is used.
28    For inputs with integer dtype, the behavior matches that for
29    :func:`numpy.cumsum`. Floating point inputs will be promoted to at least
30    double precision. The user can set `dtype` to override this behavior.
31
32    References
33    ----------
34    .. [1] F.C. Crow, "Summed-area tables for texture mapping,"
35           ACM SIGGRAPH Computer Graphics, vol. 18, 1984, pp. 207-212.
36
37    """
38    if dtype is None and image.real.dtype.kind == 'f':
39        # default to at least double precision cumsum for accuracy
40        dtype = np.promote_types(image.dtype, np.float64)
41
42    S = image
43    for i in range(image.ndim):
44        S = S.cumsum(axis=i, dtype=dtype)
45    return S
46
47
48def integrate(ii, start, end):
49    """Use an integral image to integrate over a given window.
50
51    Parameters
52    ----------
53    ii : ndarray
54        Integral image.
55    start : List of tuples, each tuple of length equal to dimension of `ii`
56        Coordinates of top left corner of window(s).
57        Each tuple in the list contains the starting row, col, ... index
58        i.e `[(row_win1, col_win1, ...), (row_win2, col_win2,...), ...]`.
59    end : List of tuples, each tuple of length equal to dimension of `ii`
60        Coordinates of bottom right corner of window(s).
61        Each tuple in the list containing the end row, col, ... index i.e
62        `[(row_win1, col_win1, ...), (row_win2, col_win2, ...), ...]`.
63
64    Returns
65    -------
66    S : scalar or ndarray
67        Integral (sum) over the given window(s).
68
69
70    Examples
71    --------
72    >>> arr = np.ones((5, 6), dtype=float)
73    >>> ii = integral_image(arr)
74    >>> integrate(ii, (1, 0), (1, 2))  # sum from (1, 0) to (1, 2)
75    array([3.])
76    >>> integrate(ii, [(3, 3)], [(4, 5)])  # sum from (3, 3) to (4, 5)
77    array([6.])
78    >>> # sum from (1, 0) to (1, 2) and from (3, 3) to (4, 5)
79    >>> integrate(ii, [(1, 0), (3, 3)], [(1, 2), (4, 5)])
80    array([3., 6.])
81    """
82    start = np.atleast_2d(np.array(start))
83    end = np.atleast_2d(np.array(end))
84    rows = start.shape[0]
85
86    total_shape = ii.shape
87    total_shape = np.tile(total_shape, [rows, 1])
88
89    # convert negative indices into equivalent positive indices
90    start_negatives = start < 0
91    end_negatives = end < 0
92    start = (start + total_shape) * start_negatives + \
93             start * ~(start_negatives)
94    end = (end + total_shape) * end_negatives + \
95           end * ~(end_negatives)
96
97    if np.any((end - start) < 0):
98        raise IndexError('end coordinates must be greater or equal to start')
99
100    # bit_perm is the total number of terms in the expression
101    # of S. For example, in the case of a 4x4 2D image
102    # sum of image from (1,1) to (2,2) is given by
103    # S = + ii[2, 2]
104    #     - ii[0, 2] - ii[2, 0]
105    #     + ii[0, 0]
106    # The total terms = 4 = 2 ** 2(dims)
107
108    S = np.zeros(rows)
109    bit_perm = 2 ** ii.ndim
110    width = len(bin(bit_perm - 1)[2:])
111
112    # Sum of a (hyper)cube, from an integral image is computed using
113    # values at the corners of the cube. The corners of cube are
114    # selected using binary numbers as described in the following example.
115    # In a 3D cube there are 8 corners. The corners are selected using
116    # binary numbers 000 to 111. Each number is called a permutation, where
117    # perm(000) means, select end corner where none of the coordinates
118    # is replaced, i.e ii[end_row, end_col, end_depth]. Similarly, perm(001)
119    # means replace last coordinate by start - 1, i.e
120    # ii[end_row, end_col, start_depth - 1], and so on.
121    # Sign of even permutations is positive, while those of odd is negative.
122    # If 'start_coord - 1' is -ve it is labeled bad and not considered in
123    # the final sum.
124
125    for i in range(bit_perm):  # for all permutations
126        # boolean permutation array eg [True, False] for '10'
127        binary = bin(i)[2:].zfill(width)
128        bool_mask = [bit == '1' for bit in binary]
129
130        sign = (-1)**sum(bool_mask)  # determine sign of permutation
131
132        bad = [np.any(((start[r] - 1) * bool_mask) < 0)
133               for r in range(rows)]  # find out bad start rows
134
135        corner_points = (end * (np.invert(bool_mask))) + \
136                         ((start - 1) * bool_mask)  # find corner for each row
137
138        S += [sign * ii[tuple(corner_points[r])] if(not bad[r]) else 0
139              for r in range(rows)]  # add only good rows
140    return S
141