JP-2256: Apply median filter to IVM weight

[mcara]

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

[schlafly]

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]

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

[codecov]

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]

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.

[schlafly]

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.

[schlafly]

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.

[mcara]

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

[mcara]

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]]

[mcara]

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...

[schlafly]

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?

[mcara]

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.

[schlafly]

Ah, right, I missed that. Thumbs up!

[mcara]

Latest commit should fix this. Please review

[nden]

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

[mcara]

@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]

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]

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

[hbushouse]

@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.