LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 from scipy.stats import norm
24 
25 from collections import defaultdict
26 
27 import lsst.pipe.base as pipeBase
29 
30 from lsstDebug import getDebugFrame
31 from lsst.afw.detection import FootprintSet, Threshold
32 from lsst.afw.display import getDisplay
33 from lsst.pex.config import Config, Field, ListField, ConfigurableField
34 from lsst.ip.isr import CrosstalkCalib, IsrProvenance
35 from lsst.pipe.tasks.getRepositoryData import DataRefListRunner
36 from lsst.cp.pipe.utils import ddict2dict
37 
38 from ._lookupStaticCalibration import lookupStaticCalibration
39 
40 __all__ = ["CrosstalkExtractConfig", "CrosstalkExtractTask",
41  "CrosstalkSolveTask", "CrosstalkSolveConfig",
42  "MeasureCrosstalkConfig", "MeasureCrosstalkTask"]
43 
44 
45 class CrosstalkExtractConnections(pipeBase.PipelineTaskConnections,
46  dimensions=("instrument", "exposure", "detector")):
47  inputExp = cT.Input(
48  name="crosstalkInputs",
49  doc="Input post-ISR processed exposure to measure crosstalk from.",
50  storageClass="Exposure",
51  dimensions=("instrument", "exposure", "detector"),
52  multiple=False,
53  )
54  # TODO: Depends on DM-21904.
55  sourceExp = cT.Input(
56  name="crosstalkSource",
57  doc="Post-ISR exposure to measure for inter-chip crosstalk onto inputExp.",
58  storageClass="Exposure",
59  dimensions=("instrument", "exposure", "detector"),
60  multiple=True,
61  deferLoad=True,
62  # lookupFunction=None,
63  )
64 
65  outputRatios = cT.Output(
66  name="crosstalkRatios",
67  doc="Extracted crosstalk pixel ratios.",
68  storageClass="StructuredDataDict",
69  dimensions=("instrument", "exposure", "detector"),
70  )
71  outputFluxes = cT.Output(
72  name="crosstalkFluxes",
73  doc="Source pixel fluxes used in ratios.",
74  storageClass="StructuredDataDict",
75  dimensions=("instrument", "exposure", "detector"),
76  )
77 
78  def __init__(self, *, config=None):
79  super().__init__(config=config)
80  # Discard sourceExp until DM-21904 allows full interchip
81  # measurements.
82  self.inputs.discard("sourceExp")
83 
84 
85 class CrosstalkExtractConfig(pipeBase.PipelineTaskConfig,
86  pipelineConnections=CrosstalkExtractConnections):
87  """Configuration for the measurement of pixel ratios.
88  """
89  doMeasureInterchip = Field(
90  dtype=bool,
91  default=False,
92  doc="Measure inter-chip crosstalk as well?",
93  )
94  threshold = Field(
95  dtype=float,
96  default=30000,
97  doc="Minimum level of source pixels for which to measure crosstalk."
98  )
99  ignoreSaturatedPixels = Field(
100  dtype=bool,
101  default=True,
102  doc="Should saturated pixels be ignored?"
103  )
104  badMask = ListField(
105  dtype=str,
106  default=["BAD", "INTRP"],
107  doc="Mask planes to ignore when identifying source pixels."
108  )
109  isTrimmed = Field(
110  dtype=bool,
111  default=True,
112  doc="Is the input exposure trimmed?"
113  )
114 
115  def validate(self):
116  super().validate()
117 
118  # Ensure the handling of the SAT mask plane is consistent
119  # with the ignoreSaturatedPixels value.
120  if self.ignoreSaturatedPixelsignoreSaturatedPixels:
121  if 'SAT' not in self.badMaskbadMask:
122  self.badMaskbadMask.append('SAT')
123  else:
124  if 'SAT' in self.badMaskbadMask:
125  self.badMaskbadMask = [mask for mask in self.badMaskbadMask if mask != 'SAT']
126 
127 
128 class CrosstalkExtractTask(pipeBase.PipelineTask,
129  pipeBase.CmdLineTask):
130  """Task to measure pixel ratios to find crosstalk.
131  """
132  ConfigClass = CrosstalkExtractConfig
133  _DefaultName = 'cpCrosstalkExtract'
134 
135  def run(self, inputExp, sourceExps=[]):
136  """Measure pixel ratios between amplifiers in inputExp.
137 
138  Extract crosstalk ratios between different amplifiers.
139 
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
144  ratio.
145 
146  Parameters
147  ----------
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
152  crosstalk.
153 
154  Returns
155  -------
156  results : `lsst.pipe.base.Struct`
157  The results struct containing:
158 
159  ``outputRatios`` : `dict` [`dict` [`dict` [`dict` [`list`]]]]
160  A catalog of ratio lists. The dictionaries are
161  indexed such that:
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
166  indexed such that:
167  outputFluxes[sourceChip][sourceAmp]
168  contains the flux list used in the outputRatios.
169 
170  Notes
171  -----
172  The lsstDebug.Info() method can be rewritten for __name__ =
173  `lsst.cp.pipe.measureCrosstalk`, and supports the parameters:
174 
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
181  for reference.
182  """
183  outputRatios = defaultdict(lambda: defaultdict(dict))
184  outputFluxes = defaultdict(lambda: defaultdict(dict))
185 
186  threshold = self.config.threshold
187  badPixels = list(self.config.badMask)
188 
189  targetDetector = inputExp.getDetector()
190  targetChip = targetDetector.getName()
191 
192  # Always look at the target chip first, then go to any other supplied exposures.
193  sourceExtractExps = [inputExp]
194  sourceExtractExps.extend(sourceExps)
195 
196  self.log.info("Measuring full detector background for target: %s", targetChip)
197  targetIm = inputExp.getMaskedImage()
198  FootprintSet(targetIm, Threshold(threshold), "DETECTED")
199  detected = targetIm.getMask().getPlaneBitMask("DETECTED")
200  bg = CrosstalkCalib.calculateBackground(targetIm, badPixels + ["DETECTED"])
201 
202  self.debugViewdebugView('extract', inputExp)
203 
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)
210 
211  if sourceExp != inputExp:
212  FootprintSet(sourceIm, Threshold(threshold), "DETECTED")
213  detected = sourceIm.getMask().getPlaneBitMask("DETECTED")
214 
215  # The dictionary of amp-to-amp ratios for this pair of source->target detectors.
216  ratioDict = defaultdict(lambda: defaultdict(list))
217  extractedCount = 0
218 
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)
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  self.log.info("Combining measurements from %d ratios and %d fluxes",
468  len(inputRatios), len(inputFluxes) if inputFluxes else 0)
469 
470  if inputFluxes is None:
471  inputFluxes = [None for exp in inputRatios]
472 
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!")
479 
480  sourceChip = targetChip
481  if sourceChip in ratioDict[targetChip]:
482  ratios = ratioDict[targetChip][sourceChip]
483 
484  for targetAmp in ratios:
485  for sourceAmp in ratios[targetAmp]:
486  combinedRatios[targetAmp][sourceAmp].extend(ratios[targetAmp][sourceAmp])
487  if fluxDict:
488  combinedFluxes[targetAmp][sourceAmp].extend(fluxDict[sourceChip][sourceAmp])
489  # TODO: DM-21904
490  # Iterating over all other entries in ratioDict[targetChip] will yield
491  # inter-chip terms.
492 
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)
500 
501  if self.config.fluxOrder == 0:
502  self.log.info("Fitting crosstalk coefficients.")
503  calib = self.measureCrosstalkCoefficientsmeasureCrosstalkCoefficients(combinedRatios,
504  self.config.rejIter, self.config.rejSigma)
505  else:
506  raise NotImplementedError("Non-linear crosstalk terms are not yet supported.")
507 
508  self.log.info("Number of valid coefficients: %d", np.sum(calib.coeffValid))
509 
510  if self.config.doFiltering:
511  # This step will apply the calculated validity values to
512  # censor poorly measured coefficients.
513  self.log.info("Filtering measured crosstalk to remove invalid solutions.")
514  calib = self.filterCrosstalkCalibfilterCrosstalkCalib(calib)
515 
516  # Populate the remainder of the calibration information.
517  calib.hasCrosstalk = True
518  calib.interChip = {}
519 
520  # calibChip is the detector dimension, which is the detector Id
521  calib._detectorId = calibChip
522  if camera:
523  calib._detectorName = camera[calibChip].getName()
524  calib._detectorSerial = camera[calibChip].getSerial()
525 
526  calib._instrument = instrument
527  calib.updateMetadata(setCalibId=True, setDate=True)
528 
529  # Make an IsrProvenance().
530  provenance = IsrProvenance(calibType="CROSSTALK")
531  provenance._detectorName = calibChip
532  if inputDims:
533  provenance.fromDataIds(inputDims)
534  provenance._instrument = instrument
535  provenance.updateMetadata()
536 
537  return pipeBase.Struct(
538  outputCrosstalk=calib,
539  outputProvenance=provenance,
540  )
541 
542  def measureCrosstalkCoefficients(self, ratios, rejIter, rejSigma):
543  """Measure crosstalk coefficients from the ratios.
544 
545  Given a list of ratios for each target/source amp combination,
546  we measure a sigma clipped mean and error.
547 
548  The coefficient errors returned are the standard deviation of
549  the final set of clipped input ratios.
550 
551  Parameters
552  ----------
553  ratios : `dict` of `dict` of `numpy.ndarray`
554  Catalog of arrays of ratios.
555  rejIter : `int`
556  Number of rejection iterations.
557  rejSigma : `float`
558  Rejection threshold (sigma).
559 
560  Returns
561  -------
562  calib : `lsst.ip.isr.CrosstalkCalib`
563  The output crosstalk calibration.
564 
565  Notes
566  -----
567  The lsstDebug.Info() method can be rewritten for __name__ =
568  `lsst.ip.isr.measureCrosstalk`, and supports the parameters:
569 
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.
574  """
575  calib = CrosstalkCalib(nAmp=len(ratios))
576 
577  # Calibration stores coefficients as a numpy ndarray.
578  ordering = list(ratios.keys())
579  for ii, jj in itertools.product(range(calib.nAmp), range(calib.nAmp)):
580  if ii == jj:
581  values = [0.0]
582  else:
583  values = np.array(ratios[ordering[ii]][ordering[jj]])
584  values = values[np.abs(values) < 1.0] # Discard unreasonable values
585 
586  calib.coeffNum[ii][jj] = len(values)
587 
588  if len(values) == 0:
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
593  else:
594  if ii != jj:
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):
600  break
601  values = values[good]
602 
603  calib.coeffs[ii][jj] = np.mean(values)
604  if calib.coeffNum[ii][jj] == 1:
605  calib.coeffErr[ii][jj] = np.nan
606  else:
607  correctionFactor = self.sigmaClipCorrectionsigmaClipCorrection(rejSigma)
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]))
611 
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])
615 
616  return calib
617 
618  @staticmethod
619  def sigmaClipCorrection(nSigClip):
620  """Correct measured sigma to account for clipping.
621 
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).
630 
631  Parameters
632  ----------
633  nSigClip : `float`
634  Number of sigma the measurement was clipped by.
635 
636  Returns
637  -------
638  scaleFactor : `float`
639  Scale factor to increase the measured sigma by.
640 
641  """
642  varFactor = 1.0 + (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))
643  return 1.0 / np.sqrt(varFactor)
644 
645  @staticmethod
646  def filterCrosstalkCalib(inCalib):
647  """Apply valid constraints to the measured values.
648 
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.
653 
654  Parameters
655  ----------
656  inCalib : `lsst.ip.isr.CrosstalkCalib`
657  Input calibration to filter.
658 
659  Returns
660  -------
661  outCalib : `lsst.ip.isr.CrosstalkCalib`
662  Filtered calibration.
663  """
664  outCalib = CrosstalkCalib()
665  outCalib.numAmps = inCalib.numAmps
666 
667  outCalib.coeffs = inCalib.coeffs
668  outCalib.coeffs[~inCalib.coeffValid] = 0.0
669 
670  outCalib.coeffErr = inCalib.coeffErr
671  outCalib.coeffErr[~inCalib.coeffValid] = np.nan
672 
673  outCalib.coeffNum = inCalib.coeffNum
674  outCalib.coeffValid = inCalib.coeffValid
675 
676  return outCalib
677 
678  def debugRatios(self, stepname, ratios, i, j, coeff=0.0, valid=False):
679  """Utility function to examine the final CT ratio set.
680 
681  Parameters
682  ----------
683  stepname : `str`
684  State of processing to view.
685  ratios : `dict` of `dict` of `np.ndarray`
686  Array of measured CT ratios, indexed by source/victim
687  amplifier.
688  i : `str`
689  Index of the source amplifier.
690  j : `str`
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.
696  """
697  frame = getDebugFrame(self._display, stepname)
698  if frame:
699  if i == j or ratios is None or len(ratios) < 1:
700  pass
701 
702  ratioList = ratios[i][j]
703  if ratioList is None or len(ratioList) < 1:
704  pass
705 
706  mean = np.mean(ratioList)
707  std = np.std(ratioList)
708  import matplotlib.pyplot as plt
709  figure = plt.figure(1)
710  figure.clear()
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}")
721  figure.show()
722 
723  prompt = "Press Enter to continue: "
724  while True:
725  ans = input(prompt).lower()
726  if ans in ("", "c",):
727  break
728  elif ans in ("pdb", "p",):
729  import pdb
730  pdb.set_trace()
731  plt.close()
732 
733 
735  extract = ConfigurableField(
736  target=CrosstalkExtractTask,
737  doc="Task to measure pixel ratios.",
738  )
740  target=CrosstalkSolveTask,
741  doc="Task to convert ratio lists to crosstalk coefficients.",
742  )
743 
744 
745 class MeasureCrosstalkTask(pipeBase.CmdLineTask):
746  """Measure intra-detector crosstalk.
747 
748  See also
749  --------
750  lsst.ip.isr.crosstalk.CrosstalkCalib
751  lsst.cp.pipe.measureCrosstalk.CrosstalkExtractTask
752  lsst.cp.pipe.measureCrosstalk.CrosstalkSolveTask
753 
754  Notes
755  -----
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.
761 
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
769  measurements.
770 
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.
776 
777  This Task simply calls the pipetask versions of the measure
778  crosstalk code.
779  """
780  ConfigClass = MeasureCrosstalkConfig
781  _DefaultName = "measureCrosstalk"
782 
783  # Let's use this instead of messing with parseAndRun.
784  RunnerClass = DataRefListRunner
785 
786  def __init__(self, **kwargs):
787  super().__init__(**kwargs)
788  self.makeSubtask("extract")
789  self.makeSubtask("solver")
790 
791  def runDataRef(self, dataRefList):
792  """Run extract task on each of inputs in the dataRef list, then pass
793  that to the solver task.
794 
795  Parameters
796  ----------
797  dataRefList : `list` [`lsst.daf.peristence.ButlerDataRef`]
798  Data references for exposures for detectors to process.
799 
800  Returns
801  -------
802  results : `lsst.pipe.base.Struct`
803  The results struct containing:
804 
805  ``outputCrosstalk`` : `lsst.ip.isr.CrosstalkCalib`
806  Final crosstalk calibration.
807  ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
808  Provenance data for the new calibration.
809 
810  Raises
811  ------
812  RuntimeError
813  Raised if multiple target detectors are supplied.
814  """
815  dataRef = dataRefList[0]
816  camera = dataRef.get("camera")
817 
818  ratios = []
819  activeChip = None
820  for dataRef in dataRefList:
821  exposure = dataRef.get("postISRCCD")
822  if activeChip:
823  if exposure.getDetector().getName() != activeChip:
824  raise RuntimeError("Too many input detectors supplied!")
825  else:
826  activeChip = exposure.getDetector().getName()
827 
828  self.extract.debugView("extract", exposure)
829  result = self.extract.run(exposure)
830  ratios.append(result.outputRatios)
831 
832  for detIter, detector in enumerate(camera):
833  if detector.getName() == activeChip:
834  detectorId = detIter
835  outputDims = {'instrument': camera.getName(),
836  'detector': detectorId,
837  }
838 
839  finalResults = self.solver.run(ratios, camera=camera, outputDims=outputDims)
840  dataRef.put(finalResults.outputCrosstalk, "crosstalk")
841 
842  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 ddict2dict(d)
Definition: utils.py:721
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95