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