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-2256: Apply median filter to IVM weight #7563

Open
wants to merge 4 commits into
base: main
Choose a base branch
from

Conversation

mcara
Copy link
Member

@mcara mcara commented Apr 21, 2023

Resolves JP-2256

This PR addresses the issue of lower flux in saturated sources in resampled images when resampling is done using the IVM weighting. By applying median filtering to the weight mask we will be able to replace low weights of saturated pixels with a median value of the weight from neighbors (suggested by @schlafly).

Checklist for maintainers

  • added entry in CHANGES.rst within the relevant release section
  • updated or added relevant tests
  • updated relevant documentation
  • added relevant milestone
  • added relevant label(s)
  • ran regression tests, post a link to the Jenkins job below.
    How to run regression tests on a PR
  • Make sure the JIRA ticket is resolved properly

@mcara mcara added the resample label Apr 21, 2023
@mcara mcara added this to the Build 9.3 milestone Apr 21, 2023
@mcara mcara requested a review from hbushouse April 21, 2023 04:46
@mcara mcara requested a review from a team as a code owner April 21, 2023 04:46
@mcara mcara self-assigned this Apr 21, 2023
@schlafly
Copy link

Thanks Mihai. I like this, but I think we only want to touch the pixels that are partially saturated. i.e., some pseudocode might be:

tmp_weight = orig_weight.copy()
msaturated = (dq & PARTIALLY_SATURATED) != 0
tmp_weight[masaturated | (orig_weight == 0)] = np.nan
tmp_weight = astropy.convolution.interpolate_replace_nans(tmp_weight, Gaussian2DKernel(x_stddev=2, y_stddev=2))
interp_weight = orig_weight.copy()
interp_weight[msaturated] = new_weight[msaturated]

i.e., roughly: make a temporary weight array, replacing bad pixels and partially saturated pixels with NaN. Interpolate over those pixels to replace them. Then replace the values of partially saturated pixels with the ones from the interpolated array.

The little bit of extra care around weight zero pixels is intended to make sure that known bad pixel weights don't bleed into good neighbors. interpolate_replace_nans does a Gaussian convolution weighting rather than the median I originally proposed; I don't think the choice of statistic there matters much, and I didn't quickly find a routine that has the NaN filtering I wanted that used medians.

@mcara
Copy link
Member Author

mcara commented Apr 21, 2023

@schlafly Yes, I started implementing alternative approach a little earlier. Will make a new commit soon.

@mcara mcara requested a review from schlafly April 22, 2023 16:04
@codecov
Copy link

codecov bot commented Apr 22, 2023

Codecov Report

Attention: Patch coverage is 91.66667% with 3 lines in your changes are missing coverage. Please review.

Project coverage is 75.67%. Comparing base (2fb073e) to head (d4c79b3).
Report is 17 commits behind head on master.

Files Patch % Lines
jwst/resample/resample_utils.py 91.66% 3 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #7563      +/-   ##
==========================================
+ Coverage   75.31%   75.67%   +0.36%     
==========================================
  Files         474      476       +2     
  Lines       38965    39465     +500     
==========================================
+ Hits        29345    29866     +521     
+ Misses       9620     9599      -21     
Flag Coverage Δ *Carryforward flag
nightly 77.63% <ø> (+0.29%) ⬆️ Carriedforward from 2fb073e

*This pull request uses carry forward flags. Click here to find out more.

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

@mcara
Copy link
Member Author

mcara commented Apr 22, 2023

Regression tests are running here: https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/678/

I expect to see differences in pixel values for tests that use IVM weighting.

Copy link

@schlafly schlafly left a comment

Choose a reason for hiding this comment

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

Looks good to me. I have one comment inline trying to make sure we handle the fully saturated pixels right; we want their weights to remain zero.

if selective_median:
# apply median filter selectively only at
# points of partially saturated sources:
jumps = np.where(model.dq & saturation)

Choose a reason for hiding this comment

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

This looks good to me. For the fully saturated case, do we set the saturation bit and the do_not_use bit, or only the do_not_use bit? If the former, we probably want this selection to be (model.dq & saturation) & ~(model.dq & do_not_use) to avoid treating fully saturated pixels as good.

Copy link
Member Author

Choose a reason for hiding this comment

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

I agree. That would be an oversight on my part. Thanks for catching this. Let me check

