From a79dbc7e2b15f995c17b06705b45554fdc980915 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 15 Aug 2011 14:40:05 +0200
Subject: [PATCH 01/14] Probabilistic hough transform

---
 scikits/image/transform/_hough_transform.pyx | 163 ++++++++++++++++++-
 scikits/image/transform/hough_transform.py   |  37 ++++-
 2 files changed, 194 insertions(+), 6 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index ec85743d..9992eee6 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -2,11 +2,12 @@ cimport cython
 
 import numpy as np
 cimport numpy as np
-
+from random import randint
 np.import_array()
 
 
 cdef extern from "math.h":
+    double fabs(double)
     double sqrt(double)
     double ceil(double)
     double floor(double)
@@ -49,17 +50,173 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     cdef np.ndarray[ndim=1, dtype=np.int_t] x_idxs, y_idxs
     y_idxs, x_idxs = np.PyArray_Nonzero(img)
 
+
     # finally, run the transform
     cdef int nidxs, nthetas, i, j, x, y, out_idx
     nidxs = y_idxs.shape[0] # x and y are the same shape
     nthetas = theta.shape[0]
     for i in range(nidxs):
         x = x_idxs[i]
-        y = y_idxs[i]
+        y = y_idxs[i]   
         for j in range(nthetas):
             out_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
             out[out_idx, j] += 1
-
     return out, theta, bins
 
 
+@cython.boundscheck(False)
+def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, int line_gap, \
+   np.ndarray[ndim=1, dtype=np.double_t] theta=None):
+    if img.ndim != 2:
+        raise ValueError('The input image must be 2D.')
+    # compute the array of angles and their sine and cosine
+    cdef np.ndarray[ndim=1, dtype=np.double_t] ctheta
+    cdef np.ndarray[ndim=1, dtype=np.double_t] stheta
+    # calculate thetas if none specified
+    if theta is None:
+        theta = np.linspace(PI_2, NEG_PI_2, 180)
+    ctheta = np.cos(theta)
+    stheta = np.sin(theta)
+    cdef int height = img.shape[0]
+    cdef int width = img.shape[1]
+    # compute the bins and allocate the output array
+    cdef np.ndarray[ndim=2, dtype=np.uint64_t] out
+    cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = np.zeros((height, width), dtype=np.uint8)
+    cdef np.ndarray[ndim=2, dtype=np.uint32_t] line_end = np.zeros((2, 2), dtype=np.uint32)
+    cdef np.ndarray[ndim=1, dtype=np.double_t] bins
+    cdef int max_distance, offset, num_indexes, index    
+    cdef double a, b
+    cdef int nidxs, nthetas, i, j, x, y, px, py, out_idx, value, max_value, max_theta
+    cdef int shift = 16
+    # maximum line number cutoff
+    cdef int lines_max = 2 ** 15
+    cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line    
+    max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
+                                       img.shape[1] * img.shape[1])))
+    out = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
+    bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance)
+    offset = max_distance / 2
+    # find the nonzero indexes
+    cdef np.ndarray[ndim=1, dtype=np.int_t] x_idxs, y_idxs
+    y_idxs, x_idxs =  np.PyArray_Nonzero(img)
+    num_indexes = y_idxs.shape[0] # x and y are the same shape
+    nthetas = theta.shape[0]
+    lines = []
+    # create mask of all non-zero indexes
+    for i in range(num_indexes):
+        mask[y_idxs[i], x_idxs[i]] = 1
+
+    for i in range(num_indexes):
+        # select random non-zero point
+        index = randint(0, num_indexes-1)
+        x = x_idxs[i]
+        y = y_idxs[i]
+        # if previously eliminated, skip
+        if not mask[y, x]:
+            continue
+        value = 0
+        max_value = 0
+        max_theta = 0
+        # apply hough transform on point
+        for j in range(nthetas):
+            out_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
+            out[out_idx, j] += 1
+            value = out[out_idx, j]
+            if value > max_value:
+                max_value = value
+                max_theta = j
+        # accumulator value of point strong enough
+        if max_value < value_threshold:
+            continue
+        # from the random point walk in opposite directions and find line beginning and end
+        a = -stheta[max_theta]
+        b = ctheta[max_theta]
+        x0 = x
+        y0 = y
+       # calculate gradient of walks using fixed point math
+        xflag = fabs(a) > fabs(b)
+        if xflag:
+            if a > 0:
+                dx0 = 1
+            else:
+                dx0 = -1
+            dy0 = <int>round(b*(1 << shift)/fabs(a) ) 
+            y0 = (y0 << shift) + (1 << (shift - 1))
+        else:
+            if b > 0:
+                dy0 = 1
+            else:
+                dy0 = -1
+            dx0 = <int>round( a*(1 << shift)/fabs(b))
+            x0 = (x0 << shift) + (1 << (shift-1))
+        # pass 1: walk the line, merging lines less than specified gap length
+        for k in range(2):
+            gap = 0
+            px = x0
+            py = y0
+            dx = dx0
+            dy = dy0
+            if k > 0:
+                dx = -dx
+                dy = -dy
+            while 1:
+                if xflag:
+                    x1 = px
+                    y1 = py >> shift
+                else:
+                    x1 = px >> shift
+                    y1 = py;
+                # check when line exits image boundary
+                if x1 < 0 or x1 >= width or y1 < 0 or y1 >= height:
+                    break
+                gap += 1
+                # if non-zero point found, continue the line
+                if mask[y1, x1]:
+                    gap = 0;
+                    line_end[k, 1] = y1
+                    line_end[k, 0] = x1
+                # if gap to this point was too large, end the line
+                elif gap > line_gap:
+                    break
+                px += dx
+                py += dy            
+        # confirm line length is sufficient
+        good_line = fabs(line_end[1, 1] - line_end[0, 1]) >= line_length or \
+                    fabs(line_end[1, 0] - line_end[0, 0]) >= line_length
+        # pass 2: walk the line again and reset accumulator and mask
+        for k in range(2):
+            px = x0
+            py = y0
+            dx = dx0
+            dy = dy0
+            if k > 0:
+                dx = -dx
+                dy = -dy
+            while 1:
+                if xflag:
+                    x1 = px
+                    y1 = py >> shift
+                else:
+                    x1 = px >> shift
+                    y1 = py
+                # if non-zero point found, continue the line
+                if mask[y1, x1]:
+                    if good_line:
+                        for j in range(nthetas):
+                            out_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
+                            out[out_idx, j] -= 1
+                    mask[y1, x1] = 0
+                # exit when the point is the line end
+                if x1 == line_end[k, 0] and y1 == line_end[k, 1]:
+                    break;
+                px += dx
+                py += dy
+
+        # add line to the result
+        if good_line:
+            lines.append(((line_end[0, 0], line_end[0, 1]), (line_end[1, 0], line_end[1, 1])))
+            if len(lines) > lines_max:
+                return lines
+    return lines
+
+
diff --git a/scikits/image/transform/hough_transform.py b/scikits/image/transform/hough_transform.py
index 5d7a28ea..95d4c46b 100644
--- a/scikits/image/transform/hough_transform.py
+++ b/scikits/image/transform/hough_transform.py
@@ -1,7 +1,8 @@
-__all__ = ['hough']
+__all__ = ['hough', 'probabilistic_hough']
 
 from itertools import izip
 import numpy as np
