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