diff --git a/doc/source/tutorials/hough_transform.txt b/doc/source/tutorials/hough_transform.txt
new file mode 100644
index 00000000..be4d7571
--- /dev/null
+++ b/doc/source/tutorials/hough_transform.txt
@@ -0,0 +1,133 @@
+***************
+Hough transform
+***************
+The Hough transform in its simplest form is a method to detect straight lines.
+
+http://en.wikipedia.org/wiki/Hough_transform
+
+As a first example we construct a line intersection.
+
+.. ipython::
+
+    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)
+    
+    @savefig hough_original.png width=4in
+    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.
+
+.. ipython::
+
+    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)
+
+    @savefig hough_transform.png width=4in
+    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 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
+line merging.
+
+.. ipython::
+
+    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]))
+       ....: 
+
+    @savefig hough_probabilistic1.png width=4in
+    In [16]: plt.show()
+
+
+The Hough transform are often used on edge detected images.
+
+.. ipython::
+
+    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)
+
+    @savefig hough_edge_detected.png width=4in
+    In [23]: plt.show()
+
+
+Apply the Probabilistic Hough Transform and find lines longer than 10 with a
+gap less than 3 pixels.
+
+.. ipython::
+
+    In [24]: plt.figure()
+
+    In [25]: plt.imshow(np.zeros(edges.shape))
+
+    In [26]: lines = probabilistic_hough(edges, threshold=10, line_length=5, line_gap=3)
+
+    In [27]: for line in lines:
+       ....:         p0, p1 = line
+       ....:     plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
+       ....: 
+
+    @savefig hough_lines.png width=4in
+    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)
diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx
index ec85743d..96c88ab1 100644
--- a/scikits/image/transform/_hough_transform.pyx
+++ b/scikits/image/transform/_hough_transform.pyx
@@ -1,12 +1,15 @@
 cimport cython
-
 import numpy as np
 cimport numpy as np
-
+from random import randint
 np.import_array()
 
-
+cdef extern from "stdlib.h":
+    int rand()
+        
 cdef extern from "math.h":
+    int abs(int)
+    double fabs(double)
     double sqrt(double)
     double ceil(double)
     double floor(double)
@@ -17,7 +20,6 @@ cdef double round(double val):
 cdef double PI_2 = 1.5707963267948966
 cdef double NEG_PI_2 = -PI_2
 
-
 @cython.boundscheck(False)
 def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     
@@ -34,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
     
@@ -49,17 +51,179 @@ 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
+    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]
+        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
+            accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
+            accum[accum_idx, j] += 1
+    return accum, theta, bins
 
-    return out, 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):
+    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(math.pi/2, -math.pi/2, 180)
+        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]
+    cdef int width = img.shape[1]
+    # 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.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, 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, count
+    max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
+                                       img.shape[1] * img.shape[1])))
+    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.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
+    while 1:
+        # select random non-zero point
+        count = len(points)
+        if count == 0:
+            break        
+        index = rand() % (count)
+        x = points[index][0]
+        y = points[index][1]
+        del points[index]        
+        # if previously eliminated, skip
+        if not mask[y, x]:
+            continue
+        value = 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]
+            if value > max_value:
+                max_value = value
+                max_theta = j
+        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 = 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
+            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:                  
+                        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]:
+                    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..31ffa565 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:
@@ -58,6 +59,39 @@ except ImportError:
     pass
 
 
+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
+    ----------
+    img : (M, N) ndarray
+        Input image with nonzero values representing edges.
+    value_threshold: int
+        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.
+        
+    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, threshold, line_length, line_gap, theta)
+
+
 def hough(img, theta=None):
     """Perform a straight line Hough transform.
 
@@ -67,7 +101,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 +140,5 @@ def hough(img, theta=None):
 
     """
     return _hough(img, theta)
+
+
diff --git a/scikits/image/transform/tests/test_hough_transform.py b/scikits/image/transform/tests/test_hough_transform.py
index 5cf10013..175c1b37 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
@@ -12,6 +13,8 @@ def append_desc(func, description):
 
     return func
 
+from scikits.image.transform import *
+
 def test_hough():
     # Generate a test image
     img = np.zeros((100, 100), dtype=int)
@@ -27,6 +30,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
@@ -43,6 +47,26 @@ 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
+    # decrease default theta sampling because similar orientations may confuse
+    # as mentioned in article of Galambos et al
+    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 = []
+    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()