+from _hough_transform import _probabilistic_hough  
 
 def _hough(img, theta=None):
     if img.ndim != 2:
@@ -53,11 +54,39 @@ _py_hough = _hough
 
 # try to import and use the faster Cython version if it exists
 try:
-    from ._hough_transform import _hough
+    from ._hough_transform import _hough  
 except ImportError:
     pass
 
 
+def probabilistic_hough(img, value_threshold=50, line_length=50, line_gap=10, theta=None):
+    """Performs a progressive probabilistic line Hough transform and returns the detected lines.
+
+    Parameters
+    ----------
+    img : (M, N) ndarray
+        Input image with nonzero values representing edges.
+    value_threshold: int
+        Threshold
+    theta :1D ndarray, dtype=double
+        Angles at which to compute the transform, in radians.
+        Defaults to -pi/2 .. pi/2
+
+    Returns
+    -------
+    lines : list
+      List of lines identified, lines in format ((x0, y0), (x1, y0)), indicating 
+      line start and end.
+
+    References
+    ----------
+    .. [1]  C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic Hough 
+            transform for line detection", in  IEEE Computer Society Conference on 
+            Computer Vision and Pattern Recognition, 1999. 
+    """    
+    return _probabilistic_hough(img, value_threshold, line_length, line_gap, theta)
+
+
 def hough(img, theta=None):
     """Perform a straight line Hough transform.
 
@@ -67,7 +96,7 @@ def hough(img, theta=None):
         Input image with nonzero values representing edges.
     theta : 1D ndarray of double
         Angles at which to compute the transform, in radians.
-        Defaults to -pi/2 - pi/2
+        Defaults to -pi/2 .. pi/2
 
     Returns
     -------
@@ -106,3 +135,5 @@ def hough(img, theta=None):
 
     """
     return _hough(img, theta)
+
+

From 24e3ee7e1b36447faeda7257f29e971a3f1a9154 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 15 Aug 2011 22:58:28 +0200
Subject: [PATCH 02/14] Testing different accumulator clearing on good lines.

---
 scikits/image/transform/_hough_transform.pyx | 79 ++++++++++++--------
 scikits/image/transform/hough_transform.py   |  2 +-
 2 files changed, 48 insertions(+), 33 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index 9992eee6..98693dbe 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -7,6 +7,7 @@ np.import_array()
 
 
 cdef extern from "math.h":
+    int abs(int)
     double fabs(double)
     double sqrt(double)
     double ceil(double)
@@ -35,14 +36,14 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     ctheta = np.cos(theta)
     stheta = np.sin(theta)
 
-    # compute the bins and allocate the output array
-    cdef np.ndarray[ndim=2, dtype=np.uint64_t] out
+    # compute the bins and allocate the accumulator array
+    cdef np.ndarray[ndim=2, dtype=np.uint64_t] accum
     cdef np.ndarray[ndim=1, dtype=np.double_t] bins
     cdef int max_distance, offset 
 
     max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
                                        img.shape[1] * img.shape[1])))
-    out = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
+    accum = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
     bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance)
     offset = max_distance / 2
     
@@ -52,16 +53,16 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
 
 
     # finally, run the transform
-    cdef int nidxs, nthetas, i, j, x, y, out_idx
+    cdef int nidxs, nthetas, i, j, x, y, accum_idx
     nidxs = y_idxs.shape[0] # x and y are the same shape
     nthetas = theta.shape[0]
     for i in range(nidxs):
         x = x_idxs[i]
         y = y_idxs[i]   
         for j in range(nthetas):
-            out_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
-            out[out_idx, j] += 1
-    return out, theta, bins
+            accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
+            accum[accum_idx, j] += 1
+    return accum, theta, bins
 
 
 @cython.boundscheck(False)
@@ -79,38 +80,45 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, i
     stheta = np.sin(theta)
     cdef int height = img.shape[0]
     cdef int width = img.shape[1]
-    # compute the bins and allocate the output array
-    cdef np.ndarray[ndim=2, dtype=np.uint64_t] out
+    # compute the bins and allocate the accumulator array
+    cdef np.ndarray[ndim=2, dtype=np.int64_t] accum
     cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = np.zeros((height, width), dtype=np.uint8)
-    cdef np.ndarray[ndim=2, dtype=np.uint32_t] line_end = np.zeros((2, 2), dtype=np.uint32)
-    cdef np.ndarray[ndim=1, dtype=np.double_t] bins
+    cdef np.ndarray[ndim=2, dtype=np.int32_t] line_end = np.zeros((2, 2), dtype=np.int32)
     cdef int max_distance, offset, num_indexes, index    
     cdef double a, b
-    cdef int nidxs, nthetas, i, j, x, y, px, py, out_idx, value, max_value, max_theta
+    cdef int nidxs, nthetas, i, j, x, y, px, py, accum_idx, value, max_value, max_theta
     cdef int shift = 16
     # maximum line number cutoff
     cdef int lines_max = 2 ** 15
-    cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line    
+    cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count
     max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
                                        img.shape[1] * img.shape[1])))
-    out = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
-    bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance)
+    accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64)
     offset = max_distance / 2
     # find the nonzero indexes
     cdef np.ndarray[ndim=1, dtype=np.int_t] x_idxs, y_idxs
-    y_idxs, x_idxs =  np.PyArray_Nonzero(img)
+    y_idxs, x_idxs =  np.nonzero(img)
     num_indexes = y_idxs.shape[0] # x and y are the same shape
     nthetas = theta.shape[0]
+    
+    points = []
+    for i in range(num_indexes):
+        points.append((x_idxs[i], y_idxs[i]))
     lines = []
     # create mask of all non-zero indexes
     for i in range(num_indexes):
         mask[y_idxs[i], x_idxs[i]] = 1
 
-    for i in range(num_indexes):
+    while 1:
+    #for i in range(num_indexes):
         # select random non-zero point
-        index = randint(0, num_indexes-1)
-        x = x_idxs[i]
-        y = y_idxs[i]
+        count = len(points)
+        if count == 0:
+            break        
+        index = randint(0, count-1)
+        x = points[index][0]
+        y = points[index][1]
+        del points[index]
         # if previously eliminated, skip
         if not mask[y, x]:
             continue
