diff --git a/doc/examples/plot_circular_elliptical_hough_transform.py b/doc/examples/plot_circular_elliptical_hough_transform.py
index 7ada848c..4c148b39 100755
--- a/doc/examples/plot_circular_elliptical_hough_transform.py
+++ b/doc/examples/plot_circular_elliptical_hough_transform.py
@@ -104,7 +104,6 @@ References
        Conference on. Vol. 2. IEEE, 2002
 """
 import matplotlib.pyplot as plt
-import heapq
 
 from skimage import data, filter, color
 from skimage.transform import hough_ellipse
@@ -120,9 +119,9 @@ edges = filter.canny(image_gray, sigma=2.0,
 # The accuracy corresponds to the bin size of a major axis.
 # The value is chosen in order to get a single high accumulator.
 # The threshold eliminates low accumulators
-accum = hough_ellipse(edges, accuracy=10, threshold=150, min_size=50)
-# Select the highest accumulator
-best = heapq.nlargest(1, accum)[0]
+accum = hough_ellipse(edges, accuracy=10, threshold=170, min_size=50)
+accum.sort(key=lambda x:x[0])
+best = accum[-1]
 # Estimated parameters for the ellipse
 center_y = int(best[1])
 center_x = int(best[2])
diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx
index b14b2d5b..47185d79 100644
--- a/skimage/transform/_hough_transform.pyx
+++ b/skimage/transform/_hough_transform.pyx
@@ -3,7 +3,6 @@
 #cython: nonecheck=False
 #cython: wraparound=False
 import numpy as np
-import heapq
 
 cimport numpy as cnp
 cimport cython
@@ -132,12 +131,8 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
     >>> img = np.zeros((25, 25), dtype=int)
     >>> rr, cc = draw.ellipse_perimeter(10, 10, 6, 8)
     >>> img[cc, rr] = 1
-    >>> result = hough_ellipse(img, threshold=4)
-    >>> # extract the highest accumulator
-    >>> heapq.nlargest(1, result)
-    [(10, 10.0, 10.0, 6.0, 8.0, 0.0)]
-    >>> # To sort them all
-    >>> results = [heappop(results) for i in range(len(results))]
+    >>> result = hough_ellipse(img, threshold=8)
+    [(10, 10.0, 8.0, 6.0, 0.0, 10.0)]
 
     Notes
     -----
@@ -221,8 +216,7 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
                             if angle > np.pi:
                                 angle = angle - np.pi / 2.
                                 a, b = b, a
-                        heapq.heappush(results,
-                                       (hist_max,  # Accumulator
+                        results.append((hist_max, # Accumulator
                                         y0,
                                         x0,
                                         a,
diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py
index 281d049d..82b538b2 100644
--- a/skimage/transform/tests/test_hough_transform.py
+++ b/skimage/transform/tests/test_hough_transform.py
@@ -1,5 +1,4 @@
 import numpy as np
-import heapq
 from numpy.testing import (assert_almost_equal,
                            assert_equal,
                            )
@@ -158,7 +157,7 @@ def test_hough_ellipse_zero_angle():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=9)
-    best = heapq.nlargest(1, result)[0]
+    best = result[-1]
     assert_equal(best[1], y0)
     assert_equal(best[2], x0)
     assert_almost_equal(best[3], ry, decimal=1)
@@ -182,7 +181,8 @@ def test_hough_ellipse_non_zero_posangle1():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     assert_almost_equal(best[1] / 100., y0 / 100., decimal=1)
     assert_almost_equal(best[2] / 100., x0 / 100., decimal=1)
     assert_almost_equal(best[3] / 10., ry / 10., decimal=1)
@@ -206,7 +206,8 @@ def test_hough_ellipse_non_zero_posangle2():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     assert_almost_equal(best[1] / 100., y0 / 100., decimal=1)
     assert_almost_equal(best[2] / 100., x0 / 100., decimal=1)
     assert_almost_equal(best[3] / 10., ry / 10., decimal=1)
@@ -230,7 +231,8 @@ def test_hough_ellipse_non_zero_posangle3():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     # Check if I re-draw the ellipse, points are the same!
     # ie check API compatibility between hough_ellipse and ellipse_perimeter
     rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), orientation=best[5])
@@ -249,7 +251,8 @@ def test_hough_ellipse_non_zero_posangle4():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     # Check if I re-draw the ellipse, points are the same!
     # ie check API compatibility between hough_ellipse and ellipse_perimeter
     rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), orientation=best[5])
@@ -268,7 +271,8 @@ def test_hough_ellipse_non_zero_negangle1():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     # Check if I re-draw the ellipse, points are the same!
     # ie check API compatibility between hough_ellipse and ellipse_perimeter
     rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), orientation=best[5])
@@ -287,7 +291,8 @@ def test_hough_ellipse_non_zero_negangle2():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     # Check if I re-draw the ellipse, points are the same!
     # ie check API compatibility between hough_ellipse and ellipse_perimeter
     rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), orientation=best[5])
@@ -306,7 +311,8 @@ def test_hough_ellipse_non_zero_negangle3():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     # Check if I re-draw the ellipse, points are the same!
     # ie check API compatibility between hough_ellipse and ellipse_perimeter
     rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), orientation=best[5])
@@ -325,7 +331,8 @@ def test_hough_ellipse_non_zero_negangle4():
     rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    best = heapq.nlargest(1, result)[0]
+    result.sort(key=lambda x:x[0])
+    best = result[-1]
     # Check if I re-draw the ellipse, points are the same!
     # ie check API compatibility between hough_ellipse and ellipse_perimeter
     rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), orientation=best[5])