-
Notifications
You must be signed in to change notification settings - Fork 20
DM-48966: Prototype an adaptive threshold detection scheme for initial PSF modeling #1170
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
Conversation
36aca54 to
2aa0488
Compare
TallJimbo
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There's a lot of technical debt here, but given how much testing has gone into this version of the branch, how badly we're doing on main in crowded fields, and the need to start adapting downstream pipeline components to these changes, I don't think we have much choice but to just accept it for now and clean it up afterwards. But we really do need to plan to devote a lot of time to that, because this is not a state we want to be in for long. I'll post more ideas on how to proceed with that on Jira later.
| ) | ||
|
|
||
| result.psf_stars_footprints, result.background, _ = self._compute_psf( | ||
| result.psf_stars_footprints, result.background, _, adaptiveDetResStruct = self._compute_psf( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This task is pretty consistently snake_case; it'd be good to continue that in at least this existing method.
| result.exposure.wcs, | ||
| sourceList=result.psf_stars_footprints, | ||
| include_covariance=self.config.do_include_astrometric_errors, | ||
| ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The "current wcs" is just raw.wcs at this point, right? if we need to initialize the catalog coordinates to something like that before running the astrometry fitter, I think doing something like this belongs in at least the _fit_astrometry method here, and maybe even AstrometryTask.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, actually we now update with the new bulterized gbdes-built camera model from Clare (up at line ~924). This gets printed to the logs:
lsst.calibrateImage INFO: Updating WCS with the provided camera model.
Given that, do you still think I should move it?
Also, this was one of those slightly mysterious cases where all of a sudden the coordinate values in the input catalogs were all NaN and the matcher was failing (but only in its measurement of the fit stats...the matcher seems to be inconsistent with when/why it gets coords via the values in the catalog vs. recomputing them with the WCS) I did try to put this in meas_astrom, but it was difficult to convince myself I wasn't going to break some logic in there (I find a lot of the code there difficult to follow/rationalize...)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, ok, that's reasonable, then. Definitely shouldn't move it to meas_astrom, but putting it in the astrometry method here might still make sense, provided it already is being passed what it needs to do it.
| if background is not None: | ||
| median_background = np.median(background.getImage().array) | ||
| else: | ||
| median_background = float("nan") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Zero seems like a better fallback value, since that's what we've actually subtracted from the image in this round. Or is this logging right before the task is about to fall over and die from a background estimation failure?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Heh, I think I was just trying to make the distinction between "no background supplied" (NaN filling in for what I really mean as N/A) vs. "measured background value had a median of 0" (and, for perhaps obvious reasons, didn't want to make any assumptions about potential future implementations flavor/order of operations in this task 😝). I could go either way...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I prefer zero. Having a NaN here is just going to guarantee that something downstream silently fails, and with zero some future intentional change might just work.
| if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: | ||
| n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ | ||
| self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated, | ||
| detected_mask_planes, bad_mask_planes) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I really don't like the idea of calibrateImage knowing about amplifier boundaries at all. Maybe we don't have any choice because there are effects we can't fix in ISR without going further in ISR than we'd like (i.e. detecting sources), but in that case we really want to wall it off into its own subclasses that different instruments could override more easily.
| if nPeakSrc == 1: | ||
| nIsolated += 1 | ||
| nPeak += nPeakSrc | ||
| minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Lots of magic numbers in this block of code that need to be turned into configuration options at some point, but not on this ticket; I don't think we want them appearing in the CalibrateImageTask configuration class directly.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Agreed...
| adaptiveDetIter = 0 | ||
| threshFactor = 0.2 | ||
| if adaptiveDetResStruct is not None: | ||
| while inAdaptiveDet and adaptiveDetIter < maxAdaptiveDetIter: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think the logic flow would be cleaner if you rewrite this loop as:
for adaptiveDetIter in range(maxAdaptiveDetIter):
and then instead of having the inAdaptiveDet variable, do a break if you've met the criteria for exiting the loop (i.e. invert the sign of your if checks at the end of the loop body).
| # Dilate the current detected mask planes and don't clear | ||
| # them in the detection step. | ||
| inDilating = True | ||
| while inDilating: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is another loop that I think would be worth rewriting as a for with break calls.
| "detection fraction checks.", exposure.detector.getId()) | ||
| return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps | ||
|
|
||
| def _remeasure_star_background(self, result): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I put a lot of work into refactoring most of this method into a task in meas_algorithms on the u/jbosch/DM-48966 branch. I don't think we want to do that on this ticket after all because of the implications for the configuration interface, but you might want to take a look at it and see if you want to instead want to apply some of the same ideas here now (entirely up to you):
- I do think what you have here is too big for one method, and the handful I split it into on my branch seems like a decent split.
- There's are a couple more loops below that would benefit from being rewritten as
forloops thatbreak, instead ofwhileloops with a "are we still looping" boolean variable.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I fixed the loop logic, but may have to punt on the refactoring...
| for mp in detected_mask_planes: | ||
| if mp not in mask.getMaskPlaneDict(): | ||
| mask.addMaskPlane(mp) | ||
| mask &= ~mask.getPlaneBitMask(detected_mask_planes) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think the detection task clears the detection mask planes when it starts already. In terms of making this code clean up after itself on principle, it might make more sense to just do
result.exposure.mask = origMask
(this is done by a context manager in my refactored branch).
tests/test_calibrateImage.py
Outdated
| # Use PCA psf fitter, as psfex fails if there are only 4 stars. | ||
| self.config.psf_measure_psf.psfDeterminer = 'pca' | ||
| # We don't have many test points, so can't match on complicated shapes. | ||
| # self.config.astrometry.sourceSelector["science"].flags.good = ["calib_psf_candidate"] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please delete this commented-out code.
| "was too large in star_background_detection = %.3f (max = %.3f). " | ||
| "Reverting to dilated mask from PSF detection...", | ||
| detected_fraction, maxDetFracForFinalBg) | ||
| star_background = self.star_background.run(exposure=result.exposure).background |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This background subtraction run and the pedestal one below do not correctly take the illumination correction account. We need to:
- pass
backgroundToPhotometricRatioto the background subtraction task in both places; - update the config validation
if self.do_illumination_correctionblock to check the configurations for these background subtasks as well.
This doesn't matter for LSSTCam yet, but I think it might for ComCam.
| minDetFracForFinalBg = 0.02 | ||
| maxDetFracForFinalBg = 0.93 | ||
| doDilatedOrigDetectionMask = True | ||
| if doDilatedOrigDetectionMask: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You're setting a variable to True here and the immediately doing an if check on it; I can guess why (A/B testing?) but it's time to remove the variable and the if block now.
2aa0488 to
3150b1e
Compare
laurenam
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not quite done yet, but I've pushed the changes I've made so far, so may as well send my replies.
| result.exposure.wcs, | ||
| sourceList=result.psf_stars_footprints, | ||
| include_covariance=self.config.do_include_astrometric_errors, | ||
| ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, actually we now update with the new bulterized gbdes-built camera model from Clare (up at line ~924). This gets printed to the logs:
lsst.calibrateImage INFO: Updating WCS with the provided camera model.
Given that, do you still think I should move it?
Also, this was one of those slightly mysterious cases where all of a sudden the coordinate values in the input catalogs were all NaN and the matcher was failing (but only in its measurement of the fit stats...the matcher seems to be inconsistent with when/why it gets coords via the values in the catalog vs. recomputing them with the WCS) I did try to put this in meas_astrom, but it was difficult to convince myself I wasn't going to break some logic in there (I find a lot of the code there difficult to follow/rationalize...)
| if background is not None: | ||
| median_background = np.median(background.getImage().array) | ||
| else: | ||
| median_background = float("nan") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Heh, I think I was just trying to make the distinction between "no background supplied" (NaN filling in for what I really mean as N/A) vs. "measured background value had a median of 0" (and, for perhaps obvious reasons, didn't want to make any assumptions about potential future implementations flavor/order of operations in this task 😝). I could go either way...
| background=background, | ||
| backgroundToPhotometricRatio=background_to_photometric_ratio, | ||
| ) | ||
| # DO WE REALLY NEED THE 2nd ROUND OF PSF FITTING?? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, I updated the comment.
| if nPeakSrc == 1: | ||
| nIsolated += 1 | ||
| nPeak += nPeakSrc | ||
| minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Agreed...
| "detection fraction checks.", exposure.detector.getId()) | ||
| return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps | ||
|
|
||
| def _remeasure_star_background(self, result): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I fixed the loop logic, but may have to punt on the refactoring...
b659745 to
d3aa468
Compare
This is to allow for proper motion corrections.
This implementation is to allow the single frame processing calibration to handle a wide range of "scenes", e.g. crowded and/or high nebulosity ields in addition to the more "typical" regions, by implementing an adaptive and iterative thresholding based on the current detection masks and tuning the thresholds to adapt to the current scene. In general, higher detection thresholds are requird for non-standard fields to avoid large contiguous detection footprints which leave too few good isolated sources for calibration. To achieve the previous behaviour, the adaptive scheme can be turned off by setting the do_adaptive_threshold_detection config to False.
d3aa468 to
2341fdc
Compare
No description provided.