1#cython: cdivision=True 2#cython: boundscheck=False 3#cython: nonecheck=False 4#cython: wraparound=False 5import numpy as np 6 7cimport numpy as cnp 8 9from cpython.mem cimport PyMem_Malloc, PyMem_Free 10from libc.stdlib cimport abs 11from libc.math cimport fabs, sqrt, ceil, atan2, M_PI 12 13from ..draw import circle_perimeter 14 15from .._shared.interpolation cimport round 16 17cnp.import_array() 18 19def _hough_circle(cnp.ndarray img, 20 cnp.ndarray[ndim=1, dtype=cnp.intp_t] radius, 21 char normalize=True, char full_output=False): 22 """Perform a circular Hough transform. 23 24 Parameters 25 ---------- 26 img : (M, N) ndarray 27 Input image with nonzero values representing edges. 28 radius : ndarray 29 Radii at which to compute the Hough transform. 30 normalize : boolean, optional (default True) 31 Normalize the accumulator with the number 32 of pixels used to draw the radius. 33 full_output : boolean, optional (default False) 34 Extend the output size by twice the largest 35 radius in order to detect centers outside the 36 input picture. 37 38 Returns 39 ------- 40 H : 3D ndarray (radius index, (M + 2R, N + 2R) ndarray) 41 Hough transform accumulator for each radius. 42 R designates the larger radius if full_output is True. 43 Otherwise, R = 0. 44 """ 45 if img.ndim != 2: 46 raise ValueError('The input image must be 2D.') 47 48 cdef Py_ssize_t xmax = img.shape[0] 49 cdef Py_ssize_t ymax = img.shape[1] 50 51 # compute the nonzero indexes 52 cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] x, y 53 x, y = np.nonzero(img) 54 55 cdef Py_ssize_t num_pixels = x.size 56 57 cdef Py_ssize_t offset = 0 58 if full_output: 59 # Offset the image 60 offset = radius.max() 61 x = x + offset 62 y = y + offset 63 64 cdef Py_ssize_t i, p, c, num_circle_pixels, tx, ty 65 cdef double incr 66 cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] circle_x, circle_y 67 68 cdef cnp.ndarray[ndim=3, dtype=cnp.double_t] acc = \ 69 np.zeros((radius.size, 70 img.shape[0] + 2 * offset, 71 img.shape[1] + 2 * offset), dtype=np.double) 72 73 for i, rad in enumerate(radius): 74 # Store in memory the circle of given radius 75 # centered at (0,0) 76 circle_x, circle_y = circle_perimeter(0, 0, rad) 77 78 num_circle_pixels = circle_x.size 79 80 with nogil: 81 82 if normalize: 83 incr = 1.0 / num_circle_pixels 84 else: 85 incr = 1 86 87 # For each non zero pixel 88 for p in range(num_pixels): 89 # Plug the circle at (px, py), 90 # its coordinates are (tx, ty) 91 for c in range(num_circle_pixels): 92 tx = circle_x[c] + x[p] 93 ty = circle_y[c] + y[p] 94 if offset: 95 acc[i, tx, ty] += incr 96 elif 0 <= tx < xmax and 0 <= ty < ymax: 97 acc[i, tx, ty] += incr 98 99 return acc 100 101 102def _hough_ellipse(cnp.ndarray img, Py_ssize_t threshold=4, double accuracy=1, 103 Py_ssize_t min_size=4, max_size=None): 104 """Perform an elliptical Hough transform. 105 106 Parameters 107 ---------- 108 img : (M, N) ndarray 109 Input image with nonzero values representing edges. 110 threshold: int, optional (default 4) 111 Accumulator threshold value. 112 accuracy : double, optional (default 1) 113 Bin size on the minor axis used in the accumulator. 114 min_size : int, optional (default 4) 115 Minimal major axis length. 116 max_size : int, optional 117 Maximal minor axis length. (default None) 118 If None, the value is set to the half of the smaller 119 image dimension. 120 121 Returns 122 ------- 123 result : ndarray with fields [(accumulator, yc, xc, a, b, orientation)] 124 Where ``(yc, xc)`` is the center, ``(a, b)`` the major and minor 125 axes, respectively. The `orientation` value follows 126 `skimage.draw.ellipse_perimeter` convention. 127 128 Examples 129 -------- 130 >>> img = np.zeros((25, 25), dtype=np.uint8) 131 >>> rr, cc = ellipse_perimeter(10, 10, 6, 8) 132 >>> img[cc, rr] = 1 133 >>> result = hough_ellipse(img, threshold=8) 134 [(10, 10.0, 8.0, 6.0, 0.0, 10.0)] 135 136 Notes 137 ----- 138 The accuracy must be chosen to produce a peak in the accumulator 139 distribution. In other words, a flat accumulator distribution with low 140 values may be caused by a too low bin size. 141 142 References 143 ---------- 144 .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection 145 method." Pattern Recognition, 2002. Proceedings. 16th International 146 Conference on. Vol. 2. IEEE, 2002 147 """ 148 if img.ndim != 2: 149 raise ValueError('The input image must be 2D.') 150 151 # The creation of the array `pixels` results in a rather nasty error 152 # when the image is empty. 153 # As discussed in GitHub #2820 and #2996, we opt to return an empty array. 154 if not np.any(img): 155 return np.zeros((0, 6)) 156 157 cdef Py_ssize_t[:, ::1] pixels = np.row_stack(np.nonzero(img)) 158 159 cdef Py_ssize_t num_pixels = pixels.shape[1] 160 cdef list acc = list() 161 cdef list results = list() 162 cdef double bin_size = accuracy * accuracy 163 164 cdef double max_b_squared 165 if max_size is None: 166 if img.shape[0] < img.shape[1]: 167 max_b_squared = np.round(0.5 * img.shape[0]) 168 else: 169 max_b_squared = np.round(0.5 * img.shape[1]) 170 max_b_squared *= max_b_squared 171 else: 172 max_b_squared = max_size * max_size 173 174 cdef Py_ssize_t p1, p2, p3, p1x, p1y, p2x, p2y, p3x, p3y 175 cdef double xc, yc, a, b, d, k, dx, dy 176 cdef double cos_tau_squared, b_squared, orientation 177 178 for p1 in range(num_pixels): 179 p1x = pixels[1, p1] 180 p1y = pixels[0, p1] 181 182 for p2 in range(p1): 183 p2x = pixels[1, p2] 184 p2y = pixels[0, p2] 185 186 # Candidate: center (xc, yc) and main axis a 187 dx = p1x - p2x 188 dy = p1y - p2y 189 a = 0.5 * sqrt(dx * dx + dy * dy) 190 if a > 0.5 * min_size: 191 xc = 0.5 * (p1x + p2x) 192 yc = 0.5 * (p1y + p2y) 193 194 for p3 in range(num_pixels): 195 p3x = pixels[1, p3] 196 p3y = pixels[0, p3] 197 dx = p3x - xc 198 dy = p3y - yc 199 d = sqrt(dx * dx + dy * dy) 200 if d > min_size: 201 dx = p3x - p1x 202 dy = p3y - p1y 203 cos_tau_squared = ((a*a + d*d - dx*dx - dy*dy) 204 / (2 * a * d)) 205 cos_tau_squared *= cos_tau_squared 206 # Consider b2 > 0 and avoid division by zero 207 k = a*a - d*d * cos_tau_squared 208 if k > 0 and cos_tau_squared < 1: 209 b_squared = a*a * d*d * (1 - cos_tau_squared) / k 210 # b2 range is limited to avoid histogram memory 211 # overflow 212 if b_squared <= max_b_squared: 213 acc.append(b_squared) 214 215 if len(acc) > 0: 216 bins = np.arange(0, np.max(acc) + bin_size, bin_size) 217 hist, bin_edges = np.histogram(acc, bins=bins) 218 hist_max = np.max(hist) 219 if hist_max > threshold: 220 orientation = atan2(p1x - p2x, p1y - p2y) 221 b = sqrt(bin_edges[hist.argmax()]) 222 # to keep ellipse_perimeter() convention 223 if orientation != 0: 224 orientation = M_PI - orientation 225 # When orientation is not in [-pi:pi] 226 # it would mean in ellipse_perimeter() 227 # that a < b. But we keep a > b. 228 if orientation > M_PI: 229 orientation = orientation - M_PI / 2. 230 a, b = b, a 231 results.append((hist_max, # Accumulator 232 yc, xc, 233 a, b, 234 orientation)) 235 acc = [] 236 237 return np.array(results, dtype=[('accumulator', np.intp), 238 ('yc', np.double), 239 ('xc', np.double), 240 ('a', np.double), 241 ('b', np.double), 242 ('orientation', np.double)]) 243 244 245def _hough_line(cnp.ndarray img, 246 cnp.ndarray[ndim=1, dtype=cnp.double_t] theta): 247 """Perform a straight line Hough transform. 248 249 Parameters 250 ---------- 251 img : (M, N) ndarray 252 Input image with nonzero values representing edges. 253 theta : 1D ndarray of double 254 Angles at which to compute the transform, in radians. 255 256 Returns 257 ------- 258 H : 2-D ndarray of uint64 259 Hough transform accumulator. 260 theta : ndarray 261 Angles at which the transform was computed, in radians. 262 distances : ndarray 263 Distance values. 264 265 Notes 266 ----- 267 The origin is the top left corner of the original image. 268 X and Y axis are horizontal and vertical edges respectively. 269 The distance is the minimal algebraic distance from the origin 270 to the detected line. 271 272 Examples 273 -------- 274 Generate a test image: 275 276 >>> img = np.zeros((100, 150), dtype=bool) 277 >>> img[30, :] = 1 278 >>> img[:, 65] = 1 279 >>> img[35:45, 35:50] = 1 280 >>> for i in range(90): 281 ... img[i, i] = 1 282 >>> rng = np.random.default_rng() 283 >>> img += rng.random(img.shape) > 0.95 284 285 Apply the Hough transform: 286 287 >>> out, angles, d = hough_line(img) 288 289 .. plot:: hough_tf.py 290 291 """ 292 # Compute the array of angles and their sine and cosine 293 cdef cnp.ndarray[ndim=1, dtype=cnp.double_t] ctheta 294 cdef cnp.ndarray[ndim=1, dtype=cnp.double_t] stheta 295 296 ctheta = np.cos(theta) 297 stheta = np.sin(theta) 298 299 # compute the bins and allocate the accumulator array 300 cdef cnp.ndarray[ndim=2, dtype=cnp.uint64_t] accum 301 cdef cnp.ndarray[ndim=1, dtype=cnp.double_t] bins 302 cdef Py_ssize_t max_distance, offset 303 304 offset = <Py_ssize_t>ceil(sqrt(img.shape[0] * img.shape[0] + 305 img.shape[1] * img.shape[1])) 306 max_distance = 2 * offset + 1 307 accum = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64) 308 bins = np.linspace(-offset, offset, max_distance) 309 310 # compute the nonzero indexes 311 cdef cnp.ndarray[ndim=1, dtype=cnp.npy_intp] x_idxs, y_idxs 312 y_idxs, x_idxs = np.nonzero(img) 313 314 # finally, run the transform 315 cdef Py_ssize_t nidxs, nthetas, i, j, x, y, accum_idx 316 317 nidxs = y_idxs.shape[0] # x and y are the same shape 318 nthetas = theta.shape[0] 319 with nogil: 320 for i in range(nidxs): 321 x = x_idxs[i] 322 y = y_idxs[i] 323 for j in range(nthetas): 324 accum_idx = round((ctheta[j] * x + stheta[j] * y)) + offset 325 accum[accum_idx, j] += 1 326 327 return accum, theta, bins 328 329 330def _probabilistic_hough_line(cnp.ndarray img, Py_ssize_t threshold, 331 Py_ssize_t line_length, Py_ssize_t line_gap, 332 cnp.ndarray[ndim=1, dtype=cnp.double_t] theta, 333 seed=None): 334 """Return lines from a progressive probabilistic line Hough transform. 335 336 Parameters 337 ---------- 338 img : (M, N) ndarray 339 Input image with nonzero values representing edges. 340 threshold : int 341 Threshold 342 line_length : int 343 Minimum accepted length of detected lines. 344 Increase the parameter to extract longer lines. 345 line_gap : int 346 Maximum gap between pixels to still form a line. 347 Increase the parameter to merge broken lines more aggressively. 348 theta : 1D ndarray, dtype=double 349 Angles at which to compute the transform, in radians. 350 seed : {None, int, `numpy.random.Generator`, optional} 351 If `seed` is None the `numpy.random.Generator` singleton is used. 352 If `seed` is an int, a new ``Generator`` instance is used, 353 seeded with `seed`. 354 If `seed` is already a ``Generator`` instance then that instance is 355 used. 356 357 Seed to initialize the random number generator. 358 359 Returns 360 ------- 361 lines : list 362 List of lines identified, lines in format ((x0, y0), (x1, y0)), 363 indicating line start and end. 364 365 References 366 ---------- 367 .. [1] C. Galamhos, J. Matas and J. Kittler, "Progressive probabilistic 368 Hough transform for line detection", in IEEE Computer Society 369 Conference on Computer Vision and Pattern Recognition, 1999. 370 """ 371 cdef Py_ssize_t height = img.shape[0] 372 cdef Py_ssize_t width = img.shape[1] 373 374 # compute the bins and allocate the accumulator array 375 cdef cnp.ndarray[ndim=2, dtype=cnp.uint8_t] mask = \ 376 np.zeros((height, width), dtype=np.uint8) 377 cdef Py_ssize_t *line_end = \ 378 <Py_ssize_t *>PyMem_Malloc(4 * sizeof(Py_ssize_t)) 379 if not line_end: 380 raise MemoryError('could not allocate line_end') 381 cdef Py_ssize_t max_distance, offset, num_indexes, index 382 cdef double a, b 383 cdef Py_ssize_t nidxs, i, j, k, x, y, px, py, accum_idx, max_theta 384 cdef Py_ssize_t xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, count 385 cdef cnp.int64_t value, max_value, 386 cdef int shift = 16 387 cdef int good_line 388 cdef Py_ssize_t nlines = 0 389 cdef Py_ssize_t lines_max = 2 ** 15 # maximum line number cutoff 390 cdef cnp.intp_t[:, :, ::1] lines = np.zeros((lines_max, 2, 2), 391 dtype=np.intp) 392 max_distance = 2 * <Py_ssize_t>ceil((sqrt(img.shape[0] * img.shape[0] + 393 img.shape[1] * img.shape[1]))) 394 cdef cnp.int64_t[:, ::1] accum = np.zeros((max_distance, theta.shape[0]), 395 dtype=np.int64) 396 offset = max_distance / 2 397 cdef Py_ssize_t nthetas = theta.shape[0] 398 399 # compute sine and cosine of angles 400 cdef cnp.double_t[::1] ctheta = np.cos(theta) 401 cdef cnp.double_t[::1] stheta = np.sin(theta) 402 403 # find the nonzero indexes 404 cdef cnp.intp_t[:] y_idxs, x_idxs 405 y_idxs, x_idxs = np.nonzero(img) 406 407 # mask all non-zero indexes 408 mask[y_idxs, x_idxs] = 1 409 410 count = len(x_idxs) 411 random_state = np.random.default_rng(seed) 412 random_ = np.arange(count, dtype=np.intp) 413 random_state.shuffle(random_) 414 cdef cnp.intp_t[::1] random = random_ 415 416 with nogil: 417 while count > 0: 418 count -= 1 419 # select random non-zero point 420 index = random[count] 421 x = x_idxs[index] 422 y = y_idxs[index] 423 424 # if previously eliminated, skip 425 if not mask[y, x]: 426 continue 427 428 value = 0 429 max_value = threshold - 1 430 max_theta = -1 431 432 # apply hough transform on point 433 for j in range(nthetas): 434 accum_idx = round((ctheta[j] * x + stheta[j] * y)) + offset 435 accum[accum_idx, j] += 1 436 value = accum[accum_idx, j] 437 if value > max_value: 438 max_value = value 439 max_theta = j 440 if max_value < threshold: 441 continue 442 443 # from the random point walk in opposite directions and find line 444 # beginning and end 445 a = -stheta[max_theta] 446 b = ctheta[max_theta] 447 x0 = x 448 y0 = y 449 # calculate gradient of walks using fixed point math 450 xflag = fabs(a) > fabs(b) 451 if xflag: 452 if a > 0: 453 dx0 = 1 454 else: 455 dx0 = -1 456 dy0 = round(b * (1 << shift) / fabs(a)) 457 y0 = (y0 << shift) + (1 << (shift - 1)) 458 else: 459 if b > 0: 460 dy0 = 1 461 else: 462 dy0 = -1 463 dx0 = round(a * (1 << shift) / fabs(b)) 464 x0 = (x0 << shift) + (1 << (shift - 1)) 465 466 # pass 1: walk the line, merging lines less than specified gap 467 # length 468 for k in range(2): 469 gap = 0 470 px = x0 471 py = y0 472 dx = dx0 473 dy = dy0 474 if k > 0: 475 dx = -dx 476 dy = -dy 477 while 1: 478 if xflag: 479 x1 = px 480 y1 = py >> shift 481 else: 482 x1 = px >> shift 483 y1 = py 484 # check when line exits image boundary 485 if x1 < 0 or x1 >= width or y1 < 0 or y1 >= height: 486 break 487 gap += 1 488 # if non-zero point found, continue the line 489 if mask[y1, x1]: 490 gap = 0 491 line_end[2*k] = x1 492 line_end[2*k + 1] = y1 493 # if gap to this point was too large, end the line 494 elif gap > line_gap: 495 break 496 px += dx 497 py += dy 498 499 # confirm line length is sufficient 500 good_line = (abs(line_end[3] - line_end[1]) >= line_length or 501 abs(line_end[2] - line_end[0]) >= line_length) 502 503 # pass 2: walk the line again and reset accumulator and mask 504 for k in range(2): 505 px = x0 506 py = y0 507 dx = dx0 508 dy = dy0 509 if k > 0: 510 dx = -dx 511 dy = -dy 512 while 1: 513 if xflag: 514 x1 = px 515 y1 = py >> shift 516 else: 517 x1 = px >> shift 518 y1 = py 519 # if non-zero point found, continue the line 520 if mask[y1, x1]: 521 if good_line: 522 accum_idx = round( 523 (ctheta[j] * x1 + stheta[j] * y1)) + offset 524 accum[accum_idx, max_theta] -= 1 525 mask[y1, x1] = 0 526 # exit when the point is the line end 527 if x1 == line_end[2*k] and y1 == line_end[2*k + 1]: 528 break 529 px += dx 530 py += dy 531 532 # add line to the result 533 if good_line: 534 lines[nlines, 0, 0] = line_end[0] 535 lines[nlines, 0, 1] = line_end[1] 536 lines[nlines, 1, 0] = line_end[2] 537 lines[nlines, 1, 1] = line_end[3] 538 nlines += 1 539 if nlines >= lines_max: 540 break 541 542 PyMem_Free(line_end) 543 return [((line[0, 0], line[0, 1]), (line[1, 0], line[1, 1])) 544 for line in lines[:nlines]] 545