23 from scipy.stats
import norm
25 from collections
import defaultdict
30 from lsstDebug
import getDebugFrame
34 from lsst.ip.isr import CrosstalkCalib, IsrProvenance
38 from ._lookupStaticCalibration
import lookupStaticCalibration
40 __all__ = [
"CrosstalkExtractConfig",
"CrosstalkExtractTask",
41 "CrosstalkSolveTask",
"CrosstalkSolveConfig",
42 "MeasureCrosstalkConfig",
"MeasureCrosstalkTask"]
46 dimensions=(
"instrument",
"exposure",
"detector")):
48 name=
"crosstalkInputs",
49 doc=
"Input post-ISR processed exposure to measure crosstalk from.",
50 storageClass=
"Exposure",
51 dimensions=(
"instrument",
"exposure",
"detector"),
56 name=
"crosstalkSource",
57 doc=
"Post-ISR exposure to measure for inter-chip crosstalk onto inputExp.",
58 storageClass=
"Exposure",
59 dimensions=(
"instrument",
"exposure",
"detector"),
65 outputRatios = cT.Output(
66 name=
"crosstalkRatios",
67 doc=
"Extracted crosstalk pixel ratios.",
68 storageClass=
"StructuredDataDict",
69 dimensions=(
"instrument",
"exposure",
"detector"),
71 outputFluxes = cT.Output(
72 name=
"crosstalkFluxes",
73 doc=
"Source pixel fluxes used in ratios.",
74 storageClass=
"StructuredDataDict",
75 dimensions=(
"instrument",
"exposure",
"detector"),
82 self.inputs.discard(
"sourceExp")
86 pipelineConnections=CrosstalkExtractConnections):
87 """Configuration for the measurement of pixel ratios.
92 doc=
"Measure inter-chip crosstalk as well?",
97 doc=
"Minimum level of source pixels for which to measure crosstalk."
102 doc=
"Should saturated pixels be ignored?"
106 default=[
"BAD",
"INTRP"],
107 doc=
"Mask planes to ignore when identifying source pixels."
112 doc=
"Is the input exposure trimmed?"
121 if 'SAT' not in self.
badMaskbadMask:
124 if 'SAT' in self.
badMaskbadMask:
125 self.
badMaskbadMask = [mask
for mask
in self.
badMaskbadMask
if mask !=
'SAT']
129 pipeBase.CmdLineTask):
130 """Task to measure pixel ratios to find crosstalk.
132 ConfigClass = CrosstalkExtractConfig
133 _DefaultName =
'cpCrosstalkExtract'
135 def run(self, inputExp, sourceExps=[]):
136 """Measure pixel ratios between amplifiers in inputExp.
138 Extract crosstalk ratios between different amplifiers.
140 For pixels above ``config.threshold``, we calculate the ratio
141 between each background-subtracted target amp and the source
142 amp. We return a list of ratios for each pixel for each
143 target/source combination, as nested dictionary containing the
148 inputExp : `lsst.afw.image.Exposure`
149 Input exposure to measure pixel ratios on.
150 sourceExp : `list` [`lsst.afw.image.Exposure`], optional
151 List of chips to use as sources to measure inter-chip
156 results : `lsst.pipe.base.Struct`
157 The results struct containing:
159 ``outputRatios`` : `dict` [`dict` [`dict` [`dict` [`list`]]]]
160 A catalog of ratio lists. The dictionaries are
162 outputRatios[targetChip][sourceChip][targetAmp][sourceAmp]
163 contains the ratio list for that combination.
164 ``outputFluxes`` : `dict` [`dict` [`list`]]
165 A catalog of flux lists. The dictionaries are
167 outputFluxes[sourceChip][sourceAmp]
168 contains the flux list used in the outputRatios.
172 The lsstDebug.Info() method can be rewritten for __name__ =
173 `lsst.cp.pipe.measureCrosstalk`, and supports the parameters:
175 debug.display['extract'] : `bool`
176 Display the exposure under consideration, with the pixels used
177 for crosstalk measurement indicated by the DETECTED mask plane.
178 debug.display['pixels'] : `bool`
179 Display a plot of the ratio calculated for each pixel used in this
180 exposure, split by amplifier pairs. The median value is listed
183 outputRatios = defaultdict(
lambda: defaultdict(dict))
184 outputFluxes = defaultdict(
lambda: defaultdict(dict))
186 threshold = self.config.threshold
187 badPixels =
list(self.config.badMask)
189 targetDetector = inputExp.getDetector()
190 targetChip = targetDetector.getName()
193 sourceExtractExps = [inputExp]
194 sourceExtractExps.extend(sourceExps)
196 self.log.
info(
"Measuring full detector background for target: %s", targetChip)
197 targetIm = inputExp.getMaskedImage()
199 detected = targetIm.getMask().getPlaneBitMask(
"DETECTED")
200 bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + [
"DETECTED"])
202 self.
debugViewdebugView(
'extract', inputExp)
204 for sourceExp
in sourceExtractExps:
205 sourceDetector = sourceExp.getDetector()
206 sourceChip = sourceDetector.getName()
207 sourceIm = sourceExp.getMaskedImage()
208 bad = sourceIm.getMask().getPlaneBitMask(badPixels)
209 self.log.
info(
"Measuring crosstalk from source: %s", sourceChip)
211 if sourceExp != inputExp:
213 detected = sourceIm.getMask().getPlaneBitMask(
"DETECTED")
216 ratioDict = defaultdict(
lambda: defaultdict(list))
219 for sourceAmp
in sourceDetector:
220 sourceAmpName = sourceAmp.getName()
221 sourceAmpImage = sourceIm[sourceAmp.getBBox()]
222 sourceMask = sourceAmpImage.mask.array
223 select = ((sourceMask & detected > 0)
224 & (sourceMask & bad == 0)
225 & np.isfinite(sourceAmpImage.image.array))
226 count = np.sum(select)
227 self.log.
debug(
" Source amplifier: %s", sourceAmpName)
229 outputFluxes[sourceChip][sourceAmpName] = sourceAmpImage.image.array[select].tolist()
231 for targetAmp
in targetDetector:
233 targetAmpName = targetAmp.getName()
234 if sourceAmpName == targetAmpName
and sourceChip == targetChip:
235 ratioDict[sourceAmpName][targetAmpName] = []
237 self.log.
debug(
" Target amplifier: %s", targetAmpName)
239 targetAmpImage = CrosstalkCalib.extractAmp(targetIm.image,
240 targetAmp, sourceAmp,
241 isTrimmed=self.config.isTrimmed)
242 ratios = (targetAmpImage.array[select] - bg)/sourceAmpImage.image.array[select]
243 ratioDict[targetAmpName][sourceAmpName] = ratios.tolist()
244 extractedCount += count
247 sourceAmpImage.image.array[select],
248 targetAmpImage.array[select] - bg,
249 sourceAmpName, targetAmpName)
251 self.log.
info(
"Extracted %d pixels from %s -> %s (targetBG: %f)",
252 extractedCount, sourceChip, targetChip, bg)
253 outputRatios[targetChip][sourceChip] = ratioDict
255 return pipeBase.Struct(
261 """Utility function to examine the image being processed.
266 State of processing to view.
267 exposure : `lsst.afw.image.Exposure`
273 display.scale(
'asinh',
'zscale')
274 display.mtv(exposure)
276 prompt =
"Press Enter to continue: "
278 ans = input(prompt).lower()
279 if ans
in (
"",
"c",):
282 def debugPixels(self, stepname, pixelsIn, pixelsOut, sourceName, targetName):
283 """Utility function to examine the CT ratio pixel values.
288 State of processing to view.
289 pixelsIn : `np.ndarray`
290 Pixel values from the potential crosstalk source.
291 pixelsOut : `np.ndarray`
292 Pixel values from the potential crosstalk target.
294 Source amplifier name
296 Target amplifier name
300 import matplotlib.pyplot
as plt
301 figure = plt.figure(1)
304 axes = figure.add_axes((0.1, 0.1, 0.8, 0.8))
305 axes.plot(pixelsIn, pixelsOut / pixelsIn,
'k+')
306 plt.xlabel(
"Source amplifier pixel value")
307 plt.ylabel(
"Measured pixel ratio")
308 plt.title(f
"(Source {sourceName} -> Target {targetName}) median ratio: "
309 f
"{(np.median(pixelsOut / pixelsIn))}")
312 prompt =
"Press Enter to continue: "
314 ans = input(prompt).lower()
315 if ans
in (
"",
"c",):
321 dimensions=(
"instrument",
"detector")):
322 inputRatios = cT.Input(
323 name=
"crosstalkRatios",
324 doc=
"Ratios measured for an input exposure.",
325 storageClass=
"StructuredDataDict",
326 dimensions=(
"instrument",
"exposure",
"detector"),
329 inputFluxes = cT.Input(
330 name=
"crosstalkFluxes",
331 doc=
"Fluxes of CT source pixels, for nonlinear fits.",
332 storageClass=
"StructuredDataDict",
333 dimensions=(
"instrument",
"exposure",
"detector"),
336 camera = cT.PrerequisiteInput(
338 doc=
"Camera the input data comes from.",
339 storageClass=
"Camera",
340 dimensions=(
"instrument",),
342 lookupFunction=lookupStaticCalibration,
345 outputCrosstalk = cT.Output(
347 doc=
"Output proposed crosstalk calibration.",
348 storageClass=
"CrosstalkCalib",
349 dimensions=(
"instrument",
"detector"),
357 if config.fluxOrder == 0:
358 self.inputs.discard(
"inputFluxes")
362 pipelineConnections=CrosstalkSolveConnections):
363 """Configuration for the solving of crosstalk from pixel ratios.
368 doc=
"Number of rejection iterations for final coefficient calculation.",
373 doc=
"Rejection threshold (sigma) for final coefficient calculation.",
378 doc=
"Polynomial order in source flux to fit crosstalk.",
383 doc=
"Filter generated crosstalk to remove marginal measurements.",
388 pipeBase.CmdLineTask):
389 """Task to solve crosstalk from pixel ratios.
391 ConfigClass = CrosstalkSolveConfig
392 _DefaultName =
'cpCrosstalkSolve'
395 """Ensure that the input and output dimensions are passed along.
399 butlerQC : `lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
400 Butler to operate on.
401 inputRefs : `lsst.pipe.base.connections.InputQuantizedConnection`
402 Input data refs to load.
403 ouptutRefs : `lsst.pipe.base.connections.OutputQuantizedConnection`
404 Output data refs to persist.
406 inputs = butlerQC.get(inputRefs)
409 inputs[
'inputDims'] = [exp.dataId.byName()
for exp
in inputRefs.inputRatios]
410 inputs[
'outputDims'] = outputRefs.outputCrosstalk.dataId.byName()
412 outputs = self.
runrun(**inputs)
413 butlerQC.put(outputs, outputRefs)
415 def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, outputDims=None):
416 """Combine ratios to produce crosstalk coefficients.
420 inputRatios : `list` [`dict` [`dict` [`dict` [`dict` [`list`]]]]]
421 A list of nested dictionaries of ratios indexed by target
422 and source chip, then by target and source amplifier.
423 inputFluxes : `list` [`dict` [`dict` [`list`]]]
424 A list of nested dictionaries of source pixel fluxes, indexed
425 by source chip and amplifier.
426 camera : `lsst.afw.cameraGeom.Camera`
428 inputDims : `list` [`lsst.daf.butler.DataCoordinate`]
429 DataIds to use to construct provenance.
430 outputDims : `list` [`lsst.daf.butler.DataCoordinate`]
431 DataIds to use to populate the output calibration.
435 results : `lsst.pipe.base.Struct`
436 The results struct containing:
438 ``outputCrosstalk`` : `lsst.ip.isr.CrosstalkCalib`
439 Final crosstalk calibration.
440 ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
441 Provenance data for the new calibration.
446 Raised if the input data contains multiple target detectors.
450 The lsstDebug.Info() method can be rewritten for __name__ =
451 `lsst.ip.isr.measureCrosstalk`, and supports the parameters:
453 debug.display['reduce'] : `bool`
454 Display a histogram of the combined ratio measurements for
455 a pair of source/target amplifiers from all input
460 calibChip = outputDims[
'detector']
461 instrument = outputDims[
'instrument']
467 self.log.
info(
"Combining measurements from %d ratios and %d fluxes",
468 len(inputRatios), len(inputFluxes)
if inputFluxes
else 0)
470 if inputFluxes
is None:
471 inputFluxes = [
None for exp
in inputRatios]
473 combinedRatios = defaultdict(
lambda: defaultdict(list))
474 combinedFluxes = defaultdict(
lambda: defaultdict(list))
475 for ratioDict, fluxDict
in zip(inputRatios, inputFluxes):
476 for targetChip
in ratioDict:
477 if calibChip
and targetChip != calibChip:
478 raise RuntimeError(
"Received multiple target chips!")
480 sourceChip = targetChip
481 if sourceChip
in ratioDict[targetChip]:
482 ratios = ratioDict[targetChip][sourceChip]
484 for targetAmp
in ratios:
485 for sourceAmp
in ratios[targetAmp]:
486 combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp])
488 combinedFluxes[targetAmp][sourceAmp].extend(fluxDict[sourceChip][sourceAmp])
493 for targetAmp
in combinedRatios:
494 for sourceAmp
in combinedRatios[targetAmp]:
495 self.log.
info(
"Read %d pixels for %s -> %s",
496 len(combinedRatios[targetAmp][sourceAmp]),
497 targetAmp, sourceAmp)
498 if len(combinedRatios[targetAmp][sourceAmp]) > 1:
499 self.
debugRatiosdebugRatios(
'reduce', combinedRatios, targetAmp, sourceAmp)
501 if self.config.fluxOrder == 0:
502 self.log.
info(
"Fitting crosstalk coefficients.")
504 self.config.rejIter, self.config.rejSigma)
506 raise NotImplementedError(
"Non-linear crosstalk terms are not yet supported.")
508 self.log.
info(
"Number of valid coefficients: %d", np.sum(calib.coeffValid))
510 if self.config.doFiltering:
513 self.log.
info(
"Filtering measured crosstalk to remove invalid solutions.")
517 calib.hasCrosstalk =
True
521 calib._detectorId = calibChip
523 calib._detectorName = camera[calibChip].
getName()
524 calib._detectorSerial = camera[calibChip].getSerial()
526 calib._instrument = instrument
527 calib.updateMetadata(setCalibId=
True, setDate=
True)
531 provenance._detectorName = calibChip
533 provenance.fromDataIds(inputDims)
534 provenance._instrument = instrument
535 provenance.updateMetadata()
537 return pipeBase.Struct(
538 outputCrosstalk=calib,
539 outputProvenance=provenance,
543 """Measure crosstalk coefficients from the ratios.
545 Given a list of ratios for each target/source amp combination,
546 we measure a sigma clipped mean and error.
548 The coefficient errors returned are the standard deviation of
549 the final set of clipped input ratios.
553 ratios : `dict` of `dict` of `numpy.ndarray`
554 Catalog of arrays of ratios.
556 Number of rejection iterations.
558 Rejection threshold (sigma).
562 calib : `lsst.ip.isr.CrosstalkCalib`
563 The output crosstalk calibration.
567 The lsstDebug.Info() method can be rewritten for __name__ =
568 `lsst.ip.isr.measureCrosstalk`, and supports the parameters:
570 debug.display['measure'] : `bool`
571 Display the CDF of the combined ratio measurements for
572 a pair of source/target amplifiers from the final set of
573 clipped input ratios.
578 ordering =
list(ratios.keys())
579 for ii, jj
in itertools.product(range(calib.nAmp), range(calib.nAmp)):
583 values = np.array(ratios[ordering[ii]][ordering[jj]])
584 values = values[np.abs(values) < 1.0]
586 calib.coeffNum[ii][jj] = len(values)
589 self.log.
warn(
"No values for matrix element %d,%d" % (ii, jj))
590 calib.coeffs[ii][jj] = np.nan
591 calib.coeffErr[ii][jj] = np.nan
592 calib.coeffValid[ii][jj] =
False
595 for rej
in range(rejIter):
596 lo, med, hi = np.percentile(values, [25.0, 50.0, 75.0])
597 sigma = 0.741*(hi - lo)
598 good = np.abs(values - med) < rejSigma*sigma
599 if good.sum() == len(good):
601 values = values[good]
603 calib.coeffs[ii][jj] = np.mean(values)
604 if calib.coeffNum[ii][jj] == 1:
605 calib.coeffErr[ii][jj] = np.nan
608 calib.coeffErr[ii][jj] = np.std(values) * correctionFactor
609 calib.coeffValid[ii][jj] = (np.abs(calib.coeffs[ii][jj])
610 > calib.coeffErr[ii][jj] / np.sqrt(calib.coeffNum[ii][jj]))
612 if calib.coeffNum[ii][jj] > 1:
613 self.
debugRatiosdebugRatios(
'measure', ratios, ordering[ii], ordering[jj],
614 calib.coeffs[ii][jj], calib.coeffValid[ii][jj])
620 """Correct measured sigma to account for clipping.
622 If we clip our input data and then measure sigma, then the
623 measured sigma is smaller than the true value because real
624 points beyond the clip threshold have been removed. This is a
625 small (1.5% at nSigClip=3) effect when nSigClip >~ 3, but the
626 default parameters for measure crosstalk use nSigClip=2.0.
627 This causes the measured sigma to be about 15% smaller than
628 real. This formula corrects the issue, for the symmetric case
629 (upper clip threshold equal to lower clip threshold).
634 Number of sigma the measurement was clipped by.
638 scaleFactor : `float`
639 Scale factor to increase the measured sigma by.
642 varFactor = 1.0 + (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))
643 return 1.0 / np.sqrt(varFactor)
647 """Apply valid constraints to the measured values.
649 Any measured coefficient that is determined to be invalid is
650 set to zero, and has the error set to nan. The validation is
651 determined by checking that the measured coefficient is larger
652 than the calculated standard error of the mean.
656 inCalib : `lsst.ip.isr.CrosstalkCalib`
657 Input calibration to filter.
661 outCalib : `lsst.ip.isr.CrosstalkCalib`
662 Filtered calibration.
665 outCalib.numAmps = inCalib.numAmps
667 outCalib.coeffs = inCalib.coeffs
668 outCalib.coeffs[~inCalib.coeffValid] = 0.0
670 outCalib.coeffErr = inCalib.coeffErr
671 outCalib.coeffErr[~inCalib.coeffValid] = np.nan
673 outCalib.coeffNum = inCalib.coeffNum
674 outCalib.coeffValid = inCalib.coeffValid
678 def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False):
679 """Utility function to examine the final CT ratio set.
684 State of processing to view.
685 ratios : `dict` of `dict` of `np.ndarray`
686 Array of measured CT ratios, indexed by source/victim
689 Index of the source amplifier.
691 Index of the target amplifier.
692 coeff : `float`, optional
693 Coefficient calculated to plot along with the simple mean.
694 valid : `bool`, optional
695 Validity to be added to the plot title.
699 if i == j
or ratios
is None or len(ratios) < 1:
702 ratioList = ratios[i][j]
703 if ratioList
is None or len(ratioList) < 1:
706 mean = np.mean(ratioList)
707 std = np.std(ratioList)
708 import matplotlib.pyplot
as plt
709 figure = plt.figure(1)
711 plt.hist(x=ratioList, bins=len(ratioList),
712 cumulative=
True, color=
'b', density=
True, histtype=
'step')
713 plt.xlabel(
"Measured pixel ratio")
714 plt.ylabel(f
"CDF: n={len(ratioList)}")
715 plt.xlim(np.percentile(ratioList, [1.0, 99]))
716 plt.axvline(x=mean, color=
"k")
717 plt.axvline(x=coeff, color=
'g')
718 plt.axvline(x=(std / np.sqrt(len(ratioList))), color=
'r')
719 plt.axvline(x=-(std / np.sqrt(len(ratioList))), color=
'r')
720 plt.title(f
"(Source {i} -> Target {j}) mean: {mean:.2g} coeff: {coeff:.2g} valid: {valid}")
723 prompt =
"Press Enter to continue: "
725 ans = input(prompt).lower()
726 if ans
in (
"",
"c",):
728 elif ans
in (
"pdb",
"p",):
736 target=CrosstalkExtractTask,
737 doc=
"Task to measure pixel ratios.",
740 target=CrosstalkSolveTask,
741 doc=
"Task to convert ratio lists to crosstalk coefficients.",
746 """Measure intra-detector crosstalk.
750 lsst.ip.isr.crosstalk.CrosstalkCalib
751 lsst.cp.pipe.measureCrosstalk.CrosstalkExtractTask
752 lsst.cp.pipe.measureCrosstalk.CrosstalkSolveTask
756 The crosstalk this method measures assumes that when a bright
757 pixel is found in one detector amplifier, all other detector
758 amplifiers may see a signal change in the same pixel location
759 (relative to the readout amplifier) as these other pixels are read
760 out at the same time.
762 After processing each input exposure through a limited set of ISR
763 stages, bright unmasked pixels above the threshold are identified.
764 The potential CT signal is found by taking the ratio of the
765 appropriate background-subtracted pixel value on the other
766 amplifiers to the input value on the source amplifier. If the
767 source amplifier has a large number of bright pixels as well, the
768 background level may be elevated, leading to poor ratio
771 The set of ratios found between each pair of amplifiers across all
772 input exposures is then gathered to produce the final CT
773 coefficients. The sigma-clipped mean and sigma are returned from
774 these sets of ratios, with the coefficient to supply to the ISR
775 CrosstalkTask() being the multiplicative inverse of these values.
777 This Task simply calls the pipetask versions of the measure
780 ConfigClass = MeasureCrosstalkConfig
781 _DefaultName =
"measureCrosstalk"
784 RunnerClass = DataRefListRunner
788 self.makeSubtask(
"extract")
789 self.makeSubtask(
"solver")
792 """Run extract task on each of inputs in the dataRef list, then pass
793 that to the solver task.
797 dataRefList : `list` [`lsst.daf.peristence.ButlerDataRef`]
798 Data references for exposures for detectors to process.
802 results : `lsst.pipe.base.Struct`
803 The results struct containing:
805 ``outputCrosstalk`` : `lsst.ip.isr.CrosstalkCalib`
806 Final crosstalk calibration.
807 ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
808 Provenance data for the new calibration.
813 Raised if multiple target detectors are supplied.
815 dataRef = dataRefList[0]
816 camera = dataRef.get(
"camera")
820 for dataRef
in dataRefList:
821 exposure = dataRef.get(
"postISRCCD")
823 if exposure.getDetector().
getName() != activeChip:
824 raise RuntimeError(
"Too many input detectors supplied!")
826 activeChip = exposure.getDetector().
getName()
828 self.extract.debugView(
"extract", exposure)
829 result = self.extract.
run(exposure)
830 ratios.append(result.outputRatios)
832 for detIter, detector
in enumerate(camera):
833 if detector.getName() == activeChip:
835 outputDims = {
'instrument': camera.getName(),
836 'detector': detectorId,
839 finalResults = self.solver.
run(ratios, camera=camera, outputDims=outputDims)
840 dataRef.put(finalResults.outputCrosstalk,
"crosstalk")
A Threshold is used to pass a threshold value to detection algorithms.
def __init__(self, *config=None)
def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, outputDims=None)
def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False)
def filterCrosstalkCalib(inCalib)
def measureCrosstalkCoefficients(self, ratios, rejIter, rejSigma)
def runQuantum(self, butlerQC, inputRefs, outputRefs)
def sigmaClipCorrection(nSigClip)
def __init__(self, **kwargs)
def runDataRef(self, dataRefList)
daf::base::PropertyList * list
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
std::string const & getName() const noexcept
Return a filter's name.
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
def getDebugFrame(debugDisplay, name)