@@ -119,9 +127,9 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, i
         max_theta = 0
         # apply hough transform on point
         for j in range(nthetas):
-            out_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
-            out[out_idx, j] += 1
-            value = out[out_idx, j]
+            accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
+            accum[accum_idx, j] += 1
+            value = accum[accum_idx, j]
             if value > max_value:
                 max_value = value
                 max_theta = j
@@ -181,8 +189,8 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, i
                 px += dx
                 py += dy            
         # confirm line length is sufficient
-        good_line = fabs(line_end[1, 1] - line_end[0, 1]) >= line_length or \
-                    fabs(line_end[1, 0] - line_end[0, 0]) >= line_length
+        good_line = abs(line_end[1, 1] - line_end[0, 1]) >= line_length or \
+                    abs(line_end[1, 0] - line_end[0, 0]) >= line_length
         # pass 2: walk the line again and reset accumulator and mask
         for k in range(2):
             px = x0
@@ -200,15 +208,22 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, i
                     x1 = px >> shift
                     y1 = py
                 # if non-zero point found, continue the line
-                if mask[y1, x1]:
-                    if good_line:
-                        for j in range(nthetas):
-                            out_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
-                            out[out_idx, j] -= 1
-                    mask[y1, x1] = 0
+                if 1:
+                    if mask[y1, x1]:  
+                        if good_line:                  
+                            accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset                    
+                            accum[accum_idx, max_theta] -= 1
+                        mask[y1, x1] = 0
+                else:
+                    if mask[y1, x1]:
+                        if good_line:
+                            for j in range(nthetas):
+                                accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
+                                accum[accum_idx, j] -= 1
+                        mask[y1, x1] = 0
                 # exit when the point is the line end
                 if x1 == line_end[k, 0] and y1 == line_end[k, 1]:
-                    break;
+                    break
                 px += dx
                 py += dy
 
diff --git a/scikits/image/transform/hough_transform.py b/scikits/image/transform/hough_transform.py
index 95d4c46b..d530fa0f 100644
--- a/scikits/image/transform/hough_transform.py
+++ b/scikits/image/transform/hough_transform.py
@@ -59,7 +59,7 @@ except ImportError:
     pass
 
 
-def probabilistic_hough(img, value_threshold=50, line_length=50, line_gap=10, theta=None):
+def probabilistic_hough(img, value_threshold=10, line_length=50, line_gap=10, theta=None):
     """Performs a progressive probabilistic line Hough transform and returns the detected lines.
 
     Parameters

From 4718e6e1d0adf016627d118df90d63d82f73b7d2 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 22 Aug 2011 15:32:44 +0200
Subject: [PATCH 03/14] Merge

---
 .../transform/tests/test_hough_transform.py    | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/scikits/image/transform/tests/test_hough_transform.py b/scikits/image/transform/tests/test_hough_transform.py
index 5cf10013..d4bd380d 100644
--- a/scikits/image/transform/tests/test_hough_transform.py
+++ b/scikits/image/transform/tests/test_hough_transform.py
@@ -3,6 +3,7 @@ from numpy.testing import *
 
 import scikits.image.transform as tf
 import scikits.image.transform.hough_transform as ht
+from scikits.image.transform import probabilistic_hough
 
 def append_desc(func, description):
     """Append the test function ``func`` and append
@@ -43,6 +44,23 @@ def test_py_hough():
 
     tf._hough = fast_hough
 
+def test_probabilistic_hough():
+    # Generate a test image
+    img = np.zeros((100, 100), dtype=int)
+    for i in range(25, 75):
+        img[100 - i, i] = 100
+        img[i, i] = 100
+    lines = probabilistic_hough(img, 10, line_length=10, line_gap=1)
+    # sort the lines according to the x-axis
+    sorted_lines = []
+    for line in lines:
+        line = list(line)
+        line.sort(lambda x,y: cmp(x[0], y[0]))
+        sorted_lines.append(line)
+    assert([(25, 75), (74, 26)] in sorted_lines)
+    assert([(25, 25), (74, 74)] in sorted_lines)
+
+
 if __name__ == "__main__":
     run_module_suite()
 

From 055358bde786cda29fbe82d822f3d3bbad8560a8 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 16:35:07 +0200
Subject: [PATCH 04/14] Hough tutorial

---
 doc/source/tutorials/hough_transform.txt | 114 +++++++++++++++++++++++
 1 file changed, 114 insertions(+)
 create mode 100644 doc/source/tutorials/hough_transform.txt

diff --git a/doc/source/tutorials/hough_transform.txt b/doc/source/tutorials/hough_transform.txt
new file mode 100644
index 00000000..be31d172
--- /dev/null
+++ b/doc/source/tutorials/hough_transform.txt
@@ -0,0 +1,114 @@
+***************
+Hough transform
+***************
+The Hough transform in its simplest form is a method to detect straight lines.
+
+http://en.wikipedia.org/wiki/Hough_transform
+
+For a first example we construct an example of intersecting lines.
+
+In [1]: import numpy as np
+
+In [2]: from scikits.image.transform import hough, probabilistic_hough
+
+In [3]: import matplotlib.pyplot as plt
+
+In [4]: from matplotlib.lines import Line2D                        
+
+In [5]: image = np.zeros((100, 100))
+
+In [6]: for i in range(25, 75):
+   ...:         image[100 - i, i] = 255
+   ...:     image[i, i] = 255
+   ...: 
+
+In [7]: plt.imshow(image)
+
+In [8]: plt.show()
+
+
+The Hough transform converts the image into a parameter space that represents 
+lines. A line can be represented by the distance r of its closest point to the 
+origin and by the angle theta of this vector.
+ 
+Every non-zero pixel of the image votes for potential line candidates, and the 
+local maxima represents the parameters of probable lines.
+
+In [9]: h, theta, d = hough(image)
+
+In [10]: plt.figure()
+
+In [10]: plt.title("hough transform")
+
+In [10]: plt.xlabel("degrees")
+
+In [10]: plt.ylabel("distance")
+
+In [11]: plt.imshow(h)
+
+In [12]: plt.show()
+
+As can be seen, the maxima occur at 45 and 135 degrees, corresponding to the 
+normal vector angles of each line.
+
+Another method is to use probabilistic_hough, an implementation
+based on the Progressive Probabilistic Hough Transform [1]. It states that a 
+random subset of voting points give good enough results, and that lines can
+be extracted during the voting process by walking along connected components.
+
+The function has three parameters: a general threshold that is applied to 
+the hough accumulator, a minimum line length and the line gap that influences
+line merging.
+
+In [13]: lines = probabilistic_hough(image, threshold=10, line_length=10, line_gap=1)
+
+In [14]: plt.figure()
+
+In [15]: for line in lines:
+   ....:         p0, p1 = line
+   ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
+   ....: 
+
+In [16]: plt.show()
+
+The Hough transform are often used on edge detected images.
+
+In [17]: from scikits.image.io import imread
+
+In [18]: from scikits.image import data_dir
+
+In [19]: from scikits.image.filter import canny
+
+In [20]: image = imread(data_dir + "/camera.png")
+
+In [21]: edges = canny(image, 2, 1, 25)
+
+In [22]: plt.imshow(edges)
+
+In [23]: plt.show()
+
+Apply the Probabilistic Hough Transform and find lines longer than 10 with a
+gap less than 3 pixels.
+
+In [24]: plt.figure()
+
+In [25]: plt.imshow(np.zeros(edges.shape))
+
+In [26]: lines = probabilistic_hough(edges, threshold=1, line_length=10, line_gap=3)
+
+In [27]: for line in lines:
+   ....:         p0, p1 = line
+   ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
+   ....: 
+
+In [28]: plt.show()
+
+
+References
+----------
+.. [1]  C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic Hough 
+        transform for line detection", in  IEEE Computer Society Conference on 
+        Computer Vision and Pattern Recognition, 1999. 
+.. [2]  Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to Detect 
+        Lines and Curves in Pictures," Comm. ACM, Vol. 15, pp. 11–15 (January, 
+        1972)