Copy link
Member Author

@mcara mcara Apr 26, 2023

Choose a reason for hiding this comment

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

Everything is fine. On line 238 below, this is handled via dqmask itself computed on line 187:

# line 187:
dqmask = build_mask(model.dq, good_bits)
# line 238:
data = ivm_copy[y1:y2, x1:x2][dqmask[y1:y2, x1:x2]]

Copy link
Member Author

Choose a reason for hiding this comment

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

However, if user says that fully saturated pixels are "good" (i.e., should not be ignored), then it will be exactly the situation that you describe but that's what the user wants...

Copy link

@schlafly schlafly May 1, 2023

Choose a reason for hiding this comment

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

Sorry for the delay. I'm missing something, though. Say there's a fully saturated pixel that has SATURATED & DO_NOT_USE set. It gets an entry in jumps, and so we loop over it. As you point out, we don't include the fully bad pixels in the median we use to replace it (because dqmask doesn't include them), but it does look to me like we would replace the fully bad pixel with the median of its neighbors?

Copy link
Member Author

Choose a reason for hiding this comment

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

Ah, I see what you mean (initially I thought you were referring to using DO_NOT_USE pixels in the computation of the median). I still do not think this is an issue because as long as dqmask is 0, the associated weight for that pixel will be set to 0 on line 260 (result = inv_variance * dqmask). However, you are right and thanks for noticing this: it doesn't make sense to compute median if it is going to be ignored.

Copy link

Choose a reason for hiding this comment

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

Ah, right, I missed that. Thumbs up!

Copy link
Member Author

Choose a reason for hiding this comment

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

Latest commit should fix this. Please review

Copy link
Collaborator

@nden nden left a comment

Choose a reason for hiding this comment

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

Just so this doesn't get merged by accident - it needs approval by the CALWG.

@mcara
Copy link
Member Author

mcara commented May 1, 2023

@schlafly @hbushouse The weighting types 'ivm-med' and 'ivm-smed' allow specifying window size for the median, for example: 'ivm-med5', or 'ivm-med7', ... Would that be of interest? (I thought it might, depending on how many adjacent pixels could saturate at once) If so, the current spec specification would not allow specifying window size and it would need to be changed from option to string. An alternative would be to introduce another weight_wndsize or so parameter to the resample step.

@schlafly
Copy link

Sorry for dropping this. I think in general going out 1 or 2 FWHM is going to be what people want. The PSF usually falls off fast enough that you should always be able to get a normal pixel within an FWHM of a partially saturated pixel? But the FWHM size in pixels is variable enough that it's good to make this configurable as you have done. I don't have a sense for what options would best be of use here, though.

@mcara
Copy link
Member Author

mcara commented May 16, 2023

@hbushouse @nden What is your preference for an option to specify median filter size? see #7563 (comment)

@hbushouse
Copy link
Collaborator

@hbushouse @nden What is your preference for an option to specify median filter size? see #7563 (comment)

Given that the whole purpose of this scheme is to reset/rebalance the weighting for pixels that are only partially saturated, I agree that it's probably only necessary to go out 1 or 2 FWHM to find pixels that have no groups saturated. The rationale being that if the core pixel(s) is only partially saturated, the PSF drops off fast enough in the neighboring pixels that the next closest 1 or 2 pix probably aren't saturated at all. I'm OK with adding a param to set this, with a suitable default in the spec. Of course the definition of "suitable" will highly filter (wavelength) dependent, because of the variation in FWHM with wavelength. To properly optimize for each instrument and filter will likely require parameter reference files in CRDS. But that would be up to the teams to create and deliver.

Comment on lines +4 to +9
resample
--------

- Apply a median filter to IVM weight array to better handle saturated
pixels. [#7563]

Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
resample
--------
- Apply a median filter to IVM weight array to better handle saturated
pixels. [#7563]

now that #8671 is merged (switching change log handling to towncrier) this change log entry should be a file in changes/ instead:

echo "Apply a median filter to IVM weight array to better handle saturated pixels." > changes/7563.resample.rst

(new PRs will include instructions on how to do this)

Copy link
Collaborator

Choose a reason for hiding this comment

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

(ignore this if this PR is no longer active or needed)

@hbushouse hbushouse removed their request for review September 23, 2024 18:37
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants