Initial import

diff --git a/line_equation.py b/line_equation.py
new file mode 100755 (executable)
index 0000000..0bce46c
--- /dev/null
+++ b/line_equation.py
@@ -0,0 +1,227 @@
+#!/usr/bin/env python3
+#
+# line_equation - some experiments about representing a line
+#
+# Copyright (C) 2016  Antonio Ospite <ao2@ao2.it>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Some experiments about representing a line passing for two points"""
+
+import math
+import matplotlib.pyplot as plt
+import numpy as np
+
+EPSILON = 10e-6
+INFINITY = 10e6
+
+# pylint: disable=invalid-name
+
+
+def distance(p1, p2):
+    """Return the distance between the two points p1 and p2"""
+    dx = p2[0] - p1[0]
+    dy = p2[1] - p1[1]
+
+    return math.sqrt(dx * dx + dy * dy)
+
+
+def line_eq_slope_intercept(p1, p2):
+    """Return the coefficients of the line equation in the slope/intercept
+    form:
+
+        y = mx + q
+    """
+    dx = p2[0] - p1[0]
+    dy = p2[1] - p1[1]
+
+    # For vertical lines the slope would be undefined in theory.
+    # For practical purposes a large value can work.
+    if dx == 0:
+        slope = INFINITY * math.copysign(1, dy)
+    else:
+        slope = dy / dx
+
+    intercept = -slope * p1[0] + p1[1]
+
+    return slope, intercept
+
+
+def line_angle(p1, p2):
+    """Return the angle between the line passing for p1 and p2, and the positive
+    x axis.
+
+    The return value is in radians.
+    """
+    dx = p2[0] - p1[0]
+    dy = p2[1] - p1[1]
+
+    angle = math.atan2(dy, dx)
+    if angle < 0:
+        angle += 2 * math.pi
+
+    return angle
+
+
+def slope_intercept_to_a_b_c(slope, intercept):
+    """Convert the coefficients from the slope/intercept form:
+
+        y = mx + q
+
+    to the coefficients in the linear equation form:
+
+        ax + by + c = 0
+    """
+    a = -slope
+    b = 1.
+    c = -intercept
+
+    return a, b, c
+
+
+def line_eq_two_points(p1, p2):
+    """Return the coefficients for the line equation, in the linear equation
+    form:
+
+        ax + by + c = 0
+    """
+    a = p2[1] - p1[1]
+    b = - (p2[0] - p1[0])
+    c = - (p1[0] * p2[1] - p2[0] * p1[1])
+
+    return a, b, c
+
+
+def line_calc_normal(p1, p2):
+    """Return the unit normal to the line passing for two points."""
+    dx = p2[0] - p1[0]
+    dy = p2[1] - p1[1]
+
+    n_length = math.sqrt(dx*dx + dy*dy)
+    if n_length < EPSILON:
+        nx = 0
+        ny = 0
+    else:
+        nx = dy / n_length
+        ny = -dx / n_length
+
+    return nx, ny
+
+
+def list_rotate(l, n):
+    return l[n:] + l[:n]
+
+
+def main():
+    points = [
+        [132.143, 941.429],
+        [480.714, 905.714],
+        [480.714, 735.714],
+        [180.0, 783.571]
+    ]
+
+    # The points are assumed to be in clockwise order
+    for p1, p2 in zip(points, list_rotate(points, 1)):
+        print("p1", p1)
+        print("p2", p2)
+
+        line_x_extremes = [p1[0], p2[0]]
+        line_y_extremes = [p1[1], p2[1]]
+
+        plt.plot(line_x_extremes, line_y_extremes, 'k')
+
+        x_samples = np.linspace(min(line_x_extremes), max(line_x_extremes), 15)
+
+        # slope/intercept form
+        m, q = line_eq_slope_intercept(p1, p2)
+
+        if abs(m) == INFINITY:
+            print("m is undefined: vertical line!")
+            y_samples = \
+                np.linspace(min(line_y_extremes), max(line_y_extremes), 15)
+        else:
+            y_samples = x_samples * m + q
+
+        plt.plot(x_samples, y_samples, 'wo', markersize=15)
+
+        # linear equation form
+        a, b, c = line_eq_two_points(p1, p2)
+
+        if b == 0:
+            print("b is 0: vertical line!")
+            y_samples = \
+                np.linspace(min(line_y_extremes), max(line_y_extremes), 15)
+        else:
+            y_samples = (-x_samples * a - c) / b
+
+        plt.plot(x_samples, y_samples, 'ws', markersize=10)
+
+        # parametric forms
+        t = np.linspace(0, 1, len(x_samples))
+
+        # point and direction form
+        dx = p2[0] - p1[0]
+        dy = p2[1] - p1[1]
+
+        x_samples = p1[0] + dx * t
+        y_samples = p1[1] + dy * t
+        plt.plot(x_samples, y_samples, 'm+', markersize=10)
+
+        # or using a point and the normal, but in this case the normal has to
+        # be rotated by 90 degrees counterclockwise
+        nx, ny = line_calc_normal(p1, p2)
+        rnx = -ny
+        rny = nx
+
+        length = distance(p1, p2)
+        x_samples = p1[0] + rnx * t * length
+        y_samples = p1[1] + rny * t * length
+        plt.plot(x_samples, y_samples, 'bx', markersize=10)
+
+        # draw the normal
+        mid_x = (p1[0] + p2[0]) / 2
+        mid_y = (p1[1] + p2[1]) / 2
+        n_length = 20
+        plt.arrow(mid_x, mid_y, nx * n_length, ny * n_length,
+                  head_width=10, head_length=5, fc='k', ec='k')
+
+        # check the other functions
+        print("angle: %g" % math.degrees(line_angle(p1, p2)))
+        a1, b1, c1 = slope_intercept_to_a_b_c(m, q)
+        print(a1, b1, c1)
+        if b != 0:
+            print(a / b, b / b, c / b)
+        else:
+            print(a, b, c)
+        print()
+
+    # draw the vertices on top of everything else
+    xvalues = [p[0] for p in points]
+    yvalues = [p[1] for p in points]
+    plt.plot(xvalues, yvalues, 'ro')
+    plt.title('Line passing for two points')
+    plt.xlabel('x')
+    plt.ylabel('y')
+
+    plt.plot([], 'k', label="line")
+    plt.plot([], 'wo', label="y = mx + q")
+    plt.plot([], 'ws', label="ax + by + c = 0")
+    plt.plot([], 'm+', label="P = P1 + D * t")
+    plt.plot([], 'bx', label="P = P1 + N * t * length")
+    plt.legend(loc='lower left', numpoints=1, fancybox=True, prop={'size': 10})
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()