LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
diffimTools.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 # python
24 import time
25 import os
26 from collections import Counter
27 import numpy as np
28 
29 # all the c++ level classes and routines
30 import diffimLib
31 
32 # all the other LSST packages
33 import lsst.afw.geom as afwGeom
34 import lsst.afw.image as afwImage
35 import lsst.afw.table as afwTable
36 import lsst.afw.detection as afwDetect
37 import lsst.afw.math.mathLib as afwMath
38 import lsst.pex.logging as pexLog
39 import lsst.pex.config as pexConfig
40 from .makeKernelBasisList import makeKernelBasisList
41 
42 # Helper functions for ipDiffim; mostly viewing of results and writing
43 # debugging info to disk.
44 
45 #######
46 # Add noise
47 #######
48 
49 def makeFlatNoiseImage(mi, seedStat = afwMath.MAX):
50  img = mi.getImage()
51  seed = int(10. * afwMath.makeStatistics(mi.getImage(), seedStat).getValue() + 1)
52  rdm = afwMath.Random(afwMath.Random.MT19937, seed)
53  rdmImage = img.Factory(img.getDimensions())
54  afwMath.randomGaussianImage(rdmImage, rdm)
55  return rdmImage
56 
58  """Return a Poisson noise image based on im
59 
60  Uses numpy.random; you may wish to call numpy.random.seed first.
61 
62  @warning This uses an undocumented numpy API (the documented API
63  uses a single float expectation value instead of an array).
64 
65  @param[in] im image; the output image has the same dimensions and shape
66  and its expectation value is the value of im at each pixel
67  """
68  import numpy.random as rand
69  imArr = im.getArray()
70  noiseIm = im.Factory(im.getBBox())
71  noiseArr = noiseIm.getArray()
72 
73  intNoiseArr = rand.poisson(imArr)
74  noiseArr[:, :] = intNoiseArr.astype(noiseArr.dtype)
75  return noiseIm
76 
77 #######
78 # Make fake images for testing; one is a delta function (or narrow
79 # gaussian) and the other is a convolution of this with a spatially
80 # varying kernel.
81 #######
82 def fakeCoeffs():
83  kCoeffs = (( 1.0, 0.0, 0.0),
84  ( 0.005, -0.000001, 0.000001),
85  ( 0.005, 0.000004, 0.000004),
86  ( -0.001, -0.000030, 0.000030),
87  ( -0.001, 0.000015, 0.000015),
88  ( -0.005, -0.000050, 0.000050))
89  return kCoeffs
90 
91 def makeFakeKernelSet(sizeCell = 128, nCell = 3,
92  deltaFunctionCounts = 1.e4, tGaussianWidth = 1.0,
93  addNoise = True, bgValue = 100., display = False):
94 
95  from . import imagePsfMatch
97  configFake.kernel.name = "AL"
98  subconfigFake = configFake.kernel.active
99  subconfigFake.alardNGauss = 1
100  subconfigFake.alardSigGauss = [2.5,]
101  subconfigFake.alardDegGauss = [2,]
102  subconfigFake.sizeCellX = sizeCell
103  subconfigFake.sizeCellY = sizeCell
104  subconfigFake.spatialKernelOrder = 1
105  subconfigFake.spatialModelType = "polynomial"
106  subconfigFake.singleKernelClipping = False # variance is a hack
107  subconfigFake.spatialKernelClipping = False # variance is a hack
108  if bgValue > 0.0:
109  subconfigFake.fitForBackground = True
110 
111  policyFake = pexConfig.makePolicy(subconfigFake)
112 
113  basisList = makeKernelBasisList(subconfigFake)
114  kSize = subconfigFake.kernelSize
115 
116  # This sets the final extent of each convolved delta function
117  gaussKernelWidth = sizeCell // 2
118 
119  # This sets the scale over which pixels are correlated in the
120  # spatial convolution; should be at least as big as the kernel you
121  # are trying to fit for
122  spatialKernelWidth = kSize
123 
124  # Number of bad pixels due to convolutions
125  border = (gaussKernelWidth + spatialKernelWidth)//2
126 
127  # Make a fake image with a matrix of delta functions
128  totalSize = nCell * sizeCell + 2 * border
129  tim = afwImage.ImageF(afwGeom.Extent2I(totalSize, totalSize))
130  for x in range(nCell):
131  for y in range(nCell):
132  tim.set(x * sizeCell + sizeCell // 2 + border - 1,
133  y * sizeCell + sizeCell // 2 + border - 1,
134  deltaFunctionCounts)
135 
136  # Turn this into stars with a narrow width; conserve counts
137  gaussFunction = afwMath.GaussianFunction2D(tGaussianWidth, tGaussianWidth)
138  gaussKernel = afwMath.AnalyticKernel(gaussKernelWidth, gaussKernelWidth, gaussFunction)
139  cim = afwImage.ImageF(tim.getDimensions())
140  afwMath.convolve(cim, tim, gaussKernel, True)
141  tim = cim
142 
143  # Trim off border pixels
144  bbox = gaussKernel.shrinkBBox(tim.getBBox(afwImage.LOCAL))
145  tim = afwImage.ImageF(tim, bbox, afwImage.LOCAL)
146 
147  # Now make a science image which is this convolved with some
148  # spatial function. Use input basis list.
149  polyFunc = afwMath.PolynomialFunction2D(1)
150  kCoeffs = fakeCoeffs()
151  nToUse = min(len(kCoeffs), len(basisList))
152 
153  # Make the full convolved science image
154  sKernel = afwMath.LinearCombinationKernel(afwMath.KernelList(basisList[:nToUse]), polyFunc)
155  sKernel.setSpatialParameters(kCoeffs[:nToUse])
156  sim = afwImage.ImageF(tim.getDimensions())
157  afwMath.convolve(sim, tim, sKernel, True)
158 
159  # Get the good subregion
160  bbox = sKernel.shrinkBBox(sim.getBBox(afwImage.LOCAL))
161 
162  # Add background
163  sim += bgValue
164 
165  # Watch out for negative values
166  tim += 2 * np.abs(np.min(tim.getArray()))
167 
168  # Add noise?
169  if addNoise:
170  sim = makePoissonNoiseImage(sim)
171  tim = makePoissonNoiseImage(tim)
172 
173  # And turn into MaskedImages
174  sim = afwImage.ImageF(sim, bbox, afwImage.LOCAL)
175  svar = afwImage.ImageF(sim, True)
176  smask = afwImage.MaskU(sim.getDimensions())
177  smask.set(0x0)
178  sMi = afwImage.MaskedImageF(sim, smask, svar)
179 
180  tim = afwImage.ImageF(tim, bbox, afwImage.LOCAL)
181  tvar = afwImage.ImageF(tim, True)
182  tmask = afwImage.MaskU(tim.getDimensions())
183  tmask.set(0x0)
184  tMi = afwImage.MaskedImageF(tim, tmask, tvar)
185 
186  if display:
187  import lsst.afw.display.ds9 as ds9
188  ds9.mtv(tMi, frame=1)
189  ds9.mtv(sMi, frame=2)
190 
191  # Finally, make a kernelSet from these 2 images
193  afwGeom.Extent2I(sizeCell * nCell,
194  sizeCell * nCell)),
195  sizeCell,
196  sizeCell)
197  stampHalfWidth = 2 * kSize
198  for x in range(nCell):
199  for y in range(nCell):
200  xCoord = x * sizeCell + sizeCell // 2
201  yCoord = y * sizeCell + sizeCell // 2
202  p0 = afwGeom.Point2I(xCoord - stampHalfWidth,
203  yCoord - stampHalfWidth)
204  p1 = afwGeom.Point2I(xCoord + stampHalfWidth,
205  yCoord + stampHalfWidth)
206  bbox = afwGeom.Box2I(p0, p1)
207  tsi = afwImage.MaskedImageF(tMi, bbox, afwImage.LOCAL)
208  ssi = afwImage.MaskedImageF(sMi, bbox, afwImage.LOCAL)
209 
210  kc = diffimLib.makeKernelCandidate(xCoord, yCoord, tsi, ssi, policyFake)
211  kernelCellSet.insertCandidate(kc)
212 
213  tMi.setXY0(0,0)
214  sMi.setXY0(0,0)
215  return tMi, sMi, sKernel, kernelCellSet, configFake
216 
217 
218 #######
219 # Background subtraction for ip_diffim
220 #######
221 
222 def backgroundSubtract(config, maskedImages):
223  backgrounds = []
224  t0 = time.time()
225  algorithm = config.algorithm
226  binsize = config.binSize
227  undersample = config.undersampleStyle
228  bctrl = afwMath.BackgroundControl(algorithm)
229  bctrl.setUndersampleStyle(undersample)
230  for maskedImage in maskedImages:
231  bctrl.setNxSample(maskedImage.getWidth()//binsize + 1)
232  bctrl.setNySample(maskedImage.getHeight()//binsize + 1)
233  image = maskedImage.getImage()
234  backobj = afwMath.makeBackground(image, bctrl)
235 
236  image -= backobj.getImageF()
237  backgrounds.append(backobj.getImageF())
238  del backobj
239 
240  t1 = time.time()
241  pexLog.Trace("lsst.ip.diffim.backgroundSubtract", 1,
242  "Total time for background subtraction : %.2f s" % (t1-t0))
243  return backgrounds
244 
245 #######
246 # More coarse debugging
247 #######
248 
249 def writeKernelCellSet(kernelCellSet, psfMatchingKernel, backgroundModel, outdir):
250  if not os.path.isdir(outdir):
251  os.makedirs(outdir)
252 
253  for cell in kernelCellSet.getCellList():
254  for cand in cell.begin(False): # False = include bad candidates
255  cand = diffimLib.cast_KernelCandidateF(cand)
256  if cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
257  xCand = int(cand.getXCenter())
258  yCand = int(cand.getYCenter())
259  idCand = cand.getId()
260  diffIm = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
261  kernel = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG)
262  diffIm.writeFits(os.path.join(outdir, 'diffim_c%d_x%d_y%d.fits' % (idCand, xCand, yCand)))
263  kernel.writeFits(os.path.join(outdir, 'kernel_c%d_x%d_y%d.fits' % (idCand, xCand, yCand)))
264 
265  # Diffim from spatial model
266  ski = afwImage.ImageD(kernel.getDimensions())
267  psfMatchingKernel.computeImage(ski, False, xCand, yCand)
268  sk = afwMath.FixedKernel(ski)
269  sbg = backgroundModel(xCand, yCand)
270  sdmi = cand.getDifferenceImage(sk, sbg)
271  sdmi.writeFits(os.path.join(outdir, 'sdiffim_c%d_x%d_y%d.fits' % (idCand, xCand, yCand)))
272 
273 #######
274 # Converting types
275 #######
276 
277 def sourceToFootprintList(candidateInList, templateExposure, scienceExposure, kernelSize, config, log):
278  """ Takes an input list of Sources that were selected to constrain
279  the Psf-matching Kernel and turns them into a List of Footprints,
280  which are used to seed a set of KernelCandidates. The function
281  checks both the template and science image for masked pixels,
282  rejecting the Source if certain Mask bits (defined in config) are
283  set within the Footprint.
284 
285  @param candidateInList: Input list of Sources
286  @param templateExposure: Template image, to be checked for Mask bits in Source Footprint
287  @param scienceExposure: Science image, to be checked for Mask bits in Source Footprint
288  @param config: Config that defines the Mask planes that indicate an invalid Source and Bbox grow radius
289  @param log: Log for output
290 
291  @return a list of dicts having a "source" and "footprint" field, to be used for Psf-matching
292  """
293 
294  candidateOutList = []
295  fsb = diffimLib.FindSetBitsU()
296  badBitMask = 0
297  for mp in config.badMaskPlanes:
298  badBitMask |= afwImage.MaskU.getPlaneBitMask(mp)
299  bbox = scienceExposure.getBBox()
300 
301  # Size to grow Sources
302  if config.scaleByFwhm:
303  fpGrowPix = int(config.fpGrowKernelScaling * kernelSize + 0.5)
304  else:
305  fpGrowPix = config.fpGrowPix
306  log.info("Growing %d kernel candidate stars by %d pixels" % (len(candidateInList), fpGrowPix))
307 
308  for kernelCandidate in candidateInList:
309  if not type(kernelCandidate) == afwTable.SourceRecord:
310  raise RuntimeError, ("Candiate not of type afwTable.SourceRecord")
311  bm1 = 0
312  bm2 = 0
313  center = afwGeom.Point2I(scienceExposure.getWcs().skyToPixel(kernelCandidate.getCoord()))
314  if center[0] < bbox.getMinX() or center[0] > bbox.getMaxX():
315  continue
316  if center[1] < bbox.getMinY() or center[1] > bbox.getMaxY():
317  continue
318 
319  xmin = center[0] - fpGrowPix
320  xmax = center[0] + fpGrowPix
321  ymin = center[1] - fpGrowPix
322  ymax = center[1] + fpGrowPix
323 
324  # Keep object centered
325  if (xmin - bbox.getMinX()) < 0:
326  xmax += (xmin - bbox.getMinX())
327  xmin -= (xmin - bbox.getMinX())
328  if (ymin - bbox.getMinY()) < 0:
329  ymax += (ymin - bbox.getMinY())
330  ymin -= (ymin - bbox.getMinY())
331  if (bbox.getMaxX() - xmax) < 0:
332  xmin -= (bbox.getMaxX() - xmax)
333  xmax += (bbox.getMaxX() - xmax)
334  if (bbox.getMaxY() - ymax) < 0:
335  ymin -= (bbox.getMaxY() - ymax)
336  ymax += (bbox.getMaxY() - ymax)
337  if xmin > xmax or ymin > ymax:
338  continue
339 
340  kbbox = afwGeom.Box2I(afwGeom.Point2I(xmin, ymin), afwGeom.Point2I(xmax, ymax))
341  try:
342  fsb.apply(afwImage.MaskedImageF(templateExposure.getMaskedImage(), kbbox, False).getMask())
343  bm1 = fsb.getBits()
344  fsb.apply(afwImage.MaskedImageF(scienceExposure.getMaskedImage(), kbbox, False).getMask())
345  bm2 = fsb.getBits()
346  except Exception:
347  pass
348  else:
349  if not((bm1 & badBitMask) or (bm2 & badBitMask)):
350  candidateOutList.append({'source':kernelCandidate, 'footprint':afwDetect.Footprint(kbbox)})
351  log.info("Selected %d / %d sources for KernelCandidacy" % (len(candidateOutList), len(candidateInList)))
352  return candidateOutList
353 
354 def sourceTableToCandidateList(sourceTable, templateExposure, scienceExposure, kConfig, dConfig, log,
355  basisList, doBuild=False):
356  """Takes an input list of Sources, and turns them into
357  KernelCandidates for fitting of the Psf-matching kernel."""
358  kernelSize = basisList[0].getWidth()
359  footprintList = sourceToFootprintList(list(sourceTable), templateExposure, scienceExposure,
360  kernelSize, dConfig, log)
361  candList = []
362 
363  if doBuild and not basisList:
364  doBuild = False
365  else:
366  policy = pexConfig.makePolicy(kConfig)
367  visitor = diffimLib.BuildSingleKernelVisitorF(basisList, policy)
368 
369  policy = pexConfig.makePolicy(kConfig)
370  for cand in footprintList:
371  bbox = cand['footprint'].getBBox()
372  tmi = afwImage.MaskedImageF(templateExposure.getMaskedImage(), bbox)
373  smi = afwImage.MaskedImageF(scienceExposure.getMaskedImage(), bbox)
374  kCand = diffimLib.makeKernelCandidate(cand['source'], tmi, smi, policy)
375  if doBuild:
376  visitor.processCandidate(kCand)
377  kCand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN)
378  candList.append(kCand)
379  return candList
380 
381 
382 #######
383 #
384 #######
385 
386 
387 class NbasisEvaluator(object):
388  """A functor to evaluate the Bayesian Information Criterion for the number of basis sets going into the kernel fitting"""
389  def __init__(self, psfMatchConfig, psfFwhmPixTc, psfFwhmPixTnc):
390  self.psfMatchConfig = psfMatchConfig
391  self.psfFwhmPixTc = psfFwhmPixTc
392  self.psfFwhmPixTnc = psfFwhmPixTnc
393  if not self.psfMatchConfig.kernelBasisSet == "alard-lupton":
394  raise RuntimeError, "BIC only implemnted for AL (alard lupton) basis"
395 
396  def __call__(self, kernelCellSet, log):
397  d1, d2, d3 = self.psfMatchConfig.alardDegGauss
398  bicArray = {}
399  for d1i in range(1, d1+1):
400  for d2i in range(1, d2+1):
401  for d3i in range(1, d3+1):
402  dList = [d1i, d2i, d3i]
403  bicConfig = type(self.psfMatchConfig)(self.psfMatchConfig, alardDegGauss=dList)
404  kList = makeKernelBasisList(bicConfig, self.psfFwhmPixTc, self.psfFwhmPixTnc)
405  k = len(kList)
406  visitor = diffimLib.BuildSingleKernelVisitorF(kList, pexConfig.makePolicy(bicConfig))
407  visitor.setSkipBuilt(False)
408  kernelCellSet.visitCandidates(visitor, bicConfig.nStarPerCell)
409 
410  for cell in kernelCellSet.getCellList():
411  for cand in cell.begin(False): # False = include bad candidates
412  cand = diffimLib.cast_KernelCandidateF(cand)
413  if cand.getStatus() != afwMath.SpatialCellCandidate.GOOD:
414  continue
415  diffIm = cand.getDifferenceImage(diffimLib.KernelCandidateF.RECENT)
416  bbox = cand.getKernel(diffimLib.KernelCandidateF.RECENT).shrinkBBox(diffIm.getBBox(afwImage.LOCAL))
417  diffIm = type(diffIm)(diffIm, bbox, True)
418  chi2 = diffIm.getImage().getArray()**2 / diffIm.getVariance().getArray()
419  n = chi2.shape[0] * chi2.shape[1]
420  bic = np.sum(chi2) + k * np.log(n)
421  if not bicArray.has_key(cand.getId()):
422  bicArray[cand.getId()] = {}
423  bicArray[cand.getId()][(d1i, d2i, d3i)] = bic
424 
425  bestConfigs = []
426  candIds = bicArray.keys()
427  for candId in candIds:
428  cconfig, cvals = bicArray[candId].keys(), bicArray[candId].values()
429  idx = np.argsort(cvals)
430  bestConfig = cconfig[idx[0]]
431  bestConfigs.append(bestConfig)
432 
433  counter = Counter(bestConfigs).most_common(3)
434  log.info("B.I.C. prefers basis complexity %s %d times; %s %d times; %s %d times" % (counter[0][0], counter[0][1],
435  counter[1][0], counter[1][1],
436  counter[2][0], counter[2][1]))
437  return counter[0][0], counter[1][0], counter[2][0]
boost::shared_ptr< Background > makeBackground(ImageT const &img, BackgroundControl const &bgCtrl)
A convenience function that uses function overloading to make the correct type of Background...
Definition: Background.h:467
def backgroundSubtract
Background subtraction for ip_diffim.
Definition: diffimTools.py:222
Configuration for image-to-image Psf matching.
limited backward compatibility to the DC2 run-time trace facilities
Definition: Trace.h:93
Pass parameters to a Background object.
Definition: Background.h:61
An integer coordinate rectangle.
Definition: Box.h:53
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:378
def writeKernelCellSet
More coarse debugging.
Definition: diffimTools.py:249
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:814
void randomGaussianImage(ImageT *image, Random &rand)
Definition: RandomImage.cc:155
A set of pixels in an Image.
Definition: Footprint.h:62
def fakeCoeffs
Make fake images for testing; one is a delta function (or narrow gaussian) and the other is a convolu...
Definition: diffimTools.py:82
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
def makeFlatNoiseImage
Add noise.
Definition: diffimTools.py:49
Record class that contains measurements made on a single exposure.
Definition: Source.h:81
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge=false)
Old, deprecated version of convolve.
std::vector< boost::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:542
A kernel described by a function.
Definition: Kernel.h:628
def sourceToFootprintList
Converting types.
Definition: diffimTools.py:277
A kernel created from an Image.
Definition: Kernel.h:551