Skip to content

Commit

Permalink
use scikit-image for ellipse functions
Browse files Browse the repository at this point in the history
  • Loading branch information
zacharyburnett committed Aug 27, 2024
1 parent bd5191c commit a1f05b1
Show file tree
Hide file tree
Showing 6 changed files with 490 additions and 133 deletions.
5 changes: 5 additions & 0 deletions CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,11 @@ Changes to API
1.7.3 (2024-07-05)
==================

Changes to API
--------------

- replace usage of ``opencv-python`` with analagous functionality from ``scikit-image`` [#138]

Bug Fixes
---------

Expand Down
1 change: 0 additions & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@ dependencies = [
"scipy >=1.7.2",
"scikit-image>=0.19",
"numpy >=1.21.2",
"opencv-python-headless >=4.6.0.66",
"asdf >=2.15.0",
"gwcs >= 0.18.1",
"tweakwcs >=0.8.8",
Expand Down
205 changes: 205 additions & 0 deletions src/stcal/jump/circle.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,205 @@
import random
from typing import Union, Optional

import numpy

RELATIVE_TOLERANCE = 1 + 1e-14


class Circle:
def __init__(self, center: tuple[float, float], radius: float):
self.center = numpy.array(center)
self.radius = radius

Check warning on line 12 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L11-L12

Added lines #L11 - L12 were not covered by tests

@classmethod
def from_points(cls, points: list[tuple[float, float]]) -> 'Circle' | None:
"""
Returns the smallest circle that encloses all the given points.
from https://www.nayuki.io/page/smallest-enclosing-circle
If 0 points are given, `None` is returned.
If 1 point is given, a circle of radius 0 is returned.
If 2 points are given, a circle with diameter between them is returned.
If 3 or more points are given, uses the algorithm described in this PDF:
https://www.cise.ufl.edu/~sitharam/COURSES/CG/kreveldnbhd.pdf
:param points: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)].
:return: the smallest circle that encloses all points, to within the relative tolerance defined by the class
"""

if len(points) == 0:
return None
elif len(points) == 1:
return cls(center=points[0], radius=0.0)
elif len(points) == 2:
points = numpy.array(points)
center = numpy.mean(points, axis=0)
radius = numpy.max(numpy.hypot(*(center - points).T))

Check warning on line 37 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L30-L37

Added lines #L30 - L37 were not covered by tests

return cls(center=center, radius=radius)

Check warning on line 39 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L39

Added line #L39 was not covered by tests
else:
# Convert to float and randomize order
shuffled_points = [(float(point[0]), float(point[1])) for point in points]
random.shuffle(shuffled_points)

Check warning on line 43 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L42-L43

Added lines #L42 - L43 were not covered by tests

# Progressively add points to circle or recompute circle
circle = None
for (index, point) in enumerate(shuffled_points):
if circle is None or point not in circle:
circle = _expand_circle_from_one_point(point, shuffled_points[:index + 1])

Check warning on line 49 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L46-L49

Added lines #L46 - L49 were not covered by tests

return circle

Check warning on line 51 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L51

Added line #L51 was not covered by tests

def __getitem__(self, index: int) -> Union[tuple, float]:
if index == 0:
return tuple(self.center)
elif index == 1:
return self.radius

Check warning on line 57 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L53-L57

Added lines #L53 - L57 were not covered by tests
else:
raise IndexError(f'{self.__class__.__name__} index out of range')

Check warning on line 59 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L59

Added line #L59 was not covered by tests

def __add__(self, delta: tuple[float, float]) -> 'Circle':
if isinstance(delta, float):
delta = (delta, delta)
return self.__class__(center=self.center + numpy.array(delta), radius=self.radius)

Check warning on line 64 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L61-L64

Added lines #L61 - L64 were not covered by tests

def __mul__(self, factor: float) -> 'Circle':
return self.__class__(center=self.center, radius=self.radius + factor)

Check warning on line 67 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L66-L67

Added lines #L66 - L67 were not covered by tests

def __contains__(self, point: tuple[float, float]) -> bool:
return numpy.hypot(*(numpy.array(point) - self.center)) <= self.radius * RELATIVE_TOLERANCE

Check warning on line 70 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L69-L70

Added lines #L69 - L70 were not covered by tests

def __eq__(self, other: 'Circle') -> bool:
return numpy.all(self.center == other.center) and self.radius == other.radius

Check warning on line 73 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L72-L73

Added lines #L72 - L73 were not covered by tests

def almost_equals(self, other: 'Circle', delta: Optional[float] = None) -> bool:
if delta is None:
delta = RELATIVE_TOLERANCE
return numpy.allclose(self.center, other.center, atol=delta) and \

Check warning on line 78 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L75-L78

Added lines #L75 - L78 were not covered by tests
numpy.allclose(self.radius, other.radius, atol=delta)

def __repr__(self) -> str:
return f'{self.__class__.__name__}({self.center}, {self.radius})'

Check warning on line 82 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L81-L82

Added lines #L81 - L82 were not covered by tests


def _expand_circle_from_one_point(

Check warning on line 85 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L85

Added line #L85 was not covered by tests
known_boundary_point: tuple[float, float],
points: list[tuple[float, float]],
) -> Circle | None:
"""
iteratively expand a circle from one known boundary point to enclose the given set of points
from https://www.nayuki.io/page/smallest-enclosing-circle
"""

circle = Circle(known_boundary_point, 0.0)
for point_index, point in enumerate(points):
if point not in circle:
if circle.radius == 0.0:
circle = Circle.from_points([known_boundary_point, point])

Check warning on line 98 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L94-L98

Added lines #L94 - L98 were not covered by tests
else:
circle = _expand_circle_from_two_points(known_boundary_point, point, points[: point_index + 1])
return circle

Check warning on line 101 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L100-L101

Added lines #L100 - L101 were not covered by tests


def _expand_circle_from_two_points(

Check warning on line 104 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L104

Added line #L104 was not covered by tests
known_boundary_point_a: tuple[float, float],
known_boundary_point_b: tuple[float, float],
points: list[tuple[float, float]],
) -> Circle | None:
"""
iteratively expand a circle from two known boundary points to enclose the given set of points
from https://www.nayuki.io/page/smallest-enclosing-circle
"""

known_boundary_points = numpy.array([known_boundary_point_a, known_boundary_point_b])

Check warning on line 114 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L114

Added line #L114 was not covered by tests

circle = Circle.from_points(known_boundary_points)
left = None
right = None

Check warning on line 118 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L116-L118

Added lines #L116 - L118 were not covered by tests

# For each point not in the two-point circle
for point in points:
if circle is None or point in circle:
continue

Check warning on line 123 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L121-L123

Added lines #L121 - L123 were not covered by tests

# Form a circumcircle and classify it on left or right side
circumcircle_cross_product = _triangle_cross_product((*known_boundary_points, point))
circumcircle = circumcircle_from_points(known_boundary_point_a, known_boundary_point_b, point)

Check warning on line 127 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L126-L127

Added lines #L126 - L127 were not covered by tests

if circumcircle is None:
continue

Check warning on line 130 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L129-L130

Added lines #L129 - L130 were not covered by tests

circumcenter_cross_product = _triangle_cross_product((*known_boundary_points, circumcircle.center))

Check warning on line 132 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L132

Added line #L132 was not covered by tests

if circumcircle_cross_product > 0.0 and \

Check warning on line 134 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L134

Added line #L134 was not covered by tests
(
left is None or
circumcenter_cross_product > _triangle_cross_product((*known_boundary_points, left.center))
):
left = circumcircle
elif circumcircle_cross_product < 0.0 and \

Check warning on line 140 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L139-L140

Added lines #L139 - L140 were not covered by tests
(
right is None or
circumcenter_cross_product < _triangle_cross_product((*known_boundary_points, right.center))
):
right = circumcircle

Check warning on line 145 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L145

Added line #L145 was not covered by tests

# Select which circle to return
if left is None and right is None:
return circle
elif left is None:
return right
elif right is None:
return left

Check warning on line 153 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L148-L153

Added lines #L148 - L153 were not covered by tests
else:
return left if (left.radius <= right.radius) else right

Check warning on line 155 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L155

Added line #L155 was not covered by tests


def circumcircle_from_points(

Check warning on line 158 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L158

Added line #L158 was not covered by tests
a: tuple[float, float],
b: tuple[float, float],
c: tuple[float, float],
) -> Circle | None:
"""
build a circumcircle from three points, using the algorithm from https://en.wikipedia.org/wiki/Circumscribed_circle
from https://www.nayuki.io/page/smallest-enclosing-circle
"""

points = numpy.array([a, b, c])
incenter = ((numpy.min(points, axis=0) + numpy.max(points, axis=0)) / 2)

Check warning on line 169 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L168-L169

Added lines #L168 - L169 were not covered by tests

relative_points = points - incenter
a, b, c = relative_points

Check warning on line 172 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L171-L172

Added lines #L171 - L172 were not covered by tests

intermediate = 2 * (a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]))
if intermediate == 0:
return None

