22__all__ = [
"backgroundSubtract",
"writeKernelCellSet",
"sourceToFootprintList",
"NbasisEvaluator"]
27from collections
import Counter
31from .
import diffimLib
41from lsst.utils.logging
import getLogger
42from .makeKernelBasisList
import makeKernelBasisList
56 rdmImage = img.Factory(img.getDimensions())
62 """Return a Poisson noise image based on im
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
73 Newly constructed image instance, same type
as ``im``.
77 - Warning: This uses an undocumented numpy API (the documented API
78 uses a single float expectation value instead of an array).
80 - Uses numpy.random; you may wish to call numpy.random.seed first.
82 import numpy.random
as rand
84 noiseIm = im.Factory(im.getBBox())
85 noiseArr = noiseIm.getArray()
87 intNoiseArr = rand.poisson(np.where(np.isfinite(imArr), imArr, 0.0))
89 noiseArr[:, :] = intNoiseArr.astype(noiseArr.dtype)
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))
110 deltaFunctionCounts=1.e4, tGaussianWidth=1.0,
111 addNoise=True, bgValue=100., display=False):
112 """Generate test template and science images with sources.
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
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
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.
150 Generated template image.
152 Generated science image.
154 The spatial kernel used to generate the sources
in the science image.
158 in the science image.
160 Config instance used
in the image generation.
162 from .
import imagePsfMatch
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
174 subconfigFake.spatialKernelClipping =
False
176 subconfigFake.fitForBackground =
True
178 psFake = pexConfig.makePropertySet(subconfigFake)
181 kSize = subconfigFake.kernelSize
184 gaussKernelWidth = sizeCell//2
189 spatialKernelWidth = kSize
192 border = (gaussKernelWidth + spatialKernelWidth)//2
195 totalSize = nCell*sizeCell + 2*border
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
204 gaussFunction = afwMath.GaussianFunction2D(tGaussianWidth, tGaussianWidth)
206 cim = afwImage.ImageF(tim.getDimensions())
208 convolutionControl.setDoNormalize(
True)
213 bbox = gaussKernel.shrinkBBox(tim.getBBox(afwImage.LOCAL))
214 tim = afwImage.ImageF(tim, bbox, afwImage.LOCAL)
218 polyFunc = afwMath.PolynomialFunction2D(1)
220 nToUse =
min(len(kCoeffs), len(basisList))
224 sKernel.setSpatialParameters(kCoeffs[:nToUse])
225 sim = afwImage.ImageF(tim.getDimensions())
227 convolutionControl.setDoNormalize(
True)
231 bbox = sKernel.shrinkBBox(sim.getBBox(afwImage.LOCAL))
237 tim += 2*np.abs(np.min(tim.getArray()))
245 sim = afwImage.ImageF(sim, bbox, afwImage.LOCAL)
246 svar = afwImage.ImageF(sim,
True)
249 sMi = afwImage.MaskedImageF(sim, smask, svar)
251 tim = afwImage.ImageF(tim, bbox, afwImage.LOCAL)
252 tvar = afwImage.ImageF(tim,
True)
255 tMi = afwImage.MaskedImageF(tim, tmask, tvar)
259 afwDisplay.Display(frame=1).
mtv(tMi)
260 afwDisplay.Display(frame=2).
mtv(sMi)
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
274 yCoord - stampHalfWidth)
276 yCoord + stampHalfWidth)
278 tsi = afwImage.MaskedImageF(tMi, bbox, origin=afwImage.LOCAL)
279 ssi = afwImage.MaskedImageF(sMi, bbox, origin=afwImage.LOCAL)
281 kc = diffimLib.makeKernelCandidate(xCoord, yCoord, tsi, ssi, psFake)
282 kernelCellSet.insertCandidate(kc)
286 return tMi, sMi, sKernel, kernelCellSet, configFake
294 """Subtract the background from masked images.
298 config : TODO: DM-17458
310 algorithm = config.algorithm
311 binsize = config.binSize
312 undersample = config.undersampleStyle
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()
321 image -= backobj.getImageF()
322 backgrounds.append(backobj.getImageF())
326 logger =
getLogger(
"lsst.ip.diffim.backgroundSubtract")
327 logger.debug(
"Total time for background subtraction : %.2f s", (t1 - t0))
340 kernelCellSet : TODO: DM-17458
342 psfMatchingKernel : TODO: DM-17458
344 backgroundModel : TODO: DM-17458
346 outdir : TODO: DM-17458
349 if not os.path.isdir(outdir):
352 for cell
in kernelCellSet.getCellList():
353 for cand
in cell.begin(
False):
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)))
364 ski = afwImage.ImageD(kernel.getDimensions())
365 psfMatchingKernel.computeImage(ski,
False, xCand, yCand)
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)))
377 """Convert a list of sources for the PSF-matching Kernel to Footprints.
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
389 config : TODO: DM-17458
390 Config that defines the Mask planes that indicate an invalid Source
and Bbox grow radius
396 candidateOutList : `list`
397 a list of dicts having a
"source" and "footprint" field, to be used
for Psf-matching
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.
414 candidateOutList = []
415 fsb = diffimLib.FindSetBitsU()
417 for mp
in config.badMaskPlanes:
418 badBitMask |= afwImage.Mask.getPlaneBitMask(mp)
419 bbox = scienceExposure.getBBox()
422 if config.scaleByFwhm:
423 fpGrowPix =
int(config.fpGrowKernelScaling*kernelSize + 0.5)
425 fpGrowPix = config.fpGrowPix
426 log.info(
"Growing %d kernel candidate stars by %d pixels", len(candidateInList), fpGrowPix)
428 for kernelCandidate
in candidateInList:
430 raise RuntimeError(
"Candiate not of type afwTable.SourceRecord")
433 center =
geom.Point2I(scienceExposure.getWcs().skyToPixel(kernelCandidate.getCoord()))
434 if center[0] < bbox.getMinX()
or center[0] > bbox.getMaxX():
436 if center[1] < bbox.getMinY()
or center[1] > bbox.getMaxY():
439 xmin = center[0] - fpGrowPix
440 xmax = center[0] + fpGrowPix
441 ymin = center[1] - fpGrowPix
442 ymax = center[1] + fpGrowPix
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:
462 fsb.apply(afwImage.MaskedImageF(templateExposure.getMaskedImage(), kbbox, deep=
False).getMask())
464 fsb.apply(afwImage.MaskedImageF(scienceExposure.getMaskedImage(), kbbox, deep=
False).getMask())
469 if not((bm1 & badBitMask)
or (bm2 & badBitMask)):
470 candidateOutList.append({
'source': kernelCandidate,
472 log.info(
"Selected %d / %d sources for KernelCandidacy", len(candidateOutList), len(candidateInList))
473 return candidateOutList
477 basisList, doBuild=False):
478 """Convert a list of Sources into KernelCandidates.
480 The KernelCandidates are used for fitting the Psf-matching kernel.
484 sourceTable : TODO: DM-17458
486 templateExposure : TODO: DM-17458
488 scienceExposure : TODO: DM-17458
490 kConfig : TODO: DM-17458
492 dConfig : TODO: DM-17458
496 basisList : TODO: DM-17458
498 doBuild : `bool`, optional
506 kernelSize = basisList[0].getWidth()
508 kernelSize, dConfig, log)
511 if doBuild
and not basisList:
514 ps = pexConfig.makePropertySet(kConfig)
515 visitor = diffimLib.BuildSingleKernelVisitorF(basisList, ps)
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)
524 visitor.processCandidate(kCand)
525 kCand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN)
526 candList.append(kCand)
536 """A functor to evaluate the Bayesian Information Criterion for the number of basis sets
537 going into the kernel fitting"""
539 def __init__(self, psfMatchConfig, psfFwhmPixTc, psfFwhmPixTnc):
543 if not self.
psfMatchConfigpsfMatchConfig.kernelBasisSet ==
"alard-lupton":
544 raise RuntimeError(
"BIC only implemnted for AL (alard lupton) basis")
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]
556 visitor = diffimLib.BuildSingleKernelVisitorF(kList,
557 pexConfig.makePropertySet(bicConfig))
558 visitor.setSkipBuilt(
False)
559 kernelCellSet.visitCandidates(visitor, bicConfig.nStarPerCell)
561 for cell
in kernelCellSet.getCellList():
562 for cand
in cell.begin(
False):
563 if cand.getStatus() != afwMath.SpatialCellCandidate.GOOD:
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
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)
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]
A compact representation of a collection of pixels.
A class to represent a 2-dimensional array of pixels.
Represent a 2-dimensional array of bitmask pixels.
A class to manipulate images, masks, and variance as a single object.
A kernel described by a function.
Pass parameters to a Background object.
Parameters to control convolution.
A kernel created from an Image.
A kernel that is a linear combination of fixed basis kernels.
A class that can be used to generate sequences of random numbers according to a number of different a...
Class to ensure constraints for spatial modeling.
A collection of SpatialCells covering an entire image.
Record class that contains measurements made on a single exposure.
An integer coordinate rectangle.
Class stored in SpatialCells for spatial Kernel fitting.
daf::base::PropertyList * list
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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.
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
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)
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.
def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, basisSigmaGauss=None, metadata=None)
def getLogger(loggername)