143 expand linexy constructors and slope information

opencsp/common/lib/geometry/LineXY.py

@@ -2,6 +2,7 @@
 from numpy.random import RandomState, SeedSequence, MT19937
 from scipy.optimize import minimize
 
+from opencsp.common.lib.geometry.angle import normalize as normalize_angle
 from opencsp.common.lib.geometry.Vxy import Vxy
 
 
@@ -36,6 +37,10 @@ def __init__(self, A: float, B: float, C: float):
  self.B = B / mag
  self.C = C / mag
 
+ # The original two points used to create this line, if created with
+ # from_two_points
+ self._original_two_points: tuple[Vxy, Vxy] | None = None
+
  def __repr__(self):
  return '2D Line: ' + self.A.__repr__() + ', ' + self.B.__repr__() + ', ' + self.C.__repr__()
 
@@ -67,11 +72,68 @@ def ABC(self) -> np.ndarray:
 
  @property
  def slope(self) -> float:
- """Get the slope of the line (could be infinity!)"""
+ """Get the slope of the line, as rise over run (could be infinity!)"""
+ # check for infinity
  if abs(self.B) < 1e-10:
- return np.inf
+ if self._original_two_points is None:
+ # can't tell the difference between pos/neg infinity
+ return np.inf
+
+ else:
+ # use the original two points to determine the positivity of the slope
+ if self._original_two_points[1].y[0] > self._original_two_points[0].y[0]:
+ return np.PINF
+ else:
+ return np.NINF
+
+ # return the slope
  return -self.A / self.B
 
