+#!/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()
projects / experiments / line_equation.git / commitdiff