From: Antonio Ospite Date: Sat, 6 Aug 2016 17:27:25 +0000 (+0200) Subject: Initial import X-Git-Url: https://git.ao2.it/experiments/line_equation.git/commitdiff_plain/HEAD Initial import --- 9a97c7bceb623e978c2f34cfbbbbc7c5ea21a5ee diff --git a/line_equation.py b/line_equation.py new file mode 100755 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 +# +# 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 . + +"""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()