LSSTApplications  18.1.0
LSSTDataManagementBasePackage
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 import numpy
23 
24 import lsst.geom
25 import lsst.afw.image as afwImage
26 import lsst.afw.math as afwMath
27 import lsst.afw.display as afwDisplay
28 
29 from lsst.pipe.base import Task, Struct, timeMethod
30 from lsst.pex.config import Config, Field, ListField, ConfigField
31 
32 afwDisplay.setDefaultMaskTransparency(75)
33 
34 
35 def getFrame():
36  """Produce a new frame number each time"""
37  getFrame.frame += 1
38  return getFrame.frame
39 
40 
41 getFrame.frame = 0
42 
43 
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)")
52 
53 
55  """Fringe subtraction options"""
56  filters = ListField(dtype=str, default=[], doc="Only fringe-subtract these filters")
57  num = Field(dtype=int, default=30000, doc="Number of fringe measurements")
58  small = Field(dtype=int, default=3, doc="Half-size of small (fringe) measurements (pixels)")
59  large = Field(dtype=int, default=30, doc="Half-size of large (background) measurements (pixels)")
60  iterations = Field(dtype=int, default=20, doc="Number of fitting iterations")
61  clip = Field(dtype=float, default=3.0, doc="Sigma clip threshold")
62  stats = ConfigField(dtype=FringeStatisticsConfig, doc="Statistics for measuring fringes")
63  pedestal = Field(dtype=bool, default=False, doc="Remove fringe pedestal?")
64 
65 
67  """Task to remove fringes from a science exposure
68 
69  We measure fringe amplitudes at random positions on the science exposure
70  and at the same positions on the (potentially multiple) fringe frames
71  and solve for the scales simultaneously.
72  """
73  ConfigClass = FringeConfig
74  _DefaultName = 'isrFringe'
75 
76  def readFringes(self, dataRef, assembler=None):
77  """Read the fringe frame(s)
78 
79  The current implementation assumes only a single fringe frame and
80  will have to be updated to support multi-mode fringe subtraction.
81 
82  This implementation could be optimised by persisting the fringe
83  positions and fluxes.
84 
85  @param dataRef Data reference for the science exposure
86  @param assembler An instance of AssembleCcdTask (for assembling fringe frames)
87  @return Struct(fringes: fringe exposure or list of fringe exposures;
88  seed: 32-bit uint derived from ccdExposureId for random number generator
89  """
90  try:
91  fringe = dataRef.get("fringe", immediate=True)
92  except Exception as e:
93  raise RuntimeError("Unable to retrieve fringe for %s: %s" % (dataRef.dataId, e))
94  if assembler is not None:
95  fringe = assembler.assembleCcd(fringe)
96 
97  seed = self.config.stats.rngSeedOffset + dataRef.get("ccdExposureId", immediate=True)
98  # Seed for numpy.random.RandomState must be convertable to a 32 bit unsigned integer
99  seed %= 2**32
100 
101  return Struct(fringes=fringe,
102  seed=seed)
103 
104  @timeMethod
105  def run(self, exposure, fringes, seed=None):
106  """Remove fringes from the provided science exposure.
107 
108  Primary method of FringeTask. Fringes are only subtracted if the
109  science exposure has a filter listed in the configuration.
110 
111  @param exposure Science exposure from which to remove fringes
112  @param fringes Exposure or list of Exposures
113  @param seed 32-bit unsigned integer for random number generator
114  """
115  import lsstDebug
116  display = lsstDebug.Info(__name__).display
117 
118  if not self.checkFilter(exposure):
119  return
120 
121  if seed is None:
122  seed = self.config.stats.rngSeedOffset
123  rng = numpy.random.RandomState(seed=seed)
124 
125  if not hasattr(fringes, '__iter__'):
126  fringes = [fringes]
127 
128  mask = exposure.getMaskedImage().getMask()
129  for fringe in fringes:
130  fringe.getMaskedImage().getMask().__ior__(mask)
131  if self.config.pedestal:
132  self.removePedestal(fringe)
133 
134  positions = self.generatePositions(fringes[0], rng)
135  fluxes = numpy.ndarray([self.config.num, len(fringes)])
136  for i, f in enumerate(fringes):
137  fluxes[:, i] = self.measureExposure(f, positions, title="Fringe frame")
138 
139  expFringes = self.measureExposure(exposure, positions, title="Science")
140  solution, rms = self.solve(expFringes, fluxes)
141  self.subtract(exposure, fringes, solution)
142  if display:
143  afwDisplay.Display(frame=getFrame()).mtv(exposure, title="Fringe subtracted")
144  return solution, rms
145 
146  @timeMethod
147  def runDataRef(self, exposure, dataRef, assembler=None):
148  """Remove fringes from the provided science exposure.
149 
150  Retrieve fringes from butler dataRef provided and remove from
151  provided science exposure.
152  Fringes are only subtracted if the science exposure has a filter
153  listed in the configuration.
154 
155  @param exposure Science exposure from which to remove fringes
156  @param dataRef Data reference for the science exposure
157  @param assembler An instance of AssembleCcdTask (for assembling fringe frames)
158  """
159  if not self.checkFilter(exposure):
160  return
161  fringeStruct = self.readFringes(dataRef, assembler=assembler)
162  return self.run(exposure, **fringeStruct.getDict())
163 
164  def checkFilter(self, exposure):
165  """Check whether we should fringe-subtract the science exposure"""
166  return exposure.getFilter().getName() in self.config.filters
167 
168  def removePedestal(self, fringe):
169  """Remove pedestal from fringe exposure"""
170  stats = afwMath.StatisticsControl()
171  stats.setNumSigmaClip(self.config.stats.clip)
172  stats.setNumIter(self.config.stats.iterations)
173  mi = fringe.getMaskedImage()
174  pedestal = afwMath.makeStatistics(mi, afwMath.MEDIAN, stats).getValue()
175  self.log.info("Removing fringe pedestal: %f", pedestal)
176  mi -= pedestal
177 
178  def generatePositions(self, exposure, rng):
179  """Generate a random distribution of positions for measuring fringe amplitudes"""
180  start = self.config.large
181  num = self.config.num
182  width = exposure.getWidth() - self.config.large
183  height = exposure.getHeight() - self.config.large
184  return numpy.array([rng.randint(start, width, size=num),
185  rng.randint(start, height, size=num)]).swapaxes(0, 1)
186 
187  @timeMethod
188  def measureExposure(self, exposure, positions, title="Fringe"):
189  """Measure fringe amplitudes for an exposure
190 
191  The fringe amplitudes are measured as the statistic within a square
192  aperture. The statistic within a larger aperture are subtracted so
193  as to remove the background.
194 
195  @param exposure Exposure to measure
196  @param positions Array of (x,y) for fringe measurement
197  @param title Title for display
198  @return Array of fringe measurements
199  """
200  stats = afwMath.StatisticsControl()
201  stats.setNumSigmaClip(self.config.stats.clip)
202  stats.setNumIter(self.config.stats.iterations)
203  stats.setAndMask(exposure.getMaskedImage().getMask().getPlaneBitMask(self.config.stats.badMaskPlanes))
204 
205  num = self.config.num
206  fringes = numpy.ndarray(num)
207 
208  for i in range(num):
209  x, y = positions[i]
210  small = measure(exposure.getMaskedImage(), x, y, self.config.small, self.config.stats.stat, stats)
211  large = measure(exposure.getMaskedImage(), x, y, self.config.large, self.config.stats.stat, stats)
212  fringes[i] = small - large
213 
214  import lsstDebug
215  display = lsstDebug.Info(__name__).display
216  if display:
217  disp = afwDisplay.Display(frame=getFrame())
218  disp.mtv(exposure, title=title)
219  if False:
220  with disp.Buffering():
221  for x, y in positions:
222  corners = numpy.array([[-1, -1], [1, -1], [1, 1], [-1, 1], [-1, -1]]) + [[x, y]]
223  disp.line(corners*self.config.small, ctype=afwDisplay.GREEN)
224  disp.line(corners*self.config.large, ctype=afwDisplay.BLUE)
225 
226  return fringes
227 
228  @timeMethod
229  def solve(self, science, fringes):
230  """Solve (with iterative clipping) for the scale factors
231 
232  @param science Array of science exposure fringe amplitudes
233  @param fringes Array of arrays of fringe frame fringe amplitudes
234  @return Array of scale factors for the fringe frames
235  """
236  import lsstDebug
237  doPlot = lsstDebug.Info(__name__).plot
238 
239  origNum = len(science)
240 
241  def emptyResult(msg=""):
242  """Generate an empty result for return to the user
243 
244  There are no good pixels; doesn't matter what we return.
245  """
246  self.log.warn("Unable to solve for fringes: no good pixels%s", msg)
247  out = [0]
248  if len(fringes) > 1:
249  out = out*len(fringes)
250  return numpy.array(out), numpy.nan
251 
252  good = numpy.where(numpy.logical_and(numpy.isfinite(science), numpy.any(numpy.isfinite(fringes), 1)))
253  science = science[good]
254  fringes = fringes[good]
255  oldNum = len(science)
256  if oldNum == 0:
257  return emptyResult()
258 
259  # Up-front rejection to get rid of extreme, potentially troublesome values
260  # (e.g., fringe apertures that fall on objects).
261  good = select(science, self.config.clip)
262  for ff in range(fringes.shape[1]):
263  good &= select(fringes[:, ff], self.config.clip)
264  science = science[good]
265  fringes = fringes[good]
266  oldNum = len(science)
267  if oldNum == 0:
268  return emptyResult(" after initial rejection")
269 
270  for i in range(self.config.iterations):
271  solution = self._solve(science, fringes)
272  resid = science - numpy.sum(solution*fringes, 1)
273  rms = stdev(resid)
274  good = numpy.logical_not(abs(resid) > self.config.clip*rms)
275  self.log.debug("Iteration %d: RMS=%f numGood=%d", i, rms, good.sum())
276  self.log.debug("Solution %d: %s", i, solution)
277  newNum = good.sum()
278  if newNum == 0:
279  return emptyResult(" after %d rejection iterations" % i)
280 
281  if doPlot:
282  import matplotlib.pyplot as plot
283  for j in range(fringes.shape[1]):
284  fig = plot.figure(j)
285  fig.clf()
286  try:
287  fig.canvas._tkcanvas._root().lift() # == Tk's raise
288  except Exception:
289  pass
290  ax = fig.add_subplot(1, 1, 1)
291  adjust = science.copy()
292  others = set(range(fringes.shape[1]))
293  others.discard(j)
294  for k in others:
295  adjust -= solution[k]*fringes[:, k]
296  ax.plot(fringes[:, j], adjust, 'r.')
297  xmin = fringes[:, j].min()
298  xmax = fringes[:, j].max()
299  ymin = solution[j]*xmin
300  ymax = solution[j]*xmax
301  ax.plot([xmin, xmax], [ymin, ymax], 'b-')
302  ax.set_title("Fringe %d: %f" % (j, solution[j]))
303  ax.set_xlabel("Fringe amplitude")
304  ax.set_ylabel("Science amplitude")
305  ax.set_autoscale_on(False)
306  ax.set_xbound(lower=xmin, upper=xmax)
307  ax.set_ybound(lower=ymin, upper=ymax)
308  fig.show()
309  while True:
310  ans = input("Enter or c to continue [chp]").lower()
311  if ans in ("", "c",):
312  break
313  if ans in ("p",):
314  import pdb
315  pdb.set_trace()
316  elif ans in ("h", ):
317  print("h[elp] c[ontinue] p[db]")
318 
319  if newNum == oldNum:
320  # Not gaining
321  break
322  oldNum = newNum
323  good = numpy.where(good)
324  science = science[good]
325  fringes = fringes[good]
326 
327  # Final solution without rejection
328  solution = self._solve(science, fringes)
329  self.log.info("Fringe solution: %s RMS: %f Good: %d/%d", solution, rms, len(science), origNum)
330  return solution, rms
331 
332  def _solve(self, science, fringes):
333  """Solve for the scale factors
334 
335  @param science Array of science exposure fringe amplitudes
336  @param fringes Array of arrays of fringe frame fringe amplitudes
337  @return Array of scale factors for the fringe frames
338  """
339  return afwMath.LeastSquares.fromDesignMatrix(fringes, science,
340  afwMath.LeastSquares.DIRECT_SVD).getSolution()
341 
342  def subtract(self, science, fringes, solution):
343  """Subtract the fringes
344 
345  @param science Science exposure
346  @param fringes List of fringe frames
347  @param solution Array of scale factors for the fringe frames
348  """
349  if len(solution) != len(fringes):
350  raise RuntimeError("Number of fringe frames (%s) != number of scale factors (%s)" %
351  (len(fringes), len(solution)))
352 
353  for s, f in zip(solution, fringes):
354  science.getMaskedImage().scaledMinus(s, f.getMaskedImage())
355 
356 
357 def measure(mi, x, y, size, statistic, stats):
358  """Measure a statistic within an aperture
359 
360  @param mi MaskedImage to measure
361  @param x, y Center for aperture
362  @param size Size of aperture
363  @param statistic Statistic to measure
364  @param stats StatisticsControl object
365  @return Value of statistic within aperture
366  """
367  bbox = lsst.geom.Box2I(lsst.geom.Point2I(int(x) - size, int(y - size)),
368  lsst.geom.Extent2I(2*size, 2*size))
369  subImage = mi.Factory(mi, bbox, afwImage.LOCAL)
370  return afwMath.makeStatistics(subImage, statistic, stats).getValue()
371 
372 
373 def stdev(vector):
374  """Calculate a robust standard deviation of an array of values
375 
376  @param vector Array of values
377  @return Standard deviation
378  """
379  q1, q3 = numpy.percentile(vector, (25, 75))
380  return 0.74*(q3 - q1)
381 
382 
383 def select(vector, clip):
384  """Select values within 'clip' standard deviations of the median
385 
386  Returns a boolean array.
387  """
388  q1, q2, q3 = numpy.percentile(vector, (25, 50, 75))
389  return numpy.abs(vector - q2) < clip*0.74*(q3 - q1)
def runDataRef(self, exposure, dataRef, assembler=None)
Definition: fringe.py:147
Angle abs(Angle const &a)
Definition: Angle.h:106
def run(self, exposure, fringes, seed=None)
Definition: fringe.py:105
def checkFilter(self, exposure)
Definition: fringe.py:164
def subtract(self, science, fringes, solution)
Definition: fringe.py:342
def _solve(self, science, fringes)
Definition: fringe.py:332
daf::base::PropertySet * set
Definition: fits.cc:884
int min
Statistics makeStatistics(lsst::afw::math::MaskedVector< EntryT > const &mv, std::vector< WeightPixel > const &vweights, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle lsst::afw::math::MaskedVector<>
Definition: Statistics.h:520
def stdev(vector)
Definition: fringe.py:373
Pass parameters to a Statistics object.
Definition: Statistics.h:93
def generatePositions(self, exposure, rng)
Definition: fringe.py:178
def measureExposure(self, exposure, positions, title="Fringe")
Definition: fringe.py:188
def getName(self)
Definition: task.py:250
int max
def removePedestal(self, fringe)
Definition: fringe.py:168
def mtv(data, frame=None, title="", wcs=None, args, kwargs)
Definition: ds9.py:93
def readFringes(self, dataRef, assembler=None)
Definition: fringe.py:76
def select(vector, clip)
Definition: fringe.py:383
def measure(mi, x, y, size, statistic, stats)
Definition: fringe.py:357
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
An integer coordinate rectangle.
Definition: Box.h:54
def solve(self, science, fringes)
Definition: fringe.py:229