diff --git a/doc/examples/plot_circular_elliptical_hough_transform.py b/doc/examples/plot_circular_elliptical_hough_transform.py
index bd036e03..7fb67046 100755
--- a/doc/examples/plot_circular_elliptical_hough_transform.py
+++ b/doc/examples/plot_circular_elliptical_hough_transform.py
@@ -74,7 +74,6 @@ for idx in np.argsort(accums)[::-1][:5]:
     image[cy, cx] = (220, 20, 20)
 
 ax.imshow(image, cmap=plt.cm.gray)
-plt.show()
 
 
 """
@@ -96,13 +95,13 @@ an ellipse passes to them. A good match corresponds to high accumulator values.
 
 A full description of the algorithm can be found in reference [1]_.
 
-
 References
 ----------
 .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection
        method." Pattern Recognition, 2002. Proceedings. 16th International
        Conference on. Vol. 2. IEEE, 2002
 """
+
 import matplotlib.pyplot as plt
 
 from skimage import data, filter, color
@@ -110,7 +109,7 @@ from skimage.transform import hough_ellipse
 from skimage.draw import ellipse_perimeter
 
 # Load picture, convert to grayscale and detect edges
-image_rgb = data.load('coffee.png')[0:220, 100:450]
+image_rgb = data.coffee()[0:220, 160:420]
 image_gray = color.rgb2gray(image_rgb)
 edges = filter.canny(image_gray, sigma=2.0,
                      low_threshold=0.55, high_threshold=0.8)
@@ -119,29 +118,31 @@ 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=170, min_size=50)
-accum.sort(key=lambda x:x[5])
+result = hough_ellipse(edges, accuracy=20, threshold=250,
+                       min_size=100, max_size=120)
+result.sort(order='accumulator')
+
 # Estimated parameters for the ellipse
-center_y = int(accum[-1][0])
-center_x = int(accum[-1][1])
-xradius = int(accum[-1][2])
-yradius = int(accum[-1][3])
-angle = np.pi - accum[-1][4]
+best = result[-1]
+yc = int(best[1])
+xc = int(best[2])
+a = int(best[3])
+b = int(best[4])
+orientation = best[5]
 
 # Draw the ellipse on the original image
-cx, cy = ellipse_perimeter(center_y, center_x,
-                           yradius, xradius, orientation=angle)
-image_rgb[cy, cx] = (0, 0, 1)
+cy, cx = ellipse_perimeter(yc, xc, a, b, orientation)
+image_rgb[cy, cx] = (0, 0, 255)
 # Draw the edge (white) and the resulting ellipse (red)
 edges = color.gray2rgb(edges)
 edges[cy, cx] = (250, 0, 0)
 
-fig = plt.subplots(figsize=(10, 6))
-plt.subplot(1, 2, 1)
-plt.title('Original picture')
-plt.imshow(image_rgb)
-plt.subplot(1, 2, 2)
-plt.title('Edge (white) and result (red)')
-plt.imshow(edges)
+fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 6))
+
+ax1.set_title('Original picture')
+ax1.imshow(image_rgb)
+
+ax2.set_title('Edge (white) and result (red)')
+ax2.imshow(edges)
 
 plt.show()
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index 3d48ce86..c22600a4 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -449,9 +449,9 @@ def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius,
     ----------
     cy, cx : int
         Centre coordinate of ellipse.
-    yradius, xradius: int
+    yradius, xradius : int
         Minor and major semi-axes. ``(x/xradius)**2 + (y/yradius)**2 = 1``.
-    orientation: double, optional (default 0)
+    orientation : double, optional (default 0)
         Major axis orientation in clockwise direction as radians.
 
     Returns
diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx
index e5fd3ed8..29344fa8 100644
--- a/skimage/transform/_hough_transform.pyx
+++ b/skimage/transform/_hough_transform.pyx
@@ -7,7 +7,7 @@ import numpy as np
 cimport numpy as cnp
 cimport cython
 