+ @property
+ def angle(self) -> float:
+ """
+ Get the angle of this vector in radians, as measured by its slope.
+
+ As in all of OpenCSP, the angle coordinate system is defined as:
+ - 0 along the positive x axis (pointing to the right)
+ - increasing counter-clockwise
+ """
+ slope = self.slope
+
+ # get the angle between -pi/2 and pi/2
+ atan = np.arctan(slope)
+
+ # is the slope positive or negative infinity?
+ if np.isinf(slope) or np.isneginf(slope):
+ return normalize_angle(atan)
+
+ # is the slope zero?
+ if slope == 0:
+ if self._original_two_points is None:
+ # can't tell if right or left pointing
+ return 0
+
+ else:
+ # use the original two points to determine if right or left pointing
+ if self._original_two_points[1].x[0] > self._original_two_points[0].x[0]:
+ return 0
+ else:
+ return np.pi
+
+ # correct for 2nd and 3rd quadrants
+ if self._original_two_points is not None:
+ vector = self._original_two_points[1] - self._original_two_points[0]
+ if atan < 0:
+ if vector.y[0] > 0:
+ # 2nd quadrant
+ return atan + np.pi
+ elif atan > 0:
+ if vector.y[0] < 0:
+ # 3rd quadrant
+ return (np.pi * 3 / 2) - atan
+
+ return normalize_angle(atan)
+
  @classmethod
  def fit_from_points(cls, Pxy: Vxy, seed: int = 1, neighbor_dist: float = 1.0):
  """
@@ -188,7 +250,70 @@ def from_two_points(cls, Pxy1: Vxy, Pxy2: Vxy):
 
  ABC = np.concatenate((AB, [C]), axis=0)
 
- return LineXY(*ABC)
+ ret = LineXY(*ABC)
+ ret._original_two_points = Pxy1, Pxy2
+ return ret
+
+ @classmethod
+ def from_rho_theta(cls, rho: float, theta: float) -> "LineXY":
+ """
+ Get a new instance of this class built from the rho + theta
+ representation of a line. This is particularly useful for representation
+ of lines found via the Hough Transform of an image
+ (https://docs.opencv.org/3.4/d9/db0/tutorial_hough_lines.html).
+
+ Parameters
+ ----------
+ rho : float
+ The right angle distance between the line and the origin (0,0). For
+ images, this will be the top-left corner of the image.
+ theta : float
+ The angle between the right angle distance vector and the X axis.
+ Units are radians on the standard graphing coordinate system (0 is
+ on the positive x-axis to the right, and the angle increases
+ counter-clockwise).
+ """
+ a = np.cos(theta)
+ b = np.sin(theta)
+ x0 = a * rho
+ y0 = b * rho
+ pt1 = (x0 + 1000 * (-b), y0 + 1000 * (a))
+ pt2 = (x0 - 1000 * (-b), y0 - 1000 * (a))
+ ret = cls.from_two_points(Vxy(pt1), Vxy(pt2))
+
+ return ret
+
+ @classmethod
+ def from_location_angle(cls, location: Vxy, angle: float) -> "LineXY":
+ """
+ Get a new instance of this class built from the xy location + angle representation of a line.
+
+ Parameters
+ ----------
+ location : Vxy
+ A point that the line travels through on the cartesion x/y grid.
+ angle : float
+ The angle of the line, in radians, for the standard graphing
+ coordinate system. 0 is on the positive x-axis to the right, and the
+ angle increases counter-clockwise.
+ """
+ pt1 = location
+
+ # normalize input
+ angle = normalize_angle(angle)
+
+ # # values to help with small angle issues
+ # onepi_angle: float = angle if angle < np.pi else angle-np.pi
+ # small_angle: float = np.deg2rad(3)
+
+ # sohcahtoa
+ hypotenuse = 1000.0
+ x = hypotenuse * np.cos(angle)
+ y = hypotenuse * np.sin(angle)
+ pt2 = Vxy(np.array([[x + pt1.x[0]], [y + pt1.y[0]]]))
+
+ # build the line
+ return cls.from_two_points(pt1, pt2)
 
  def y_from_x(self, xs: np.ndarray | float) -> np.ndarray | float:
  """
@@ -287,7 +412,7 @@ def intersect_with(self, Lxy: 'LineXY') -> Vxy | None:
  v = np.cross(self.ABC, Lxy.ABC)
  return Vxy(v[:2] / v[2])
 
- def flip(self):
+ def flip(self) -> "LineXY":
  """
  Flips the orientation of a line. Returns a flipped copy of the LineXY.
 
@@ -298,6 +423,7 @@ def flip(self):
 
  """
  L_out = LineXY(*self.ABC)
+ L_out._original_two_points = self._original_two_points
  L_out.flip_in_place()
  return L_out
 
@@ -310,3 +436,6 @@ def flip_in_place(self) -> None:
  self.A *= -1
  self.B *= -1
  self.C *= -1
+
+ if self._original_two_points is not None:
+ self._original_two_points = self._original_two_points[1], self._original_two_points[0]

opencsp/common/lib/geometry/test/test_LineXY.py

@@ -55,6 +55,38 @@ def test_from_two_points_edge_case(self):
  self.assertAlmostEqual(-linev.C / linev.A, 1)
  self.assertAlmostEqual(linev.B, 0)
 
+ def test_from_rho_theta(self):
+ # vertical line
+ l1 = LineXY.from_rho_theta(1, 0)
+ self.assertAlmostEqual(l1.x_from_y(-1), 1)
+ self.assertAlmostEqual(l1.x_from_y(1), 1)
+
+ # 45-degree downward slope
+ l2 = LineXY.from_rho_theta(np.sqrt(2) / 2, np.pi / 4)
+ self.assertAlmostEqual(l2.y_from_x(0), 1)
+ self.assertAlmostEqual(l2.y_from_x(1), 0)
+
+ # horizontal line
+ l3 = LineXY.from_rho_theta(1, np.pi / 2)
+ self.assertAlmostEqual(l3.y_from_x(-1), 1)
+ self.assertAlmostEqual(l3.y_from_x(1), 1)
+
+ def test_from_location_angle(self):
+ # horizontal line
+ l1 = LineXY.from_location_angle(Vxy([1, 1]), 0)
+ self.assertAlmostEqual(l1.y_from_x(-1), 1)
+ self.assertAlmostEqual(l1.y_from_x(2), 1)
+
+ # 45-degree upward slope
+ l2 = LineXY.from_location_angle(Vxy([1, 1]), np.pi / 4)
+ self.assertAlmostEqual(l2.y_from_x(0), 0)
+ self.assertAlmostEqual(l2.y_from_x(2), 2)
+
+ # vertical line
+ l3 = LineXY.from_location_angle(Vxy([1, 1]), np.pi / 2)
+ self.assertAlmostEqual(l3.x_from_y(-1), 1)
+ self.assertAlmostEqual(l3.x_from_y(2), 1)
+
  def test_y_from_x(self):
  # Line y = -x
  line = LineXY(1, 1, 0)
@@ -173,6 +205,72 @@ def flip_in_place(self):
  line.flip()
  np.testing.assert_almost_equal(line.ABC, -ABC)
 
+ def test_slope(self):
+ # flat line to the right
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((1, 0)))
+ self.assertAlmostEqual(line.slope, 0)
+
+ # flat line to the left
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((-1, 0)))
+ self.assertAlmostEqual(line.slope, 0)
+
+ # vertical line up
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((0, 1)))
+ self.assertTrue(np.isinf(line.slope) and not np.isneginf(line.slope))
+
+ # vertical line down
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((0, -1)))
+ self.assertTrue(np.isneginf(line.slope))
+
+ # 45-degree, up and to the right
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((1, 1)))
+ self.assertAlmostEqual(line.slope, 1)
+
+ # 135-degree, up and to the left
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((-1, 1)))
+ self.assertAlmostEqual(line.slope, -1)
+
+ # 225-degree, down and to the left
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((-1, -1)))
+ self.assertAlmostEqual(line.slope, 1)
+
+ # 315-degree, down and to the right
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((1, -1)))
+ self.assertAlmostEqual(line.slope, -1)
+
+ def test_angle(self):
+ # flat line to the right
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((1, 0)))
+ self.assertAlmostEqual(line.angle, 0)
+
+ # flat line to the left
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((-1, 0)))
+ self.assertAlmostEqual(line.angle, np.pi)
+
+ # vertical line up
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((0, 1)))
+ self.assertAlmostEqual(line.angle, np.pi / 2)
+
+ # vertical line down
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((0, -1)))
+ self.assertAlmostEqual(line.angle, np.pi * 3 / 2)
+
+ # 45-degree, up and to the right
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((1, 1)))
+ self.assertAlmostEqual(line.angle, np.pi / 4)
+
+ # 135-degree, up and to the left
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((-1, 1)))
+ self.assertAlmostEqual(line.angle, np.pi * 3 / 4)
+
+ # 225-degree, down and to the left
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((-1, -1)))
+ self.assertAlmostEqual(line.angle, np.pi * 5 / 4)
+
+ # 315-degree, down and to the right
+ line = LineXY.from_two_points(Vxy((0, 0)), Vxy((1, -1)))
+ self.assertAlmostEqual(line.angle, np.pi * 7 / 4)
+
 
 if __name__ == '__main__':
  unittest.main()