From 90335efdda23570fdb2a3c6f54b83c526fa5110f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= <fboulogne sciunto org>
Date: Tue, 27 Aug 2013 16:36:19 +0200
Subject: [PATCH] Implement line_aa

---
 skimage/draw/__init__.py        |  6 ++-
 skimage/draw/_draw.pyx          | 81 +++++++++++++++++++++++++++++++++
 skimage/draw/tests/test_draw.py | 45 ++++++++++++++++--
 3 files changed, 127 insertions(+), 5 deletions(-)

diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index 285ba49a..5ca5afb1 100644
--- a/skimage/draw/__init__.py
+++ b/skimage/draw/__init__.py
@@ -1,9 +1,11 @@
 from .draw import circle, ellipse, set_color
-from ._draw import (line, polygon, ellipse_perimeter, circle_perimeter,
-                    circle_perimeter_aa, bezier_segment)
 from .draw3d import ellipsoid, ellipsoid_stats
+from ._draw import (line, line_aa, polygon, ellipse_perimeter,
+                    circle_perimeter, circle_perimeter_aa,
+                    bezier_segment)
 
 __all__ = ['line',
+           'line_aa',
            'polygon',
            'ellipse',
            'ellipse_perimeter',
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index a8716983..3e835d93 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -27,6 +27,10 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2):
         May be used to directly index into an array, e.g.
         ``img[rr, cc] = 1``.
 
+    Notes
+    -----
+    Anti-aliased line generator is available with `line_aa`.
+
     Examples
     --------
     >>> from skimage.draw import line
@@ -89,6 +93,83 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2):
     return np.asarray(rr), np.asarray(cc)
 
 
+def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2):
+    """Generate line pixel coordinates.
+
+    Parameters
+    ----------
+    y1, x1 : int
+        Starting position (row, column).
+    y2, x2 : int
+        End position (row, column).
+
+    Returns
+    -------
+    rr, cc, val : (N,) ndarray of int
+        Indices of pixels that belong to the line.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = val``.
+
+    References
+    ----------
+    .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012
+           http://members.chello.at/easyfilter/Bresenham.pdf
+
+    """
+    cdef list rr = list()
+    cdef list cc = list()
+    cdef list val = list()
+
+    cdef int dx = abs(x1 - x2)
+    cdef int dy = abs(y1 - y2)
+    cdef int err = dx - dy
+    cdef int x, y, e, ed, sign_x, sign_y
+
+    if x1 < x2:
+        sign_x = 1
+    else:
+        sign_x = -1
+
+    if y1 < y2:
+        sign_y = 1
+    else:
+        sign_y = -1
+
+    if dx + dy == 0:
+        ed = 1
+    else:
+        ed = <int>(sqrt(dx*dx + dy*dy))
+
+    x, y = x1, y1
+    while True:
+        cc.append(x)
+        rr.append(y)
+        val.append(255 * abs(err - dx + dy) / ed)
+        e = err
+        if 2 * e >= -dx:
+            if x == x2:
+                break
+            if e + dy < ed:
+                cc.append(x)
+                rr.append(y + sign_y)
+                val.append(255 * abs(e + dy) / ed)
+            err -= dy
+            x += sign_x
+        if 2 * e <= dy:
+            if y == y2:
+                break
+            if dx - e < ed:
+                cc.append(x)
+                rr.append(y)
+                val.append(255 * abs(dx - e) / ed)
+            err += dx
+            y += sign_y
+
+    return (np.array(rr, dtype=np.intp),
+            np.array(cc, dtype=np.intp),
+            255 - np.array(val, dtype=np.intp))
+
+
 def polygon(y, x, shape=None):
     """Generate coordinates of pixels within polygon.
 
diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py
index 152797a7..23bbff4e 100644
--- a/skimage/draw/tests/test_draw.py
+++ b/skimage/draw/tests/test_draw.py
@@ -1,8 +1,10 @@
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_equal, assert_equal
 import numpy as np
 
-from skimage.draw import (line, polygon, circle, circle_perimeter,
-                          circle_perimeter_aa, ellipse,
+from skimage.draw import (line, line_aa,
+                          polygon, circle,
+                          circle_perimeter, circle_perimeter_aa,
+                          ellipse,
                           ellipse_perimeter, bezier_segment,
                           )
 
@@ -54,6 +56,43 @@ def test_line_diag():
     assert_array_equal(img, img_)
 
 
+def test_line_aa_horizontal():
+    img = np.zeros((10, 10))
+
+    rr, cc, val = line_aa(0, 0, 0, 9)
+    img[rr, cc] = val
+
+    img_ = np.zeros((10, 10))
+    img_[0, :] = 255
+
+    assert_array_equal(img, img_)
+
+
+def test_line_aa_vertical():
+    img = np.zeros((10, 10))
+
+    rr, cc, val = line_aa(0, 0, 9, 0)
+    img[rr, cc] = val
+
+    img_ = np.zeros((10, 10))
+    img_[:, 0] = 255
+
+    assert_array_equal(img, img_)
+
+
+def test_line_aa_diagonal():
+    img = np.zeros((10, 10))
+
+    rr, cc, val = line_aa(0, 0, 9, 6)
+    img[rr, cc] = 1
+
+    # Check that each pixel belonging to line,
+    # also belongs to line_aa
+    r, c = line(0, 0, 9, 6)
+    for x, y in zip(r, c):
+        assert_equal(img[r, c], 1)
+
+
 def test_polygon_rectangle():
     img = np.zeros((10, 10), 'uint8')
     poly = np.array(((1, 1), (4, 1), (4, 4), (1, 4), (1, 1)))