LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
diffimTools.py
Go to the documentation of this file.
1 # This file is part of ip_diffim.
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__ = ["backgroundSubtract", "writeKernelCellSet", "sourceToFootprintList", "NbasisEvaluator"]
23 
24 # python
25 import time
26 import os
27 from collections import Counter
28 import numpy as np
29 
30 # all the c++ level classes and routines
31 from . import diffimLib
32 
33 # all the other LSST packages
34 import lsst.afw.geom as afwGeom
35 import lsst.afw.image as afwImage
36 import lsst.afw.table as afwTable
37 import lsst.afw.detection as afwDetect
38 import lsst.afw.math as afwMath
39 import lsst.geom as geom
40 from lsst.log import Log
41 import lsst.pex.config as pexConfig
42 from .makeKernelBasisList import makeKernelBasisList
43 
44 # Helper functions for ipDiffim; mostly viewing of results and writing
45 # debugging info to disk.
46 
47 
50 
51 
52 def makeFlatNoiseImage(mi, seedStat=afwMath.MAX):
53  img = mi.getImage()
54  seed = int(10.*afwMath.makeStatistics(mi.getImage(), seedStat).getValue() + 1)
55  rdm = afwMath.Random(afwMath.Random.MT19937, seed)
56  rdmImage = img.Factory(img.getDimensions())
57  afwMath.randomGaussianImage(rdmImage, rdm)
58  return rdmImage
59 
60 
62  """Return a Poisson noise image based on im
63 
64  Parameters
65  ----------
66  im : `lsst.afw.image.Image`
67  image; the output image has the same dtype, dimensions, and shape
68  and its expectation value is the value of ``im`` at each pixel
69 
70  Returns
71  -------
72  noiseIm : `lsst.afw.image.Image`
73  Newly constructed image instance, same type as ``im``.
74 
75  Notes
76  -----
77  - Warning: This uses an undocumented numpy API (the documented API
78  uses a single float expectation value instead of an array).
79 
80  - Uses numpy.random; you may wish to call numpy.random.seed first.
81  """
82  import numpy.random as rand
83  imArr = im.getArray()
84  noiseIm = im.Factory(im.getBBox())
85  noiseArr = noiseIm.getArray()
86 
87  intNoiseArr = rand.poisson(np.where(np.isfinite(imArr), imArr, 0.0))
88 
89  noiseArr[:, :] = intNoiseArr.astype(noiseArr.dtype)
90  return noiseIm
91 
92 
97 
98 
99 def fakeCoeffs():
100  kCoeffs = ((1.0, 0.0, 0.0),
101  (0.005, -0.000001, 0.000001),
102  (0.005, 0.000004, 0.000004),
103  (-0.001, -0.000030, 0.000030),
104  (-0.001, 0.000015, 0.000015),
105  (-0.005, -0.000050, 0.000050))
106  return kCoeffs
107 
108 
109 def makeFakeKernelSet(sizeCell=128, nCell=3,
110  deltaFunctionCounts=1.e4, tGaussianWidth=1.0,
111  addNoise=True, bgValue=100., display=False):
112  """Generate test template and science images with sources.
113 
114  Parameters
115  ----------
116  sizeCell : `int`, optional
117  Size of the square spatial cells in pixels.
118  nCell : `int`, optional
119  Number of adjacent spatial cells in both direction in both images.
120  deltaFunctionCounts : `float`, optional
121  Flux value for the template image sources.
122  tGaussianWidth : `float`, optional
123  Sigma of the generated Gaussian PSF sources in the template image.
124  addNoise : `bool`, optional
125  If `True`, Poisson noise is added to both the generated template
126  and science images.
127  bgValue : `float`, optional
128  Background level to be added to the generated science image.
129  display : `bool`, optional
130  If `True` displays the generated template and science images by
131  `lsst.afw.display.Display`.
132 
133  Notes
134  -----
135  - The generated images consist of adjacent ``nCell x nCell`` cells, each
136  of pixel size ``sizeCell x sizeCell``.
137  - The sources in the science image are generated by convolving the
138  template by ``sKernel``. ``sKernel`` is a spatial `LinearCombinationKernel`
139  of hard wired kernel bases functions. The linear combination has first
140  order polynomial spatial dependence with polynomial parameters from ``fakeCoeffs()``.
141  - The template image sources are generated in the center of each spatial
142  cell from one pixel, set to `deltaFunctionCounts` counts, then convolved
143  by a 2D Gaussian with sigma of `tGaussianWidth` along each axis.
144  - The sources are also returned in ``kernelCellSet`` each source is "detected"
145  exactly at the center of a cell.
146 
147  Returns
148  -------
149  tMi : `lsst.afw.image.MaskedImage`
150  Generated template image.
151  sMi : `lsst.afw.image.MaskedImage`
152  Generated science image.
153  sKernel : `lsst.afw.math.LinearCombinationKernel`
154  The spatial kernel used to generate the sources in the science image.
155  kernelCellSet : `lsst.afw.math.SpatialCellSet`
156  Cell grid of `lsst.afw.math.SpatialCell` instances, containing
157  `lsst.ip.diffim.KernelCandidate` instances around all the generated sources
158  in the science image.
159  configFake : `lsst.ip.diffim.ImagePsfMatchConfig`
160  Config instance used in the image generation.
161  """
162  from . import imagePsfMatch
163  configFake = imagePsfMatch.ImagePsfMatchConfig()
164  configFake.kernel.name = "AL"
165  subconfigFake = configFake.kernel.active
166  subconfigFake.alardNGauss = 1
167  subconfigFake.alardSigGauss = [2.5, ]
168  subconfigFake.alardDegGauss = [2, ]
169  subconfigFake.sizeCellX = sizeCell
170  subconfigFake.sizeCellY = sizeCell
171  subconfigFake.spatialKernelOrder = 1
172  subconfigFake.spatialModelType = "polynomial"
173  subconfigFake.singleKernelClipping = False # variance is a hack
174  subconfigFake.spatialKernelClipping = False # variance is a hack
175  if bgValue > 0.0:
176  subconfigFake.fitForBackground = True
177 
178  psFake = pexConfig.makePropertySet(subconfigFake)
179 
180  basisList = makeKernelBasisList(subconfigFake)
181  kSize = subconfigFake.kernelSize
182 
183  # This sets the final extent of each convolved delta function
184  gaussKernelWidth = sizeCell//2
185 
186  # This sets the scale over which pixels are correlated in the
187  # spatial convolution; should be at least as big as the kernel you
188  # are trying to fit for
189  spatialKernelWidth = kSize
190 
191  # Number of bad pixels due to convolutions
192  border = (gaussKernelWidth + spatialKernelWidth)//2
193 
194  # Make a fake image with a matrix of delta functions
195  totalSize = nCell*sizeCell + 2*border
196  tim = afwImage.ImageF(geom.Extent2I(totalSize, totalSize))
197  for x in range(nCell):
198  for y in range(nCell):
199  tim[x*sizeCell + sizeCell//2 + border - 1,
200  y*sizeCell + sizeCell//2 + border - 1,
201  afwImage.LOCAL] = deltaFunctionCounts
202 
203  # Turn this into stars with a narrow width; conserve counts
204  gaussFunction = afwMath.GaussianFunction2D(tGaussianWidth, tGaussianWidth)
205  gaussKernel = afwMath.AnalyticKernel(gaussKernelWidth, gaussKernelWidth, gaussFunction)
206  cim = afwImage.ImageF(tim.getDimensions())
207  convolutionControl = afwMath.ConvolutionControl()
208  convolutionControl.setDoNormalize(True)
209  afwMath.convolve(cim, tim, gaussKernel, convolutionControl)
210  tim = cim
211 
212  # Trim off border pixels
213  bbox = gaussKernel.shrinkBBox(tim.getBBox(afwImage.LOCAL))
214  tim = afwImage.ImageF(tim, bbox, afwImage.LOCAL)
215 
216  # Now make a science image which is this convolved with some
217  # spatial function. Use input basis list.
218  polyFunc = afwMath.PolynomialFunction2D(1)
219  kCoeffs = fakeCoeffs()
220  nToUse = min(len(kCoeffs), len(basisList))
221 
222  # Make the full convolved science image
223  sKernel = afwMath.LinearCombinationKernel(basisList[:nToUse], polyFunc)
224  sKernel.setSpatialParameters(kCoeffs[:nToUse])
225  sim = afwImage.ImageF(tim.getDimensions())
226  convolutionControl = afwMath.ConvolutionControl()
227  convolutionControl.setDoNormalize(True)
228  afwMath.convolve(sim, tim, sKernel, convolutionControl)
229 
230  # Get the good subregion
231  bbox = sKernel.shrinkBBox(sim.getBBox(afwImage.LOCAL))
232 
233  # Add background
234  sim += bgValue
235 
236  # Watch out for negative values
237  tim += 2*np.abs(np.min(tim.getArray()))
238 
239  # Add noise?
240  if addNoise:
241  sim = makePoissonNoiseImage(sim)
242  tim = makePoissonNoiseImage(tim)
243 
244  # And turn into MaskedImages
245  sim = afwImage.ImageF(sim, bbox, afwImage.LOCAL)
246  svar = afwImage.ImageF(sim, True)
247  smask = afwImage.Mask(sim.getDimensions())
248  smask.set(0x0)
249  sMi = afwImage.MaskedImageF(sim, smask, svar)
250 
251  tim = afwImage.ImageF(tim, bbox, afwImage.LOCAL)
252  tvar = afwImage.ImageF(tim, True)
253  tmask = afwImage.Mask(tim.getDimensions())
254  tmask.set(0x0)
255  tMi = afwImage.MaskedImageF(tim, tmask, tvar)
256 
257  if display:
258  import lsst.afw.display as afwDisplay
259  afwDisplay.Display(frame=1).mtv(tMi)
260  afwDisplay.Display(frame=2).mtv(sMi)
261 
262  # Finally, make a kernelSet from these 2 images
263  kernelCellSet = afwMath.SpatialCellSet(geom.Box2I(geom.Point2I(0, 0),
264  geom.Extent2I(sizeCell*nCell,
265  sizeCell*nCell)),
266  sizeCell,
267  sizeCell)
268  stampHalfWidth = 2*kSize
269  for x in range(nCell):
270  for y in range(nCell):
271  xCoord = x*sizeCell + sizeCell//2
272  yCoord = y*sizeCell + sizeCell//2
273  p0 = geom.Point2I(xCoord - stampHalfWidth,
274  yCoord - stampHalfWidth)
275  p1 = geom.Point2I(xCoord + stampHalfWidth,
276  yCoord + stampHalfWidth)
277  bbox = geom.Box2I(p0, p1)
278  tsi = afwImage.MaskedImageF(tMi, bbox, origin=afwImage.LOCAL)
279  ssi = afwImage.MaskedImageF(sMi, bbox, origin=afwImage.LOCAL)
280 
281  kc = diffimLib.makeKernelCandidate(xCoord, yCoord, tsi, ssi, psFake)
282  kernelCellSet.insertCandidate(kc)
283 
284  tMi.setXY0(0, 0)
285  sMi.setXY0(0, 0)
286  return tMi, sMi, sKernel, kernelCellSet, configFake
287 
288 
289 
292 
293 def backgroundSubtract(config, maskedImages):
294  """Subtract the background from masked images.
295 
296  Parameters
297  ----------
298  config : TODO: DM-17458
299  TODO: DM-17458
300  maskedImages : `list` of `lsst.afw.image.MaskedImage`
301  TODO: DM-17458
302 
303  Returns
304  -------
305  TODO: DM-17458
306  TODO: DM-17458
307  """
308  backgrounds = []
309  t0 = time.time()
310  algorithm = config.algorithm
311  binsize = config.binSize
312  undersample = config.undersampleStyle
313  bctrl = afwMath.BackgroundControl(algorithm)
314  bctrl.setUndersampleStyle(undersample)
315  for maskedImage in maskedImages:
316  bctrl.setNxSample(maskedImage.getWidth()//binsize + 1)
317  bctrl.setNySample(maskedImage.getHeight()//binsize + 1)
318  image = maskedImage.getImage()
319  backobj = afwMath.makeBackground(image, bctrl)
320 
321  image -= backobj.getImageF()
322  backgrounds.append(backobj.getImageF())
323  del backobj
324 
325  t1 = time.time()
326  logger = Log.getLogger("ip.diffim.backgroundSubtract")
327  logger.debug("Total time for background subtraction : %.2f s", (t1 - t0))
328  return backgrounds
329 
330 
333 
334 
335 def writeKernelCellSet(kernelCellSet, psfMatchingKernel, backgroundModel, outdir):
336  """TODO: DM-17458
337 
338  Parameters
339  ----------
340  kernelCellSet : TODO: DM-17458
341  TODO: DM-17458
342  psfMatchingKernel : TODO: DM-17458
343  TODO: DM-17458
344  backgroundModel : TODO: DM-17458
345  TODO: DM-17458
346  outdir : TODO: DM-17458
347  TODO: DM-17458
348  """
349  if not os.path.isdir(outdir):
350  os.makedirs(outdir)
351 
352  for cell in kernelCellSet.getCellList():
353  for cand in cell.begin(False): # False = include bad candidates
354  if cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
355  xCand = int(cand.getXCenter())
356  yCand = int(cand.getYCenter())
357  idCand = cand.getId()
358  diffIm = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
359  kernel = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG)
360  diffIm.writeFits(os.path.join(outdir, 'diffim_c%d_x%d_y%d.fits' % (idCand, xCand, yCand)))
361  kernel.writeFits(os.path.join(outdir, 'kernel_c%d_x%d_y%d.fits' % (idCand, xCand, yCand)))
362 
363  # Diffim from spatial model
364  ski = afwImage.ImageD(kernel.getDimensions())
365  psfMatchingKernel.computeImage(ski, False, xCand, yCand)
366  sk = afwMath.FixedKernel(ski)
367  sbg = backgroundModel(xCand, yCand)
368  sdmi = cand.getDifferenceImage(sk, sbg)
369  sdmi.writeFits(os.path.join(outdir, 'sdiffim_c%d_x%d_y%d.fits' % (idCand, xCand, yCand)))
370 
371 
374 
375 
376 def sourceToFootprintList(candidateInList, templateExposure, scienceExposure, kernelSize, config, log):
377  """Convert a list of sources for the PSF-matching Kernel to Footprints.
378 
379  Parameters
380  ----------
381  candidateInList : TODO: DM-17458
382  Input list of Sources
383  templateExposure : TODO: DM-17458
384  Template image, to be checked for Mask bits in Source Footprint
385  scienceExposure : TODO: DM-17458
386  Science image, to be checked for Mask bits in Source Footprint
387  kernelSize : TODO: DM-17458
388  TODO: DM-17458
389  config : TODO: DM-17458
390  Config that defines the Mask planes that indicate an invalid Source and Bbox grow radius
391  log : TODO: DM-17458
392  Log for output
393 
394  Returns
395  -------
396  candidateOutList : `list`
397  a list of dicts having a "source" and "footprint" field, to be used for Psf-matching
398 
399  Raises
400  ------
401  RuntimeError
402  TODO: DM-17458
403 
404  Notes
405  -----
406  Takes an input list of Sources that were selected to constrain
407  the Psf-matching Kernel and turns them into a List of Footprints,
408  which are used to seed a set of KernelCandidates. The function
409  checks both the template and science image for masked pixels,
410  rejecting the Source if certain Mask bits (defined in config) are
411  set within the Footprint.
412  """
413 
414  candidateOutList = []
415  fsb = diffimLib.FindSetBitsU()
416  badBitMask = 0
417  for mp in config.badMaskPlanes:
418  badBitMask |= afwImage.Mask.getPlaneBitMask(mp)
419  bbox = scienceExposure.getBBox()
420 
421  # Size to grow Sources
422  if config.scaleByFwhm:
423  fpGrowPix = int(config.fpGrowKernelScaling*kernelSize + 0.5)
424  else:
425  fpGrowPix = config.fpGrowPix
426  log.info("Growing %d kernel candidate stars by %d pixels", len(candidateInList), fpGrowPix)
427 
428  for kernelCandidate in candidateInList:
429  if not type(kernelCandidate) == afwTable.SourceRecord:
430  raise RuntimeError("Candiate not of type afwTable.SourceRecord")
431  bm1 = 0
432  bm2 = 0
433  center = geom.Point2I(scienceExposure.getWcs().skyToPixel(kernelCandidate.getCoord()))
434  if center[0] < bbox.getMinX() or center[0] > bbox.getMaxX():
435  continue
436  if center[1] < bbox.getMinY() or center[1] > bbox.getMaxY():
437  continue
438 
439  xmin = center[0] - fpGrowPix
440  xmax = center[0] + fpGrowPix
441  ymin = center[1] - fpGrowPix
442  ymax = center[1] + fpGrowPix
443 
444  # Keep object centered
445  if (xmin - bbox.getMinX()) < 0:
446  xmax += (xmin - bbox.getMinX())
447  xmin -= (xmin - bbox.getMinX())
448  if (ymin - bbox.getMinY()) < 0:
449  ymax += (ymin - bbox.getMinY())
450  ymin -= (ymin - bbox.getMinY())
451  if (bbox.getMaxX() - xmax) < 0:
452  xmin -= (bbox.getMaxX() - xmax)
453  xmax += (bbox.getMaxX() - xmax)
454  if (bbox.getMaxY() - ymax) < 0:
455  ymin -= (bbox.getMaxY() - ymax)
456  ymax += (bbox.getMaxY() - ymax)
457  if xmin > xmax or ymin > ymax:
458  continue
459 
460  kbbox = geom.Box2I(geom.Point2I(xmin, ymin), geom.Point2I(xmax, ymax))
461  try:
462  fsb.apply(afwImage.MaskedImageF(templateExposure.getMaskedImage(), kbbox, deep=False).getMask())
463  bm1 = fsb.getBits()
464  fsb.apply(afwImage.MaskedImageF(scienceExposure.getMaskedImage(), kbbox, deep=False).getMask())
465  bm2 = fsb.getBits()
466  except Exception:
467  pass
468  else:
469  if not((bm1 & badBitMask) or (bm2 & badBitMask)):
470  candidateOutList.append({'source': kernelCandidate,
471  'footprint': afwDetect.Footprint(afwGeom.SpanSet(kbbox))})
472  log.info("Selected %d / %d sources for KernelCandidacy", len(candidateOutList), len(candidateInList))
473  return candidateOutList
474 
475 
476 def sourceTableToCandidateList(sourceTable, templateExposure, scienceExposure, kConfig, dConfig, log,
477  basisList, doBuild=False):
478  """Convert a list of Sources into KernelCandidates.
479 
480  The KernelCandidates are used for fitting the Psf-matching kernel.
481 
482  Parameters
483  ----------
484  sourceTable : TODO: DM-17458
485  TODO: DM-17458
486  templateExposure : TODO: DM-17458
487  TODO: DM-17458
488  scienceExposure : TODO: DM-17458
489  TODO: DM-17458
490  kConfig : TODO: DM-17458
491  TODO: DM-17458
492  dConfig : TODO: DM-17458
493  TODO: DM-17458
494  log : TODO: DM-17458
495  TODO: DM-17458
496  basisList : TODO: DM-17458
497  TODO: DM-17458
498  doBuild : `bool`, optional
499  TODO: DM-17458
500 
501  Returns
502  -------
503  TODO: DM-17458
504  TODO: DM-17458
505  """
506  kernelSize = basisList[0].getWidth()
507  footprintList = sourceToFootprintList(list(sourceTable), templateExposure, scienceExposure,
508  kernelSize, dConfig, log)
509  candList = []
510 
511  if doBuild and not basisList:
512  doBuild = False
513  else:
514  ps = pexConfig.makePropertySet(kConfig)
515  visitor = diffimLib.BuildSingleKernelVisitorF(basisList, ps)
516 
517  ps = pexConfig.makePropertySet(kConfig)
518  for cand in footprintList:
519  bbox = cand['footprint'].getBBox()
520  tmi = afwImage.MaskedImageF(templateExposure.getMaskedImage(), bbox)
521  smi = afwImage.MaskedImageF(scienceExposure.getMaskedImage(), bbox)
522  kCand = diffimLib.makeKernelCandidate(cand['source'], tmi, smi, ps)
523  if doBuild:
524  visitor.processCandidate(kCand)
525  kCand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN)
526  candList.append(kCand)
527  return candList
528 
529 
530 
533 
534 
536  """A functor to evaluate the Bayesian Information Criterion for the number of basis sets
537  going into the kernel fitting"""
538 
539  def __init__(self, psfMatchConfig, psfFwhmPixTc, psfFwhmPixTnc):
540  self.psfMatchConfigpsfMatchConfig = psfMatchConfig
541  self.psfFwhmPixTcpsfFwhmPixTc = psfFwhmPixTc
542  self.psfFwhmPixTncpsfFwhmPixTnc = psfFwhmPixTnc
543  if not self.psfMatchConfigpsfMatchConfig.kernelBasisSet == "alard-lupton":
544  raise RuntimeError("BIC only implemnted for AL (alard lupton) basis")
545 
546  def __call__(self, kernelCellSet, log):
547  d1, d2, d3 = self.psfMatchConfigpsfMatchConfig.alardDegGauss
548  bicArray = {}
549  for d1i in range(1, d1 + 1):
550  for d2i in range(1, d2 + 1):
551  for d3i in range(1, d3 + 1):
552  dList = [d1i, d2i, d3i]
553  bicConfig = type(self.psfMatchConfigpsfMatchConfig)(self.psfMatchConfigpsfMatchConfig, alardDegGauss=dList)
554  kList = makeKernelBasisList(bicConfig, self.psfFwhmPixTcpsfFwhmPixTc, self.psfFwhmPixTncpsfFwhmPixTnc)
555  k = len(kList)
556  visitor = diffimLib.BuildSingleKernelVisitorF(kList,
557  pexConfig.makePropertySet(bicConfig))
558  visitor.setSkipBuilt(False)
559  kernelCellSet.visitCandidates(visitor, bicConfig.nStarPerCell)
560 
561  for cell in kernelCellSet.getCellList():
562  for cand in cell.begin(False): # False = include bad candidates
563  if cand.getStatus() != afwMath.SpatialCellCandidate.GOOD:
564  continue
565  diffIm = cand.getDifferenceImage(diffimLib.KernelCandidateF.RECENT)
566  bbox = cand.getKernel(diffimLib.KernelCandidateF.RECENT).shrinkBBox(
567  diffIm.getBBox(afwImage.LOCAL))
568  diffIm = type(diffIm)(diffIm, bbox, True)
569  chi2 = diffIm.getImage().getArray()**2/diffIm.getVariance().getArray()
570  n = chi2.shape[0]*chi2.shape[1]
571  bic = np.sum(chi2) + k*np.log(n)
572  if cand.getId() not in bicArray:
573  bicArray[cand.getId()] = {}
574  bicArray[cand.getId()][(d1i, d2i, d3i)] = bic
575 
576  bestConfigs = []
577  for candId in bicArray:
578  cconfig, cvals = list(bicArray[candId].keys()), list(bicArray[candId].values())
579  idx = np.argsort(cvals)
580  bestConfig = cconfig[idx[0]]
581  bestConfigs.append(bestConfig)
582 
583  counter = Counter(bestConfigs).most_common(3)
584  log.info("B.I.C. prefers basis complexity %s %d times; %s %d times; %s %d times",
585  counter[0][0], counter[0][1],
586  counter[1][0], counter[1][1],
587  counter[2][0], counter[2][1])
588  return counter[0][0], counter[1][0], counter[2][0]
int min
table::Key< int > type
Definition: Detector.cc:163
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
A compact representation of a collection of pixels.
Definition: SpanSet.h:78
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
A kernel described by a function.
Definition: Kernel.h:535
Pass parameters to a Background object.
Definition: Background.h:56
Parameters to control convolution.
Definition: ConvolveImage.h:50
A kernel created from an Image.
Definition: Kernel.h:471
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:704
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition: Random.h:57
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:383
Record class that contains measurements made on a single exposure.
Definition: Source.h:78
An integer coordinate rectangle.
Definition: Box.h:55
def __init__(self, psfMatchConfig, psfFwhmPixTc, psfFwhmPixTnc)
Definition: diffimTools.py:539
def __call__(self, kernelCellSet, log)
Definition: diffimTools.py:546
daf::base::PropertyList * list
Definition: fits.cc:913
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Definition: ds9.py:92
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
Definition: RandomImage.cc: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:359
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
std::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:526
def makeFlatNoiseImage(mi, seedStat=afwMath.MAX)
Add noise.
Definition: diffimTools.py:52
def makeFakeKernelSet(sizeCell=128, nCell=3, deltaFunctionCounts=1.e4, tGaussianWidth=1.0, addNoise=True, bgValue=100., display=False)
Definition: diffimTools.py:111
def fakeCoeffs()
Make fake images for testing; one is a delta function (or narrow gaussian) and the other is a convolu...
Definition: diffimTools.py:99
def writeKernelCellSet(kernelCellSet, psfMatchingKernel, backgroundModel, outdir)
More coarse debugging.
Definition: diffimTools.py:335
def sourceTableToCandidateList(sourceTable, templateExposure, scienceExposure, kConfig, dConfig, log, basisList, doBuild=False)
Definition: diffimTools.py:477
def backgroundSubtract(config, maskedImages)
Background subtraction for ip_diffim.
Definition: diffimTools.py:293
def sourceToFootprintList(candidateInList, templateExposure, scienceExposure, kernelSize, config, log)
Converting types.
Definition: diffimTools.py:376
def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, basisSigmaGauss=None, metadata=None)
Definition: Log.h:717