Check warning on line 176 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L174-L176

Added lines #L174 - L176 were not covered by tests

relative_circumcenter = numpy.array([

Check warning on line 178 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L178

Added line #L178 was not covered by tests
(a[0] ** 2 + a[1] ** 2) * (b[1] - c[1]) +
(b[0] ** 2 + b[1] ** 2) * (c[1] - a[1]) +
(c[0] ** 2 + c[1] ** 2) * (a[1] - b[1]),
(a[0] ** 2 + a[1] ** 2) * (c[0] - b[0]) +
(b[0] ** 2 + b[1] ** 2) * (a[0] - c[0]) +
(c[0] ** 2 + c[1] ** 2) * (b[0] - a[0]),
]) / intermediate

return Circle(

Check warning on line 187 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L187

Added line #L187 was not covered by tests
center=relative_circumcenter + incenter,
radius=numpy.max(numpy.hypot(*(relative_circumcenter - relative_points).T)),
)


def _triangle_cross_product(triangle: tuple[tuple[float, float], tuple[float, float], tuple[float, float]]) -> float:

Check warning on line 193 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L193

Added line #L193 was not covered by tests
"""
calculates twice the signed area of the provided triangle
from https://www.nayuki.io/page/smallest-enclosing-circle
:param triangle: three points defining a triangle
:return: twice the signed area of triangle
"""

return (triangle[1][0] - triangle[0][0]) \

Check warning on line 202 in src/stcal/jump/circle.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/circle.py#L202

Added line #L202 was not covered by tests
* (triangle[2][1] - triangle[0][1]) \
- (triangle[1][1] - triangle[0][1]) \
* (triangle[2][0] - triangle[0][0])
Loading

0 comments on commit a1f05b1

Please sign in to comment.