LSST Applications  21.0.0-131-g8cabc107+528f53ee53,22.0.0+00495a2688,22.0.0+0ef2527977,22.0.0+11a2aa21cd,22.0.0+269b7e55e3,22.0.0+2c6b6677a3,22.0.0+64c1bc5aa5,22.0.0+7b3a3f865e,22.0.0+e1b6d2281c,22.0.0+ff3c34362c,22.0.1-1-g1b65d06+c95cbdf3df,22.0.1-1-g7058be7+1cf78af69b,22.0.1-1-g7dab645+2a65e40b06,22.0.1-1-g8760c09+64c1bc5aa5,22.0.1-1-g949febb+64c1bc5aa5,22.0.1-1-ga324b9c+269b7e55e3,22.0.1-1-gf9d8b05+ff3c34362c,22.0.1-10-g781e53d+9b51d1cd24,22.0.1-10-gba590ab+b9624b875d,22.0.1-13-g76f9b8d+2c6b6677a3,22.0.1-14-g22236948+57af756299,22.0.1-18-g3db9cf4b+9b7092c56c,22.0.1-18-gb17765a+2264247a6b,22.0.1-2-g8ef0a89+2c6b6677a3,22.0.1-2-gcb770ba+c99495d3c6,22.0.1-24-g2e899d296+4206820b0d,22.0.1-3-g7aa11f2+2c6b6677a3,22.0.1-3-g8c1d971+f253ffa91f,22.0.1-3-g997b569+ff3b2f8649,22.0.1-4-g1930a60+6871d0c7f6,22.0.1-4-g5b7b756+6b209d634c,22.0.1-6-ga02864e+6871d0c7f6,22.0.1-7-g3402376+a1a2182ac4,22.0.1-7-g65f59fa+54b92689ce,master-gcc5351303a+e1b6d2281c,w.2021.32
LSST Data Management Base Package
measureCrosstalk.py
Go to the documentation of this file.
1 # This file is part of cp_pipe
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <http://www.gnu.org/licenses/>.
21 import itertools
22 import numpy as np
23 
24 from collections import defaultdict
25 
26 import lsst.pipe.base as pipeBase
27 import lsst.pipe.base.connectionTypes as cT
28 
29 from lsstDebug import getDebugFrame
30 from lsst.afw.detection import FootprintSet, Threshold
31 from lsst.afw.display import getDisplay
32 from lsst.pex.config import Config, Field, ListField, ConfigurableField
33 from lsst.ip.isr import CrosstalkCalib, IsrProvenance
34 from lsst.pipe.tasks.getRepositoryData import DataRefListRunner
35 from lsst.cp.pipe.utils import (ddict2dict, sigmaClipCorrection)
36 
37 from ._lookupStaticCalibration import lookupStaticCalibration
38 
39 __all__ = ["CrosstalkExtractConfig", "CrosstalkExtractTask",
40  "CrosstalkSolveTask", "CrosstalkSolveConfig",
41  "MeasureCrosstalkConfig", "MeasureCrosstalkTask"]
42 
43 
44 class CrosstalkExtractConnections(pipeBase.PipelineTaskConnections,
45  dimensions=("instrument", "exposure", "detector")):
46  inputExp = cT.Input(
47  name="crosstalkInputs",
48  doc="Input post-ISR processed exposure to measure crosstalk from.",
49  storageClass="Exposure",
50  dimensions=("instrument", "exposure", "detector"),
51  multiple=False,
52  )
53  # TODO: Depends on DM-21904.
54  sourceExp = cT.Input(
55  name="crosstalkSource",
56  doc="Post-ISR exposure to measure for inter-chip crosstalk onto inputExp.",
57  storageClass="Exposure",
58  dimensions=("instrument", "exposure", "detector"),
59  multiple=True,
60  deferLoad=True,
61  # lookupFunction=None,
62  )
63 
64  outputRatios = cT.Output(
65  name="crosstalkRatios",
66  doc="Extracted crosstalk pixel ratios.",
67  storageClass="StructuredDataDict",
68  dimensions=("instrument", "exposure", "detector"),
69  )
70  outputFluxes = cT.Output(
71  name="crosstalkFluxes",
72  doc="Source pixel fluxes used in ratios.",
73  storageClass="StructuredDataDict",
74  dimensions=("instrument", "exposure", "detector"),
75  )
76 
77  def __init__(self, *, config=None):
78  super().__init__(config=config)
79  # Discard sourceExp until DM-21904 allows full interchip
80  # measurements.
81  self.inputs.discard("sourceExp")
82 
83 
84 class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig,
85  pipelineConnections=CrosstalkExtractConnections):
86  """Configuration for the measurement of pixel ratios.
87  """
88  doMeasureInterchip = Field(
89  dtype=bool,
90  default=False,
91  doc="Measure inter-chip crosstalk as well?",
92  )
93  threshold = Field(
94  dtype=float,
95  default=30000,
96  doc="Minimum level of source pixels for which to measure crosstalk."
97  )
98  ignoreSaturatedPixels = Field(
99  dtype=bool,
100  default=True,
101  doc="Should saturated pixels be ignored?"
102  )
103  badMask = ListField(
104  dtype=str,
105  default=["BAD", "INTRP"],
106  doc="Mask planes to ignore when identifying source pixels."
107  )
108  isTrimmed = Field(
109  dtype=bool,
110  default=True,
111  doc="Is the input exposure trimmed?"
112  )
113 
114  def validate(self):
115  super().validate()
116 
117  # Ensure the handling of the SAT mask plane is consistent
118  # with the ignoreSaturatedPixels value.
119  if self.ignoreSaturatedPixelsignoreSaturatedPixels:
120  if 'SAT' not in self.badMaskbadMask:
121  self.badMaskbadMask.append('SAT')
122  else:
123  if 'SAT' in self.badMaskbadMask:
124  self.badMaskbadMask = [mask for mask in self.badMaskbadMask if mask != 'SAT']
125 
126 
127 class CrosstalkExtractTask(pipeBase.PipelineTask,
128  pipeBase.CmdLineTask):
129  """Task to measure pixel ratios to find crosstalk.
130  """
131  ConfigClass = CrosstalkExtractConfig
132  _DefaultName = 'cpCrosstalkExtract'
133 
134  def run(self, inputExp, sourceExps=[]):
135  """Measure pixel ratios between amplifiers in inputExp.
136 
137  Extract crosstalk ratios between different amplifiers.
138 
139  For pixels above ``config.threshold``, we calculate the ratio
140  between each background-subtracted target amp and the source
141  amp. We return a list of ratios for each pixel for each
142  target/source combination, as nested dictionary containing the
143  ratio.
144 
145  Parameters
146  ----------
147  inputExp : `lsst.afw.image.Exposure`
148  Input exposure to measure pixel ratios on.
149  sourceExp : `list` [`lsst.afw.image.Exposure`], optional
150  List of chips to use as sources to measure inter-chip
151  crosstalk.
152 
153  Returns
154  -------
155  results : `lsst.pipe.base.Struct`
156  The results struct containing:
157 
158  ``outputRatios`` : `dict` [`dict` [`dict` [`dict` [`list`]]]]
159  A catalog of ratio lists. The dictionaries are
160  indexed such that:
161  outputRatios[targetChip][sourceChip][targetAmp][sourceAmp]
162  contains the ratio list for that combination.
163  ``outputFluxes`` : `dict` [`dict` [`list`]]
164  A catalog of flux lists. The dictionaries are
165  indexed such that:
166  outputFluxes[sourceChip][sourceAmp]
167  contains the flux list used in the outputRatios.
168 
169  Notes
170  -----
171  The lsstDebug.Info() method can be rewritten for __name__ =
172  `lsst.cp.pipe.measureCrosstalk`, and supports the parameters:
173 
174  debug.display['extract'] : `bool`
175  Display the exposure under consideration, with the pixels used
176  for crosstalk measurement indicated by the DETECTED mask plane.
177  debug.display['pixels'] : `bool`
178  Display a plot of the ratio calculated for each pixel used in this
179  exposure, split by amplifier pairs. The median value is listed
180  for reference.
181  """
182  outputRatios = defaultdict(lambda: defaultdict(dict))
183  outputFluxes = defaultdict(lambda: defaultdict(dict))
184 
185  threshold = self.config.threshold
186  badPixels = list(self.config.badMask)
187 
188  targetDetector = inputExp.getDetector()
189  targetChip = targetDetector.getName()
190 
191  # Always look at the target chip first, then go to any other supplied exposures.
192  sourceExtractExps = [inputExp]
193  sourceExtractExps.extend(sourceExps)
194 
195  self.log.info("Measuring full detector background for target: %s", targetChip)
196  targetIm = inputExp.getMaskedImage()
197  FootprintSet(targetIm, Threshold(threshold), "DETECTED")
198  detected = targetIm.getMask().getPlaneBitMask("DETECTED")
199  bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"])
200 
201  self.debugViewdebugView('extract', inputExp)
202 
203  for sourceExp in sourceExtractExps:
204  sourceDetector = sourceExp.getDetector()
205  sourceChip = sourceDetector.getName()
206  sourceIm = sourceExp.getMaskedImage()
207  bad = sourceIm.getMask().getPlaneBitMask(badPixels)
208  self.log.info("Measuring crosstalk from source: %s", sourceChip)
209 
210  if sourceExp != inputExp:
211  FootprintSet(sourceIm, Threshold(threshold), "DETECTED")
212  detected = sourceIm.getMask().getPlaneBitMask("DETECTED")
213 
214  # The dictionary of amp-to-amp ratios for this pair of source->target detectors.
215  ratioDict = defaultdict(lambda: defaultdict(list))
216  extractedCount = 0
217 
218  for sourceAmp in sourceDetector:
219  sourceAmpName = sourceAmp.getName()
220  sourceAmpBBox = sourceAmp.getBBox() if self.config.isTrimmed else sourceAmp.getRawDataBBox()
221  sourceAmpImage = sourceIm[sourceAmpBBox]
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)
228 
229  outputFluxes[sourceChip][sourceAmpName] = sourceAmpImage.image.array[select].tolist()
230 
231  for targetAmp in targetDetector:
232  # iterate over targetExposure
233  targetAmpName = targetAmp.getName()
234  if sourceAmpName == targetAmpName and sourceChip == targetChip:
235  ratioDict[sourceAmpName][targetAmpName] = []
236  continue
237  self.log.debug(" Target amplifier: %s", targetAmpName)
238 
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
245 
246  self.debugPixelsdebugPixels('pixels',
247  sourceAmpImage.image.array[select],
248  targetAmpImage.array[select] - bg,
249  sourceAmpName, targetAmpName)
250 
251  self.log.info("Extracted %d pixels from %s -> %s (targetBG: %f)",
252  extractedCount, sourceChip, targetChip, bg)
253  outputRatios[targetChip][sourceChip] = ratioDict
254 
255  return pipeBase.Struct(
256  outputRatios=ddict2dict(outputRatios),
257  outputFluxes=ddict2dict(outputFluxes)
258  )
259 
260  def debugView(self, stepname, exposure):
261  """Utility function to examine the image being processed.
262 
263  Parameters
264  ----------
265  stepname : `str`
266  State of processing to view.
267  exposure : `lsst.afw.image.Exposure`
268  Exposure to view.
269  """
270  frame = getDebugFrame(self._display, stepname)
271  if frame:
272  display = getDisplay(frame)
273  display.scale('asinh', 'zscale')
274  display.mtv(exposure)
275 
276  prompt = "Press Enter to continue: "
277  while True:
278  ans = input(prompt).lower()
279  if ans in ("", "c",):
280  break
281 
282  def debugPixels(self, stepname, pixelsIn, pixelsOut, sourceName, targetName):
283  """Utility function to examine the CT ratio pixel values.
284 
285  Parameters
286  ----------
287  stepname : `str`
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.
293  sourceName : `str`
294  Source amplifier name
295  targetName : `str`
296  Target amplifier name
297  """
298  frame = getDebugFrame(self._display, stepname)
299  if frame:
300  import matplotlib.pyplot as plt
301  figure = plt.figure(1)
302  figure.clear()
303 
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))}")
310  figure.show()
311 
312  prompt = "Press Enter to continue: "
313  while True:
314  ans = input(prompt).lower()
315  if ans in ("", "c",):
316  break
317  plt.close()
318 
319 
320 class CrosstalkSolveConnections(pipeBase.PipelineTaskConnections,
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"),
327  multiple=True,
328  )
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"),
334  multiple=True,
335  )
336  camera = cT.PrerequisiteInput(
337  name="camera",
338  doc="Camera the input data comes from.",
339  storageClass="Camera",
340  dimensions=("instrument",),
341  isCalibration=True,
342  lookupFunction=lookupStaticCalibration,
343  )
344 
345  outputCrosstalk = cT.Output(
346  name="crosstalk",
347  doc="Output proposed crosstalk calibration.",
348  storageClass="CrosstalkCalib",
349  dimensions=("instrument", "detector"),
350  multiple=False,
351  isCalibration=True,
352  )
353 
354  def __init__(self, *, config=None):
355  super().__init__(config=config)
356 
357  if config.fluxOrder == 0:
358  self.inputs.discard("inputFluxes")
359 
360 
361 class CrosstalkSolveConfig(pipeBase.PipelineTaskConfig,
362  pipelineConnections=CrosstalkSolveConnections):
363  """Configuration for the solving of crosstalk from pixel ratios.
364  """
365  rejIter = Field(
366  dtype=int,
367  default=3,
368  doc="Number of rejection iterations for final coefficient calculation.",
369  )
370  rejSigma = Field(
371  dtype=float,
372  default=2.0,
373  doc="Rejection threshold (sigma) for final coefficient calculation.",
374  )
375  fluxOrder = Field(
376  dtype=int,
377  default=0,
378  doc="Polynomial order in source flux to fit crosstalk.",
379  )
380  doFiltering = Field(
381  dtype=bool,
382  default=False,
383  doc="Filter generated crosstalk to remove marginal measurements.",
384  )
385 
386 
387 class CrosstalkSolveTask(pipeBase.PipelineTask,
388  pipeBase.CmdLineTask):
389  """Task to solve crosstalk from pixel ratios.
390  """
391  ConfigClass = CrosstalkSolveConfig
392  _DefaultName = 'cpCrosstalkSolve'
393 
394  def runQuantum(self, butlerQC, inputRefs, outputRefs):
395  """Ensure that the input and output dimensions are passed along.
396 
397  Parameters
398  ----------
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.
405  """
406  inputs = butlerQC.get(inputRefs)
407 
408  # Use the dimensions to set calib/provenance information.
409  inputs['inputDims'] = [exp.dataId.byName() for exp in inputRefs.inputRatios]
410  inputs['outputDims'] = outputRefs.outputCrosstalk.dataId.byName()
411 
412  outputs = self.runrun(**inputs)
413  butlerQC.put(outputs, outputRefs)
414 
415  def run(self, inputRatios, inputFluxes=None, camera=None, inputDims=None, outputDims=None):
416  """Combine ratios to produce crosstalk coefficients.
417 
418  Parameters
419  ----------
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`
427  Input 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.
432 
433  Returns
434  -------
435  results : `lsst.pipe.base.Struct`
436  The results struct containing:
437 
438  ``outputCrosstalk`` : `lsst.ip.isr.CrosstalkCalib`
439  Final crosstalk calibration.
440  ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
441  Provenance data for the new calibration.
442 
443  Raises
444  ------
445  RuntimeError
446  Raised if the input data contains multiple target detectors.
447 
448  Notes
449  -----
450  The lsstDebug.Info() method can be rewritten for __name__ =
451  `lsst.ip.isr.measureCrosstalk`, and supports the parameters:
452 
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
456  exposures/detectors.
457 
458  """
459  if outputDims:
460  calibChip = outputDims['detector']
461  instrument = outputDims['instrument']
462  else:
463  # calibChip needs to be set manually in Gen2.
464  calibChip = None
465  instrument = None
466 
467  if camera and calibChip:
468  calibDetector = camera[calibChip]
469  else:
470  calibDetector = None
471 
472  self.log.info("Combining measurements from %d ratios and %d fluxes",
473  len(inputRatios), len(inputFluxes) if inputFluxes else 0)
474 
475  if inputFluxes is None:
476  inputFluxes = [None for exp in inputRatios]
477 
478  combinedRatios = defaultdict(lambda: defaultdict(list))
479  combinedFluxes = defaultdict(lambda: defaultdict(list))
480  for ratioDict, fluxDict in zip(inputRatios, inputFluxes):
481  for targetChip in ratioDict:
482  if calibChip and targetChip != calibChip and targetChip != calibDetector.getName():
483  raise RuntimeError(f"Target chip: {targetChip} does not match calibration dimension: "
484  f"{calibChip}, {calibDetector.getName()}!")
485 
486  sourceChip = targetChip
487  if sourceChip in ratioDict[targetChip]:
488  ratios = ratioDict[targetChip][sourceChip]
489 
490  for targetAmp in ratios:
491  for sourceAmp in ratios[targetAmp]:
492  combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp])
493  if fluxDict:
494  combinedFluxes[targetAmp][sourceAmp].extend(fluxDict[sourceChip][sourceAmp])
495  # TODO: DM-21904
496  # Iterating over all other entries in ratioDict[targetChip] will yield
497  # inter-chip terms.
498 
499  for targetAmp in combinedRatios:
500  for sourceAmp in combinedRatios[targetAmp]:
501  self.log.info("Read %d pixels for %s -> %s",
502  len(combinedRatios[targetAmp][sourceAmp]),
503  targetAmp, sourceAmp)
504  if len(combinedRatios[targetAmp][sourceAmp]) > 1:
505  self.debugRatiosdebugRatios('reduce', combinedRatios, targetAmp, sourceAmp)
506 
507  if self.config.fluxOrder == 0:
508  self.log.info("Fitting crosstalk coefficients.")
509  calib = self.measureCrosstalkCoefficientsmeasureCrosstalkCoefficients(combinedRatios,
510  self.config.rejIter, self.config.rejSigma)
511  else:
512  raise NotImplementedError("Non-linear crosstalk terms are not yet supported.")
513 
514  self.log.info("Number of valid coefficients: %d", np.sum(calib.coeffValid))
515 
516  if self.config.doFiltering:
517  # This step will apply the calculated validity values to
518  # censor poorly measured coefficients.
519  self.log.info("Filtering measured crosstalk to remove invalid solutions.")
520  calib = self.filterCrosstalkCalibfilterCrosstalkCalib(calib)
521 
522  # Populate the remainder of the calibration information.
523  calib.hasCrosstalk = True
524  calib.interChip = {}
525 
526  # calibChip is the detector dimension, which is the detector Id
527  calib._detectorId = calibChip
528  if calibDetector:
529  calib._detectorName = calibDetector.getName()
530  calib._detectorSerial = calibDetector.getSerial()
531 
532  calib._instrument = instrument
533  calib.updateMetadata(setCalibId=True, setDate=True)
534 
535  # Make an IsrProvenance().
536  provenance = IsrProvenance(calibType="CROSSTALK")
537  provenance._detectorName = calibChip
538  if inputDims:
539  provenance.fromDataIds(inputDims)
540  provenance._instrument = instrument
541  provenance.updateMetadata()
542 
543  return pipeBase.Struct(
544  outputCrosstalk=calib,
545  outputProvenance=provenance,
546  )
547 
548  def measureCrosstalkCoefficients(self, ratios, rejIter, rejSigma):
549  """Measure crosstalk coefficients from the ratios.
550 
551  Given a list of ratios for each target/source amp combination,
552  we measure a sigma clipped mean and error.
553 
554  The coefficient errors returned are the standard deviation of
555  the final set of clipped input ratios.
556 
557  Parameters
558  ----------
559  ratios : `dict` of `dict` of `numpy.ndarray`
560  Catalog of arrays of ratios.
561  rejIter : `int`
562  Number of rejection iterations.
563  rejSigma : `float`
564  Rejection threshold (sigma).
565 
566  Returns
567  -------
568  calib : `lsst.ip.isr.CrosstalkCalib`
569  The output crosstalk calibration.
570 
571  Notes
572  -----
573  The lsstDebug.Info() method can be rewritten for __name__ =
574  `lsst.ip.isr.measureCrosstalk`, and supports the parameters:
575 
576  debug.display['measure'] : `bool`
577  Display the CDF of the combined ratio measurements for
578  a pair of source/target amplifiers from the final set of
579  clipped input ratios.
580  """
581  calib = CrosstalkCalib(nAmp=len(ratios))
582 
583  # Calibration stores coefficients as a numpy ndarray.
584  ordering = list(ratios.keys())
585  for ii, jj in itertools.product(range(calib.nAmp), range(calib.nAmp)):
586  if ii == jj:
587  values = [0.0]
588  else:
589  values = np.array(ratios[ordering[ii]][ordering[jj]])
590  values = values[np.abs(values) < 1.0] # Discard unreasonable values
591 
592  calib.coeffNum[ii][jj] = len(values)
593 
594  if len(values) == 0:
595  self.log.warn("No values for matrix element %d,%d" % (ii, jj))
596  calib.coeffs[ii][jj] = np.nan
597  calib.coeffErr[ii][jj] = np.nan
598  calib.coeffValid[ii][jj] = False
599  else:
600  if ii != jj:
601  for rej in range(rejIter):
602  lo, med, hi = np.percentile(values, [25.0, 50.0, 75.0])
603  sigma = 0.741*(hi - lo)
604  good = np.abs(values - med) < rejSigma*sigma
605  if good.sum() == len(good):
606  break
607  values = values[good]
608 
609  calib.coeffs[ii][jj] = np.mean(values)
610  if calib.coeffNum[ii][jj] == 1:
611  calib.coeffErr[ii][jj] = np.nan
612  else:
613  correctionFactor = sigmaClipCorrection(rejSigma)
614  calib.coeffErr[ii][jj] = np.std(values) * correctionFactor
615  calib.coeffValid[ii][jj] = (np.abs(calib.coeffs[ii][jj])
616  > calib.coeffErr[ii][jj] / np.sqrt(calib.coeffNum[ii][jj]))
617 
618  if calib.coeffNum[ii][jj] > 1:
619  self.debugRatiosdebugRatios('measure', ratios, ordering[ii], ordering[jj],
620  calib.coeffs[ii][jj], calib.coeffValid[ii][jj])
621 
622  return calib
623 
624  @staticmethod
625  def filterCrosstalkCalib(inCalib):
626  """Apply valid constraints to the measured values.
627 
628  Any measured coefficient that is determined to be invalid is
629  set to zero, and has the error set to nan. The validation is
630  determined by checking that the measured coefficient is larger
631  than the calculated standard error of the mean.
632 
633  Parameters
634  ----------
635  inCalib : `lsst.ip.isr.CrosstalkCalib`
636  Input calibration to filter.
637 
638  Returns
639  -------
640  outCalib : `lsst.ip.isr.CrosstalkCalib`
641  Filtered calibration.
642  """
643  outCalib = CrosstalkCalib()
644  outCalib.numAmps = inCalib.numAmps
645 
646  outCalib.coeffs = inCalib.coeffs
647  outCalib.coeffs[~inCalib.coeffValid] = 0.0
648 
649  outCalib.coeffErr = inCalib.coeffErr
650  outCalib.coeffErr[~inCalib.coeffValid] = np.nan
651 
652  outCalib.coeffNum = inCalib.coeffNum
653  outCalib.coeffValid = inCalib.coeffValid
654 
655  return outCalib
656 
657  def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False):
658  """Utility function to examine the final CT ratio set.
659 
660  Parameters
661  ----------
662  stepname : `str`
663  State of processing to view.
664  ratios : `dict` of `dict` of `np.ndarray`
665  Array of measured CT ratios, indexed by source/victim
666  amplifier.
667  i : `str`
668  Index of the source amplifier.
669  j : `str`
670  Index of the target amplifier.
671  coeff : `float`, optional
672  Coefficient calculated to plot along with the simple mean.
673  valid : `bool`, optional
674  Validity to be added to the plot title.
675  """
676  frame = getDebugFrame(self._display, stepname)
677  if frame:
678  if i == j or ratios is None or len(ratios) < 1:
679  pass
680 
681  ratioList = ratios[i][j]
682  if ratioList is None or len(ratioList) < 1:
683  pass
684 
685  mean = np.mean(ratioList)
686  std = np.std(ratioList)
687  import matplotlib.pyplot as plt
688  figure = plt.figure(1)
689  figure.clear()
690  plt.hist(x=ratioList, bins=len(ratioList),
691  cumulative=True, color='b', density=True, histtype='step')
692  plt.xlabel("Measured pixel ratio")
693  plt.ylabel(f"CDF: n={len(ratioList)}")
694  plt.xlim(np.percentile(ratioList, [1.0, 99]))
695  plt.axvline(x=mean, color="k")
696  plt.axvline(x=coeff, color='g')
697  plt.axvline(x=(std / np.sqrt(len(ratioList))), color='r')
698  plt.axvline(x=-(std / np.sqrt(len(ratioList))), color='r')
699  plt.title(f"(Source {i} -> Target {j}) mean: {mean:.2g} coeff: {coeff:.2g} valid: {valid}")
700  figure.show()
701 
702  prompt = "Press Enter to continue: "
703  while True:
704  ans = input(prompt).lower()
705  if ans in ("", "c",):
706  break
707  elif ans in ("pdb", "p",):
708  import pdb
709  pdb.set_trace()
710  plt.close()
711 
712 
714  extract = ConfigurableField(
715  target=CrosstalkExtractTask,
716  doc="Task to measure pixel ratios.",
717  )
719  target=CrosstalkSolveTask,
720  doc="Task to convert ratio lists to crosstalk coefficients.",
721  )
722 
723 
724 class MeasureCrosstalkTask(pipeBase.CmdLineTask):
725  """Measure intra-detector crosstalk.
726 
727  See also
728  --------
729  lsst.ip.isr.crosstalk.CrosstalkCalib
730  lsst.cp.pipe.measureCrosstalk.CrosstalkExtractTask
731  lsst.cp.pipe.measureCrosstalk.CrosstalkSolveTask
732 
733  Notes
734  -----
735  The crosstalk this method measures assumes that when a bright
736  pixel is found in one detector amplifier, all other detector
737  amplifiers may see a signal change in the same pixel location
738  (relative to the readout amplifier) as these other pixels are read
739  out at the same time.
740 
741  After processing each input exposure through a limited set of ISR
742  stages, bright unmasked pixels above the threshold are identified.
743  The potential CT signal is found by taking the ratio of the
744  appropriate background-subtracted pixel value on the other
745  amplifiers to the input value on the source amplifier. If the
746  source amplifier has a large number of bright pixels as well, the
747  background level may be elevated, leading to poor ratio
748  measurements.
749 
750  The set of ratios found between each pair of amplifiers across all
751  input exposures is then gathered to produce the final CT
752  coefficients. The sigma-clipped mean and sigma are returned from
753  these sets of ratios, with the coefficient to supply to the ISR
754  CrosstalkTask() being the multiplicative inverse of these values.
755 
756  This Task simply calls the pipetask versions of the measure
757  crosstalk code.
758  """
759  ConfigClass = MeasureCrosstalkConfig
760  _DefaultName = "measureCrosstalk"
761 
762  # Let's use this instead of messing with parseAndRun.
763  RunnerClass = DataRefListRunner
764 
765  def __init__(self, **kwargs):
766  super().__init__(**kwargs)
767  self.makeSubtask("extract")
768  self.makeSubtask("solver")
769 
770  def runDataRef(self, dataRefList):
771  """Run extract task on each of inputs in the dataRef list, then pass
772  that to the solver task.
773 
774  Parameters
775  ----------
776  dataRefList : `list` [`lsst.daf.peristence.ButlerDataRef`]
777  Data references for exposures for detectors to process.
778 
779  Returns
780  -------
781  results : `lsst.pipe.base.Struct`
782  The results struct containing:
783 
784  ``outputCrosstalk`` : `lsst.ip.isr.CrosstalkCalib`
785  Final crosstalk calibration.
786  ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
787  Provenance data for the new calibration.
788 
789  Raises
790  ------
791  RuntimeError
792  Raised if multiple target detectors are supplied.
793  """
794  dataRef = dataRefList[0]
795  camera = dataRef.get("camera")
796 
797  ratios = []
798  activeChip = None
799  for dataRef in dataRefList:
800  exposure = dataRef.get("postISRCCD")
801  if activeChip:
802  if exposure.getDetector().getName() != activeChip:
803  raise RuntimeError("Too many input detectors supplied!")
804  else:
805  activeChip = exposure.getDetector().getName()
806 
807  self.extract.debugView("extract", exposure)
808  result = self.extract.run(exposure)
809  ratios.append(result.outputRatios)
810 
811  for detIter, detector in enumerate(camera):
812  if detector.getName() == activeChip:
813  detectorId = detIter
814  outputDims = {'instrument': camera.getName(),
815  'detector': detectorId,
816  }
817 
818  finalResults = self.solver.run(ratios, camera=camera, outputDims=outputDims)
819  dataRef.put(finalResults.outputCrosstalk, "crosstalk")
820 
821  return finalResults
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
def debugPixels(self, stepname, pixelsIn, pixelsOut, sourceName, targetName)
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 measureCrosstalkCoefficients(self, ratios, rejIter, rejSigma)
def runQuantum(self, butlerQC, inputRefs, outputRefs)
daf::base::PropertyList * list
Definition: fits.cc:913
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
std::string const & getName() const noexcept
Return a filter's name.
Definition: Filter.h:78
def sigmaClipCorrection(nSigClip)
Definition: utils.py:42
def ddict2dict(d)
Definition: utils.py:806
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95