-from libc.math cimport abs, fabs, sqrt, ceil
+from libc.math cimport abs, fabs, sqrt, ceil, atan2, M_PI
 from libc.stdlib cimport rand
 
 from skimage.draw import circle_perimeter
@@ -122,17 +122,18 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
 
     Returns
     -------
-    res : list of tuples [(x0, y0, a, b, angle, accumulator)]
-          Where (x0, y0) is the center, (a, b) major and minor axis.
-          The angle value follows `draw.ellipse_perimeter()` convention.
+    result : ndarray with fields [(accumulator, y0, x0, a, b, orientation)]
+          Where ``(yc, xc)`` is the center, ``(a, b)`` the major and minor
+          axes, respectively. The `orientation` value follows
+          `skimage.draw.ellipse_perimeter` convention.
 
     Examples
     --------
-    >>> img = np.zeros((25, 25), dtype=int)
-    >>> rr, cc = draw.ellipse_perimeter(10, 10, 6, 8)
+    >>> img = np.zeros((25, 25), dtype=np.uint8)
+    >>> rr, cc = ellipse_perimeter(10, 10, 6, 8)
     >>> img[cc, rr] = 1
     >>> result = hough_ellipse(img, threshold=8)
-    [(10.0, 10.0, 8.0, 6.0, 0.0, 10)]
+    [(10, 10.0, 8.0, 6.0, 0.0, 10.0)]
 
     Notes
     -----
@@ -149,47 +150,47 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
     if img.ndim != 2:
             raise ValueError('The input image must be 2D.')
 
-    cdef Py_ssize_t[:, :] pixels = np.transpose(np.nonzero(img))
-    cdef Py_ssize_t num_pixels = pixels.shape[0]
+    cdef Py_ssize_t[:, ::1] pixels = np.row_stack(np.nonzero(img))
+    cdef Py_ssize_t num_pixels = pixels.shape[1]
     cdef list acc = list()
     cdef list results = list()
-    cdef bin_size = accuracy**2
+    cdef double bin_size = accuracy ** 2
 
     cdef int max_b_squared
     if max_size is None:
         if img.shape[0] < img.shape[1]:
-            max_b_squared = np.round(0.5 * img.shape[0])**2
+            max_b_squared = np.round(0.5 * img.shape[0]) ** 2
         else:
-            max_b_squared = np.round(0.5 * img.shape[1])**2
+            max_b_squared = np.round(0.5 * img.shape[1]) ** 2
     else:
         max_b_squared = max_size**2
 
     cdef Py_ssize_t p1, p2, p3, p1x, p1y, p2x, p2y, p3x, p3y
-    cdef double x0, y0, a, b, d, k
-    cdef double cos_tau_squared, b_squared, f_squared, angle
+    cdef double xc, yc, a, b, d, k
+    cdef double cos_tau_squared, b_squared, f_squared, orientation
 
     for p1 in range(num_pixels):
-        p1x = pixels[p1, 1]
-        p1y = pixels[p1, 0]
+        p1x = pixels[1, p1]
+        p1y = pixels[0, p1]
 
         for p2 in range(p1):
-            p2x = pixels[p2, 1]
-            p2y = pixels[p2, 0]
+            p2x = pixels[1, p2]
+            p2y = pixels[0, p2]
 
-            # Candidate: center (x0, y0) and main axis a
+            # Candidate: center (xc, yc) and main axis a
             a = 0.5 * sqrt((p1x - p2x)**2 + (p1y - p2y)**2)
             if a > 0.5 * min_size:
-                x0 = 0.5 * (p1x + p2x)
-                y0 = 0.5 * (p1y + p2y)
+                xc = 0.5 * (p1x + p2x)
+                yc = 0.5 * (p1y + p2y)
 
                 for p3 in range(num_pixels):
-                    p3x = pixels[p3, 1]
-                    p3y = pixels[p3, 0]
+                    p3x = pixels[1, p3]
+                    p3y = pixels[0, p3]
 
