diff --git a/changes/138.apichange.rst b/changes/138.apichange.rst new file mode 100644 index 00000000..774052b8 --- /dev/null +++ b/changes/138.apichange.rst @@ -0,0 +1 @@ +replace usage of ``opencv-python`` with analagous functionality from ``scikit-image`` diff --git a/pyproject.toml b/pyproject.toml index 50fdbba5..3d0c1bc2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -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", diff --git a/src/stcal/jump/circle.py b/src/stcal/jump/circle.py new file mode 100644 index 00000000..89a40a23 --- /dev/null +++ b/src/stcal/jump/circle.py @@ -0,0 +1,207 @@ +import random +from typing import Optional, Union + +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 + + @classmethod + def from_points(cls, points: list[tuple[float, float]]) -> 'Circle': + """ + 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)) + + return cls(center=center, radius=radius) + else: + # Convert to float and randomize order + shuffled_points = [(float(point[0]), float(point[1])) for point in points] + random.shuffle(shuffled_points) + + # 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]) + + return circle + + def __getitem__(self, index: int) -> Union[tuple, float]: + if index == 0: + return tuple(self.center) + elif index == 1: + return self.radius + else: + raise IndexError(f'{self.__class__.__name__} index out of range') + + 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) + + def __mul__(self, factor: float) -> 'Circle': + return self.__class__(center=self.center, radius=self.radius + factor) + + def __contains__(self, point: tuple[float, float]) -> bool: + return numpy.hypot(*(numpy.array(point) - self.center)) <= self.radius * RELATIVE_TOLERANCE + + def __eq__(self, other: object) -> bool: + if not isinstance(other, Circle): + return False + return numpy.all(self.center == other.center) and self.radius == other.radius + + 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 \ + numpy.allclose(self.radius, other.radius, atol=delta) + + def __repr__(self) -> str: + return f'{self.__class__.__name__}({self.center}, {self.radius})' + + +def _expand_circle_from_one_point( + known_boundary_point: tuple[float, float], + points: list[tuple[float, float]], +) -> Circle: + """ + 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]) + else: + circle = _expand_circle_from_two_points(known_boundary_point, point, points[: point_index + 1]) + return circle + + +def _expand_circle_from_two_points( + known_boundary_point_a: tuple[float, float], + known_boundary_point_b: tuple[float, float], + points: list[tuple[float, float]], +) -> Circle: + """ + 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]) + + circle = Circle.from_points(known_boundary_points) + left = None + right = None + + # For each point not in the two-point circle + for point in points: + if circle is None or point in circle: + continue + + # 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) + + if circumcircle is None: + continue + + circumcenter_cross_product = _triangle_cross_product((*known_boundary_points, circumcircle.center)) + + if circumcircle_cross_product > 0.0 and \ + ( + left is None or + circumcenter_cross_product > _triangle_cross_product((*known_boundary_points, left.center)) + ): + left = circumcircle + elif circumcircle_cross_product < 0.0 and \ + ( + right is None or + circumcenter_cross_product < _triangle_cross_product((*known_boundary_points, right.center)) + ): + right = circumcircle + + # Select which circle to return + if left is None and right is None: + return circle + elif left is None: + return right if right is not None else Circle(center=(0.0, 0.0), radius=0.0) + elif right is None: + return left + else: + return left if (left.radius <= right.radius) else right + + +def circumcircle_from_points( + 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) + + relative_points = points - incenter + a, b, c = relative_points + + 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 + + relative_circumcenter = numpy.array([ + (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( + 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: + """ + 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]) \ + * (triangle[2][1] - triangle[0][1]) \ + - (triangle[1][1] - triangle[0][1]) \ + * (triangle[2][0] - triangle[0][0]) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 9eaecf66..34896807 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -4,14 +4,15 @@ import warnings import numpy as np -import cv2 as cv -import astropy.stats as stats - -from astropy.convolution import Ring2DKernel -from astropy.convolution import convolve +import skimage.draw +import skimage.measure +from astropy import stats +from astropy.convolution import Ring2DKernel, convolve +from scipy.spatial import ConvexHull from . import constants from . import twopoint_difference as twopt +from .circle import Circle log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) @@ -652,8 +653,6 @@ def extend_saturation( expansion=2, max_extended_radius=200 ): ngroups, nrows, ncols = cube.shape - image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) - persist_image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) outcube = cube.copy() for ellipse in sat_ellipses: ceny = ellipse[0][0] @@ -665,34 +664,26 @@ def extend_saturation( alpha = ellipse[2] axis1 = min(axis1, max_extended_radius) axis2 = min(axis2, max_extended_radius) - image = cv.ellipse( - image, - (round(ceny), round(cenx)), - (round(axis1 / 2), round(axis2 / 2)), - alpha, - 0, - 360, - (0, 0, 22), # in the RGB cube, set blue plane pixels of the ellipse to 22 - -1, + saty, satx = skimage.draw.ellipse( + r=round(ceny), + c=round(cenx), + r_radius=round(axis1/2), + c_radius=round(axis2/2), + shape=(nrows, ncols, 3), + rotation=alpha, ) # Create another non-extended ellipse that is used to create the # persist_jumps for this integration. This will be used to mask groups # in subsequent integrations. - sat_ellipse = image[:, :, 2] # extract the Blue plane of the image - saty, satx = np.where(sat_ellipse == 22) # find all the ellipse pixels in the ellipse outcube[grp:, saty, satx] = sat_flag - persist_image = cv.ellipse( - persist_image, - (round(ceny), round(cenx)), - (round(ellipse[1][0] / 2), round(ellipse[1][1] / 2)), - alpha, - 0, - 360, - (0, 0, 22), - -1, + persist_saty, persist_satx = skimage.draw.ellipse( + r=round(ceny), + c=round(cenx), + r_radius=round(ellipse[1][0] / 2), + c_radius=round(ellipse[1][1] / 2), + shape=(nrows, ncols, 3), + rotation=alpha, ) - persist_ellipse = persist_image[:, :, 2] - persist_saty, persist_satx = np.where(persist_ellipse == 22) persist_jumps[persist_saty, persist_satx] = jump_flag return outcube, persist_jumps @@ -741,16 +732,11 @@ def extend_ellipses( axis1 = min(axis1, max_extended_radius) axis2 = min(axis2, max_extended_radius) alpha = ellipse[2] - image = cv.ellipse( - image, - (round(ceny), round(cenx)), - (round(axis1 / 2), round(axis2 / 2)), - alpha, - 0, - 360, - (0, 0, jump_flag), - -1, - ) + center = (round(ceny), round(cenx)) + axes = (round(axis1 / 2), round(axis2 / 2)) + color = (0, 0, jump_flag) + ellipse_rr, ellipse_cc = skimage.draw.ellipse(*center, *axes, shape=image.shape, rotation=alpha) + image[ellipse_rr, ellipse_cc] = color jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] last_grp = find_last_grp(grp, ngrps, num_grps_masked) @@ -760,7 +746,7 @@ def extend_ellipses( saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) - diff_cube = out_gdq_cube - gdq_cube + # diff_cube = out_gdq_cube - gdq_cube return out_gdq_cube, num_ellipses def find_last_grp(grp, ngrps, num_grps_masked): @@ -786,26 +772,110 @@ def find_last_grp(grp, ngrps, num_grps_masked): last_grp = min(grp + num_grps_masked, ngrps) return last_grp -def find_circles(dqplane, bitmask, min_area): + +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)) + rval[0] = np.dot([x1, y2], r) + rval[1] = np.dot([x2, y2], r) + rval[2] = np.dot([x2, y1], r) + rval[3] = np.dot([x1, y1], r) + + return rval + + +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 float(0.5 * np.abs( + np.dot(xy[:, 0], np.roll(xy[:, 1], 1)) - np.dot(xy[:, 1], np.roll(xy[:, 0], 1)) + )) + + +def find_circles(dqplane: np.ndarray, bitmask: np.ndarray, min_area: float) -> list[Circle]: # 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) - contours, hierarchy = cv.findContours(pixels, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) - bigcontours = [con for con in contours if cv.contourArea(con) >= min_area] - return [cv.minEnclosingCircle(con) for con in bigcontours] + pixels = np.bitwise_and(dqplane, bitmask) if bitmask is not None else dqplane + contours = skimage.measure.find_contours(pixels) + bigcontours = [con for con in contours if area_of_polygon(con) > min_area] + return [circle for circle in (Circle.from_points(con) for con in bigcontours) if circle is not None] -def find_ellipses(dqplane, bitmask, min_area): +def find_ellipses(dqplane: np.ndarray, bitmask: np.ndarray, min_area: float) -> list[tuple[tuple[float, float], tuple[float, float], float]]: # 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 ellipse parameters. - pixels = np.bitwise_and(dqplane, bitmask) - contours, hierarchy = cv.findContours(pixels, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) - bigcontours = [con for con in contours if cv.contourArea(con) > min_area] - # minAreaRect is used because fitEllipse requires 5 points and it is - # possible to have a contour - # with just 4 points. - return [cv.minAreaRect(con) for con in bigcontours] + pixels = np.bitwise_and(dqplane, bitmask) if bitmask is not None else dqplane + + contours = skimage.measure.find_contours(pixels) + bigcontours = [con for con in contours if area_of_polygon(con) > min_area] + rectangles = [ + minimum_bounding_rectangle(con) for con in bigcontours + ] + return [ + ( + tuple(np.flip(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 + ] def make_snowballs( @@ -1035,12 +1105,22 @@ def find_faint_extended( extended_emission[exty, extx] = 1 # find the contours of the extended emission - contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) + contours = skimage.measure.find_contours(extended_emission) # get the contours that are above the minimum size - bigcontours = [con for con in contours if cv.contourArea(con) > min_shower_area] + bigcontours = [con for con in contours if area_of_polygon(con) > min_shower_area] # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse - ellipses = [cv.minAreaRect(con) for con in bigcontours] + rectangles = [ + minimum_bounding_rectangle(con) for con in bigcontours + ] + 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 + ] expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] @@ -1048,7 +1128,14 @@ def find_faint_extended( ncols = plane.shape[1] image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) image2 = np.zeros_like(image) - cv.drawContours(image2, bigcontours, -1, (0, 0, jump_flag), -1) + for contour in bigcontours: + contour_mask = skimage.draw.polygon( + r=contour[:, 1], + c=contour[:, 0], + shape=image2.shape, + ) + image2[*contour_mask, 2] = jump_flag + for ellipse in ellipses: ceny = ellipse[0][0] cenx = ellipse[0][1] @@ -1072,16 +1159,15 @@ def find_faint_extended( axis1 = min(axis1, max_extended_radius) axis2 = min(axis2, max_extended_radius) alpha = ellipse[2] - image = cv.ellipse( - image, - (round(ceny), round(cenx)), - (round(axis1 / 2), round(axis2 / 2)), - alpha, - 0, - 360, - (0, 0, jump_flag), - -1, + ellipse_rr, ellipse_cc = 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_rr, ellipse_cc, 2] = jump_flag if len(ellipses) > 0: # add all the showers for this integration to the list all_ellipses.append([intg, grp, ellipses]) diff --git a/tests/test_circle.py b/tests/test_circle.py new file mode 100644 index 00000000..c78157d2 --- /dev/null +++ b/tests/test_circle.py @@ -0,0 +1,98 @@ +# from https://www.nayuki.io/page/smallest-enclosing-circle + +import random +from typing import Tuple, List + +import pytest + +from stcal.jump.circle import Circle + +RELATIVE_TOLERANCE = 1e-12 + + +@pytest.mark.parametrize('trial', range(10)) +def test_circle_matching_naive_algorithm(trial): + points = _random_points(random.randint(1, 30)) + + reference_circle = _smallest_enclosing_circle_naive(points) + test_circle = Circle.from_points(points) + + assert test_circle.almost_equals(reference_circle, delta=RELATIVE_TOLERANCE) + + +@pytest.mark.parametrize('trial', range(10)) +def test_circle_translation(trial): + points = _random_points(random.randint(1, 300)) + + test_circle = Circle.from_points(points) + + dx = random.gauss(0, 1) + dy = random.gauss(0, 1) + translated_points = [(x + dx, y + dy) for (x, y) in points] + + translated_circle = Circle.from_points(translated_points) + reference_circle = test_circle + (dx, dy) + + assert translated_circle.almost_equals(reference_circle, delta=RELATIVE_TOLERANCE) + + +@pytest.mark.parametrize('trial', range(10)) +def test_circle_scaling(trial): + points = _random_points(random.randint(1, 300)) + + test_circle = Circle.from_points(points) + + scale = random.gauss(0, 1) + scaled_points = [(x * scale, y * scale) for (x, y) in points] + + scaled_circle = Circle.from_points(scaled_points) + reference_circle = Circle((test_circle.center[0] * scale, test_circle.center[1] * scale), + test_circle.radius * abs(scale)) + + assert scaled_circle.almost_equals(reference_circle, delta=RELATIVE_TOLERANCE) + + +def _random_points(n: int) -> List[Tuple[float, float]]: + if random.random() < 0.2: # Discrete lattice (to have a chance of duplicated points) + return [(random.randrange(10), random.randrange(10)) for _ in range(n)] + else: # Gaussian distribution + return [(random.gauss(0, 1), random.gauss(0, 1)) for _ in range(n)] + + +def _smallest_enclosing_circle_naive(points: List[Tuple[float, float]]) -> Circle: + """ + Returns the smallest enclosing circle in O(n^4) time using the naive algorithm. + """ + + # Degenerate cases + if len(points) == 0: + return None + elif len(points) == 1: + return Circle(points[0], 0.0) + + # Try all unique pairs + result = None + for i in range(len(points)): + p = points[i] + for j in range(i + 1, len(points)): + q = points[j] + c = Circle.from_points([p, q]) + if (result is None or c.radius < result.radius) and all(r in c for r in points): + result = c + if result is not None: + return result # This optimization is not mathematically proven + + # Try all unique triples + for i in range(len(points)): + p = points[i] + for j in range(i + 1, len(points)): + q = points[j] + for k in range(j + 1, len(points)): + r = points[k] + c = Circle.from_points([p, q, r]) + if c is not None and (result is None or c.radius < result.radius) and all(s in c for s in points): + result = c + + if result is None: + raise AssertionError() + return result diff --git a/tests/test_jump.py b/tests/test_jump.py index e0bf2b62..2c33cd53 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -1,16 +1,18 @@ import numpy as np import pytest from astropy.io import fits + from stcal.jump.jump import ( + area_of_polygon, calc_num_slices, + detect_jumps, extend_saturation, + find_circles, find_ellipses, find_faint_extended, + find_last_grp, flag_large_events, point_inside_ellipse, - find_first_good_group, - detect_jumps, - find_last_grp ) DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8, @@ -189,6 +191,26 @@ def test_find_simple_ellipse(): assert ellipse[0][0] == pytest.approx((2.5, 2.0)) # center +def test_area_of_polygon(): + x_1 = np.arange(0, 1, 0.001) + polygon_1 = np.array([x_1, np.sqrt(1 - x_1**2)]).T + polygon_2 = np.array([[-100, 0], [100, 0], [100, 150], [-100, 150], [-100, 0]]) + assert area_of_polygon(polygon_1) == pytest.approx(0.26353377782163534, 1e-4) + assert area_of_polygon(polygon_2) == 30000.0 + + +def test_find_simple_circle(): + plane = np.zeros(shape=(5, 5), dtype=np.uint8) + plane[2, 2] = DQFLAGS["JUMP_DET"] + plane[3, 2] = DQFLAGS["JUMP_DET"] + plane[1, 2] = DQFLAGS["JUMP_DET"] + plane[2, 3] = DQFLAGS["JUMP_DET"] + plane[2, 1] = DQFLAGS["JUMP_DET"] + circles = find_circles(plane, DQFLAGS["JUMP_DET"], 1) + assert circles[0][0] == pytest.approx((2, 2)) + assert circles[0][1] == pytest.approx(1.5, 1e-3) + + def test_find_ellipse2(): plane = np.zeros(shape=(5, 5), dtype=np.uint8) plane[1, :] = [0, DQFLAGS["JUMP_DET"], DQFLAGS["JUMP_DET"], DQFLAGS["JUMP_DET"], 0] @@ -198,8 +220,8 @@ def test_find_ellipse2(): ellipse = ellipses[0] assert ellipse[0][0] == 2 assert ellipse[0][1] == 2 - assert ellipse[1][0] == 2 - assert ellipse[1][1] == 2 + assert ellipse[1][0] == 3 + assert ellipse[1][1] == 3 assert ellipse[2] == 90.0 @@ -286,7 +308,7 @@ def test_flag_large_events_withsnowball(): ) assert cube[0, 1, 2, 2] == 0 assert cube[0, 1, 3, 5] == 0 - assert cube[0, 2, 0, 0] == 0 + assert cube[0, 2, 0, 0] == DQFLAGS["JUMP_DET"] assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended assert cube[0, 2, 3, 6] == DQFLAGS["JUMP_DET"] @@ -390,13 +412,13 @@ def test_find_faint_extended(): # assert num_showers == 1 assert np.all(gdq[0, 1, 22, 14:23] == 0) assert gdq[0, 1, 16, 18] == DQFLAGS['JUMP_DET'] - assert np.all(gdq[0, 1, 11:22, 16:19] == DQFLAGS["JUMP_DET"]) + assert np.all(gdq[0, 1, 12:21, 16:19] == DQFLAGS["JUMP_DET"]) assert np.all(gdq[0, 1, 22, 16:19] == 0) assert np.all(gdq[0, 1, 10, 16:19] == 0) # Check that the same area is flagged in the first group after the event assert np.all(gdq[0, 2, 22, 14:23] == 0) assert gdq[0, 2, 16, 18] == DQFLAGS['JUMP_DET'] - assert np.all(gdq[0, 2, 11:22, 16:19] == DQFLAGS["JUMP_DET"]) + assert np.all(gdq[0, 2, 12:21, 16:19] == DQFLAGS["JUMP_DET"]) assert np.all(gdq[0, 2, 22, 16:19] == 0) assert np.all(gdq[0, 2, 10, 16:19] == 0) @@ -405,85 +427,6 @@ def test_find_faint_extended(): # Check that the flags are not applied in the 3rd group after the event assert np.all(gdq[0, 4, 12:22, 14:23]) == 0 - def test_find_faint_extended(): - nint, ngrps, ncols, nrows = 1, 66, 5, 5 - data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) - gdq = np.zeros_like(data, dtype=np.uint32) - pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) - pdq[0, 0] = 1 - pdq[1, 1] = 2147483648 - # pdq = np.zeros(shape=(data.shape[2], data.shape[3]), dtype=np.uint8) - gain = 4 - readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain - rng = np.random.default_rng(12345) - data[0, 1:, 14:20, 15:20] = 6 * gain * 6.0 * np.sqrt(2) - data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise - gdq, num_showers = find_faint_extended( - data, - gdq, - pdq, - readnoise * np.sqrt(2), - 1, - 100, - snr_threshold=3, - min_shower_area=10, - inner=1, - outer=2.6, - sat_flag=2, - jump_flag=4, - ellipse_expand=1.1, - num_grps_masked=0, - ) - - -# No shower is found because the event is identical in all ints -def test_find_faint_extended_sigclip(): - nint, ngrps, ncols, nrows = 101, 6, 30, 30 - data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) - gdq = np.zeros_like(data, dtype=np.uint8) - pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint8) - gain = 4 - readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain - rng = np.random.default_rng(12345) - data[0, 1:, 14:20, 15:20] = 6 * gain * 1.7 - data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise - min_shower_area=20 - gdq, num_showers = find_faint_extended( - data, - gdq, - pdq, - readnoise, - 1, - 100, - DQFLAGS, - snr_threshold=1.3, - min_shower_area=min_shower_area, - inner=1, - outer=2, - sat_flag=2, - jump_flag=4, - ellipse_expand=1.1, - num_grps_masked=3, - ) - # Check that all the expected samples in group 2 are flagged as jump and - # that they are not flagged outside - assert num_showers == 0 - assert np.all(gdq[0, 1, 22, 14:23] == 0) - assert np.all(gdq[0, 1, 21, 16:20] == 0) - assert np.all(gdq[0, 1, 20, 15:22] == 0) - assert np.all(gdq[0, 1, 19, 15:23] == 0) - assert np.all(gdq[0, 1, 18, 14:23] == 0) - assert np.all(gdq[0, 1, 17, 14:23] == 0) - assert np.all(gdq[0, 1, 16, 14:23] == 0) - assert np.all(gdq[0, 1, 15, 14:22] == 0) - assert np.all(gdq[0, 1, 14, 16:22] == 0) - assert np.all(gdq[0, 1, 13, 17:21] == 0) - assert np.all(gdq[0, 1, 12, 14:23] == 0) - assert np.all(gdq[0, 1, 12:23, 24] == 0) - assert np.all(gdq[0, 1, 12:23, 13] == 0) - - # Check that the flags are not applied in the 3rd group after the event - assert np.all(gdq[0, 4, 12:22, 14:23]) == 0 # No shower is found because the event is identical in all ints def test_find_faint_extended_sigclip():