LSST Applications 24.1.6,g063fba187b+56b85ce14a,g0f08755f38+df8a265115,g12f32b3c4e+891a09f10d,g1524ad2192+7a5d7b3fbd,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g28da252d5a+07cb1400be,g2bbee38e9b+ae03bbfc84,g2bc492864f+ae03bbfc84,g3156d2b45e+6e55a43351,g347aa1857d+ae03bbfc84,g35bb328faa+a8ce1bb630,g3a166c0a6a+ae03bbfc84,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+8f257c4c0b,g781aacb6e4+a8ce1bb630,g7af13505b9+7137b3b17d,g80478fca09+6df6903293,g82479be7b0+091ce1d07f,g858d7b2824+df8a265115,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,g9726552aa6+414189b318,ga5288a1d22+4a2bca08d7,gacef1a1666+c9a8ff65f4,gb58c049af0+d64f4d3760,gbcfae0f0a0+de1d42d831,gc28159a63d+ae03bbfc84,gcf0d15dbbd+72117bf34e,gda6a2b7d83+72117bf34e,gdaeeff99f8+1711a396fd,ge500cccec5+c8c9c9af63,ge79ae78c31+ae03bbfc84,gf0baf85859+c1f95f4921,gfa517265be+df8a265115,gfa999e8aa5+17cd334064,gfb92a5be7c+df8a265115
LSST Data Management Base Package
Loading...
Searching...
No Matches
fringe.py
Go to the documentation of this file.
1# This file is part of ip_isr.
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 <https://www.gnu.org/licenses/>.
21
22__all__ = ["FringeStatisticsConfig", "FringeConfig", "FringeTask"]
23
24import numpy
25
26import lsst.geom
27import lsst.afw.image as afwImage
28import lsst.afw.math as afwMath
29import lsst.afw.display as afwDisplay
30
31from lsst.pipe.base import Task, Struct
32from lsst.pex.config import Config, Field, ListField, ConfigField
33from lsst.utils.timer import timeMethod
34from .isrFunctions import checkFilter
35
36afwDisplay.setDefaultMaskTransparency(75)
37
38
40 """Produce a new frame number each time"""
41 getFrame.frame += 1
42 return getFrame.frame
43
44
45getFrame.frame = 0
46
47
49 """Options for measuring fringes on an exposure"""
50 badMaskPlanes = ListField(dtype=str, default=["SAT"], doc="Ignore pixels with these masks")
51 stat = Field(dtype=int, default=int(afwMath.MEDIAN), doc="Statistic to use")
52 clip = Field(dtype=float, default=3.0, doc="Sigma clip threshold")
53 iterations = Field(dtype=int, default=3, doc="Number of fitting iterations")
54 rngSeedOffset = Field(dtype=int, default=0,
55 doc="Offset to the random number generator seed (full seed includes exposure ID)")
56
57
59 """Fringe subtraction options"""
60 # TODO DM-28093: change the doc to specify that these are physical labels
61 filters = ListField(dtype=str, default=[], doc="Only fringe-subtract these filters")
62 # TODO: remove in DM-27177
63 useFilterAliases = Field(dtype=bool, default=False, doc="Search filter aliases during check.",
64 deprecated=("Removed with no replacement (FilterLabel has no aliases)."
65 "Will be removed after v22."))
66 num = Field(dtype=int, default=30000, doc="Number of fringe measurements")
67 small = Field(dtype=int, default=3, doc="Half-size of small (fringe) measurements (pixels)")
68 large = Field(dtype=int, default=30, doc="Half-size of large (background) measurements (pixels)")
69 iterations = Field(dtype=int, default=20, doc="Number of fitting iterations")
70 clip = Field(dtype=float, default=3.0, doc="Sigma clip threshold")
71 stats = ConfigField(dtype=FringeStatisticsConfig, doc="Statistics for measuring fringes")
72 pedestal = Field(dtype=bool, default=False, doc="Remove fringe pedestal?")
73
74
75class FringeTask(Task):
76 """Task to remove fringes from a science exposure
77
78 We measure fringe amplitudes at random positions on the science exposure
79 and at the same positions on the (potentially multiple) fringe frames
80 and solve for the scales simultaneously.
81 """
82 ConfigClass = FringeConfig
83 _DefaultName = 'isrFringe'
84
85 def loadFringes(self, fringeExp, expId=None, assembler=None):
86 """Pack the fringe data into a Struct.
87
88 This method moves the struct parsing code into a butler
89 generation agnostic handler.
90
91 Parameters
92 ----------
93 fringeExp : `lsst.afw.exposure.Exposure`
94 The exposure containing the fringe data.
95 expId : `int`, optional
96 Exposure id to be fringe corrected, used to set RNG seed.
97 assembler : `lsst.ip.isr.AssembleCcdTask`, optional
98 An instance of AssembleCcdTask (for assembling fringe
99 frames).
100
101 Returns
102 -------
103 fringeData : `pipeBase.Struct`
104 Struct containing fringe data:
105
106 ``fringes``
107 Calibration fringe files containing master fringe frames.
108 ( : `lsst.afw.image.Exposure` or `list` thereof)
109 ``seed``
110 Seed for random number generation. (`int`, optional)
111 """
112 if assembler is not None:
113 fringeExp = assembler.assembleCcd(fringeExp)
114
115 if expId is None:
116 seed = self.config.stats.rngSeedOffset
117 else:
118 self.log.debug("Seeding with offset %d and ccdExposureId %d.",
119 self.config.stats.rngSeedOffset, expId)
120 seed = self.config.stats.rngSeedOffset + expId
121
122 # Seed for numpy.random.RandomState must be convertable to a 32 bit
123 # unsigned integer.
124 seed %= 2**32
125
126 return Struct(fringes=fringeExp,
127 seed=seed)
128
129 @timeMethod
130 def run(self, exposure, fringes, seed=None):
131 """Remove fringes from the provided science exposure.
132
133 Primary method of FringeTask. Fringes are only subtracted if the
134 science exposure has a filter listed in the configuration.
135
136 Parameters
137 ----------
138 exposure : `lsst.afw.image.Exposure`
139 Science exposure from which to remove fringes.
140 fringes : `lsst.afw.image.Exposure` or `list` thereof
141 Calibration fringe files containing master fringe frames.
142 seed : `int`, optional
143 Seed for random number generation.
144
145 Returns
146 -------
147 solution : `np.array`
148 Fringe solution amplitudes for each input fringe frame.
149 rms : `float`
150 RMS error for the fit solution for this exposure.
151 """
152 import lsstDebug
153 display = lsstDebug.Info(__name__).display
154
155 if not self.checkFilter(exposure):
156 self.log.info("Filter not found in FringeTaskConfig.filters. Skipping fringe correction.")
157 return
158
159 if seed is None:
160 seed = self.config.stats.rngSeedOffset
161 rng = numpy.random.RandomState(seed=seed)
162
163 if not hasattr(fringes, '__iter__'):
164 fringes = [fringes]
165
166 mask = exposure.getMaskedImage().getMask()
167 for fringe in fringes:
168 fringe.getMaskedImage().getMask().__ior__(mask)
169 if self.config.pedestal:
170 self.removePedestal(fringe)
171
172 positions = self.generatePositions(fringes[0], rng)
173 fluxes = numpy.ndarray([self.config.num, len(fringes)])
174 for i, f in enumerate(fringes):
175 fluxes[:, i] = self.measureExposure(f, positions, title="Fringe frame")
176
177 expFringes = self.measureExposure(exposure, positions, title="Science")
178 solution, rms = self.solve(expFringes, fluxes)
179 self.subtract(exposure, fringes, solution)
180 if display:
181 afwDisplay.Display(frame=getFrame()).mtv(exposure, title="Fringe subtracted")
182 return solution, rms
183
184 def checkFilter(self, exposure):
185 """Check whether we should fringe-subtract the science exposure.
186
187 Parameters
188 ----------
189 exposure : `lsst.afw.image.Exposure`
190 Exposure to check the filter of.
191
192 Returns
193 -------
194 needsFringe : `bool`
195 If True, then the exposure has a filter listed in the
196 configuration, and should have the fringe applied.
197 """
198 return checkFilter(exposure, self.config.filters, log=self.log)
199
200 def removePedestal(self, fringe):
201 """Remove pedestal from fringe exposure.
202
203 Parameters
204 ----------
205 fringe : `lsst.afw.image.Exposure`
206 Fringe data to subtract the pedestal value from.
207 """
209 stats.setNumSigmaClip(self.config.stats.clip)
210 stats.setNumIter(self.config.stats.iterations)
211 mi = fringe.getMaskedImage()
212 pedestal = afwMath.makeStatistics(mi, afwMath.MEDIAN, stats).getValue()
213 self.log.info("Removing fringe pedestal: %f", pedestal)
214 mi -= pedestal
215
216 def generatePositions(self, exposure, rng):
217 """Generate a random distribution of positions for measuring fringe
218 amplitudes.
219
220 Parameters
221 ----------
222 exposure : `lsst.afw.image.Exposure`
223 Exposure to measure the positions on.
224 rng : `numpy.random.RandomState`
225 Random number generator to use.
226
227 Returns
228 -------
229 positions : `numpy.array`
230 Two-dimensional array containing the positions to sample
231 for fringe amplitudes.
232 """
233 start = self.config.large
234 num = self.config.num
235 width = exposure.getWidth() - self.config.large
236 height = exposure.getHeight() - self.config.large
237 return numpy.array([rng.randint(start, width, size=num),
238 rng.randint(start, height, size=num)]).swapaxes(0, 1)
239
240 @timeMethod
241 def measureExposure(self, exposure, positions, title="Fringe"):
242 """Measure fringe amplitudes for an exposure
243
244 The fringe amplitudes are measured as the statistic within a square
245 aperture. The statistic within a larger aperture are subtracted so
246 as to remove the background.
247
248 Parameters
249 ----------
250 exposure : `lsst.afw.image.Exposure`
251 Exposure to measure the positions on.
252 positions : `numpy.array`
253 Two-dimensional array containing the positions to sample
254 for fringe amplitudes.
255 title : `str`, optional
256 Title used for debug out plots.
257
258 Returns
259 -------
260 fringes : `numpy.array`
261 Array of measured exposure values at each of the positions
262 supplied.
263 """
265 stats.setNumSigmaClip(self.config.stats.clip)
266 stats.setNumIter(self.config.stats.iterations)
267 stats.setAndMask(exposure.getMaskedImage().getMask().getPlaneBitMask(self.config.stats.badMaskPlanes))
268
269 num = self.config.num
270 fringes = numpy.ndarray(num)
271
272 for i in range(num):
273 x, y = positions[i]
274 small = measure(exposure.getMaskedImage(), x, y, self.config.small, self.config.stats.stat, stats)
275 large = measure(exposure.getMaskedImage(), x, y, self.config.large, self.config.stats.stat, stats)
276 fringes[i] = small - large
277
278 import lsstDebug
279 display = lsstDebug.Info(__name__).display
280 if display:
281 disp = afwDisplay.Display(frame=getFrame())
282 disp.mtv(exposure, title=title)
283 if False:
284 with disp.Buffering():
285 for x, y in positions:
286 corners = numpy.array([[-1, -1], [1, -1], [1, 1], [-1, 1], [-1, -1]]) + [[x, y]]
287 disp.line(corners*self.config.small, ctype=afwDisplay.GREEN)
288 disp.line(corners*self.config.large, ctype=afwDisplay.BLUE)
289
290 return fringes
291
292 @timeMethod
293 def solve(self, science, fringes):
294 """Solve for the scale factors with iterative clipping.
295
296 Parameters
297 ----------
298 science : `numpy.array`
299 Array of measured science image values at each of the
300 positions supplied.
301 fringes : `numpy.array`
302 Array of measured fringe values at each of the positions
303 supplied.
304
305 Returns
306 -------
307 solution : `np.array`
308 Fringe solution amplitudes for each input fringe frame.
309 rms : `float`
310 RMS error for the fit solution for this exposure.
311 """
312 import lsstDebug
313 doPlot = lsstDebug.Info(__name__).plot
314
315 origNum = len(science)
316
317 def emptyResult(msg=""):
318 """Generate an empty result for return to the user
319
320 There are no good pixels; doesn't matter what we return.
321 """
322 self.log.warning("Unable to solve for fringes: no good pixels%s", msg)
323 out = [0]
324 if len(fringes) > 1:
325 out = out*len(fringes)
326 return numpy.array(out), numpy.nan
327
328 good = numpy.where(numpy.logical_and(numpy.isfinite(science), numpy.any(numpy.isfinite(fringes), 1)))
329 science = science[good]
330 fringes = fringes[good]
331 oldNum = len(science)
332 if oldNum == 0:
333 return emptyResult()
334
335 # Up-front rejection to get rid of extreme, potentially troublesome
336 # values (e.g., fringe apertures that fall on objects).
337 good = select(science, self.config.clip)
338 for ff in range(fringes.shape[1]):
339 good &= select(fringes[:, ff], self.config.clip)
340 science = science[good]
341 fringes = fringes[good]
342 oldNum = len(science)
343 if oldNum == 0:
344 return emptyResult(" after initial rejection")
345
346 for i in range(self.config.iterations):
347 solution = self._solve(science, fringes)
348 resid = science - numpy.sum(solution*fringes, 1)
349 rms = stdev(resid)
350 good = numpy.logical_not(abs(resid) > self.config.clip*rms)
351 self.log.debug("Iteration %d: RMS=%f numGood=%d", i, rms, good.sum())
352 self.log.debug("Solution %d: %s", i, solution)
353 newNum = good.sum()
354 if newNum == 0:
355 return emptyResult(" after %d rejection iterations" % i)
356
357 if doPlot:
358 import matplotlib.pyplot as plot
359 for j in range(fringes.shape[1]):
360 fig = plot.figure(j)
361 fig.clf()
362 try:
363 fig.canvas._tkcanvas._root().lift() # == Tk's raise
364 except Exception:
365 pass
366 ax = fig.add_subplot(1, 1, 1)
367 adjust = science.copy()
368 others = set(range(fringes.shape[1]))
369 others.discard(j)
370 for k in others:
371 adjust -= solution[k]*fringes[:, k]
372 ax.plot(fringes[:, j], adjust, 'r.')
373 xmin = fringes[:, j].min()
374 xmax = fringes[:, j].max()
375 ymin = solution[j]*xmin
376 ymax = solution[j]*xmax
377 ax.plot([xmin, xmax], [ymin, ymax], 'b-')
378 ax.set_title("Fringe %d: %f" % (j, solution[j]))
379 ax.set_xlabel("Fringe amplitude")
380 ax.set_ylabel("Science amplitude")
381 ax.set_autoscale_on(False)
382 ax.set_xbound(lower=xmin, upper=xmax)
383 ax.set_ybound(lower=ymin, upper=ymax)
384 fig.show()
385 while True:
386 ans = input("Enter or c to continue [chp]").lower()
387 if ans in ("", "c",):
388 break
389 if ans in ("p",):
390 import pdb
391 pdb.set_trace()
392 elif ans in ("h", ):
393 print("h[elp] c[ontinue] p[db]")
394
395 if newNum == oldNum:
396 # Not gaining
397 break
398 oldNum = newNum
399 good = numpy.where(good)
400 science = science[good]
401 fringes = fringes[good]
402
403 # Final solution without rejection
404 solution = self._solve(science, fringes)
405 self.log.info("Fringe solution: %s RMS: %f Good: %d/%d", solution, rms, len(science), origNum)
406 return solution, rms
407
408 def _solve(self, science, fringes):
409 """Solve for the scale factors.
410
411 Parameters
412 ----------
413 science : `numpy.array`
414 Array of measured science image values at each of the
415 positions supplied.
416 fringes : `numpy.array`
417 Array of measured fringe values at each of the positions
418 supplied.
419
420 Returns
421 -------
422 solution : `np.array`
423 Fringe solution amplitudes for each input fringe frame.
424 """
425 return afwMath.LeastSquares.fromDesignMatrix(fringes, science,
426 afwMath.LeastSquares.DIRECT_SVD).getSolution()
427
428 def subtract(self, science, fringes, solution):
429 """Subtract the fringes.
430
431 Parameters
432 ----------
433 science : `lsst.afw.image.Exposure`
434 Science exposure from which to remove fringes.
435 fringes : `lsst.afw.image.Exposure` or `list` thereof
436 Calibration fringe files containing master fringe frames.
437 solution : `np.array`
438 Fringe solution amplitudes for each input fringe frame.
439
440 Raises
441 ------
442 RuntimeError
443 Raised if the number of fringe frames does not match the
444 number of measured amplitudes.
445 """
446 if len(solution) != len(fringes):
447 raise RuntimeError("Number of fringe frames (%s) != number of scale factors (%s)" %
448 (len(fringes), len(solution)))
449
450 for s, f in zip(solution, fringes):
451 # We do not want to add the mask from the fringe to the image.
452 f.getMaskedImage().getMask().getArray()[:] = 0
453 science.getMaskedImage().scaledMinus(s, f.getMaskedImage())
454
455
456def measure(mi, x, y, size, statistic, stats):
457 """Measure a statistic within an aperture
458
459 @param mi MaskedImage to measure
460 @param x, y Center for aperture
461 @param size Size of aperture
462 @param statistic Statistic to measure
463 @param stats StatisticsControl object
464 @return Value of statistic within aperture
465 """
466 bbox = lsst.geom.Box2I(lsst.geom.Point2I(int(x) - size, int(y - size)),
467 lsst.geom.Extent2I(2*size, 2*size))
468 subImage = mi.Factory(mi, bbox, afwImage.LOCAL)
469 return afwMath.makeStatistics(subImage, statistic, stats).getValue()
470
471
472def stdev(vector):
473 """Calculate a robust standard deviation of an array of values
474
475 @param vector Array of values
476 @return Standard deviation
477 """
478 q1, q3 = numpy.percentile(vector, (25, 75))
479 return 0.74*(q3 - q1)
480
481
482def select(vector, clip):
483 """Select values within 'clip' standard deviations of the median
484
485 Returns a boolean array.
486 """
487 q1, q2, q3 = numpy.percentile(vector, (25, 50, 75))
488 return numpy.abs(vector - q2) < clip*0.74*(q3 - q1)
int min
int max
Pass parameters to a Statistics object.
Definition Statistics.h:83
An integer coordinate rectangle.
Definition Box.h:55
loadFringes(self, fringeExp, expId=None, assembler=None)
Definition fringe.py:85
checkFilter(self, exposure)
Definition fringe.py:184
generatePositions(self, exposure, rng)
Definition fringe.py:216
solve(self, science, fringes)
Definition fringe.py:293
_solve(self, science, fringes)
Definition fringe.py:408
subtract(self, science, fringes, solution)
Definition fringe.py:428
measureExposure(self, exposure, positions, title="Fringe")
Definition fringe.py:241
run(self, exposure, fringes, seed=None)
Definition fringe.py:130
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361
select(vector, clip)
Definition fringe.py:482
measure(mi, x, y, size, statistic, stats)
Definition fringe.py:456