histogram.py 10.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271
  1. # LICENSE HEADER MANAGED BY add-license-header
  2. #
  3. # Copyright 2018 Kornia Team
  4. #
  5. # Licensed under the Apache License, Version 2.0 (the "License");
  6. # you may not use this file except in compliance with the License.
  7. # You may obtain a copy of the License at
  8. #
  9. # http://www.apache.org/licenses/LICENSE-2.0
  10. #
  11. # Unless required by applicable law or agreed to in writing, software
  12. # distributed under the License is distributed on an "AS IS" BASIS,
  13. # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14. # See the License for the specific language governing permissions and
  15. # limitations under the License.
  16. #
  17. from typing import Optional, Tuple
  18. import torch
  19. from kornia.core import Tensor
  20. def marginal_pdf(values: Tensor, bins: Tensor, sigma: Tensor, epsilon: float = 1e-10) -> Tuple[Tensor, Tensor]:
  21. """Calculate the marginal probability distribution function of the input based on the number of histogram bins.
  22. Args:
  23. values: shape [BxNx1].
  24. bins: shape [NUM_BINS].
  25. sigma: shape [1], gaussian smoothing factor.
  26. epsilon: scalar, for numerical stability.
  27. Returns:
  28. Tuple[Tensor, Tensor]:
  29. - Tensor: shape [BxN].
  30. - Tensor: shape [BxNxNUM_BINS].
  31. """
  32. if not isinstance(values, Tensor):
  33. raise TypeError(f"Input values type is not a Tensor. Got {type(values)}")
  34. if not isinstance(bins, Tensor):
  35. raise TypeError(f"Input bins type is not a Tensor. Got {type(bins)}")
  36. if not isinstance(sigma, Tensor):
  37. raise TypeError(f"Input sigma type is not a Tensor. Got {type(sigma)}")
  38. if not values.dim() == 3:
  39. raise ValueError(f"Input values must be a of the shape BxNx1. Got {values.shape}")
  40. if not bins.dim() == 1:
  41. raise ValueError(f"Input bins must be a of the shape NUM_BINS. Got {bins.shape}")
  42. if not sigma.dim() == 0:
  43. raise ValueError(f"Input sigma must be a of the shape 1. Got {sigma.shape}")
  44. residuals = values - bins.unsqueeze(0).unsqueeze(0)
  45. kernel_values = torch.exp(-0.5 * (residuals / sigma).pow(2))
  46. pdf = torch.mean(kernel_values, dim=1)
  47. normalization = torch.sum(pdf, dim=1).unsqueeze(1) + epsilon
  48. pdf = pdf / normalization
  49. return pdf, kernel_values
  50. def joint_pdf(kernel_values1: Tensor, kernel_values2: Tensor, epsilon: float = 1e-10) -> Tensor:
  51. """Calculate the joint probability distribution function of the input tensors based on the number of histogram bins.
  52. Args:
  53. kernel_values1: shape [BxNxNUM_BINS].
  54. kernel_values2: shape [BxNxNUM_BINS].
  55. epsilon: scalar, for numerical stability.
  56. Returns:
  57. shape [BxNUM_BINSxNUM_BINS].
  58. """
  59. if not isinstance(kernel_values1, Tensor):
  60. raise TypeError(f"Input kernel_values1 type is not a Tensor. Got {type(kernel_values1)}")
  61. if not isinstance(kernel_values2, Tensor):
  62. raise TypeError(f"Input kernel_values2 type is not a Tensor. Got {type(kernel_values2)}")
  63. if not kernel_values1.dim() == 3:
  64. raise ValueError(f"Input kernel_values1 must be a of the shape BxN. Got {kernel_values1.shape}")
  65. if not kernel_values2.dim() == 3:
  66. raise ValueError(f"Input kernel_values2 must be a of the shape BxN. Got {kernel_values2.shape}")
  67. if kernel_values1.shape != kernel_values2.shape:
  68. raise ValueError(
  69. "Inputs kernel_values1 and kernel_values2 must have the same shape."
  70. f" Got {kernel_values1.shape} and {kernel_values2.shape}"
  71. )
  72. joint_kernel_values = torch.matmul(kernel_values1.transpose(1, 2), kernel_values2)
  73. normalization = torch.sum(joint_kernel_values, dim=(1, 2)).view(-1, 1, 1) + epsilon
  74. pdf = joint_kernel_values / normalization
  75. return pdf
  76. def histogram(x: Tensor, bins: Tensor, bandwidth: Tensor, epsilon: float = 1e-10) -> Tensor:
  77. """Estimate the histogram of the input tensor.
  78. The calculation uses kernel density estimation which requires a bandwidth (smoothing) parameter.
  79. Args:
  80. x: Input tensor to compute the histogram with shape :math:`(B, D)`.
  81. bins: The number of bins to use the histogram :math:`(N_{bins})`.
  82. bandwidth: Gaussian smoothing factor with shape shape [1].
  83. epsilon: A scalar, for numerical stability.
  84. Returns:
  85. Computed histogram of shape :math:`(B, N_{bins})`.
  86. Examples:
  87. >>> x = torch.rand(1, 10)
  88. >>> bins = torch.torch.linspace(0, 255, 128)
  89. >>> hist = histogram(x, bins, bandwidth=torch.tensor(0.9))
  90. >>> hist.shape
  91. torch.Size([1, 128])
  92. """
  93. pdf, _ = marginal_pdf(x.unsqueeze(2), bins, bandwidth, epsilon)
  94. return pdf
  95. def histogram2d(x1: Tensor, x2: Tensor, bins: Tensor, bandwidth: Tensor, epsilon: float = 1e-10) -> Tensor:
  96. """Estimate the 2d histogram of the input tensor.
  97. The calculation uses kernel density estimation which requires a bandwidth (smoothing) parameter.
  98. Args:
  99. x1: Input tensor to compute the histogram with shape :math:`(B, D1)`.
  100. x2: Input tensor to compute the histogram with shape :math:`(B, D2)`.
  101. bins: The number of bins to use the histogram :math:`(N_{bins})`.
  102. bandwidth: Gaussian smoothing factor with shape shape [1].
  103. epsilon: A scalar, for numerical stability. Default: 1e-10.
  104. Returns:
  105. Computed histogram of shape :math:`(B, N_{bins}), N_{bins})`.
  106. Examples:
  107. >>> x1 = torch.rand(2, 32)
  108. >>> x2 = torch.rand(2, 32)
  109. >>> bins = torch.torch.linspace(0, 255, 128)
  110. >>> hist = histogram2d(x1, x2, bins, bandwidth=torch.tensor(0.9))
  111. >>> hist.shape
  112. torch.Size([2, 128, 128])
  113. """
  114. _, kernel_values1 = marginal_pdf(x1.unsqueeze(2), bins, bandwidth, epsilon)
  115. _, kernel_values2 = marginal_pdf(x2.unsqueeze(2), bins, bandwidth, epsilon)
  116. pdf = joint_pdf(kernel_values1, kernel_values2)
  117. return pdf
  118. def image_histogram2d(
  119. image: Tensor,
  120. min: float = 0.0,
  121. max: float = 255.0,
  122. n_bins: int = 256,
  123. bandwidth: Optional[float] = None,
  124. centers: Optional[Tensor] = None,
  125. return_pdf: bool = False,
  126. kernel: str = "triangular",
  127. eps: float = 1e-10,
  128. ) -> Tuple[Tensor, Tensor]:
  129. """Estimate the histogram of the input image(s).
  130. The calculation uses triangular kernel density estimation.
  131. Args:
  132. image: Input tensor to compute the histogram with shape
  133. :math:`(H, W)`, :math:`(C, H, W)` or :math:`(B, C, H, W)`.
  134. min: Lower end of the interval (inclusive).
  135. max: Upper end of the interval (inclusive). Ignored when
  136. :attr:`centers` is specified.
  137. n_bins: The number of histogram bins. Ignored when
  138. :attr:`centers` is specified.
  139. bandwidth: Smoothing factor. If not specified or equal to -1,
  140. :math:`(bandwidth = (max - min) / n_bins)`.
  141. centers: Centers of the bins with shape :math:`(n_bins,)`.
  142. If not specified or empty, it is calculated as centers of
  143. equal width bins of [min, max] range.
  144. return_pdf: If True, also return probability densities for
  145. each bin.
  146. kernel: kernel to perform kernel density estimation
  147. ``(`triangular`, `gaussian`, `uniform`, `epanechnikov`)``.
  148. eps: epsilon for numerical stability.
  149. Returns:
  150. Computed histogram of shape :math:`(bins)`, :math:`(C, bins)`,
  151. :math:`(B, C, bins)`.
  152. Computed probability densities of shape :math:`(bins)`, :math:`(C, bins)`,
  153. :math:`(B, C, bins)`, if return_pdf is ``True``. Tensor of zeros with shape
  154. of the histogram otherwise.
  155. """
  156. if image is not None and not isinstance(image, Tensor):
  157. raise TypeError(f"Input image type is not a Tensor. Got {type(image)}.")
  158. if centers is not None and not isinstance(centers, Tensor):
  159. raise TypeError(f"Bins' centers type is not a Tensor. Got {type(centers)}.")
  160. if centers is not None and len(centers.shape) > 0 and centers.dim() != 1:
  161. raise ValueError(f"Bins' centers must be a Tensor of the shape (n_bins,). Got {centers.shape}.")
  162. if not isinstance(min, float):
  163. raise TypeError(f"Type of lower end of the range is not a float. Got {type(min)}.")
  164. if not isinstance(max, float):
  165. raise TypeError(f"Type of upper end of the range is not a float. Got {type(min)}.")
  166. if not isinstance(n_bins, int):
  167. raise TypeError(f"Type of number of bins is not an int. Got {type(n_bins)}.")
  168. if bandwidth is not None and not isinstance(bandwidth, float):
  169. raise TypeError(f"Bandwidth type is not a float. Got {type(bandwidth)}.")
  170. if not isinstance(return_pdf, bool):
  171. raise TypeError(f"Return_pdf type is not a bool. Got {type(return_pdf)}.")
  172. if bandwidth is None:
  173. bandwidth = (max - min) / n_bins
  174. if centers is None:
  175. centers = min + bandwidth * (torch.arange(n_bins, device=image.device, dtype=image.dtype) + 0.5)
  176. centers = centers.reshape(-1, 1, 1, 1, 1)
  177. u = torch.abs(image.unsqueeze(0) - centers) / bandwidth
  178. if kernel == "gaussian":
  179. kernel_values = torch.exp(-0.5 * u**2)
  180. elif kernel in ("triangular", "uniform", "epanechnikov"):
  181. # compute the mask and cast to floating point
  182. mask = (u <= 1).to(u.dtype)
  183. if kernel == "triangular":
  184. kernel_values = (1.0 - u) * mask
  185. elif kernel == "uniform":
  186. kernel_values = mask
  187. else: # kernel == "epanechnikov"
  188. kernel_values = (1.0 - u**2) * mask
  189. else:
  190. raise ValueError(f"Kernel must be 'triangular', 'gaussian', 'uniform' or 'epanechnikov'. Got {kernel}.")
  191. hist = torch.sum(kernel_values, dim=(-2, -1)).permute(1, 2, 0)
  192. if return_pdf:
  193. normalization = torch.sum(hist, dim=-1, keepdim=True) + eps
  194. pdf = hist / normalization
  195. if image.dim() == 2:
  196. hist = hist.squeeze()
  197. pdf = pdf.squeeze()
  198. elif image.dim() == 3:
  199. hist = hist.squeeze(0)
  200. pdf = pdf.squeeze(0)
  201. return hist, pdf
  202. if image.dim() == 2:
  203. hist = hist.squeeze()
  204. elif image.dim() == 3:
  205. hist = hist.squeeze(0)
  206. return hist, torch.zeros_like(hist)