32 afwDisplay.setDefaultMaskTransparency(75)
36 """Produce a new frame number each time"""
45 """Options for measuring fringes on an exposure"""
46 badMaskPlanes =
ListField(dtype=str, default=[
"SAT"], doc=
"Ignore pixels with these masks")
47 stat =
Field(dtype=int, default=int(afwMath.MEDIAN), doc=
"Statistic to use")
48 clip =
Field(dtype=float, default=3.0, doc=
"Sigma clip threshold")
49 iterations =
Field(dtype=int, default=3, doc=
"Number of fitting iterations")
50 rngSeedOffset =
Field(dtype=int, default=0,
51 doc=
"Offset to the random number generator seed (full seed includes exposure ID)")
55 """Fringe subtraction options"""
56 filters =
ListField(dtype=str, default=[], doc=
"Only fringe-subtract these filters")
57 useFilterAliases =
Field(dtype=bool, default=
False, doc=
"Search filter aliases during check.")
58 num =
Field(dtype=int, default=30000, doc=
"Number of fringe measurements")
59 small =
Field(dtype=int, default=3, doc=
"Half-size of small (fringe) measurements (pixels)")
60 large =
Field(dtype=int, default=30, doc=
"Half-size of large (background) measurements (pixels)")
61 iterations =
Field(dtype=int, default=20, doc=
"Number of fitting iterations")
62 clip =
Field(dtype=float, default=3.0, doc=
"Sigma clip threshold")
63 stats =
ConfigField(dtype=FringeStatisticsConfig, doc=
"Statistics for measuring fringes")
64 pedestal =
Field(dtype=bool, default=
False, doc=
"Remove fringe pedestal?")
68 """Task to remove fringes from a science exposure
70 We measure fringe amplitudes at random positions on the science exposure
71 and at the same positions on the (potentially multiple) fringe frames
72 and solve for the scales simultaneously.
74 ConfigClass = FringeConfig
75 _DefaultName =
'isrFringe'
78 """Read the fringe frame(s), and pack data into a Struct
80 The current implementation assumes only a single fringe frame and
81 will have to be updated to support multi-mode fringe subtraction.
83 This implementation could be optimised by persisting the fringe
88 dataRef : `daf.butler.butlerSubset.ButlerDataRef`
89 Butler reference for the exposure that will have fringing
91 assembler : `lsst.ip.isr.AssembleCcdTask`, optional
92 An instance of AssembleCcdTask (for assembling fringe
97 fringeData : `pipeBase.Struct`
98 Struct containing fringe data:
99 - ``fringes`` : `lsst.afw.image.Exposure` or `list` thereof
100 Calibration fringe files containing master fringe frames.
101 - ``seed`` : `int`, optional
102 Seed for random number generation.
105 fringe = dataRef.get(
"fringe", immediate=
True)
106 except Exception
as e:
107 raise RuntimeError(
"Unable to retrieve fringe for %s: %s" % (dataRef.dataId, e))
112 """Pack the fringe data into a Struct.
114 This method moves the struct parsing code into a butler
115 generation agnostic handler.
119 fringeExp : `lsst.afw.exposure.Exposure`
120 The exposure containing the fringe data.
121 expId : `int`, optional
122 Exposure id to be fringe corrected, used to set RNG seed.
123 assembler : `lsst.ip.isr.AssembleCcdTask`, optional
124 An instance of AssembleCcdTask (for assembling fringe
129 fringeData : `pipeBase.Struct`
130 Struct containing fringe data:
131 - ``fringes`` : `lsst.afw.image.Exposure` or `list` thereof
132 Calibration fringe files containing master fringe frames.
133 - ``seed`` : `int`, optional
134 Seed for random number generation.
136 if assembler
is not None:
137 fringeExp = assembler.assembleCcd(fringeExp)
140 seed = self.
config.stats.rngSeedOffset
142 print(f
"{self.config.stats.rngSeedOffset} {expId}")
143 seed = self.
config.stats.rngSeedOffset + expId
148 return Struct(fringes=fringeExp,
152 def run(self, exposure, fringes, seed=None):
153 """Remove fringes from the provided science exposure.
155 Primary method of FringeTask. Fringes are only subtracted if the
156 science exposure has a filter listed in the configuration.
160 exposure : `lsst.afw.image.Exposure`
161 Science exposure from which to remove fringes.
162 fringes : `lsst.afw.image.Exposure` or `list` thereof
163 Calibration fringe files containing master fringe frames.
164 seed : `int`, optional
165 Seed for random number generation.
169 solution : `np.array`
170 Fringe solution amplitudes for each input fringe frame.
172 RMS error for the fit solution for this exposure.
178 self.
log.
info(
"Filter not found in FringeTaskConfig.filters. Skipping fringe correction.")
182 seed = self.
config.stats.rngSeedOffset
183 rng = numpy.random.RandomState(seed=seed)
185 if not hasattr(fringes,
'__iter__'):
188 mask = exposure.getMaskedImage().getMask()
189 for fringe
in fringes:
190 fringe.getMaskedImage().getMask().__ior__(mask)
195 fluxes = numpy.ndarray([self.
config.num, len(fringes)])
196 for i, f
in enumerate(fringes):
197 fluxes[:, i] = self.
measureExposure(f, positions, title=
"Fringe frame")
199 expFringes = self.
measureExposure(exposure, positions, title=
"Science")
200 solution, rms = self.
solve(expFringes, fluxes)
201 self.
subtract(exposure, fringes, solution)
203 afwDisplay.Display(frame=
getFrame()).
mtv(exposure, title=
"Fringe subtracted")
208 """Remove fringes from the provided science exposure.
210 Retrieve fringes from butler dataRef provided and remove from
211 provided science exposure. Fringes are only subtracted if the
212 science exposure has a filter listed in the configuration.
216 exposure : `lsst.afw.image.Exposure`
217 Science exposure from which to remove fringes.
218 dataRef : `daf.persistence.butlerSubset.ButlerDataRef`
219 Butler reference to the exposure. Used to find
220 appropriate fringe data.
221 assembler : `lsst.ip.isr.AssembleCcdTask`, optional
222 An instance of AssembleCcdTask (for assembling fringe
227 solution : `np.array`
228 Fringe solution amplitudes for each input fringe frame.
230 RMS error for the fit solution for this exposure.
233 self.
log.
info(
"Filter not found in FringeTaskConfig.filters. Skipping fringe correction.")
235 fringeStruct = self.
readFringes(dataRef, assembler=assembler)
236 return self.
run(exposure, **fringeStruct.getDict())
239 """Check whether we should fringe-subtract the science exposure.
243 exposure : `lsst.afw.image.Exposure`
244 Exposure to check the filter of.
249 If True, then the exposure has a filter listed in the
250 configuration, and should have the fringe applied.
253 if self.
config.useFilterAliases:
254 filterNameSet =
set(filterObj.getAliases() + [filterObj.getName()])
256 filterNameSet =
set([filterObj.getName(), ])
257 return bool(len(filterNameSet.intersection(self.
config.filters)))
260 """Remove pedestal from fringe exposure.
264 fringe : `lsst.afw.image.Exposure`
265 Fringe data to subtract the pedestal value from.
268 stats.setNumSigmaClip(self.
config.stats.clip)
269 stats.setNumIter(self.
config.stats.iterations)
270 mi = fringe.getMaskedImage()
272 self.
log.
info(
"Removing fringe pedestal: %f", pedestal)
276 """Generate a random distribution of positions for measuring fringe amplitudes.
280 exposure : `lsst.afw.image.Exposure`
281 Exposure to measure the positions on.
282 rng : `numpy.random.RandomState`
283 Random number generator to use.
287 positions : `numpy.array`
288 Two-dimensional array containing the positions to sample
289 for fringe amplitudes.
293 width = exposure.getWidth() - self.
config.large
294 height = exposure.getHeight() - self.
config.large
295 return numpy.array([rng.randint(start, width, size=num),
296 rng.randint(start, height, size=num)]).swapaxes(0, 1)
300 """Measure fringe amplitudes for an exposure
302 The fringe amplitudes are measured as the statistic within a square
303 aperture. The statistic within a larger aperture are subtracted so
304 as to remove the background.
308 exposure : `lsst.afw.image.Exposure`
309 Exposure to measure the positions on.
310 positions : `numpy.array`
311 Two-dimensional array containing the positions to sample
312 for fringe amplitudes.
313 title : `str`, optional
314 Title used for debug out plots.
318 fringes : `numpy.array`
319 Array of measured exposure values at each of the positions
323 stats.setNumSigmaClip(self.
config.stats.clip)
324 stats.setNumIter(self.
config.stats.iterations)
325 stats.setAndMask(exposure.getMaskedImage().getMask().getPlaneBitMask(self.
config.stats.badMaskPlanes))
328 fringes = numpy.ndarray(num)
332 small =
measure(exposure.getMaskedImage(), x, y, self.
config.small, self.
config.stats.stat, stats)
333 large =
measure(exposure.getMaskedImage(), x, y, self.
config.large, self.
config.stats.stat, stats)
334 fringes[i] = small - large
339 disp = afwDisplay.Display(frame=
getFrame())
340 disp.mtv(exposure, title=title)
342 with disp.Buffering():
343 for x, y
in positions:
344 corners = numpy.array([[-1, -1], [1, -1], [1, 1], [-1, 1], [-1, -1]]) + [[x, y]]
345 disp.line(corners*self.
config.small, ctype=afwDisplay.GREEN)
346 disp.line(corners*self.
config.large, ctype=afwDisplay.BLUE)
352 """Solve for the scale factors with iterative clipping.
356 science : `numpy.array`
357 Array of measured science image values at each of the
359 fringes : `numpy.array`
360 Array of measured fringe values at each of the positions
365 solution : `np.array`
366 Fringe solution amplitudes for each input fringe frame.
368 RMS error for the fit solution for this exposure.
373 origNum = len(science)
375 def emptyResult(msg=""):
376 """Generate an empty result for return to the user
378 There are no good pixels; doesn't matter what we return.
380 self.
log.
warn(
"Unable to solve for fringes: no good pixels%s", msg)
383 out = out*len(fringes)
384 return numpy.array(out), numpy.nan
386 good = numpy.where(numpy.logical_and(numpy.isfinite(science), numpy.any(numpy.isfinite(fringes), 1)))
387 science = science[good]
388 fringes = fringes[good]
389 oldNum = len(science)
396 for ff
in range(fringes.shape[1]):
398 science = science[good]
399 fringes = fringes[good]
400 oldNum = len(science)
402 return emptyResult(
" after initial rejection")
404 for i
in range(self.
config.iterations):
405 solution = self.
_solve(science, fringes)
406 resid = science - numpy.sum(solution*fringes, 1)
408 good = numpy.logical_not(
abs(resid) > self.
config.clip*rms)
409 self.
log.
debug(
"Iteration %d: RMS=%f numGood=%d", i, rms, good.sum())
410 self.
log.
debug(
"Solution %d: %s", i, solution)
413 return emptyResult(
" after %d rejection iterations" % i)
416 import matplotlib.pyplot
as plot
417 for j
in range(fringes.shape[1]):
421 fig.canvas._tkcanvas._root().lift()
424 ax = fig.add_subplot(1, 1, 1)
425 adjust = science.copy()
426 others =
set(range(fringes.shape[1]))
429 adjust -= solution[k]*fringes[:, k]
430 ax.plot(fringes[:, j], adjust,
'r.')
431 xmin = fringes[:, j].
min()
432 xmax = fringes[:, j].
max()
433 ymin = solution[j]*xmin
434 ymax = solution[j]*xmax
435 ax.plot([xmin, xmax], [ymin, ymax],
'b-')
436 ax.set_title(
"Fringe %d: %f" % (j, solution[j]))
437 ax.set_xlabel(
"Fringe amplitude")
438 ax.set_ylabel(
"Science amplitude")
439 ax.set_autoscale_on(
False)
440 ax.set_xbound(lower=xmin, upper=xmax)
441 ax.set_ybound(lower=ymin, upper=ymax)
444 ans = input(
"Enter or c to continue [chp]").lower()
445 if ans
in (
"",
"c",):
451 print(
"h[elp] c[ontinue] p[db]")
457 good = numpy.where(good)
458 science = science[good]
459 fringes = fringes[good]
462 solution = self.
_solve(science, fringes)
463 self.
log.
info(
"Fringe solution: %s RMS: %f Good: %d/%d", solution, rms, len(science), origNum)
466 def _solve(self, science, fringes):
467 """Solve for the scale factors.
471 science : `numpy.array`
472 Array of measured science image values at each of the
474 fringes : `numpy.array`
475 Array of measured fringe values at each of the positions
480 solution : `np.array`
481 Fringe solution amplitudes for each input fringe frame.
483 return afwMath.LeastSquares.fromDesignMatrix(fringes, science,
484 afwMath.LeastSquares.DIRECT_SVD).getSolution()
487 """Subtract the fringes.
491 science : `lsst.afw.image.Exposure`
492 Science exposure from which to remove fringes.
493 fringes : `lsst.afw.image.Exposure` or `list` thereof
494 Calibration fringe files containing master fringe frames.
495 solution : `np.array`
496 Fringe solution amplitudes for each input fringe frame.
501 Raised if the number of fringe frames does not match the
502 number of measured amplitudes.
504 if len(solution) != len(fringes):
505 raise RuntimeError(
"Number of fringe frames (%s) != number of scale factors (%s)" %
506 (len(fringes), len(solution)))
508 for s, f
in zip(solution, fringes):
509 science.getMaskedImage().scaledMinus(s, f.getMaskedImage())
512 def measure(mi, x, y, size, statistic, stats):
513 """Measure a statistic within an aperture
515 @param mi MaskedImage to measure
516 @param x, y Center for aperture
517 @param size Size of aperture
518 @param statistic Statistic to measure
519 @param stats StatisticsControl object
520 @return Value of statistic within aperture
524 subImage = mi.Factory(mi, bbox, afwImage.LOCAL)
529 """Calculate a robust standard deviation of an array of values
531 @param vector Array of values
532 @return Standard deviation
534 q1, q3 = numpy.percentile(vector, (25, 75))
535 return 0.74*(q3 - q1)
539 """Select values within 'clip' standard deviations of the median
541 Returns a boolean array.
543 q1, q2, q3 = numpy.percentile(vector, (25, 50, 75))
544 return numpy.abs(vector - q2) < clip*0.74*(q3 - q1)