-                    d = sqrt((p3x - x0)**2 + (p3y - y0)**2)
+                    d = sqrt((p3x - xc)**2 + (p3y - yc)**2)
                     if d > min_size:
                         f_squared = (p3x - p1x)**2 + (p3y - p1y)**2
-                        cos_tau_squared = ((a**2 + d**2 - f_squared) \
+                        cos_tau_squared = ((a**2 + d**2 - f_squared)
                                            / (2 * a * d))**2
                         # Consider b2 > 0 and avoid division by zero
                         k = a**2 - d**2 * cos_tau_squared
@@ -205,21 +206,29 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
                     hist, bin_edges = np.histogram(acc, bins=bins)
                     hist_max = np.max(hist)
                     if hist_max > threshold:
-                        angle = np.arctan2(p1x - p2x, p1y - p2y)
-                        # pi - angle to keep ellipse_perimeter() convention
-                        if angle != 0:
-                            angle = np.pi - angle
+                        orientation = atan2(p1x - p2x, p1y - p2y)
                         b = sqrt(bin_edges[hist.argmax()])
-                        results.append((x0,
-                                        y0,
-                                        a,
-                                        b,
-                                        angle,
-                                        hist_max,  # Accumulator
-                                        ))
+                        # to keep ellipse_perimeter() convention
+                        if orientation != 0:
+                            orientation = M_PI - orientation
+                            # When orientation is not in [-pi:pi]
+                            # it would mean in ellipse_perimeter()
+                            # that a < b. But we keep a > b.
+                            if orientation > M_PI:
+                                orientation = orientation - M_PI / 2.
+                                a, b = b, a
+                        results.append((hist_max, # Accumulator
+                                        yc, xc,
+                                        a, b,
+                                        orientation))
                     acc = []
 
-    return results
+    return np.array(results, dtype=[('accumulator', np.intp),
+                                    ('yc', np.double),
+                                    ('xc', np.double),
+                                    ('a', np.double),
+                                    ('b', np.double),
+                                    ('orientation', np.double)])
 
 
 def hough_line(cnp.ndarray img,
diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py
index 651861f3..fb19d8c1 100644
--- a/skimage/transform/tests/test_hough_transform.py
+++ b/skimage/transform/tests/test_hough_transform.py
@@ -1,7 +1,5 @@
 import numpy as np
-from numpy.testing import (assert_almost_equal,
-                           assert_equal,
-                           )
+from numpy.testing import assert_almost_equal, assert_equal
 
 import skimage.transform as tf
 from skimage.draw import line, circle_perimeter, ellipse_perimeter
@@ -81,8 +79,10 @@ def test_hough_line_peaks_dist():
     img[:, 30] = True
     img[:, 40] = True
     hspace, angles, dists = tf.hough_line(img)
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_distance=5)[0]) == 2
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_distance=15)[0]) == 1
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_distance=5)[0]) == 2
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_distance=15)[0]) == 1
 
 
 def test_hough_line_peaks_angle():
@@ -91,18 +91,24 @@ def test_hough_line_peaks_angle():
     img[0, :] = True
 
     hspace, angles, dists = tf.hough_line(img)
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_angle=45)[0]) == 2
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_angle=90)[0]) == 1
 
     theta = np.linspace(0, np.pi, 100)
     hspace, angles, dists = tf.hough_line(img, theta)
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_angle=45)[0]) == 2
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_angle=90)[0]) == 1
 
     theta = np.linspace(np.pi / 3, 4. / 3 * np.pi, 100)
     hspace, angles, dists = tf.hough_line(img, theta)
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2
-    assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_angle=45)[0]) == 2
+    assert len(tf.hough_line_peaks(hspace, angles, dists,
+                                   min_angle=90)[0]) == 1
 
 
 def test_hough_line_peaks_num():
@@ -149,36 +155,204 @@ def test_hough_circle_extended():
 
 def test_hough_ellipse_zero_angle():
     img = np.zeros((25, 25), dtype=int)
