line_equation.py

   1 #!/usr/bin/env python3
   2 #
   3 # line_equation - some experiments about representing a line
   4 #
   5 # Copyright (C) 2016  Antonio Ospite <ao2@ao2.it>
   6 #
   7 # This program is free software: you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation, either version 3 of the License, or
  10 # (at your option) any later version.
  11 #
  12 # This program is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  19
  20 """Some experiments about representing a line passing for two points"""
  21
  22 import math
  23 import matplotlib.pyplot as plt
  24 import numpy as np
  25
  26 EPSILON = 10e-6
  27 INFINITY = 10e6
  28
  29 # pylint: disable=invalid-name
  30
  31
  32 def distance(p1, p2):
  33     """Return the distance between the two points p1 and p2"""
  34     dx = p2[0] - p1[0]
  35     dy = p2[1] - p1[1]
  36
  37     return math.sqrt(dx * dx + dy * dy)
  38
  39
  40 def line_eq_slope_intercept(p1, p2):
  41     """Return the coefficients of the line equation in the slope/intercept
  42     form:
  43
  44         y = mx + q
  45     """
  46     dx = p2[0] - p1[0]
  47     dy = p2[1] - p1[1]
  48
  49     # For vertical lines the slope would be undefined in theory.
  50     # For practical purposes a large value can work.
  51     if dx == 0:
  52         slope = INFINITY * math.copysign(1, dy)
  53     else:
  54         slope = dy / dx
  55
  56     intercept = -slope * p1[0] + p1[1]
  57
  58     return slope, intercept
  59
  60
  61 def line_angle(p1, p2):
  62     """Return the angle between the line passing for p1 and p2, and the positive
  63     x axis.
  64
  65     The return value is in radians.
  66     """
  67     dx = p2[0] - p1[0]
  68     dy = p2[1] - p1[1]
  69
  70     angle = math.atan2(dy, dx)
  71     if angle < 0:
  72         angle += 2 * math.pi
  73
  74     return angle
  75
  76
  77 def slope_intercept_to_a_b_c(slope, intercept):
  78     """Convert the coefficients from the slope/intercept form:
  79
  80         y = mx + q
  81
  82     to the coefficients in the linear equation form:
  83
  84         ax + by + c = 0
  85     """
  86     a = -slope
  87     b = 1.
  88     c = -intercept
  89
  90     return a, b, c
  91
  92
  93 def line_eq_two_points(p1, p2):
  94     """Return the coefficients for the line equation, in the linear equation
  95     form:
  96
  97         ax + by + c = 0
  98     """
  99     a = p2[1] - p1[1]
 100     b = - (p2[0] - p1[0])
 101     c = - (p1[0] * p2[1] - p2[0] * p1[1])
 102
 103     return a, b, c
 104
 105
 106 def line_calc_normal(p1, p2):
 107     """Return the unit normal to the line passing for two points."""
 108     dx = p2[0] - p1[0]
 109     dy = p2[1] - p1[1]
 110
 111     n_length = math.sqrt(dx*dx + dy*dy)
 112     if n_length < EPSILON:
 113         nx = 0
 114         ny = 0
 115     else:
 116         nx = dy / n_length
 117         ny = -dx / n_length
 118
 119     return nx, ny
 120
 121
 122 def list_rotate(l, n):
 123     return l[n:] + l[:n]
 124
 125
 126 def main():
 127     points = [
 128         [132.143, 941.429],
 129         [480.714, 905.714],
 130         [480.714, 735.714],
 131         [180.0, 783.571]
 132     ]
 133
 134     # The points are assumed to be in clockwise order
 135     for p1, p2 in zip(points, list_rotate(points, 1)):
 136         print("p1", p1)
 137         print("p2", p2)
 138
 139         line_x_extremes = [p1[0], p2[0]]
 140         line_y_extremes = [p1[1], p2[1]]
 141
 142         plt.plot(line_x_extremes, line_y_extremes, 'k')
 143
 144         x_samples = np.linspace(min(line_x_extremes), max(line_x_extremes), 15)
 145
 146         # slope/intercept form
 147         m, q = line_eq_slope_intercept(p1, p2)
 148
 149         if abs(m) == INFINITY:
 150             print("m is undefined: vertical line!")
 151             y_samples = \
 152                 np.linspace(min(line_y_extremes), max(line_y_extremes), 15)
 153         else:
 154             y_samples = x_samples * m + q
 155
 156         plt.plot(x_samples, y_samples, 'wo', markersize=15)
 157
 158         # linear equation form
 159         a, b, c = line_eq_two_points(p1, p2)
 160
 161         if b == 0:
 162             print("b is 0: vertical line!")
 163             y_samples = \
 164                 np.linspace(min(line_y_extremes), max(line_y_extremes), 15)
 165         else:
 166             y_samples = (-x_samples * a - c) / b
 167
 168         plt.plot(x_samples, y_samples, 'ws', markersize=10)
 169
 170         # parametric forms
 171         t = np.linspace(0, 1, len(x_samples))
 172
 173         # point and direction form
 174         dx = p2[0] - p1[0]
 175         dy = p2[1] - p1[1]
 176
 177         x_samples = p1[0] + dx * t
 178         y_samples = p1[1] + dy * t
 179         plt.plot(x_samples, y_samples, 'm+', markersize=10)
 180
 181         # or using a point and the normal, but in this case the normal has to
 182         # be rotated by 90 degrees counterclockwise
 183         nx, ny = line_calc_normal(p1, p2)
 184         rnx = -ny
 185         rny = nx
 186
 187         length = distance(p1, p2)
 188         x_samples = p1[0] + rnx * t * length
 189         y_samples = p1[1] + rny * t * length
 190         plt.plot(x_samples, y_samples, 'bx', markersize=10)
 191
 192         # draw the normal
 193         mid_x = (p1[0] + p2[0]) / 2
 194         mid_y = (p1[1] + p2[1]) / 2
 195         n_length = 20
 196         plt.arrow(mid_x, mid_y, nx * n_length, ny * n_length,
 197                   head_width=10, head_length=5, fc='k', ec='k')
 198
 199         # check the other functions
 200         print("angle: %g" % math.degrees(line_angle(p1, p2)))
 201         a1, b1, c1 = slope_intercept_to_a_b_c(m, q)
 202         print(a1, b1, c1)
 203         if b != 0:
 204             print(a / b, b / b, c / b)
 205         else:
 206             print(a, b, c)
 207         print()
 208
 209     # draw the vertices on top of everything else
 210     xvalues = [p[0] for p in points]
 211     yvalues = [p[1] for p in points]
 212     plt.plot(xvalues, yvalues, 'ro')
 213     plt.title('Line passing for two points')
 214     plt.xlabel('x')
 215     plt.ylabel('y')
 216
 217     plt.plot([], 'k', label="line")
 218     plt.plot([], 'wo', label="y = mx + q")
 219     plt.plot([], 'ws', label="ax + by + c = 0")
 220     plt.plot([], 'm+', label="P = P1 + D * t")
 221     plt.plot([], 'bx', label="P = P1 + N * t * length")
 222     plt.legend(loc='lower left', numpoints=1, fancybox=True, prop={'size': 10})
 223     plt.show()
 224
 225
 226 if __name__ == "__main__":
 227     main()