From 54a9deae476c2472b564016bd64d63b31f16173c Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 16:36:01 +0200
Subject: [PATCH 05/14] Formatting

---
 scikits/image/transform/_hough_transform.pyx | 9 ++++-----
 scikits/image/transform/hough_transform.py   | 4 ++--
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index 98693dbe..9dcd7205 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -19,7 +19,7 @@ cdef double round(double val):
 cdef double PI_2 = 1.5707963267948966
 cdef double NEG_PI_2 = -PI_2
 
-
+@cython.cdivision(True)
 @cython.boundscheck(False)
 def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     
@@ -66,8 +66,8 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
 
 
 @cython.boundscheck(False)
-def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, int line_gap, \
-   np.ndarray[ndim=1, dtype=np.double_t] theta=None):
+def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
+        int line_gap, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     if img.ndim != 2:
         raise ValueError('The input image must be 2D.')
     # compute the array of angles and their sine and cosine
@@ -110,7 +110,6 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, i
         mask[y_idxs[i], x_idxs[i]] = 1
 
     while 1:
-    #for i in range(num_indexes):
         # select random non-zero point
         count = len(points)
         if count == 0:
@@ -141,7 +140,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, i
         b = ctheta[max_theta]
         x0 = x
         y0 = y
-       # calculate gradient of walks using fixed point math
+        # calculate gradient of walks using fixed point math
         xflag = fabs(a) > fabs(b)
         if xflag:
             if a > 0:
diff --git a/scikits/image/transform/hough_transform.py b/scikits/image/transform/hough_transform.py
index d530fa0f..45d44599 100644
--- a/scikits/image/transform/hough_transform.py
+++ b/scikits/image/transform/hough_transform.py
@@ -59,7 +59,7 @@ except ImportError:
     pass
 
 
-def probabilistic_hough(img, value_threshold=10, line_length=50, line_gap=10, theta=None):
+def probabilistic_hough(img, threshold=10, line_length=50, line_gap=10, theta=None):
     """Performs a progressive probabilistic line Hough transform and returns the detected lines.
 
     Parameters
@@ -84,7 +84,7 @@ def probabilistic_hough(img, value_threshold=10, line_length=50, line_gap=10, th
             transform for line detection", in  IEEE Computer Society Conference on 
             Computer Vision and Pattern Recognition, 1999. 
     """    
-    return _probabilistic_hough(img, value_threshold, line_length, line_gap, theta)
+    return _probabilistic_hough(img, threshold, line_length, line_gap, theta)
 
 
 def hough(img, theta=None):

From fb45730c785608ed79056edaba84bc3e0f976e70 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 16:41:27 +0200
Subject: [PATCH 06/14] Tutorial work

---
 doc/source/tutorials/hough_transform.txt | 104 +++++++++++++----------
 1 file changed, 59 insertions(+), 45 deletions(-)

diff --git a/doc/source/tutorials/hough_transform.txt b/doc/source/tutorials/hough_transform.txt
index be31d172..1c215ae0 100644
--- a/doc/source/tutorials/hough_transform.txt
+++ b/doc/source/tutorials/hough_transform.txt
@@ -5,26 +5,28 @@ The Hough transform in its simplest form is a method to detect straight lines.
 
 http://en.wikipedia.org/wiki/Hough_transform
 
-For a first example we construct an example of intersecting lines.
+As a first example we construct a line intersection.
 
-In [1]: import numpy as np
+.. ipython::
 
-In [2]: from scikits.image.transform import hough, probabilistic_hough
+    In [1]: import numpy as np
 
-In [3]: import matplotlib.pyplot as plt
+    In [2]: from scikits.image.transform import hough, probabilistic_hough
 
-In [4]: from matplotlib.lines import Line2D                        
+    In [3]: import matplotlib.pyplot as plt
 
-In [5]: image = np.zeros((100, 100))
+    In [4]: from matplotlib.lines import Line2D                        
 
-In [6]: for i in range(25, 75):
-   ...:         image[100 - i, i] = 255
-   ...:     image[i, i] = 255
-   ...: 
+    In [5]: image = np.zeros((100, 100))
 
-In [7]: plt.imshow(image)
+    In [6]: for i in range(25, 75):
+       ...:         image[100 - i, i] = 255
+       ...:     image[i, i] = 255
+       ...: 
 
-In [8]: plt.show()
+    In [7]: plt.imshow(image)
+
+    In [8]: plt.show()
 
 
 The Hough transform converts the image into a parameter space that represents 
@@ -34,81 +36,93 @@ origin and by the angle theta of this vector.
 Every non-zero pixel of the image votes for potential line candidates, and the 
 local maxima represents the parameters of probable lines.
 
-In [9]: h, theta, d = hough(image)
+.. ipython::
 
-In [10]: plt.figure()
+    In [9]: h, theta, d = hough(image)
 
-In [10]: plt.title("hough transform")
+    In [10]: plt.figure()
 
-In [10]: plt.xlabel("degrees")
+    In [10]: plt.title("hough transform")
 
-In [10]: plt.ylabel("distance")
+    In [10]: plt.xlabel("degrees")
 
-In [11]: plt.imshow(h)
+    In [10]: plt.ylabel("distance")
+
+    In [11]: plt.imshow(h)
+
+    In [12]: plt.show()
 
