replace `opencv-python` with `scikit-image` for snowball detection #138

zacharyburnett · 2022-12-29T16:21:56Z

This PR ~~adds fallback functions for~~ replaces ellipse construction using scikit-image ~~and shapely~~, removing opencv-python entirely. The scikit-image ~~+ shapely~~ methods produce slightly different contours, ellipses, and circles than opencv-python.

Checklist

added entry in CHANGES.rst (either in Bug Fixes or Changes to API)
updated relevant tests
updated relevant documentation
updated relevant milestone(s)
added relevant label(s)

zacharyburnett · 2022-12-29T16:22:35Z

Contour-finding also needs to be ported from opencv-python to scikit-image:

ModuleNotFoundError: `opencv-python` required for this functionality

def find_circles(dqplane, bitmask, min_area):
        # Using an input DQ plane this routine will find the groups of pixels with at least the minimum
        # area and return a list of the minimum enclosing circle parameters.
        pixels = np.bitwise_and(dqplane, bitmask)
        if ELLIPSE_PACKAGE == 'opencv-python':
            contours, hierarchy = cv.findContours(pixels, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
            bigcontours = [con for con in contours if cv.contourArea(con) >= min_area]
            circles = [cv.minEnclosingCircle(con) for con in bigcontours]
        else:
            # raise ModuleNotFoundError(ELLIPSE_PACKAGE_WARNING)
           raise ModuleNotFoundError('`opencv-python` required for this functionality')

codecov · 2022-12-29T16:25:34Z

Codecov Report

Attention: Patch coverage is 88.58025% with 37 lines in your changes missing coverage. Please review.

Project coverage is 83.99%. Comparing base (af5aefb) to head (9fca583).
Report is 10 commits behind head on main.

Files	Patch %	Lines
src/stcal/jump/jump.py	77.41%	21 Missing ⚠️
src/stcal/jump/circle.py	88.04%	11 Missing ⚠️
tests/test_circle.py	95.08%	3 Missing ⚠️
tests/test_jump.py	97.43%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #138      +/-   ##
==========================================
- Coverage   85.18%   83.99%   -1.19%     
==========================================
  Files          35       37       +2     
  Lines        6797     7211     +414     
==========================================
+ Hits         5790     6057     +267     
- Misses       1007     1154     +147

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

zacharyburnett · 2022-12-29T16:30:59Z

As described in #126 (comment), there are some differences in the methods used to construct ellipses between the two packages which result in different image masks. Perhaps a warning should make this more clear.

looking into replacing shape construction with scikit-image, it appears there is a crucial difference in the drawing method:

import cv2
import numpy
from matplotlib import pyplot
from skimage import draw

image_0 = numpy.full((10, 10), fill_value=0, dtype=numpy.uint8)

center = (4, 4)
radius = 4

disk = draw.disk(center, radius)
image_1 = image_0.copy()
image_1[disk] = 1

image_2 = cv2.circle(image_0.copy(), center, radius, 1, thickness=-1)

figure, axes = pyplot.subplots(nrows=1, ncols=2)

axes[0].imshow(image_1)
axes[0].set_title('scikit-image')

axes[1].imshow(image_2)
axes[1].set_title('opencv')

pyplot.show()

zacharyburnett · 2022-12-30T16:44:38Z

It looks like the contours found by scikit-image are consistently larger, by about half a unit, than those found by opencv:

>       assert circle[0][1] == pytest.approx(1.0, 1e-3)
E       assert 1.5 == 1.0 ± 1.0e-03

I assume this is a difference in methodology in thinking about pixels as blocks vs as points.

zacharyburnett · 2022-12-30T18:22:09Z

expanding the ellipse test case to a more oblong shape reveals a few edge cases in the scikit-image method; I'll convert this to draft for now and work on that

zacharyburnett · 2022-12-30T21:00:13Z

added some more complexity to the test shape and fixing issues in the ellipse fitting:

zacharyburnett · 2022-12-30T21:19:43Z

tests/test_jump.py

@@ -40,23 +45,25 @@ def test_find_simple_circle():
    plane[2, 3] = DQFLAGS['JUMP_DET']
    plane[2, 1] = DQFLAGS['JUMP_DET']
    circle = find_circles(plane, DQFLAGS['JUMP_DET'], 1)
+    assert circle[0][0] == pytest.approx((2, 2))
    assert circle[0][1] == pytest.approx(1.0, 1e-3)


here is where the fit circle from the scikit-image contour differs from the opencv-python circle; the radius is 1.5 instead of 1

zacharyburnett · 2023-01-03T13:58:23Z

Currently, the .minAreaRect() is replaced by .minimum_rotated_rectangle from Shapely; this algorithm could be written without Shapely, like in circle.py, to remove the need for the additional dependency

src/stcal/jump/circle.py

zacharyburnett · 2023-01-04T16:20:47Z

we should also get input from INS on whether this alternative is acceptable

zacharyburnett · 2024-05-01T14:47:54Z

I removed opencv-python entirely in this PR. We can also look into alternative methods using scipy and binary closure / dilation + erosion; however, this change will solve our installation issues with opencv for the short term, until we can figure that out.

zacharyburnett

the formatter made some irrelevant edits, so here are the relevant changes to this PR

zacharyburnett · 2024-06-06T15:40:53Z

src/stcal/jump/jump.py

+            ellipse_mask = skimage.draw.ellipse(
+                r=round(ceny),
+                c=round(cenx),
+                r_radius=round(axis1/2),
+                c_radius=round(axis2/2),
+                shape=image.shape,
+                rotation=alpha,
            )
+            image[*ellipse_mask, 2] = 22


replaces cv.ellipse with a mask derived from skimage.draw.ellipse, then applies Blue 22 to the image array (why is the blue channel set to 22 if it's discarded right after? can this just be a simple mask?)

No real reason. It's hack to create a mask of the ellipse. Sorry that's it is confusing.

zacharyburnett · 2024-06-06T15:41:44Z

src/stcal/jump/jump.py

+            ellipse_mask = skimage.draw.ellipse(
+                r=round(ceny),
+                c=round(cenx),
+                r_radius=round(ellipse[1][0] / 2),
+                c_radius=round(ellipse[1][1] / 2),
+                shape=persist_image.shape,
+                rotation=alpha,
            )
+            persist_image[*ellipse_mask, 2] = 22


same as above; why is this blue?

Same reason, a hack.

zacharyburnett · 2024-06-06T15:42:38Z

src/stcal/jump/jump.py

+        center = (round(ceny), round(cenx))
+        axes = (round(axis1 / 2), round(axis2 / 2))
+        rr, cc = skimage.draw.ellipse(*center, *axes, shape=image.shape, rotation=alpha)
+        image[rr, cc, 2] = 4


is there a reason why this is rounded to the nearest integer?

If I remember right the routine is expecting integers. My guess was this is just an interface to C code.
Eventually we have to more to integers to mask the pixels anyway.

zacharyburnett · 2024-06-06T15:43:55Z

src/stcal/jump/jump.py

+def minimum_bounding_rectangle(points: np.ndarray) -> np.ndarray:
+    """
+    Find the smallest bounding rectangle for a set of points.
+    Returns a set of points representing the corners of the bounding box.
+    https://stackoverflow.com/questions/13542855/algorithm-to-find-the-minimum-area-rectangle-for-given-points-in-order-to-comput
+
+    :param points: an nx2 matrix of coordinates
+    :rval: an nx2 matrix of coordinates
+    """
+    pi2 = np.pi / 2.0
+
+    # get the convex hull for the points
+    hull_points = points[ConvexHull(points).vertices]
+
+    # calculate edge angles
+    edges = np.zeros((len(hull_points) - 1, 2))
+    edges = hull_points[1:] - hull_points[:-1]
+
+    angles = np.zeros(len(edges))
+    angles = np.arctan2(edges[:, 1], edges[:, 0])
+
+    angles = np.abs(np.mod(angles, pi2))
+    angles = np.unique(angles)
+
+    # find rotation matrices
+    # XXX both work
+    rotations = np.vstack(
+        [np.cos(angles), np.cos(angles - pi2), np.cos(angles + pi2), np.cos(angles)]
+    ).T
+    #     rotations = np.vstack([
+    #         np.cos(angles),
+    #         -np.sin(angles),
+    #         np.sin(angles),
+    #         np.cos(angles)]).T
+    rotations = rotations.reshape((-1, 2, 2))
+
+    # apply rotations to the hull
+    rot_points = np.dot(rotations, hull_points.T)
+
+    # find the bounding points
+    min_x = np.nanmin(rot_points[:, 0], axis=1)
+    max_x = np.nanmax(rot_points[:, 0], axis=1)
+    min_y = np.nanmin(rot_points[:, 1], axis=1)
+    max_y = np.nanmax(rot_points[:, 1], axis=1)
+
+    # find the box with the best area
+    areas = (max_x - min_x) * (max_y - min_y)
+    best_idx = np.argmin(areas)
+
+    # return the best box
+    x1 = max_x[best_idx]
+    x2 = min_x[best_idx]
+    y1 = max_y[best_idx]
+    y2 = min_y[best_idx]
+    r = rotations[best_idx]
+
+    rval = np.zeros((4, 2))


this function replaces shapely.geometry.Polygon.minimum_bounding_box

Remember that the goal is to generate an ellipse. The ellipse can be specified by the minimum bounding box. I have to say that the format that CV uses for the bounding ellipse is very confusing. Hopefully, your version will be easier to understand.

zacharyburnett · 2024-06-06T15:44:12Z

src/stcal/jump/jump.py

+def area_of_polygon(xy: np.ndarray) -> float:
+    """
+    apply shoelace algorithm on collection of xy vertex pairs
+    https://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates
+    """
+    return 0.5 * np.abs(
+        np.dot(xy[:, 0], np.roll(xy[:, 1], 1)) - np.dot(xy[:, 1], np.roll(xy[:, 0], 1))
+    )


this function replaces shapely.geometry.Polygon.area

zacharyburnett · 2024-06-06T15:44:34Z

src/stcal/jump/jump.py

+    )
+
+
+def find_circles(dqplane: np.ndarray, bitmask: np.ndarray, min_area: float) -> list[Circle]:


added type hints here

zacharyburnett · 2024-06-06T15:47:16Z

src/stcal/jump/jump.py

+            ellipses = [
+                (
+                    tuple(np.mean(rectangle[[0, 2], :], axis=0)),
+                    tuple(np.hypot(*np.diff(rectangle[[0, 1, 2], :], axis=0))),
+                    np.degrees(np.arctan2(*np.flip(np.diff(rectangle[[3, 0], :], axis=0)[0]))),
+                )
+                for rectangle in rectangles
+            ]


this is where the axis differences comes from; the minor and major radii are larger than those from opencv by 1.5 (?) units

Well, that's odd. The major and minor axis are what they are. Have you tried an input ellipse and see if you get the same ellipse out?

zacharyburnett · 2024-06-06T16:20:56Z

differences in the cube values in test_flag_large_events_withsnowball (opencv on the left, skimage on the right)

zacharyburnett · 2024-06-06T16:51:27Z

comparison of the input data (left) and the data quality layer (right) from test_find_faint_extended (opencv on the left, skimage on the right):

it appears that skimage flags less pixels than the opencv method

zacharyburnett · 2024-06-06T16:53:19Z

to summarize, two tests still show differences between the two methods; test_flag_large_events_withsnowball and test_find_faint_extended. I've added comparison images, above.

mwregan2 · 2024-06-06T16:55:02Z

I don’t understand what I’m looking at From: Zach Burnett ***@***.***> Reply-To: spacetelescope/stcal ***@***.***> Date: Thursday, June 6, 2024 at 12:51 PM To: spacetelescope/stcal ***@***.***> Cc: Michael Regan ***@***.***>, Review requested ***@***.***> Subject: Re: [spacetelescope/stcal] replace `opencv-python` with `scikit-image` for snowball detection (PR #138) comparison of results from test_find_faint_extended (opencv on the left, skimage on the right): opencv_test_find_faint_extended.png (view on web)<https://urldefense.com/v3/__https:/github.com/spacetelescope/stcal/assets/16024299/e20bc052-a1f2-454f-b761-103781f487cc__;!!CrWY41Z8OgsX0i-WU-0LuAcUu2o!zoeP76EdhjR4CrgsLzrhmnUTWGdhHtDYy_1OXPuqbmu5Stz2qu4d-044BXnXDpUu618wwMy2DxTWaDEM612-mKU$> skimage_test_find_faint_extended.png (view on web)<https://urldefense.com/v3/__https:/github.com/spacetelescope/stcal/assets/16024299/2cc4a402-1471-4d38-859c-4739cd95022b__;!!CrWY41Z8OgsX0i-WU-0LuAcUu2o!zoeP76EdhjR4CrgsLzrhmnUTWGdhHtDYy_1OXPuqbmu5Stz2qu4d-044BXnXDpUu618wwMy2DxTWaDEM9oAMH2g$> — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/spacetelescope/stcal/pull/138*issuecomment-2152981460__;Iw!!CrWY41Z8OgsX0i-WU-0LuAcUu2o!zoeP76EdhjR4CrgsLzrhmnUTWGdhHtDYy_1OXPuqbmu5Stz2qu4d-044BXnXDpUu618wwMy2DxTWaDEM2pUrQKc$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AF6KJLWWIFMVZ4IGFJWRRYTZGCHSNAVCNFSM6AAAAAATMHDDFCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJSHE4DCNBWGA__;!!CrWY41Z8OgsX0i-WU-0LuAcUu2o!zoeP76EdhjR4CrgsLzrhmnUTWGdhHtDYy_1OXPuqbmu5Stz2qu4d-044BXnXDpUu618wwMy2DxTWaDEML_d3XiM$>. You are receiving this because your review was requested.Message ID: ***@***.***>

zacharyburnett · 2024-06-06T16:57:26Z

I don’t understand what I’m looking at

Sorry about that! I'll update the comments to be more descriptive

…x errors

perrygreenfield · 2024-06-06T17:53:54Z

To use your first simple example. I find that if I give the radius for draw.disk as 4.0001 I get the same image as cv2. So it wouldn't surprise me that in one the test is an < test and the other <=. Or something like that. Adding a small epsilon may do it.

zacharyburnett · 2024-06-07T18:32:02Z

a comparison of the find_contours functions between opencv and skimage:

and the minimum_bounding_box analogues:

zacharyburnett self-assigned this Dec 29, 2022

github-actions bot added automation installation jump testing labels Dec 29, 2022

zacharyburnett mentioned this pull request Dec 29, 2022

Opencv dependency in stcal being optional causes uncaught error in jump step spacetelescope/jwst#7409

Closed

zacharyburnett force-pushed the scikit-image_ellipses branch from 349d600 to ca9758c Compare December 30, 2022 15:58

zacharyburnett marked this pull request as ready for review December 30, 2022 16:06

zacharyburnett requested review from mwregan2 and hbushouse December 30, 2022 16:06

zacharyburnett marked this pull request as draft December 30, 2022 16:32

zacharyburnett force-pushed the scikit-image_ellipses branch from b73fd2d to 0c3d1ad Compare December 30, 2022 16:48

zacharyburnett marked this pull request as ready for review December 30, 2022 16:56

zacharyburnett requested a review from kmacdonald-stsci December 30, 2022 17:12

zacharyburnett marked this pull request as draft December 30, 2022 18:22

zacharyburnett force-pushed the scikit-image_ellipses branch from e7e7235 to 8793d4a Compare December 30, 2022 21:16

zacharyburnett commented Dec 30, 2022

View reviewed changes

kmacdonald-stsci approved these changes Jan 3, 2023

View reviewed changes

zacharyburnett force-pushed the scikit-image_ellipses branch from f5b745e to 6775cf3 Compare January 4, 2023 17:27

zacharyburnett requested a review from nden January 5, 2023 18:42

zacharyburnett force-pushed the scikit-image_ellipses branch from 6775cf3 to bbc4608 Compare January 25, 2023 13:53

zacharyburnett requested a review from jdavies-st February 3, 2023 18:02

update change log entry

7499170

zacharyburnett force-pushed the scikit-image_ellipses branch from cd1e6d2 to 7499170 Compare May 1, 2024 14:45

zacharyburnett marked this pull request as ready for review May 1, 2024 14:48

zacharyburnett requested a review from a team as a code owner May 1, 2024 14:48

zacharyburnett changed the title ~~add fallback to scikit-image for ellipse functions~~ replace opencv-python with scikit-image for ellipse functions and snowball detection May 1, 2024

zacharyburnett force-pushed the scikit-image_ellipses branch from b651c39 to 4d811c5 Compare June 4, 2024 18:04

fixes

6f77bc0

zacharyburnett force-pushed the scikit-image_ellipses branch from 4d811c5 to 6f77bc0 Compare June 5, 2024 13:37

zacharyburnett added 2 commits June 5, 2024 13:38

remove usage of shapely

74bf216

pin scikit-image

2e89556

zacharyburnett force-pushed the scikit-image_ellipses branch from 2ceb67d to 2e89556 Compare June 5, 2024 19:05

replace other usages of cv.ellipse

0a7ca88

zacharyburnett changed the title ~~replace opencv-python with scikit-image for ellipse functions and snowball detection~~ replace opencv-python with scikit-image for snowball detection Jun 6, 2024

zacharyburnett added 2 commits June 6, 2024 10:33

convert to tuples to fix equality comparison

c8a2dd5

Merge branch 'main' into scikit-image_ellipses

6500702

zacharyburnett commented Jun 6, 2024

View reviewed changes

resolve negative angle issue

f7786d9

zacharyburnett force-pushed the scikit-image_ellipses branch from a264bf9 to f7786d9 Compare June 6, 2024 15:51

remove extranneous function definition

dfee466

avoid partial index tuple expansion to satisfy Python 3.10 test synta…

94dee71

…x errors

adjust expected test results to match skimage ellipses

9fca583

zacharyburnett mentioned this pull request Jun 7, 2024

remove OpenCV from dependencies spacetelescope/jwst#8544

Draft

8 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

replace `opencv-python` with `scikit-image` for snowball detection #138

replace `opencv-python` with `scikit-image` for snowball detection #138

zacharyburnett commented Dec 29, 2022 •

edited

zacharyburnett commented Dec 29, 2022 •

edited

codecov bot commented Dec 29, 2022 •

edited

zacharyburnett commented Dec 29, 2022 •

edited

zacharyburnett commented Dec 30, 2022 •

edited

zacharyburnett commented Dec 30, 2022

zacharyburnett commented Dec 30, 2022 •

edited

zacharyburnett Dec 30, 2022

zacharyburnett commented Jan 3, 2023 •

edited

zacharyburnett commented Jan 4, 2023 •

edited

zacharyburnett commented May 1, 2024 •

edited

zacharyburnett left a comment

zacharyburnett Jun 6, 2024

mwregan2 Jun 6, 2024

zacharyburnett Jun 6, 2024

mwregan2 Jun 6, 2024

zacharyburnett Jun 6, 2024

mwregan2 Jun 6, 2024

zacharyburnett Jun 6, 2024

mwregan2 Jun 6, 2024

zacharyburnett Jun 6, 2024

zacharyburnett Jun 6, 2024

zacharyburnett Jun 6, 2024

mwregan2 Jun 6, 2024

zacharyburnett commented Jun 6, 2024 •

edited

zacharyburnett commented Jun 6, 2024 •

edited

zacharyburnett commented Jun 6, 2024

mwregan2 commented Jun 6, 2024 via email

zacharyburnett commented Jun 6, 2024

perrygreenfield commented Jun 6, 2024

zacharyburnett commented Jun 7, 2024 •

edited

		)


		def find_circles(dqplane: np.ndarray, bitmask: np.ndarray, min_area: float) -> list[Circle]:

replace opencv-python with scikit-image for snowball detection #138

Are you sure you want to change the base?

replace opencv-python with scikit-image for snowball detection #138

Conversation

zacharyburnett commented Dec 29, 2022 • edited

zacharyburnett commented Dec 29, 2022 • edited

codecov bot commented Dec 29, 2022 • edited

Codecov Report

zacharyburnett commented Dec 29, 2022 • edited

zacharyburnett commented Dec 30, 2022 • edited

zacharyburnett commented Dec 30, 2022

zacharyburnett commented Dec 30, 2022 • edited

Choose a reason for hiding this comment

zacharyburnett commented Jan 3, 2023 • edited

zacharyburnett commented Jan 4, 2023 • edited

zacharyburnett commented May 1, 2024 • edited

zacharyburnett left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

zacharyburnett commented Jun 6, 2024 • edited

zacharyburnett commented Jun 6, 2024 • edited

zacharyburnett commented Jun 6, 2024

mwregan2 commented Jun 6, 2024 via email

zacharyburnett commented Jun 6, 2024

perrygreenfield commented Jun 6, 2024

zacharyburnett commented Jun 7, 2024 • edited

replace `opencv-python` with `scikit-image` for snowball detection #138

replace `opencv-python` with `scikit-image` for snowball detection #138

zacharyburnett commented Dec 29, 2022 •

edited

zacharyburnett commented Dec 29, 2022 •

edited

codecov bot commented Dec 29, 2022 •

edited

zacharyburnett commented Dec 29, 2022 •

edited

zacharyburnett commented Dec 30, 2022 •

edited

zacharyburnett commented Dec 30, 2022 •

edited

zacharyburnett commented Jan 3, 2023 •

edited

zacharyburnett commented Jan 4, 2023 •

edited

zacharyburnett commented May 1, 2024 •

edited

zacharyburnett commented Jun 6, 2024 •

edited

zacharyburnett commented Jun 6, 2024 •

edited

zacharyburnett commented Jun 7, 2024 •

edited