|
30 | 30 |
|
31 | 31 | import numpy as np |
32 | 32 |
|
33 | | -from lsst.pex.config import Field, Config, ConfigField, DictField, FieldValidationError, ListField |
| 33 | +from lsst.pex.config import Field, Config, ConfigField, DictField, FieldValidationError, ListField, RangeField |
34 | 34 | from lsst.pipe.base import Task |
35 | 35 |
|
36 | 36 | from lsst.afw.geom import SpanSet |
@@ -327,17 +327,42 @@ def run(self, table, exposure, **kwargs): |
327 | 327 | class AdaptiveThresholdBackgroundConfig(SubtractBackgroundConfig): |
328 | 328 | detectedFractionBadMaskPlanes = ListField( |
329 | 329 | "Mask planes to ignore when computing the detected fraction.", dtype=str, |
330 | | - default=["BAD", "EDGE", "NO_DATA"] |
| 330 | + default=["BAD", "EDGE", "NO_DATA"], |
331 | 331 | ) |
332 | 332 | minDetFracForFinalBg = Field( |
333 | 333 | "Minimum detected fraction for the final background.", |
334 | | - dtype=float, default=0.02 |
| 334 | + dtype=float, default=0.02, |
335 | 335 | ) |
336 | 336 | maxDetFracForFinalBg = Field( |
337 | 337 | "Maximum detected fraction for the final background.", |
338 | | - dtype=float, default=0.93 |
| 338 | + dtype=float, default=0.93, |
| 339 | + ) |
| 340 | + initialDilateRadius = RangeField( |
| 341 | + "Number of pixels to dilate the original mask by.", |
| 342 | + dtype=int, default=10, min=1, |
| 343 | + ) |
| 344 | + doCheckPerAmpDetectionFraction = Field( |
| 345 | + "Whether to check per-amplifier detection fractions.", |
| 346 | + dtype=bool, default=True, |
| 347 | + ) |
| 348 | + maxDetectionIterationCount = Field( |
| 349 | + "Maximum number of adaptive detection iterations.", |
| 350 | + dtype=int, default=39, # Original code had 40, but iterated from [0, n] inclusive. |
| 351 | + ) |
| 352 | + detection = ConfigField( |
| 353 | + "Baseline configuration for SourceDetectionTask prior to iteration.", |
| 354 | + SourceDetectionConfig, |
339 | 355 | ) |
340 | 356 |
|
| 357 | + def setDefaults(self): |
| 358 | + super().setDefaults() |
| 359 | + self.detection.doTempLocalBackground = False |
| 360 | + self.detection.nSigmaToGrow = 70.0 |
| 361 | + self.detection.reEstimateBackground = False |
| 362 | + self.detection.includeThresholdMultiplier = 1.0 |
| 363 | + self.detection.thresholdValue = 2.0 # Actually used only as a floor. |
| 364 | + self.detection.thresholdType = "pixel_stdev" |
| 365 | + |
341 | 366 |
|
342 | 367 | class AdaptiveThresholdBackgroundTask(SubtractBackgroundTask): |
343 | 368 | """A background subtraction task that does its own masking of detected |
@@ -369,123 +394,107 @@ def run(self, exposure, background=None, stats=True, statsKeys=None, backgroundT |
369 | 394 | exposure.image.array += background.getImage().array |
370 | 395 |
|
371 | 396 | with self._restore_mask_when_done(exposure) as original_mask: |
372 | | - self._dilate_original_mask(exposure, original_mask) |
373 | | - self._set_adaptive_detection_mask(exposure, median_background) |
| 397 | + dilated_mask = self._dilate_original_mask(exposure, original_mask) |
| 398 | + self._set_adaptive_detection_mask(exposure, median_background, dilated_mask) |
374 | 399 | # Do not pass the original background in, since we want to wholly |
375 | 400 | # replace it. |
376 | 401 | return super().run(exposure=exposure, stats=stats, statsKeys=statsKeys, |
377 | 402 | backgroundToPhotometricRatio=backgroundToPhotometricRatio) |
378 | 403 |
|
379 | 404 | def _dilate_original_mask(self, exposure, original_mask): |
380 | | - nPixToDilate = 10 |
381 | 405 | detected_fraction_orig = self._compute_mask_fraction(exposure.mask) |
382 | 406 | # Dilate the current detected mask planes and don't clear |
383 | 407 | # them in the detection step. |
384 | | - inDilating = True |
385 | | - while inDilating: |
386 | | - dilatedMask = original_mask.clone() |
387 | | - for maskName in self._DETECTED_MASK_PLANES: |
| 408 | + for n_pix_to_dilate in range(self.config.initialDilateRadius, 0, -1): |
| 409 | + dilated_mask = original_mask.clone() |
| 410 | + for mask_name in self._DETECTED_MASK_PLANES: |
388 | 411 | # Compute the grown detection mask plane using SpanSet |
389 | | - detectedMaskBit = dilatedMask.getPlaneBitMask(maskName) |
390 | | - detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit) |
391 | | - detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate) |
392 | | - detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox()) |
| 412 | + detected_mask_bit = dilated_mask.getPlaneBitMask(mask_name) |
| 413 | + detected_mask_span_set = SpanSet.fromMask(dilated_mask, detected_mask_bit) |
| 414 | + detected_mask_span_set = detected_mask_span_set.dilated(n_pix_to_dilate) |
| 415 | + detected_mask_span_set = detected_mask_span_set.clippedTo(dilated_mask.getBBox()) |
393 | 416 | # Clear the detected mask plane |
394 | | - detectedMask = dilatedMask.getMaskPlane(maskName) |
395 | | - dilatedMask.clearMaskPlane(detectedMask) |
| 417 | + detectedMask = dilated_mask.getMaskPlane(mask_name) |
| 418 | + dilated_mask.clearMaskPlane(detectedMask) |
396 | 419 | # Set the mask plane to the dilated one |
397 | | - detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit) |
| 420 | + detected_mask_span_set.setMask(dilated_mask, detected_mask_bit) |
398 | 421 |
|
399 | | - detected_fraction_dilated = self._compute_mask_fraction(dilatedMask) |
400 | | - if detected_fraction_dilated < self.config.maxDetFracForFinalBg or nPixToDilate == 1: |
401 | | - inDilating = False |
402 | | - else: |
403 | | - nPixToDilate -= 1 |
404 | | - exposure.mask = dilatedMask |
| 422 | + detected_fraction_dilated = self._compute_mask_fraction(dilated_mask) |
| 423 | + if detected_fraction_dilated < self.config.maxDetFracForFinalBg or n_pix_to_dilate == 1: |
| 424 | + break |
| 425 | + exposure.mask = dilated_mask |
405 | 426 | self.log.warning("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f", |
406 | 427 | detected_fraction_orig, detected_fraction_dilated) |
407 | 428 | n_above_max_per_amp = -99 |
408 | 429 | highest_detected_fraction_per_amp = float("nan") |
409 | | - doCheckPerAmpDetFraction = True |
410 | | - if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: |
| 430 | + if self.config.doCheckPerAmpDetectionFraction: |
411 | 431 | n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ |
412 | 432 | self._compute_per_amp_fraction(exposure, detected_fraction_dilated) |
413 | 433 | self.log.warning("Dilated mask: n_above_max_per_amp = %d, " |
414 | 434 | "highest_detected_fraction_per_amp = %.3f", |
415 | 435 | n_above_max_per_amp, highest_detected_fraction_per_amp) |
| 436 | + return dilated_mask |
416 | 437 |
|
417 | | - def _set_adaptive_detection_mask(self, exposure, median_background): |
418 | | - inBackgroundDet = True |
| 438 | + def _set_adaptive_detection_mask(self, exposure, median_background, dilated_mask): |
419 | 439 | detected_fraction = 1.0 |
420 | | - maxIter = 40 |
421 | | - nIter = 0 |
422 | | - nFootprintTemp = 1e12 |
423 | | - starBackgroundDetectionConfig = SourceDetectionConfig() |
424 | | - starBackgroundDetectionConfig.doTempLocalBackground = False |
425 | | - starBackgroundDetectionConfig.nSigmaToGrow = 70.0 |
426 | | - starBackgroundDetectionConfig.reEstimateBackground = False |
427 | | - starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0 |
428 | | - starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background) |
429 | | - starBackgroundDetectionConfig.thresholdType = "pixel_stdev" # "stdev" |
| 440 | + n_footprint_tmp = 1e12 |
| 441 | + detection_config = self.config.detection.copy() |
| 442 | + detection_config.thresholdValue = max(detection_config.thresholdValue, 0.2*median_background) |
| 443 | + detection_config.thresholdType = "pixel_stdev" # "stdev" |
430 | 444 |
|
| 445 | + no_zero_det_amps = 0 |
431 | 446 | n_above_max_per_amp = -99 |
432 | 447 | highest_detected_fraction_per_amp = float("nan") |
433 | | - doCheckPerAmpDetFraction = True |
434 | 448 |
|
435 | | - while inBackgroundDet: |
436 | | - currentThresh = starBackgroundDetectionConfig.thresholdValue |
| 449 | + for n_iter in range(0, self.config.maxDetectionIterationCount): |
| 450 | + current_threshold = detection_config.thresholdValue |
437 | 451 | if detected_fraction > self.config.maxDetFracForFinalBg: |
438 | | - starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh |
439 | | - if nFootprintTemp < 3 and detected_fraction > 0.9*self.config.maxDetFracForFinalBg: |
440 | | - starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh |
| 452 | + detection_config.thresholdValue = 1.07*current_threshold |
| 453 | + if n_footprint_tmp < 3 and detected_fraction > 0.9*self.config.maxDetFracForFinalBg: |
| 454 | + detection_config.thresholdValue = 1.2*current_threshold |
441 | 455 | if n_above_max_per_amp > 1: |
442 | | - starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh |
| 456 | + detection_config.thresholdValue = 1.1*current_threshold |
443 | 457 | if detected_fraction < self.config.minDetFracForFinalBg: |
444 | | - starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh |
445 | | - starBackgroundDetectionTask = SourceDetectionTask( |
446 | | - config=starBackgroundDetectionConfig) |
447 | | - tempDetections = starBackgroundDetectionTask.detectFootprints( |
448 | | - exposure=exposure, clearMask=True) |
449 | | - exposure.mask |= dilatedMask |
450 | | - nFootprintTemp = ( |
| 458 | + detection_config.thresholdValue = 0.8*current_threshold |
| 459 | + detection_task = SourceDetectionTask(config=detection_config) |
| 460 | + tempDetections = detection_task.detectFootprints(exposure=exposure, clearMask=True) |
| 461 | + exposure.mask |= dilated_mask |
| 462 | + n_footprint_tmp = ( |
451 | 463 | (len(tempDetections.positive.getFootprints()) if tempDetections is not None else 0) |
452 | 464 | + (len(tempDetections.negative.getFootprints()) if tempDetections.negative is not None else 0) |
453 | 465 | ) |
454 | 466 | detected_fraction = self._compute_mask_fraction(exposure.mask) |
455 | 467 | self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or " |
456 | 468 | "DETECTED_NEGATIVE in star_background_detection = %.3f " |
457 | 469 | "(max is %.3f; min is %.3f)", |
458 | | - nIter, starBackgroundDetectionConfig.thresholdValue, |
| 470 | + n_iter, detection_config.thresholdValue, |
459 | 471 | detected_fraction, self.config.maxDetFracForFinalBg, self.config.minDetFracForFinalBg) |
460 | 472 |
|
461 | 473 | n_amp = len(exposure.detector.getAmplifiers()) |
462 | | - if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: |
| 474 | + if self.config.doCheckPerAmpDetectionFraction: |
463 | 475 | n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ |
464 | 476 | self._compute_per_amp_fraction(exposure, detected_fraction) |
465 | 477 |
|
466 | 478 | if not no_zero_det_amps: |
467 | | - starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh |
468 | | - nIter += 1 |
469 | | - if nIter > maxIter: |
470 | | - inBackgroundDet = False |
| 479 | + detection_config.thresholdValue = 0.95*current_threshold |
471 | 480 |
|
472 | 481 | if (detected_fraction < self.config.maxDetFracForFinalBg and detected_fraction > self.config.minDetFracForFinalBg |
473 | 482 | and n_above_max_per_amp < int(0.75*n_amp) |
474 | 483 | and no_zero_det_amps): |
475 | 484 | if (n_above_max_per_amp < max(1, int(0.15*n_amp)) |
476 | 485 | or detected_fraction < 0.85*self.config.maxDetFracForFinalBg): |
477 | | - inBackgroundDet = False |
| 486 | + break |
478 | 487 | else: |
479 | 488 | self.log.warning("Making small tweak....") |
480 | | - starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh |
| 489 | + detection_config.thresholdValue = 1.05*current_threshold |
481 | 490 | self.log.warning("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp)) |
482 | 491 |
|
483 | 492 | self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f " |
484 | 493 | "(highest per amp section = %.5f)", |
485 | 494 | detected_fraction, highest_detected_fraction_per_amp) |
486 | 495 |
|
487 | 496 | if detected_fraction > self.config.maxDetFracForFinalBg: |
488 | | - exposure.mask = dilatedMask |
| 497 | + exposure.mask = dilated_mask |
489 | 498 | self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE " |
490 | 499 | "was too large in star_background_detection = %.3f (max = %.3f). " |
491 | 500 | "Reverting to dilated mask from PSF detection...", |
|
0 commit comments