-In [12]: plt.show()
 
 As can be seen, the maxima occur at 45 and 135 degrees, corresponding to the 
 normal vector angles of each line.
 
-Another method is to use probabilistic_hough, an implementation
+Another method is to use the function probabilistic_hough, an implementation
 based on the Progressive Probabilistic Hough Transform [1]. It states that a 
 random subset of voting points give good enough results, and that lines can
 be extracted during the voting process by walking along connected components.
+This returns the beginning and end of line segments, which are useful.
 
 The function has three parameters: a general threshold that is applied to 
-the hough accumulator, a minimum line length and the line gap that influences
+the Hough accumulator, a minimum line length and the line gap that influences
 line merging.
 
-In [13]: lines = probabilistic_hough(image, threshold=10, line_length=10, line_gap=1)
+.. ipython::
 
-In [14]: plt.figure()
+    In [13]: lines = probabilistic_hough(image, threshold=10, line_length=10, line_gap=1)
 
-In [15]: for line in lines:
-   ....:         p0, p1 = line
-   ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
-   ....: 
+    In [14]: plt.figure()
+
+    In [15]: for line in lines:
+       ....:         p0, p1 = line
+       ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
+       ....: 
+
+    In [16]: plt.show()
 
-In [16]: plt.show()
 
 The Hough transform are often used on edge detected images.
 
-In [17]: from scikits.image.io import imread
+.. ipython::
 
-In [18]: from scikits.image import data_dir
+    In [17]: from scikits.image.io import imread
 
-In [19]: from scikits.image.filter import canny
+    In [18]: from scikits.image import data_dir
 
-In [20]: image = imread(data_dir + "/camera.png")
+    In [19]: from scikits.image.filter import canny
 
-In [21]: edges = canny(image, 2, 1, 25)
+    In [20]: image = imread(data_dir + "/camera.png")
 
-In [22]: plt.imshow(edges)
+    In [21]: edges = canny(image, 2, 1, 25)
+
+    In [22]: plt.imshow(edges)
+
+    In [23]: plt.show()
 
-In [23]: plt.show()
 
 Apply the Probabilistic Hough Transform and find lines longer than 10 with a
 gap less than 3 pixels.
 
-In [24]: plt.figure()
+.. ipython::
 
-In [25]: plt.imshow(np.zeros(edges.shape))
+    In [24]: plt.figure()
 
-In [26]: lines = probabilistic_hough(edges, threshold=1, line_length=10, line_gap=3)
+    In [25]: plt.imshow(np.zeros(edges.shape))
 
-In [27]: for line in lines:
-   ....:         p0, p1 = line
-   ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
-   ....: 
+    In [26]: lines = probabilistic_hough(edges, threshold=1, line_length=10, line_gap=3)
 
-In [28]: plt.show()
+    In [27]: for line in lines:
+       ....:         p0, p1 = line
+       ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
+       ....: 
+
+    In [28]: plt.show()
 
 
 References
 ----------
-.. [1]  C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic Hough 
+..  [1]  C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic Hough 
         transform for line detection", in  IEEE Computer Society Conference on 
         Computer Vision and Pattern Recognition, 1999. 
-.. [2]  Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to Detect 
+    [2]  Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to Detect 
         Lines and Curves in Pictures," Comm. ACM, Vol. 15, pp. 11–15 (January, 
         1972)

From 65be497b0a80bf462181f21b25f496c14dc0868f Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Sat, 20 Aug 2011 01:49:27 +0200
Subject: [PATCH 07/14] Tutorial, adjusted threshold

---
 doc/source/tutorials/hough_transform.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/source/tutorials/hough_transform.txt b/doc/source/tutorials/hough_transform.txt
index 1c215ae0..9825510e 100644
--- a/doc/source/tutorials/hough_transform.txt
+++ b/doc/source/tutorials/hough_transform.txt
@@ -108,7 +108,7 @@ gap less than 3 pixels.
 
     In [25]: plt.imshow(np.zeros(edges.shape))
 
-    In [26]: lines = probabilistic_hough(edges, threshold=1, line_length=10, line_gap=3)
+    In [26]: lines = probabilistic_hough(edges, threshold=10, line_length=5, line_gap=3)
 
     In [27]: for line in lines:
        ....:         p0, p1 = line

From 06f91af36a781e369de94359914b41a6f4ae7ae2 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Sat, 20 Aug 2011 01:49:43 +0200
Subject: [PATCH 08/14] Improved line detection

---
 scikits/image/transform/_hough_transform.pyx | 44 +++++++++-----------
 1 file changed, 20 insertions(+), 24 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index 9dcd7205..e229008c 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -19,7 +19,6 @@ cdef double round(double val):
 cdef double PI_2 = 1.5707963267948966
 cdef double NEG_PI_2 = -PI_2
 
-@cython.cdivision(True)
 @cython.boundscheck(False)
 def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     
@@ -64,7 +63,9 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
             accum[accum_idx, j] += 1
     return accum, theta, bins
 
+import math
 
+@cython.cdivision(True)
 @cython.boundscheck(False)
 def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         int line_gap, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
@@ -75,7 +76,9 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     cdef np.ndarray[ndim=1, dtype=np.double_t] stheta
     # calculate thetas if none specified
     if theta is None:
-        theta = np.linspace(PI_2, NEG_PI_2, 180)
+        theta = np.linspace(math.pi/2, -math.pi/2, 180)
+        #p_2 = math.pi/2
+        #theta = p_2-np.arange(180)/180.0*p_2*2
     ctheta = np.cos(theta)
     stheta = np.sin(theta)
     cdef int height = img.shape[0]
@@ -93,6 +96,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count
     max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
                                        img.shape[1] * img.shape[1])))
+    #max_distance = (img.shape[0] + img.shape[1]) * 2 
     accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64)
     offset = max_distance / 2
     # find the nonzero indexes
@@ -108,7 +112,6 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     # create mask of all non-zero indexes
     for i in range(num_indexes):
         mask[y_idxs[i], x_idxs[i]] = 1
-
     while 1:
         # select random non-zero point
         count = len(points)
@@ -117,22 +120,22 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         index = randint(0, count-1)
         x = points[index][0]
         y = points[index][1]
-        del points[index]
+        del points[index]        
         # if previously eliminated, skip
         if not mask[y, x]:
             continue
         value = 0
-        max_value = 0
-        max_theta = 0
+        max_value = value_threshold-1
+        max_theta = -1
+        
         # apply hough transform on point
         for j in range(nthetas):
             accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
             accum[accum_idx, j] += 1
-            value = accum[accum_idx, j]
+            value = accum[accum_idx, j]            
             if value > max_value:
                 max_value = value
                 max_theta = j
