From 745490884e08e5158501b19b9575cb4190cdb2d6 Mon Sep 17 00:00:00 2001 From: bbean Date: Fri, 2 Aug 2024 11:56:05 -0600 Subject: [PATCH 1/4] add "_original_two_points" property in order to be able to tell the difference between positive and negative infinity in the slope --- opencsp/common/lib/geometry/LineXY.py | 30 ++++++++++++++--- .../common/lib/geometry/test/test_LineXY.py | 32 +++++++++++++++++++ 2 files changed, 58 insertions(+), 4 deletions(-) diff --git a/opencsp/common/lib/geometry/LineXY.py b/opencsp/common/lib/geometry/LineXY.py index 90d44bcf..1462c366 100644 --- a/opencsp/common/lib/geometry/LineXY.py +++ b/opencsp/common/lib/geometry/LineXY.py @@ -36,6 +36,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,9 +71,21 @@ 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 @classmethod @@ -188,7 +204,9 @@ 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 def y_from_x(self, xs: np.ndarray | float) -> np.ndarray | float: """ @@ -287,7 +305,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 +316,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 +329,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] diff --git a/opencsp/common/lib/geometry/test/test_LineXY.py b/opencsp/common/lib/geometry/test/test_LineXY.py index c8e4b4ea..005602a4 100644 --- a/opencsp/common/lib/geometry/test/test_LineXY.py +++ b/opencsp/common/lib/geometry/test/test_LineXY.py @@ -173,6 +173,38 @@ 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) if __name__ == '__main__': unittest.main() From 04713dccdd04b9bba2dc4de40cf6294bca1de5ed Mon Sep 17 00:00:00 2001 From: bbean Date: Mon, 5 Aug 2024 09:38:09 -0600 Subject: [PATCH 2/4] add angle()/from_rho_theta()/from_location_angle() methods to LineXY --- opencsp/common/lib/geometry/LineXY.py | 102 ++++++++++++++++++ .../common/lib/geometry/test/test_LineXY.py | 34 ++++++ 2 files changed, 136 insertions(+) diff --git a/opencsp/common/lib/geometry/LineXY.py b/opencsp/common/lib/geometry/LineXY.py index 1462c366..f31496f2 100644 --- a/opencsp/common/lib/geometry/LineXY.py +++ b/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 @@ -88,6 +89,51 @@ def slope(self) -> float: # 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): """ @@ -208,6 +254,62 @@ def from_two_points(cls, Pxy1: Vxy, Pxy2: Vxy): 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. + + 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 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. + """ + a = np.cos(theta - np.pi / 2) + b = np.sin(theta - np.pi / 2) + x0 = a * rho + y0 = b * rho + pt1 = (int(x0 + 1000 * (-b)), int(y0 + 1000 * (a))) + pt2 = (int(x0 - 1000 * (-b)), int(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: """ Returns y points that lie on line given corresponding x points. diff --git a/opencsp/common/lib/geometry/test/test_LineXY.py b/opencsp/common/lib/geometry/test/test_LineXY.py index 005602a4..cb67dc41 100644 --- a/opencsp/common/lib/geometry/test/test_LineXY.py +++ b/opencsp/common/lib/geometry/test/test_LineXY.py @@ -206,5 +206,39 @@ def test_slope(self): 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() From 8ec7a3621a7ac496d2720146ec59fa1a99638370 Mon Sep 17 00:00:00 2001 From: bbean Date: Fri, 9 Aug 2024 09:04:31 -0600 Subject: [PATCH 3/4] fix LineXY.from_rho_theta(), update to match definition from OpenCV Hough Transforms --- opencsp/common/lib/geometry/LineXY.py | 21 ++++++++++++------- .../common/lib/geometry/test/test_LineXY.py | 16 ++++++++++++++ 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/opencsp/common/lib/geometry/LineXY.py b/opencsp/common/lib/geometry/LineXY.py index f31496f2..cfb982fa 100644 --- a/opencsp/common/lib/geometry/LineXY.py +++ b/opencsp/common/lib/geometry/LineXY.py @@ -257,7 +257,10 @@ def from_two_points(cls, Pxy1: Vxy, Pxy2: Vxy): @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. + 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 ---------- @@ -265,17 +268,19 @@ def from_rho_theta(cls, rho: float, theta: float) -> "LineXY": 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 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. + 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 - np.pi / 2) - b = np.sin(theta - np.pi / 2) + a = np.cos(theta) + b = np.sin(theta) x0 = a * rho y0 = b * rho - pt1 = (int(x0 + 1000 * (-b)), int(y0 + 1000 * (a))) - pt2 = (int(x0 - 1000 * (-b)), int(y0 - 1000 * (a))) + 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 diff --git a/opencsp/common/lib/geometry/test/test_LineXY.py b/opencsp/common/lib/geometry/test/test_LineXY.py index cb67dc41..c676f269 100644 --- a/opencsp/common/lib/geometry/test/test_LineXY.py +++ b/opencsp/common/lib/geometry/test/test_LineXY.py @@ -55,6 +55,22 @@ 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_y_from_x(self): # Line y = -x line = LineXY(1, 1, 0) From 2c42308c3cde5a3fd1863d6f116b95f286d67fa6 Mon Sep 17 00:00:00 2001 From: bbean Date: Fri, 9 Aug 2024 09:10:14 -0600 Subject: [PATCH 4/4] add test for LineXY.from_location_angle() --- opencsp/common/lib/geometry/test/test_LineXY.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/opencsp/common/lib/geometry/test/test_LineXY.py b/opencsp/common/lib/geometry/test/test_LineXY.py index c676f269..98d54c92 100644 --- a/opencsp/common/lib/geometry/test/test_LineXY.py +++ b/opencsp/common/lib/geometry/test/test_LineXY.py @@ -71,6 +71,22 @@ def test_from_rho_theta(self): 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)