Scattered light model #1704
Scattered light model #1704
Conversation
|
Partially addresses issue #1684 |
Codecov ReportAttention:
❗ Your organization needs to install the Codecov GitHub app to enable full functionality. Additional details and impacted files@@ Coverage Diff @@
## develop #1704 +/- ##
===========================================
- Coverage 40.97% 40.90% -0.08%
===========================================
Files 190 193 +3
Lines 43914 44300 +386
===========================================
+ Hits 17994 18119 +125
- Misses 25920 26181 +261 ☔ View full report in Codecov by Sentry. |
There was a problem hiding this comment.
Thanks, @rcooke-ast ! I've got lots of comments for you. The most significant ones (I think) are:
- Update the docs, and more specifically KCWI docs with a description of the new approach
- Add some tests (will there be an associated dev-suite PR?)
- I'm confused about the control flow
- It feels like we can make the scattered-light modeling code more general
I'm fine if you want to push off the latter two points to future PRs.
pypeit/spectrographs/keck_kcwi.py
Outdated
| spl = interpolate.RectBivariateSpline(specvec, spatvec, scale_img, kx=1, ky=1) | ||
| return spl(zoom * (specvec + shft_spec), zoom * (spatvec + shft_spat)) | ||
|
|
||
| def scattered_light(self, frame, offslitmask, tbl, detpad=300): |
There was a problem hiding this comment.
Similar comment: I would generalize and move much/most of this to scattlight.py.
There was a problem hiding this comment.
See comment above. Perhaps we can discuss this further on Slack, but the reason I've implemented it like this is because I suspect the scattered light model requires spectrograph specific code.
pypeit/images/rawimage.py
Outdated
| if not success: | ||
| msgs.warn("Scattered light model failed - using predefined model parameters") | ||
| # Use predefined model parameters | ||
| scatt_img = self.spectrograph.scattered_light_model(scattlight.scattlight_param, self.image[ii, ...]) |
There was a problem hiding this comment.
My understanding of the current flow is:
- If scattered-light frames are provided, fit a scattered light model to them. If they are not, scattered light is never subtracted.
- For the science frames, the code will try to independently model the scattered light in them (only if scattered light frames are provided). If that fails, use the scattered light model from the scattered-light frames.
Is this right? If not, can you describe the flow?
If so, can we enable the scattered-light model to be subtracted from the science data regardless of whether or not you provide scattered-light frames?
Also, if I understand the source of the scattered light for KCWI, at least, (i.e., its intensity is proportional to the light flux through the system), then is it appropriate to subtract the model from the science frames given the significant difference in flux?
There was a problem hiding this comment.
Your understanding of the control flow is correct. But, I'm still deciding what is the best approach. I think it's best to general the model parameters and always use those, unless the user explicitly wants to calculate the scattered light model from each individual frame. I will change the flow such that it's like that, and document it here so the flow is clearer.
If so, can we enable the scattered-light model to be subtracted from the science data regardless of whether or not you provide scattered-light frames?
Yes, that is possible, and I think that may even be better in some cases... still testing...
Also, if I understand the source of the scattered light for KCWI, at least, (i.e., its intensity is proportional to the light flux through the system), then is it appropriate to subtract the model from the science frames given the significant difference in flux?
That is correct. This is precisely what is currently being done. For each setup there's a set of model parameters that:
(1) take the light distribution of the current frame (e.g. a science frame)
(2) based on this light distribution on the detector (and the pre-defined model parameters), calculate the scattered light image
(3) subtract this scattered light image from the current frame.
|
Head's up that I am going to branch off this and try it on ESI. We can then compare to the old IDL code (by having @jhennawi run it) |
# Conflicts: # doc/releases/1.14.1dev.rst
| model : `numpy.ndarray`_ | ||
| Model of the scattered light for the input | ||
| """ | ||
| # Extract the parameters into more conveniently named variables |
There was a problem hiding this comment.
Is there a reason why you don't just use the simpler:
sigmx_g, sigmy_g, sigmx_l, sigmy_l, shft_spec, shft_spat, zoom_spec, zoom_spat, constant, kern_angle, kern_scale, polyterms_spat, polyterms_spec = param
There was a problem hiding this comment.
Not really, I just didn't want a very long line of code, so I put these into sensible groups. Also, polyterms_spat is two parameters, while polyterms_spec can have any number of free parameters, so they have to be separated.
| Raw 2D data frame to be used to compute the scattered light. | ||
| bpm : `numpy.ndarray`_ | ||
| 2D boolean array indicating the bad pixels (True=bad) | ||
| offslitmask : `numpy.ndarray`_ |
There was a problem hiding this comment.
Can you specify image shapes for frame and then state that bpm and offslitmask are the same size.
There was a problem hiding this comment.
Good idea - done!
| else: # Everything else in between | ||
| indx = (spat[None, :] > centrace[:, ii-1, None]) \ | ||
| & (spat[None, :] < centrace[:, ii, None]) | ||
| # Exclude these pixels |
There was a problem hiding this comment.
Since offslitmask is boolean, I think the logic would be clearer if you did something like
good_mask = offslitmask & np.logical_not(indx)
I would also name indx to indicate what it means, i.e. bad_pixels or something like that.
There was a problem hiding this comment.
Good idea - I've changed the logic to make it clearer, as you suggest, and changed the name.
pypeit/core/scattlight.py
Outdated
| model_med = ndimage.median_filter(model, size=(50, 1)) # Median filter to get rid of CRs | ||
| scatt_light_fine = ndimage.gaussian_filter(model_med, sigma=10) # Gaussian filter to smooth median filter | ||
| if debug: | ||
| embed() |
|
|
||
| Parameters | ||
| ---------- | ||
| frame : `numpy.ndarray`_ |
There was a problem hiding this comment.
Yep, good call 👍
pypeit/images/rawimage.py
Outdated
| # Use predefined model parameters | ||
| scatt_img = scattlight.scattered_light_model(this_modpar, _img) | ||
| if debug: | ||
| embed() |
pypeit/images/rawimage.py
Outdated
| scatt_img = scattlight.scattered_light_model(this_modpar, _img) | ||
| if debug: | ||
| embed() | ||
| # import astropy.io.fits as fits |
| f" {tmp[10]}, # Relative kernel scale (>1 means the kernel is more Gaussian, >0 but <1 makes the profile more lorentzian)\n" + \ | ||
| f" {tmp[11]}, {tmp[12]}, # Polynomial terms (coefficients of \"spat\" and \"spat*spec\")\n" + \ | ||
| f" {tmp[13]}, {tmp[14]}, {tmp[15]}]) # Polynomial terms (coefficients of spec**index)\n" | ||
| print(strprint) |
There was a problem hiding this comment.
This printing should be via msgs right?
There was a problem hiding this comment.
This is in the debug statement. It's only used for developers that want a pretty printout of the best-fit parameters to include in the spectrograph file. I've added a more detailed comment now to make this a little more obvious!
| defaults['finecorr_mask'] = None | ||
| dtypes['finecorr_mask'] = [int, list] | ||
| descr['finecorr_mask'] = 'A list containing the inter-slit regions that the user wishes to mask during ' \ | ||
| 'the fine correction to the scattered light. Each integer corresponds to an ' \ |
There was a problem hiding this comment.
Just mention that the defaults mean no additional masking.
There was a problem hiding this comment.
Hi @rcooke-ast excellent work on this PR. The code looks great. I have minor comments throughout that I hope you can look at, but I'm approving anyway. Only one substantive question which is why do you opt to hard-code the parameters in the spectrograph file versus fitting them from flats every time? Is it that the fitting is too time consuming? Also, I'm not sure I followed how things actually work. Do you fit for a template and then determine its amplitude from flats to account for the fact that the scattered light depends on the total amount of light? Sorry if these questions are naive.
There was a problem hiding this comment.
Thanks @jhennawi for the comments! Let me respond to your other queries:
Only one substantive question which is why do you opt to hard-code the parameters in the spectrograph file versus fitting them from flats every time? Is it that the fitting is too time consuming?
The hard-coded parameters in the spectrograph files just provide the initial guesses for the fitting. The scattered light model parameters are always derived using the user's data (unless they request to use the archived model parameters). It's a little time-consuming, so starting somewhere near the expected solution speeds things up and makes sure the solution doesn't get stuck in model space far from the expected "best" solution.
Do you fit for a template and then determine its amplitude from flats to account for the fact that the scattered light depends on the total amount of light?
The model currently assumes that the scattered light is a reprocessed version of the detector image (this seems to work well for both ESI and KCWI). So, with this approach, the model parameters are assumed to be the same for all frames, the only thing that changes is the image of the detector (i.e. the detector image + any fixed [pre-calculated] model parameters determine the total scattered light). This way, we can derive the model parameters using a ScatteredLight calibration frame, and then apply these model parameters (without any adjustment) to any other frame.
I've also updated the calibrations/sacttlight.rst docs to include all of these details, and hopefully that makes things a little clearer.
I hope that makes sense!
| model : `numpy.ndarray`_ | ||
| Model of the scattered light for the input | ||
| """ | ||
| # Extract the parameters into more conveniently named variables |
There was a problem hiding this comment.
Not really, I just didn't want a very long line of code, so I put these into sensible groups. Also, polyterms_spat is two parameters, while polyterms_spec can have any number of free parameters, so they have to be separated.
| Raw 2D data frame to be used to compute the scattered light. | ||
| bpm : `numpy.ndarray`_ | ||
| 2D boolean array indicating the bad pixels (True=bad) | ||
| offslitmask : `numpy.ndarray`_ |
There was a problem hiding this comment.
Good idea - done!
| else: # Everything else in between | ||
| indx = (spat[None, :] > centrace[:, ii-1, None]) \ | ||
| & (spat[None, :] < centrace[:, ii, None]) | ||
| # Exclude these pixels |
There was a problem hiding this comment.
Good idea - I've changed the logic to make it clearer, as you suggest, and changed the name.
|
|
||
| Parameters | ||
| ---------- | ||
| frame : `numpy.ndarray`_ |
There was a problem hiding this comment.
Yep, good call 👍
pypeit/core/scattlight.py
Outdated
| model_med = ndimage.median_filter(model, size=(50, 1)) # Median filter to get rid of CRs | ||
| scatt_light_fine = ndimage.gaussian_filter(model_med, sigma=10) # Gaussian filter to smooth median filter | ||
| if debug: | ||
| embed() |
pypeit/images/rawimage.py
Outdated
| # Use predefined model parameters | ||
| scatt_img = scattlight.scattered_light_model(this_modpar, _img) | ||
| if debug: | ||
| embed() |
pypeit/images/rawimage.py
Outdated
| scatt_img = scattlight.scattered_light_model(this_modpar, _img) | ||
| if debug: | ||
| embed() | ||
| # import astropy.io.fits as fits |
| f" {tmp[10]}, # Relative kernel scale (>1 means the kernel is more Gaussian, >0 but <1 makes the profile more lorentzian)\n" + \ | ||
| f" {tmp[11]}, {tmp[12]}, # Polynomial terms (coefficients of \"spat\" and \"spat*spec\")\n" + \ | ||
| f" {tmp[13]}, {tmp[14]}, {tmp[15]}]) # Polynomial terms (coefficients of spec**index)\n" | ||
| print(strprint) |
There was a problem hiding this comment.
This is in the debug statement. It's only used for developers that want a pretty printout of the best-fit parameters to include in the spectrograph file. I've added a more detailed comment now to make this a little more obvious!
| defaults['finecorr_mask'] = None | ||
| dtypes['finecorr_mask'] = [int, list] | ||
| descr['finecorr_mask'] = 'A list containing the inter-slit regions that the user wishes to mask during ' \ | ||
| 'the fine correction to the scattered light. Each integer corresponds to an ' \ |
|
All comments addressed. Running the dev suite again to make sure nothing is broken 🤞 |
|
Tests pass... merging |
This PR introduces a "proper" scattered light model using the inter-slit pixels (with some padding to reject light near the slit edges). This has been primarily tested with KCWI data (some tweaks are still being implemented), but this PR contains the main base to model the scattered light.