-        # accumulator value of point strong enough
         if max_value < value_threshold:
             continue
         # from the random point walk in opposite directions and find line beginning and end
@@ -147,15 +150,16 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
                 dx0 = 1
             else:
                 dx0 = -1
-            dy0 = <int>round(b*(1 << shift)/fabs(a) ) 
+            dy0 = <int>round(b * (1 << shift) / fabs(a)) 
             y0 = (y0 << shift) + (1 << (shift - 1))
         else:
             if b > 0:
                 dy0 = 1
             else:
                 dy0 = -1
-            dx0 = <int>round( a*(1 << shift)/fabs(b))
-            x0 = (x0 << shift) + (1 << (shift-1))
+            dx0 = <int>round(a * (1 << shift) / fabs(b))
+            x0 = (x0 << shift) + (1 << (shift - 1))
+  
         # pass 1: walk the line, merging lines less than specified gap length
         for k in range(2):
             gap = 0
@@ -186,7 +190,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
                 elif gap > line_gap:
                     break
                 px += dx
-                py += dy            
+                py += dy
         # confirm line length is sufficient
         good_line = abs(line_end[1, 1] - line_end[0, 1]) >= line_length or \
                     abs(line_end[1, 0] - line_end[0, 0]) >= line_length
@@ -207,18 +211,10 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
                     x1 = px >> shift
                     y1 = py
                 # if non-zero point found, continue the line
-                if 1:
-                    if mask[y1, x1]:  
-                        if good_line:                  
-                            accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset                    
-                            accum[accum_idx, max_theta] -= 1
-                        mask[y1, x1] = 0
-                else:
-                    if mask[y1, x1]:
-                        if good_line:
-                            for j in range(nthetas):
-                                accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
-                                accum[accum_idx, j] -= 1
+                if mask[y1, x1]:  
+                    if good_line:                  
+                        accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset                    
+                        accum[accum_idx, max_theta] -= 1
                         mask[y1, x1] = 0
                 # exit when the point is the line end
                 if x1 == line_end[k, 0] and y1 == line_end[k, 1]:

From 9ed82321c02ef3c2a5cb9fff3e3c1e0aa0a0dc22 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Sat, 20 Aug 2011 02:34:56 +0200
Subject: [PATCH 09/14] Added a settle count for similar angles to resolve
 themselves

---
 scikits/image/transform/_hough_transform.pyx  | 25 +++++++++++++------
 .../transform/tests/test_hough_transform.py   | 20 ++++++++-------
 2 files changed, 28 insertions(+), 17 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index e229008c..d1737f30 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -77,8 +77,8 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     # calculate thetas if none specified
     if theta is None:
         theta = np.linspace(math.pi/2, -math.pi/2, 180)
-        #p_2 = math.pi/2
-        #theta = p_2-np.arange(180)/180.0*p_2*2
+        p_2 = math.pi/2
+        theta = p_2-np.arange(180)/180.0*p_2*2
     ctheta = np.cos(theta)
     stheta = np.sin(theta)
     cdef int height = img.shape[0]
@@ -87,7 +87,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     cdef np.ndarray[ndim=2, dtype=np.int64_t] accum
     cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = np.zeros((height, width), dtype=np.uint8)
     cdef np.ndarray[ndim=2, dtype=np.int32_t] line_end = np.zeros((2, 2), dtype=np.int32)
-    cdef int max_distance, offset, num_indexes, index    
+    cdef int max_distance, offset, num_indexes, index, other_max, settle
     cdef double a, b
     cdef int nidxs, nthetas, i, j, x, y, px, py, accum_idx, value, max_value, max_theta
     cdef int shift = 16
@@ -112,6 +112,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     # create mask of all non-zero indexes
     for i in range(num_indexes):
         mask[y_idxs[i], x_idxs[i]] = 1
+    settle = 0
     while 1:
         # select random non-zero point
         count = len(points)
@@ -127,17 +128,25 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         value = 0
         max_value = value_threshold-1
         max_theta = -1
-        
+        other_max = 0
         # apply hough transform on point
         for j in range(nthetas):
             accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
             accum[accum_idx, j] += 1
-            value = accum[accum_idx, j]            
+            value = accum[accum_idx, j]
             if value > max_value:
                 max_value = value
                 max_theta = j
-        if max_value < value_threshold:
+            # if two max angles are very similar, extend the threshold for a while
+            if j != max_theta and value == max_value and abs(j - max_theta) == 1:
+                settle += 1
+                if settle > 100:
+                    other_max = 0
+                else:
+                    other_max = 1
+        if max_value < value_threshold or other_max:
             continue
+        settle = 0
         # from the random point walk in opposite directions and find line beginning and end
         a = -stheta[max_theta]
         b = ctheta[max_theta]
@@ -151,14 +160,14 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
             else:
                 dx0 = -1
             dy0 = <int>round(b * (1 << shift) / fabs(a)) 
-            y0 = (y0 << shift) + (1 << (shift - 1))
+            y0 = (y0 << shift) #+ (1 << (shift - 1))
         else:
             if b > 0:
                 dy0 = 1
             else:
                 dy0 = -1
             dx0 = <int>round(a * (1 << shift) / fabs(b))
-            x0 = (x0 << shift) + (1 << (shift - 1))
+            x0 = (x0 << shift) #+ (1 << (shift - 1))
   
         # pass 1: walk the line, merging lines less than specified gap length
         for k in range(2):
diff --git a/scikits/image/transform/tests/test_hough_transform.py b/scikits/image/transform/tests/test_hough_transform.py
index d4bd380d..7582035e 100644
--- a/scikits/image/transform/tests/test_hough_transform.py
+++ b/scikits/image/transform/tests/test_hough_transform.py
@@ -50,15 +50,17 @@ def test_probabilistic_hough():
     for i in range(25, 75):
         img[100 - i, i] = 100
         img[i, i] = 100
-    lines = probabilistic_hough(img, 10, line_length=10, line_gap=1)
-    # sort the lines according to the x-axis
-    sorted_lines = []
-    for line in lines:
-        line = list(line)
-        line.sort(lambda x,y: cmp(x[0], y[0]))
-        sorted_lines.append(line)
-    assert([(25, 75), (74, 26)] in sorted_lines)
-    assert([(25, 25), (74, 74)] in sorted_lines)
+    # test the line extraction a few times
+    for i in range(100):
+        lines = probabilistic_hough(img, threshold=10, line_length=10, line_gap=1)
+        # sort the lines according to the x-axis
+        sorted_lines = []
+        for line in lines:
+            line = list(line)
+            line.sort(lambda x,y: cmp(x[0], y[0]))
+            sorted_lines.append(line)
+        assert([(25, 75), (74, 26)] in sorted_lines)
+        assert([(25, 25), (74, 74)] in sorted_lines)
 
 
 if __name__ == "__main__":

