integral.py 5.2 KB

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