Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

JP-3502: Extend snowball flags to next integration #238

Merged
merged 21 commits into from
Feb 15, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,13 @@
Changes to API
--------------

jump
~~~~

- Add in the flagging of groups in the integration after a snowball
occurs. The saturated core of the snowball gets flagged as jump
for a number of groups passed in as a parameter [#238]

Bug Fixes
---------

Expand Down
76 changes: 60 additions & 16 deletions src/stcal/jump/jump.py
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,8 @@
minimum_groups=3,
minimum_sigclip_groups=100,
only_use_ints=True,
mask_persist_grps_next_int=True,
persist_grps_flagged=25,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function now has roughly 30 arguments. Would it be better served by using a class to contain the parameters, rather than continue to add parameters to this function?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that it has a lot of parameters but users need to be able to tweak the parameters. Is there any cost to the large number of parameters?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The cost is readability and maintainability. Grouping various parameters into a class or classes makes the code easier to read, as well as easier to maintain, especially when adding a parameter that needs to be passed to subroutines. If you simply pass parameters, the function signature needs to be changed for all functions that need that parameter. If you use classes, those classes will be passed to the subroutines, so when you add a parameter to a class, no subroutine signature needs to be edited, since the parameter is contained in the class.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While it is true that passing objects instead of primitives has advantages there are also disadvantages. APIs with only primitive parameters are easier to understand where the primitives are coming from. At the extreme, a function that has a single class as input in a black box. Someone browsing the code has more work to do than with primitive parameters.
Also, If we want to consolidate the parameters passed to detect_jump.py, I feel that the change should be in a separate PR. That would give us a working master to compare against.
What do you think?

):
"""
This is the high-level controlling routine for the jump detection process.
Expand Down Expand Up @@ -299,6 +301,8 @@
edge_size=edge_size,
sat_expand=sat_expand,
max_extended_radius=max_extended_radius,
mask_persist_grps_next_int=mask_persist_grps_next_int,
persist_grps_flagged=persist_grps_flagged,
)
log.info("Total snowballs = %i", total_snowballs)
number_extended_events = total_snowballs
Expand Down Expand Up @@ -457,6 +461,8 @@
edge_size=edge_size,
sat_expand=sat_expand,
max_extended_radius=max_extended_radius,
mask_persist_grps_next_int=mask_persist_grps_next_int,
persist_grps_flagged=persist_grps_flagged,
)
log.info("Total snowballs = %i", total_snowballs)
number_extended_events = total_snowballs
Expand Down Expand Up @@ -503,6 +509,8 @@
sat_expand=2,
edge_size=25,
max_extended_radius=200,
mask_persist_grps_next_int=True,
persist_grps_flagged=5,
):
"""
This routine controls the creation of expanded regions that are flagged as
Expand Down Expand Up @@ -540,7 +548,12 @@
required for a snowball to be created
max_extended_radius : int
The largest radius that a snowball or shower can be extended

mask_persist_grps_next_int : bool
The flag to turn on the extension of the flagging of the saturated cores of
snowballs.
persist_grps_flagged : int
How many groups to be flagged when the saturated cores are extended into
subsequent integrations
Returns
-------
Nothing, gdq array is modified.
Expand All @@ -550,8 +563,8 @@

n_showers_grp = []
total_snowballs = 0
nints = gdq.shape[0]
ngrps = gdq.shape[1]
nints, ngrps, nrows, ncols = gdq.shape
persist_jumps = np.zeros(shape=(nints, gdq.shape[2], gdq.shape[3]), dtype=np.uint8)
for integration in range(nints):
for group in range(1, ngrps):
current_gdq = gdq[integration, group, :, :]
Expand All @@ -565,10 +578,8 @@
jump_ellipses = find_ellipses(gdq[integration, group, :, :], jump_flag, min_jump_area)
if sat_required_snowball:
low_threshold = edge_size
nrows = gdq.shape[2]
high_threshold = max(0, nrows - edge_size)

gdq, snowballs = make_snowballs(
gdq, snowballs, persist_jumps = make_snowballs(
gdq,
integration,
group,
Expand All @@ -579,7 +590,9 @@
min_sat_radius_extend,
sat_expand,
sat_flag,
jump_flag,
max_extended_radius,
persist_jumps,
)
else:
snowballs = jump_ellipses
Expand All @@ -596,15 +609,26 @@
num_grps_masked=0,
max_extended_radius=max_extended_radius,
)
return gdq, total_snowballs

# Test to see if the flagging of the saturated cores will be extended into the
# subsequent integrations. Persist_jumps contains all the pixels that were saturated
# in the cores of snowballs.
if mask_persist_grps_next_int:
for intg in range(1, nints):
if persist_grps_flagged >= 1:
last_grp_flagged = min(persist_grps_flagged, ngrps)
gdq[intg, 1:last_grp_flagged, :, :] = np.bitwise_or(gdq[intg, 1:last_grp_flagged, :, :],

Check warning on line 620 in src/stcal/jump/jump.py

View check run for this annotation

Codecov / codecov/patch

src/stcal/jump/jump.py#L618-L620

Added lines #L618 - L620 were not covered by tests
np.repeat(persist_jumps[intg - 1, np.newaxis, :, :],
last_grp_flagged - 1, axis=0))
return gdq, total_snowballs