From e5ddf8590743b67fcafa273a68e59b830c4dca7e Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Sat, 20 Aug 2011 11:15:50 +0200
Subject: [PATCH 10/14] Debugging similar angles

---
 scikits/image/transform/_hough_transform.pyx | 27 ++++++++++++--------
 1 file changed, 16 insertions(+), 11 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index d1737f30..5d0d8b92 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -77,8 +77,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     # calculate thetas if none specified
     if theta is None:
         theta = np.linspace(math.pi/2, -math.pi/2, 180)
-        p_2 = math.pi/2
-        theta = p_2-np.arange(180)/180.0*p_2*2
+        theta =  math.pi/2-np.arange(180)/180.0* math.pi
     ctheta = np.cos(theta)
     stheta = np.sin(theta)
     cdef int height = img.shape[0]
@@ -113,6 +112,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     for i in range(num_indexes):
         mask[y_idxs[i], x_idxs[i]] = 1
     settle = 0
+    print
     while 1:
         # select random non-zero point
         count = len(points)
@@ -125,6 +125,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         # if previously eliminated, skip
         if not mask[y, x]:
             continue
+        print "select", y, x
         value = 0
         max_value = value_threshold-1
         max_theta = -1
@@ -135,18 +136,22 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
             accum[accum_idx, j] += 1
             value = accum[accum_idx, j]
             if value > max_value:
+                print "max", j, value
                 max_value = value
                 max_theta = j
-            # if two max angles are very similar, extend the threshold for a while
-            if j != max_theta and value == max_value and abs(j - max_theta) == 1:
-                settle += 1
-                if settle > 100:
-                    other_max = 0
-                else:
-                    other_max = 1
+            # if two bordering angles have max values, extend the threshold for a while
+            if j != max_theta and value == max_value:
+                if abs(j - max_theta) <= 2:
+                    if settle > 100:
+                        other_max = 0
+                    else:
+                        other_max = 1
+            if j == 134 or j == 135:
+                print j, value, accum_idx, ctheta[j], stheta[j]
         if max_value < value_threshold or other_max:
             continue
         settle = 0
+        print "THETA", max_theta, max_value
         # from the random point walk in opposite directions and find line beginning and end
         a = -stheta[max_theta]
         b = ctheta[max_theta]
@@ -160,14 +165,14 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
             else:
                 dx0 = -1
             dy0 = <int>round(b * (1 << shift) / fabs(a)) 
-            y0 = (y0 << shift) #+ (1 << (shift - 1))
+            y0 = (y0 << shift) + (1 << (shift - 1))
         else:
             if b > 0:
                 dy0 = 1
             else:
                 dy0 = -1
             dx0 = <int>round(a * (1 << shift) / fabs(b))
-            x0 = (x0 << shift) #+ (1 << (shift - 1))
+            x0 = (x0 << shift) + (1 << (shift - 1))
   
         # pass 1: walk the line, merging lines less than specified gap length
         for k in range(2):

From 9cd072d92502818c4d2ee5fef70e69bf215a98e9 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 22 Aug 2011 15:35:57 +0200
Subject: [PATCH 11/14] Merge

---
 scikits/image/transform/_hough_transform.pyx  | 20 ++------------
 .../transform/tests/test_hough_transform.py   | 27 +++++++++++--------
 2 files changed, 18 insertions(+), 29 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index 5d0d8b92..eeb4bb0d 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -86,7 +86,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     cdef np.ndarray[ndim=2, dtype=np.int64_t] accum
     cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = np.zeros((height, width), dtype=np.uint8)
     cdef np.ndarray[ndim=2, dtype=np.int32_t] line_end = np.zeros((2, 2), dtype=np.int32)
-    cdef int max_distance, offset, num_indexes, index, other_max, settle
+    cdef int max_distance, offset, num_indexes, index
     cdef double a, b
     cdef int nidxs, nthetas, i, j, x, y, px, py, accum_idx, value, max_value, max_theta
     cdef int shift = 16
@@ -111,8 +111,6 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     # create mask of all non-zero indexes
     for i in range(num_indexes):
         mask[y_idxs[i], x_idxs[i]] = 1
-    settle = 0
-    print
     while 1:
         # select random non-zero point
         count = len(points)
@@ -125,33 +123,19 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         # if previously eliminated, skip
         if not mask[y, x]:
             continue
-        print "select", y, x
         value = 0
         max_value = value_threshold-1
         max_theta = -1
-        other_max = 0
         # apply hough transform on point
         for j in range(nthetas):
             accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
             accum[accum_idx, j] += 1
             value = accum[accum_idx, j]
             if value > max_value:
-                print "max", j, value
                 max_value = value
                 max_theta = j
-            # if two bordering angles have max values, extend the threshold for a while
-            if j != max_theta and value == max_value:
-                if abs(j - max_theta) <= 2:
-                    if settle > 100:
-                        other_max = 0
-                    else:
-                        other_max = 1
-            if j == 134 or j == 135:
-                print j, value, accum_idx, ctheta[j], stheta[j]
-        if max_value < value_threshold or other_max:
+        if max_value < value_threshold:
             continue
-        settle = 0
-        print "THETA", max_theta, max_value
         # from the random point walk in opposite directions and find line beginning and end
         a = -stheta[max_theta]
         b = ctheta[max_theta]
diff --git a/scikits/image/transform/tests/test_hough_transform.py b/scikits/image/transform/tests/test_hough_transform.py
index 7582035e..3d185521 100644
--- a/scikits/image/transform/tests/test_hough_transform.py
+++ b/scikits/image/transform/tests/test_hough_transform.py
@@ -13,6 +13,9 @@ def append_desc(func, description):
 
     return func
 
+from scikits.image.transform import *
+import math
+
 def test_hough():
     # Generate a test image
     img = np.zeros((100, 100), dtype=int)
@@ -28,6 +31,7 @@ def test_hough():
     assert_equal(dist > 70, dist < 72)
     assert_equal(theta > 0.78, theta < 0.79)
 
+
 def test_hough_angles():
     img = np.zeros((10, 10))
     img[0, 0] = 1
@@ -50,17 +54,18 @@ def test_probabilistic_hough():
     for i in range(25, 75):
         img[100 - i, i] = 100
         img[i, i] = 100