-    a = 6
-    b = 8
+    rx = 6
+    ry = 8
     x0 = 12
-    y0 = 12
+    y0 = 15
     angle = 0
-    rr, cc = ellipse_perimeter(x0, x0, b, a)
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=9)
-    assert_equal(result[0][0], x0)
-    assert_equal(result[0][1], y0)
-    assert_almost_equal(result[0][2], b, decimal=1)
-    assert_almost_equal(result[0][3], a, decimal=1)
-    assert_equal(result[0][4], angle)
+    best = result[-1]
+    assert_equal(best[1], y0)
+    assert_equal(best[2], x0)
+    assert_almost_equal(best[3], ry, decimal=1)
+    assert_almost_equal(best[4], rx, decimal=1)
+    assert_equal(best[5], angle)
+    # 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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
 
 
-def test_hough_ellipse_non_zero_angle():
-    img = np.zeros((20, 20), dtype=int)
-    a = 6
-    b = 9
+def test_hough_ellipse_non_zero_posangle1():
+    # ry > rx, angle in [0:pi/2]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 6
+    ry = 12
     x0 = 10
-    y0 = 10
+    y0 = 15
     angle = np.pi / 1.35
-    rr, cc = ellipse_perimeter(x0, x0, b, a, orientation=angle)
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
     img[rr, cc] = 1
     result = tf.hough_ellipse(img, threshold=15, accuracy=3)
-    assert_almost_equal(result[0][0] / 100., x0 / 100., decimal=1)
-    assert_almost_equal(result[0][1] / 100., y0 / 100., decimal=1)
-    assert_almost_equal(result[0][2] / 100., b / 100., decimal=1)
-    assert_almost_equal(result[0][3] / 100., a / 100., decimal=1)
-    assert_almost_equal(result[0][4], angle, decimal=1)
+    result.sort(order='accumulator')
+    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)
+    assert_almost_equal(best[4] / 100., rx / 100., decimal=1)
+    assert_almost_equal(best[5], angle, decimal=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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_posangle2():
+    # ry < rx, angle in [0:pi/2]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 12
+    ry = 6
+    x0 = 10
+    y0 = 15
+    angle = np.pi / 1.35
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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)
+    assert_almost_equal(best[4] / 100., rx / 100., decimal=1)
+    assert_almost_equal(best[5], angle, decimal=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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_posangle3():
+    # ry < rx, angle in [pi/2:pi]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 12
+    ry = 6
+    x0 = 10
+    y0 = 15
+    angle = np.pi / 1.35 + np.pi / 2.
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_posangle4():
+    # ry < rx, angle in [pi:3pi/4]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 12
+    ry = 6
+    x0 = 10
+    y0 = 15
+    angle = np.pi / 1.35 + np.pi
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_negangle1():
+    # ry > rx, angle in [0:-pi/2]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 6
+    ry = 12
+    x0 = 10
+    y0 = 15
+    angle = - np.pi / 1.35
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_negangle2():
+    # ry < rx, angle in [0:-pi/2]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 12
+    ry = 6
+    x0 = 10
+    y0 = 15
+    angle = - np.pi / 1.35
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_negangle3():
+    # ry < rx, angle in [-pi/2:-pi]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 12
+    ry = 6
+    x0 = 10
+    y0 = 15
+    angle = - np.pi / 1.35 - np.pi / 2.
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
+
+
+def test_hough_ellipse_non_zero_negangle4():
+    # ry < rx, angle in [-pi:-3pi/4]
+    img = np.zeros((30, 24), dtype=int)
+    rx = 12
+    ry = 6
+    x0 = 10
+    y0 = 15
+    angle = - np.pi / 1.35 - np.pi
+    rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
+    img[rr, cc] = 1
+    result = tf.hough_ellipse(img, threshold=15, accuracy=3)
+    result.sort(order='accumulator')
+    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])
+    assert_equal(rr, rr2)
+    assert_equal(cc, cc2)
 
 
 if __name__ == "__main__":