hough_transform.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473
  1. import numpy as np
  2. from scipy.spatial import cKDTree
  3. from ._hough_transform import _hough_circle, _hough_ellipse, _hough_line
  4. from ._hough_transform import _probabilistic_hough_line as _prob_hough_line
  5. def hough_line_peaks(
  6. hspace,
  7. angles,
  8. dists,
  9. min_distance=9,
  10. min_angle=10,
  11. threshold=None,
  12. num_peaks=np.inf,
  13. ):
  14. """Return peaks in a straight line Hough transform.
  15. Identifies most prominent lines separated by a certain angle and distance
  16. in a Hough transform. Non-maximum suppression with different sizes is
  17. applied separately in the first (distances) and second (angles) dimension
  18. of the Hough space to identify peaks.
  19. Parameters
  20. ----------
  21. hspace : ndarray, shape (M, N)
  22. Hough space returned by the `hough_line` function.
  23. angles : array, shape (N,)
  24. Angles returned by the `hough_line` function. Assumed to be continuous.
  25. (`angles[-1] - angles[0] == PI`).
  26. dists : array, shape (M,)
  27. Distances returned by the `hough_line` function.
  28. min_distance : int, optional
  29. Minimum distance separating lines (maximum filter size for first
  30. dimension of hough space).
  31. min_angle : int, optional
  32. Minimum angle separating lines (maximum filter size for second
  33. dimension of hough space).
  34. threshold : float, optional
  35. Minimum intensity of peaks. Default is `0.5 * max(hspace)`.
  36. num_peaks : int, optional
  37. Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
  38. return `num_peaks` coordinates based on peak intensity.
  39. Returns
  40. -------
  41. accum, angles, dists : tuple of array
  42. Peak values in Hough space, angles and distances.
  43. Examples
  44. --------
  45. >>> from skimage.transform import hough_line, hough_line_peaks
  46. >>> from skimage.draw import line
  47. >>> img = np.zeros((15, 15), dtype=bool)
  48. >>> rr, cc = line(0, 0, 14, 14)
  49. >>> img[rr, cc] = 1
  50. >>> rr, cc = line(0, 14, 14, 0)
  51. >>> img[cc, rr] = 1
  52. >>> hspace, angles, dists = hough_line(img)
  53. >>> hspace, angles, dists = hough_line_peaks(hspace, angles, dists)
  54. >>> len(angles)
  55. 2
  56. """
  57. from ..feature.peak import _prominent_peaks
  58. min_angle = min(min_angle, hspace.shape[1])
  59. h, a, d = _prominent_peaks(
  60. hspace,
  61. min_xdistance=min_angle,
  62. min_ydistance=min_distance,
  63. threshold=threshold,
  64. num_peaks=num_peaks,
  65. )
  66. if a.size > 0:
  67. return (h, angles[a], dists[d])
  68. else:
  69. return (h, np.array([]), np.array([]))
  70. def hough_circle(image, radius, normalize=True, full_output=False):
  71. """Perform a circular Hough transform.
  72. Parameters
  73. ----------
  74. image : ndarray, shape (M, N)
  75. Input image with nonzero values representing edges.
  76. radius : scalar or sequence of scalars
  77. Radii at which to compute the Hough transform.
  78. Floats are converted to integers.
  79. normalize : bool, optional
  80. Normalize the accumulator with the number
  81. of pixels used to draw the radius.
  82. full_output : bool, optional
  83. Extend the output size by twice the largest
  84. radius in order to detect centers outside the
  85. input picture.
  86. Returns
  87. -------
  88. H : ndarray, shape (radius index, M + 2R, N + 2R)
  89. Hough transform accumulator for each radius.
  90. R designates the larger radius if full_output is True.
  91. Otherwise, R = 0.
  92. Examples
  93. --------
  94. >>> from skimage.transform import hough_circle
  95. >>> from skimage.draw import circle_perimeter
  96. >>> img = np.zeros((100, 100), dtype=bool)
  97. >>> rr, cc = circle_perimeter(25, 35, 23)
  98. >>> img[rr, cc] = 1
  99. >>> try_radii = np.arange(5, 50)
  100. >>> res = hough_circle(img, try_radii)
  101. >>> ridx, r, c = np.unravel_index(np.argmax(res), res.shape)
  102. >>> r, c, try_radii[ridx]
  103. (25, 35, 23)
  104. """
  105. radius = np.atleast_1d(np.asarray(radius))
  106. return _hough_circle(
  107. image, radius.astype(np.intp), normalize=normalize, full_output=full_output
  108. )
  109. def hough_ellipse(image, threshold=4, accuracy=1, min_size=4, max_size=None):
  110. """Perform an elliptical Hough transform.
  111. Parameters
  112. ----------
  113. image : (M, N) ndarray
  114. Input image with nonzero values representing edges.
  115. threshold : int, optional
  116. Accumulator threshold value. A lower value will return more ellipses.
  117. accuracy : double, optional
  118. Bin size on the minor axis used in the accumulator. A higher value
  119. will return more ellipses, but lead to a less precise estimation of
  120. the minor axis lengths.
  121. min_size : int, optional
  122. Minimal major axis length.
  123. max_size : int, optional
  124. Maximal minor axis length.
  125. If None, the value is set to half of the smaller
  126. image dimension.
  127. Returns
  128. -------
  129. result : ndarray with fields [(accumulator, yc, xc, a, b, orientation)].
  130. Where ``(yc, xc)`` is the center, ``(a, b)`` the major and minor
  131. axes, respectively. The `orientation` value follows the
  132. `skimage.draw.ellipse_perimeter` convention.
  133. Examples
  134. --------
  135. >>> from skimage.transform import hough_ellipse
  136. >>> from skimage.draw import ellipse_perimeter
  137. >>> img = np.zeros((25, 25), dtype=np.uint8)
  138. >>> rr, cc = ellipse_perimeter(10, 10, 6, 8)
  139. >>> img[cc, rr] = 1
  140. >>> result = hough_ellipse(img, threshold=8)
  141. >>> result.tolist()
  142. [(10, 10.0, 10.0, 8.0, 6.0, 0.0)]
  143. Notes
  144. -----
  145. Potential ellipses in the image are characterized by their major and
  146. minor axis lengths. For any pair of nonzero pixels in the image that
  147. are at least half of `min_size` apart, an accumulator keeps track of
  148. the minor axis lengths of potential ellipses formed with all the
  149. other nonzero pixels. If any bin (with `bin_size = accuracy * accuracy`)
  150. in the histogram of those accumulated minor axis lengths is above
  151. `threshold`, the corresponding ellipse is added to the results.
  152. A higher `accuracy` will therefore lead to more ellipses being found
  153. in the image, at the cost of a less precise estimation of the minor
  154. axis length.
  155. References
  156. ----------
  157. .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection
  158. method." Pattern Recognition, 2002. Proceedings. 16th International
  159. Conference on. Vol. 2. IEEE, 2002
  160. """
  161. return _hough_ellipse(
  162. image,
  163. threshold=threshold,
  164. accuracy=accuracy,
  165. min_size=min_size,
  166. max_size=max_size,
  167. )
  168. def hough_line(image, theta=None):
  169. """Perform a straight line Hough transform.
  170. Parameters
  171. ----------
  172. image : (M, N) ndarray
  173. Input image with nonzero values representing edges.
  174. theta : ndarray of double, shape (K,), optional
  175. Angles at which to compute the transform, in radians.
  176. Defaults to a vector of 180 angles evenly spaced in the
  177. range [-pi/2, pi/2).
  178. Returns
  179. -------
  180. hspace : ndarray of uint64, shape (P, Q)
  181. Hough transform accumulator.
  182. angles : ndarray
  183. Angles at which the transform is computed, in radians.
  184. distances : ndarray
  185. Distance values.
  186. Notes
  187. -----
  188. The origin is the top left corner of the original image.
  189. X and Y axis are horizontal and vertical edges respectively.
  190. The distance is the minimal algebraic distance from the origin
  191. to the detected line.
  192. The angle accuracy can be improved by decreasing the step size in
  193. the `theta` array.
  194. Examples
  195. --------
  196. Generate a test image:
  197. >>> img = np.zeros((100, 150), dtype=bool)
  198. >>> img[30, :] = 1
  199. >>> img[:, 65] = 1
  200. >>> img[35:45, 35:50] = 1
  201. >>> for i in range(90):
  202. ... img[i, i] = 1
  203. >>> rng = np.random.default_rng()
  204. >>> img += rng.random(img.shape) > 0.95
  205. Apply the Hough transform:
  206. >>> out, angles, d = hough_line(img)
  207. """
  208. if image.ndim != 2:
  209. raise ValueError('The input image `image` must be 2D.')
  210. if theta is None:
  211. # These values are approximations of pi/2
  212. theta = np.linspace(-np.pi / 2, np.pi / 2, 180, endpoint=False)
  213. return _hough_line(image, theta=theta)
  214. def probabilistic_hough_line(
  215. image, threshold=10, line_length=50, line_gap=10, theta=None, rng=None
  216. ):
  217. """Return lines from a progressive probabilistic line Hough transform.
  218. Parameters
  219. ----------
  220. image : ndarray, shape (M, N)
  221. Input image with nonzero values representing edges.
  222. threshold : int, optional
  223. Threshold
  224. line_length : int, optional
  225. Minimum accepted length of detected lines.
  226. Increase the parameter to extract longer lines.
  227. line_gap : int, optional
  228. Maximum gap between pixels to still form a line.
  229. Increase the parameter to merge broken lines more aggressively.
  230. theta : ndarray of dtype, shape (K,), optional
  231. Angles at which to compute the transform, in radians.
  232. Defaults to a vector of 180 angles evenly spaced in the
  233. range [-pi/2, pi/2).
  234. rng : {`numpy.random.Generator`, int}, optional
  235. Pseudo-random number generator.
  236. By default, a PCG64 generator is used (see :func:`numpy.random.default_rng`).
  237. If `rng` is an int, it is used to seed the generator.
  238. Returns
  239. -------
  240. lines : list
  241. List of lines identified, lines in format ((x0, y0), (x1, y1)),
  242. indicating line start and end.
  243. References
  244. ----------
  245. .. [1] C. Galamhos, J. Matas and J. Kittler, "Progressive probabilistic
  246. Hough transform for line detection", in IEEE Computer Society
  247. Conference on Computer Vision and Pattern Recognition, 1999.
  248. """
  249. if image.ndim != 2:
  250. raise ValueError('The input image `image` must be 2D.')
  251. if theta is None:
  252. theta = np.linspace(-np.pi / 2, np.pi / 2, 180, endpoint=False)
  253. return _prob_hough_line(
  254. image,
  255. threshold=threshold,
  256. line_length=line_length,
  257. line_gap=line_gap,
  258. theta=theta,
  259. rng=rng,
  260. )
  261. def hough_circle_peaks(
  262. hspaces,
  263. radii,
  264. min_xdistance=1,
  265. min_ydistance=1,
  266. threshold=None,
  267. num_peaks=np.inf,
  268. total_num_peaks=np.inf,
  269. normalize=False,
  270. ):
  271. """Return peaks in a circle Hough transform.
  272. Identifies most prominent circles separated by certain distances in given
  273. Hough spaces. Non-maximum suppression with different sizes is applied
  274. separately in the first and second dimension of the Hough space to
  275. identify peaks. For circles with different radius but close in distance,
  276. only the one with highest peak is kept.
  277. Parameters
  278. ----------
  279. hspaces : (M, N, P) array
  280. Hough spaces returned by the `hough_circle` function.
  281. radii : (M,) array
  282. Radii corresponding to Hough spaces.
  283. min_xdistance : int, optional
  284. Minimum distance separating centers in the x dimension.
  285. min_ydistance : int, optional
  286. Minimum distance separating centers in the y dimension.
  287. threshold : float, optional
  288. Minimum intensity of peaks in each Hough space.
  289. Default is `0.5 * max(hspace)`.
  290. num_peaks : int, optional
  291. Maximum number of peaks in each Hough space. When the
  292. number of peaks exceeds `num_peaks`, only `num_peaks`
  293. coordinates based on peak intensity are considered for the
  294. corresponding radius.
  295. total_num_peaks : int, optional
  296. Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
  297. return `num_peaks` coordinates based on peak intensity.
  298. normalize : bool, optional
  299. If True, normalize the accumulator by the radius to sort the prominent
  300. peaks.
  301. Returns
  302. -------
  303. accum, cx, cy, rad : tuple of array
  304. Peak values in Hough space, x and y center coordinates and radii.
  305. Examples
  306. --------
  307. >>> from skimage import transform, draw
  308. >>> img = np.zeros((120, 100), dtype=int)
  309. >>> radius, x_0, y_0 = (20, 99, 50)
  310. >>> y, x = draw.circle_perimeter(y_0, x_0, radius)
  311. >>> img[x, y] = 1
  312. >>> hspaces = transform.hough_circle(img, radius)
  313. >>> accum, cx, cy, rad = hough_circle_peaks(hspaces, [radius,])
  314. Notes
  315. -----
  316. Circles with bigger radius have higher peaks in Hough space. If larger
  317. circles are preferred over smaller ones, `normalize` should be False.
  318. Otherwise, circles will be returned in the order of decreasing voting
  319. number.
  320. """
  321. from ..feature.peak import _prominent_peaks
  322. r = []
  323. cx = []
  324. cy = []
  325. accum = []
  326. for rad, hp in zip(radii, hspaces):
  327. h_p, x_p, y_p = _prominent_peaks(
  328. hp,
  329. min_xdistance=min_xdistance,
  330. min_ydistance=min_ydistance,
  331. threshold=threshold,
  332. num_peaks=num_peaks,
  333. )
  334. r.extend((rad,) * len(h_p))
  335. cx.extend(x_p)
  336. cy.extend(y_p)
  337. accum.extend(h_p)
  338. r = np.array(r)
  339. cx = np.array(cx)
  340. cy = np.array(cy)
  341. accum = np.array(accum)
  342. if normalize:
  343. s = np.argsort(accum / r)
  344. else:
  345. s = np.argsort(accum)
  346. accum_sorted, cx_sorted, cy_sorted, r_sorted = (
  347. accum[s][::-1],
  348. cx[s][::-1],
  349. cy[s][::-1],
  350. r[s][::-1],
  351. )
  352. tnp = len(accum_sorted) if total_num_peaks == np.inf else total_num_peaks
  353. # Skip searching for neighboring circles
  354. # if default min_xdistance and min_ydistance are used
  355. # or if no peak was detected
  356. if (min_xdistance == 1 and min_ydistance == 1) or len(accum_sorted) == 0:
  357. return (accum_sorted[:tnp], cx_sorted[:tnp], cy_sorted[:tnp], r_sorted[:tnp])
  358. # For circles with centers too close, only keep the one with
  359. # the highest peak
  360. should_keep = label_distant_points(
  361. cx_sorted, cy_sorted, min_xdistance, min_ydistance, tnp
  362. )
  363. return (
  364. accum_sorted[should_keep],
  365. cx_sorted[should_keep],
  366. cy_sorted[should_keep],
  367. r_sorted[should_keep],
  368. )
  369. def label_distant_points(xs, ys, min_xdistance, min_ydistance, max_points):
  370. """Keep points that are separated by certain distance in each dimension.
  371. The first point is always accepted and all subsequent points are selected
  372. so that they are distant from all their preceding ones.
  373. Parameters
  374. ----------
  375. xs : array, shape (M,)
  376. X coordinates of points.
  377. ys : array, shape (M,)
  378. Y coordinates of points.
  379. min_xdistance : int
  380. Minimum distance separating points in the x dimension.
  381. min_ydistance : int
  382. Minimum distance separating points in the y dimension.
  383. max_points : int
  384. Max number of distant points to keep.
  385. Returns
  386. -------
  387. should_keep : array of bool
  388. A mask array for distant points to keep.
  389. """
  390. is_neighbor = np.zeros(len(xs), dtype=bool)
  391. coordinates = np.stack([xs, ys], axis=1)
  392. # Use a KDTree to search for neighboring points effectively
  393. kd_tree = cKDTree(coordinates)
  394. n_pts = 0
  395. for i in range(len(xs)):
  396. if n_pts >= max_points:
  397. # Ignore the point if points to keep reaches maximum
  398. is_neighbor[i] = True
  399. elif not is_neighbor[i]:
  400. # Find a short list of candidates to remove
  401. # by searching within a circle
  402. neighbors_i = kd_tree.query_ball_point(
  403. (xs[i], ys[i]), np.hypot(min_xdistance, min_ydistance)
  404. )
  405. # Check distance in both dimensions and mark if close
  406. for ni in neighbors_i:
  407. x_close = abs(xs[ni] - xs[i]) <= min_xdistance
  408. y_close = abs(ys[ni] - ys[i]) <= min_ydistance
  409. if x_close and y_close and ni > i:
  410. is_neighbor[ni] = True
  411. n_pts += 1
  412. should_keep = ~is_neighbor
  413. return should_keep