-    # test the line extraction a few times
-    for i in range(100):
-        lines = probabilistic_hough(img, threshold=10, line_length=10, line_gap=1)
-        # sort the lines according to the x-axis
-        sorted_lines = []
-        for line in lines:
-            line = list(line)
-            line.sort(lambda x,y: cmp(x[0], y[0]))
-            sorted_lines.append(line)
-        assert([(25, 75), (74, 26)] in sorted_lines)
-        assert([(25, 25), (74, 74)] in sorted_lines)
+    # decrease default theta sampling because similar orientations may confuse
+    # as mentioned in article of Galambos et al
+    theta=np.linspace(0, math.pi, 45)
+    lines = probabilistic_hough(img, theta=theta, threshold=10, line_length=10, line_gap=1)
+    # sort the lines according to the x-axis
+    sorted_lines = []
+    for line in lines:
+        line = list(line)
+        line.sort(lambda x,y: cmp(x[0], y[0]))
+        sorted_lines.append(line)
+    assert([(25, 75), (74, 26)] in sorted_lines)
+    assert([(25, 25), (74, 74)] in sorted_lines)
 
 
 if __name__ == "__main__":

From 6df94d2fdd1bf05cfa209f2d9478ccedfd4b217b Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Sat, 20 Aug 2011 11:45:07 +0200
Subject: [PATCH 12/14] Tutorial save figs

---
 doc/source/tutorials/hough_transform.txt | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/doc/source/tutorials/hough_transform.txt b/doc/source/tutorials/hough_transform.txt
index 9825510e..be4d7571 100644
--- a/doc/source/tutorials/hough_transform.txt
+++ b/doc/source/tutorials/hough_transform.txt
@@ -25,7 +25,8 @@ As a first example we construct a line intersection.
        ...: 
 
     In [7]: plt.imshow(image)
-
+    
+    @savefig hough_original.png width=4in
     In [8]: plt.show()
 
 
@@ -50,6 +51,7 @@ local maxima represents the parameters of probable lines.
 
     In [11]: plt.imshow(h)
 
+    @savefig hough_transform.png width=4in
     In [12]: plt.show()
 
 
@@ -77,6 +79,7 @@ line merging.
        ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
        ....: 
 
+    @savefig hough_probabilistic1.png width=4in
     In [16]: plt.show()
 
 
@@ -96,6 +99,7 @@ The Hough transform are often used on edge detected images.
 
     In [22]: plt.imshow(edges)
 
+    @savefig hough_edge_detected.png width=4in
     In [23]: plt.show()
 
 
@@ -115,6 +119,7 @@ gap less than 3 pixels.
        ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
        ....: 
 
+    @savefig hough_lines.png width=4in
     In [28]: plt.show()
 
 

From f5225d2369cb2ea5b60ab83e89c303b80c223de0 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 22 Aug 2011 15:20:41 +0200
Subject: [PATCH 13/14] Fixes

---
 scikits/image/transform/_hough_transform.pyx      | 10 +++++-----
 scikits/image/transform/hough_transform.py        | 15 ++++++++++-----
 .../image/transform/tests/test_hough_transform.py |  3 +--
 3 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index eeb4bb0d..ddcf91f6 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -1,11 +1,13 @@
 cimport cython
-
 import numpy as np
 cimport numpy as np
 from random import randint
 np.import_array()
 
-
+cdef extern from "stdlib.h":
+    int rand()
+    void randomize()
+        
 cdef extern from "math.h":
     int abs(int)
     double fabs(double)
@@ -95,7 +97,6 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count
     max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
                                        img.shape[1] * img.shape[1])))
-    #max_distance = (img.shape[0] + img.shape[1]) * 2 
     accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64)
     offset = max_distance / 2
     # find the nonzero indexes
@@ -103,7 +104,6 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     y_idxs, x_idxs =  np.nonzero(img)
     num_indexes = y_idxs.shape[0] # x and y are the same shape
     nthetas = theta.shape[0]
-    
     points = []
     for i in range(num_indexes):
         points.append((x_idxs[i], y_idxs[i]))
@@ -116,7 +116,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         count = len(points)
         if count == 0:
             break        
-        index = randint(0, count-1)
+        index = rand() % (count)
         x = points[index][0]
         y = points[index][1]
         del points[index]        
diff --git a/scikits/image/transform/hough_transform.py b/scikits/image/transform/hough_transform.py
index 45d44599..31ffa565 100644
--- a/scikits/image/transform/hough_transform.py
+++ b/scikits/image/transform/hough_transform.py
@@ -54,7 +54,7 @@ _py_hough = _hough
 
 # try to import and use the faster Cython version if it exists
 try:
-    from ._hough_transform import _hough  
+    from ._hough_transform import _hough
 except ImportError:
     pass
 
@@ -67,11 +67,16 @@ def probabilistic_hough(img, threshold=10, line_length=50, line_gap=10, theta=No
     img : (M, N) ndarray
         Input image with nonzero values representing edges.
     value_threshold: int
-        Threshold
-    theta :1D ndarray, dtype=double
+        Threshold        
+    line_length: int, optional (default 50)
+        Minimum accepted length of detected lines.
+        Increase the parameter to extract longer lines.
+    line_gap: int, optional, (default 10)
+        Maximum gap between pixels to still form a line. 
+        Increase the parameter to merge broken lines more aggresively.
+    theta :1D ndarray, dtype=double, optional, default (-pi/2 .. pi/2)
         Angles at which to compute the transform, in radians.
-        Defaults to -pi/2 .. pi/2
-
+        
     Returns
     -------
     lines : list
diff --git a/scikits/image/transform/tests/test_hough_transform.py b/scikits/image/transform/tests/test_hough_transform.py
index 3d185521..175c1b37 100644
--- a/scikits/image/transform/tests/test_hough_transform.py
+++ b/scikits/image/transform/tests/test_hough_transform.py
@@ -14,7 +14,6 @@ def append_desc(func, description):
     return func
 
 from scikits.image.transform import *
-import math
 
 def test_hough():
     # Generate a test image
@@ -56,7 +55,7 @@ def test_probabilistic_hough():
         img[i, i] = 100
     # decrease default theta sampling because similar orientations may confuse
     # as mentioned in article of Galambos et al
-    theta=np.linspace(0, math.pi, 45)
+    theta=np.linspace(0, np.pi, 45)
     lines = probabilistic_hough(img, theta=theta, threshold=10, line_length=10, line_gap=1)
     # sort the lines according to the x-axis
     sorted_lines = []

From ed011f95d1456f28ec4289a9be75a499e3efc7e3 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Wed, 24 Aug 2011 17:25:33 +0200
Subject: [PATCH 14/14] Removed unused definition

---
 scikits/image/transform/_hough_transform.pyx | 1 -
 1 file changed, 1 deletion(-)

diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index ddcf91f6..96c88ab1 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -6,7 +6,6 @@ np.import_array()
 
 cdef extern from "stdlib.h":
     int rand()
-    void randomize()
         
 cdef extern from "math.h":
     int abs(int)