def extend_saturation(
cube, grp, sat_ellipses, sat_flag, min_sat_radius_extend, expansion=2, max_extended_radius=200
cube, grp, sat_ellipses, sat_flag, jump_flag, min_sat_radius_extend, persist_jumps,
expansion=2, max_extended_radius=200
):
ncols = cube.shape[2]
nrows = cube.shape[1]
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]
Expand All @@ -623,13 +647,29 @@
alpha,
0,
360,
(0, 0, 22),
(0, 0, 22), # in the RGB cube, set blue plane pixels of the ellipse to 22
-1,
)
sat_ellipse = image[:, :, 2]
saty, satx = np.where(sat_ellipse == 22)
# 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
return outcube
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_ellipse = persist_image[:, :, 2]
persist_saty, persist_satx = np.where(persist_ellipse == 22)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is special about 22 here? A comment or declaring a descriptive variable with the value 22 would be good here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I updated the comments in this section. Let me know, if you think it needs more details.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still don't know why you use 22 instead of say 21 or 23.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No real reason that I can remember. There needs to be a number that I can use np.where to find the pixels that are affected. I'm sure there's a more elegant way to do it.

Copy link
Collaborator

@kmacdonald-stsci kmacdonald-stsci Feb 1, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be fine to make a comment saying this number could change and that this number was determined though trial and error and could be changed if needed.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The value of 22 is arbitrary. It just needs to be something that can be tested for to create the mask of pixels within the ellipse.

persist_jumps[persist_saty, persist_satx] = jump_flag
return outcube, persist_jumps


def extend_ellipses(
Expand Down Expand Up @@ -754,7 +794,9 @@
min_sat_radius,
expansion,
sat_flag,
jump_flag,
max_extended_radius,
persist_jumps,
):
# This routine will create a list of snowballs (ellipses) that have the
# center
Expand Down Expand Up @@ -786,17 +828,19 @@
):
snowballs.append(jump)
# extend the saturated ellipses that are larger than the min_sat_radius
gdq[integration, :, :, :] = extend_saturation(
gdq[integration, :, :, :], persist_jumps[integration, :, :] = extend_saturation(
gdq[integration, :, :, :],
group,
sat_ellipses,
sat_flag,
jump_flag,
min_sat_radius,
persist_jumps[integration, :, :],
expansion=expansion,
max_extended_radius=max_extended_radius,
)

return gdq, snowballs
return gdq, snowballs, persist_jumps


def point_inside_ellipse(point, ellipse):
Expand Down
27 changes: 24 additions & 3 deletions tests/test_jump.py
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,7 @@

def test_extend_saturation_simple():
cube = np.zeros(shape=(5, 7, 7), dtype=np.uint8)
persist_jumps = np.zeros(shape=(7, 7), dtype=np.uint8)
grp = 1
min_sat_radius_extend = 1
cube[1, 3, 3] = DQFLAGS["SATURATED"]
Expand All @@ -165,8 +166,9 @@
cube[1, 3, 2] = DQFLAGS["SATURATED"]
cube[1, 2, 2] = DQFLAGS["JUMP_DET"]
sat_circles = find_ellipses(cube[grp, :, :], DQFLAGS["SATURATED"], 1)
new_cube = extend_saturation(
cube, grp, sat_circles, DQFLAGS["SATURATED"], min_sat_radius_extend, expansion=1.1
new_cube, persist_jumps = extend_saturation(
cube, grp, sat_circles, DQFLAGS["SATURATED"], DQFLAGS["JUMP_DET"],
1.1, persist_jumps,
)

assert new_cube[grp, 2, 2] == DQFLAGS["SATURATED"]
Expand Down Expand Up @@ -430,7 +432,26 @@
result = point_inside_ellipse(point, ellipse)
assert result


@pytest.mark.skip(" used for local testing")
def test_flag_persist_groups():
# gdq = fits.getdata("persistgdq.fits")
gdq = np.zeros(shape=(2,2,2,2))
print(gdq.shape[0])
gdq = gdq[:, 0:10, :, :]
total_snowballs = flag_large_events(

Check warning on line 441 in tests/test_jump.py

View check run for this annotation

Codecov / codecov/patch

tests/test_jump.py#L438-L441

Added lines #L438 - L441 were not covered by tests
gdq,
DQFLAGS["JUMP_DET"],
DQFLAGS["SATURATED"],
min_sat_area=1,
min_jump_area=6,
expand_factor=1.9,
edge_size=0,
sat_required_snowball=True,
min_sat_radius_extend=2.5,
sat_expand=1.1,
mask_persist_grps_next_int=True,
persist_grps_flagged=0)
# fits.writeto("persitflaggedgdq.fits", gdq, overwrite=True)
def test_calc_num_slices():
n_rows = 20
max_available_cores = 10
Expand Down
Loading