LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
Namespaces | Classes | Typedefs | Functions
lsst::meas::algorithms Namespace Reference

Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations. More...

Namespaces

namespace  accumulator_mean_stack
 
namespace  astrometrySourceSelector
 
namespace  brightStarStamps
 
namespace  coaddPsf
 
namespace  convertRefcatManager
 
namespace  convertReferenceCatalog
 
namespace  detection
 
namespace  dynamicDetection
 
namespace  findCosmicRaysConfig
 
namespace  flaggedSourceSelector
 
namespace  gaussianPsfFactory
 
namespace  htmIndexer
 
namespace  indexerRegistry
 
namespace  ingestIndexReferenceTask
 
namespace  installGaussianPsf
 
namespace  interp
 
namespace  loadIndexedReferenceObjects
 
namespace  loadReferenceObjects
 
namespace  makeCoaddApCorrMap
 
namespace  makePsfCandidates
 
namespace  matcherSourceSelector
 
namespace  measureApCorr
 
namespace  objectSizeStarSelector
 
namespace  pcaPsfDeterminer
 
namespace  psfCandidate
 
namespace  psfDeterminer
 
namespace  psfSelectionFromMatchList
 
namespace  readFitsCatalogTask
 
namespace  readTextCatalogTask
 
namespace  reserveSourcesTask
 
namespace  scaleVariance
 
namespace  simple_curve
 
namespace  skyObjects
 
namespace  sourceSelector
 
namespace  stamps
 
namespace  starSelector
 
namespace  subtractBackground
 
namespace  testUtils
 
namespace  utils
 
namespace  version
 

Classes

class  CoaddBoundedField
 
struct  CoaddBoundedFieldElement
 Struct used to hold one Exposure's data in a CoaddBoundedField. More...
 
class  CoaddPsf
 CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images. More...
 
class  CoaddPsfControl
 
class  Defect
 Encapsulate information about a bad portion of a detector. More...
 
class  DoubleGaussianPsf
 Represent a Psf as a circularly symmetrical double Gaussian. More...
 
class  evalChi2Visitor
 A class to pass around to all our PsfCandidates to evaluate the PSF fit's X^2. More...
 
class  ExposurePatch
 A convenience container for the exposure, peak and footprint that will be measured. More...
 
class  ImagePsf
 An intermediate base class for Psfs that use an image representation. More...
 
class  ImagePsfTrampoline
 "Trampoline" for ImagePsf to let it be used as a base class in Python. More...
 
class  KernelPsf
 A Psf defined by a Kernel. More...
 
class  KernelPsfFactory
 A PersistableFactory for KernelPsf and its subclasses. More...
 
struct  KernelPsfPersistenceHelper
 A read-only singleton struct containing the schema and key used in persistence for KernelPsf. More...
 
class  MinimizeChi2
 
class  PcaPsf
 Represent a PSF as a linear combination of PCA (== Karhunen-Loeve) basis functions. More...
 
class  PsfCandidate
 Class stored in SpatialCells for spatial Psf fitting. More...
 
class  PsfImagePca
 
class  SingleGaussianPsf
 Represent a PSF as a circularly symmetrical Gaussian. More...
 
class  WarpedPsf
 A Psf class that maps an arbitrary Psf through a coordinate transformation. More...
 

Typedefs

typedef std::vector< Defect::Ptr >::const_iterator DefectCIter
 

Functions

std::shared_ptr< afw::image::TransmissionCurve const > makeCoaddTransmissionCurve (std::shared_ptr< afw::geom::SkyWcs const > coaddWcs, afw::table::ExposureCatalog const &inputSensors)
 Create a TransmissionCurve that represents the effective throughput on a coadd. More...
 
template<typename MaskedImageT >
std::vector< std::shared_ptr< afw::detection::Footprint > > findCosmicRays (MaskedImageT &mimage, afw::detection::Psf const &psf, double const bkgd, daf::base::PropertySet const &ps, bool const keep)
 Find cosmic rays in an Image, and mask and remove them. More...
 
template<typename ExposureT >
std::shared_ptr< ExposurePatch< ExposureT > > makeExposurePatch (std::shared_ptr< ExposureT const > exp, std::shared_ptr< afw::detection::Footprint const > foot, geom::Point2D const &center)
 Factory function for ExposurePatch. More...
 
template<typename ExposureT >
std::shared_ptr< ExposurePatch< ExposureT > > makeExposurePatch (std::shared_ptr< ExposureT const > exp, afw::detection::Footprint const &standardFoot, geom::Point2D const &standardCenter, afw::geom::SkyWcs const &standardWcs)
 
template<typename MaskedImageT >
void interpolateOverDefects (MaskedImageT &mimage, afw::detection::Psf const &, std::vector< Defect::Ptr > &_badList, double fallbackValue, bool useFallbackValueAtEdge)
 Process a set of known bad pixels in an image. More...
 
template<typename PixelT >
std::shared_ptr< PsfCandidate< PixelT > > makePsfCandidate (std::shared_ptr< afw::table::SourceRecord > const &source, std::shared_ptr< afw::image::Exposure< PixelT > > image)
 Return a PsfCandidate of the right sort. More...
 
template<typename PixelT >
std::pair< std::shared_ptr< afw::math::LinearCombinationKernel >, std::vector< double > > createKernelFromPsfCandidates (afw::math::SpatialCellSet const &psfCells, lsst::geom::Extent2I const &dims, lsst::geom::Point2I const &xy0, int const nEigenComponents, int const spatialOrder, int const ksize, int const nStarPerCell, bool const constantWeight, int const border)
 Return a Kernel pointer and a list of eigenvalues resulting from analysing the provided SpatialCellSet. More...
 
template<typename PixelT >
int countPsfCandidates (afw::math::SpatialCellSet const &psfCells, int const nStarPerCell)
 Count the number of candidates in use. More...
 
template<typename PixelT >
std::pair< bool, double > fitSpatialKernelFromPsfCandidates (afw::math::Kernel *kernel, afw::math::SpatialCellSet const &psfCells, int const nStarPerCell, double const tolerance, double const lambda)
 Fit spatial kernel using full-nonlinear optimization to estimate candidate amplitudes. More...
 
template<typename PixelT >
std::pair< bool, double > fitSpatialKernelFromPsfCandidates (afw::math::Kernel *kernel, afw::math::SpatialCellSet const &psfCells, bool const doNonLinearFit, int const nStarPerCell=-1, double const tolerance=1e-5, double const lambda=0.0)
 
template<typename ImageT >
double subtractPsf (afw::detection::Psf const &psf, ImageT *data, double x, double y, double psfFlux=std::numeric_limits< double >::quiet_NaN())
 
template<typename Image >
std::pair< std::vector< double >, afw::math::KernelListfitKernelParamsToImage (afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
 Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary. More...
 
template<typename Image >
std::pair< std::shared_ptr< afw::math::Kernel >, std::pair< double, double > > fitKernelToImage (afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
 Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary. More...
 
 PYBIND11_MODULE (psfCandidate, mod)
 
geom::Box2I getOverallBBox (std::vector< std::shared_ptr< afw::image::Image< double > > > const &imgVector)
 
void setSpatialParameters (afw::math::Kernel *kernel, std::vector< double > const &coeffs)
 Fit a Kernel's spatial variability from a set of stars. More...
 
void setSpatialParameters (afw::math::Kernel *kernel, Eigen::VectorXd const &vec)
 Fit a Kernel's spatial variability from a set of stars. More...
 
template<typename MaskedImageT >
double subtractPsf (afw::detection::Psf const &psf, MaskedImageT *data, double x, double y, double psfFlux)
 Subtract a PSF from an image at a given position. More...
 

Detailed Description

Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.

Typedef Documentation

◆ DefectCIter

Definition at line 49 of file Interp.cc.

Function Documentation

◆ countPsfCandidates()

template<typename PixelT >
int lsst::meas::algorithms::countPsfCandidates ( afw::math::SpatialCellSet const &  psfCells,
int const  nStarPerCell = -1 
)

Count the number of candidates in use.

Definition at line 337 of file SpatialModelPsf.cc.

337 {
338 countVisitor<PixelT> counter;
339 psfCells.visitCandidates(&counter, nStarPerCell);
340
341 return counter.getN();
342}

◆ createKernelFromPsfCandidates()

template<typename PixelT >
std::pair< std::shared_ptr< afw::math::LinearCombinationKernel >, std::vector< double > > lsst::meas::algorithms::createKernelFromPsfCandidates ( afw::math::SpatialCellSet const &  psfCells,
lsst::geom::Extent2I const &  dims,
lsst::geom::Point2I const &  xy0,
int const  nEigenComponents,
int const  spatialOrder,
int const  ksize,
int const  nStarPerCell,
bool const  constantWeight,
int const  border 
)

Return a Kernel pointer and a list of eigenvalues resulting from analysing the provided SpatialCellSet.

The Kernel is a LinearCombinationKernel of the first nEigenComponents eigenImages

N.b. This is templated over the Pixel type of the science image

Parameters
psfCellsA SpatialCellSet containing PsfCandidates
dimsDimensions of image
xy0Origin of image
nEigenComponentsnumber of eigen components to keep; <= 0 => infty
spatialOrderOrder of spatial variation (cf. afw::math::PolynomialFunction2)
ksizeSize of generated Kernel images
nStarPerCellmax no. of stars per cell; <= 0 => infty
constantWeightshould each star have equal weight in the fit?
borderBorder size for background subtraction

Definition at line 193 of file SpatialModelPsf.cc.

203 {
204 typedef typename afw::image::Image<PixelT> ImageT;
205 typedef typename afw::image::MaskedImage<PixelT> MaskedImageT;
206
207 //
208 // Set the sizes for PsfCandidates made from either Images or MaskedImages
209 //
210 // lsst::meas::algorithms::PsfCandidate<ImageT>::setWidth(ksize);
211 // lsst::meas::algorithms::PsfCandidate<ImageT>::setHeight(ksize);
212 // lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setWidth(ksize);
213 // lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setHeight(ksize);
216
217 PsfImagePca<MaskedImageT> imagePca(constantWeight, border); // Here's the set of images we'll analyze
218
219 {
220 SetPcaImageVisitor<PixelT> importStarVisitor(&imagePca);
221 bool const ignoreExceptions = true;
222 psfCells.visitCandidates(&importStarVisitor, nStarPerCell, ignoreExceptions);
223 }
224
225 //
226 // Do a PCA decomposition of those PSF candidates.
227 //
228 // We have "gappy" data; in other words we don't want to include any pixels with INTRP set
229 //
230 int niter = 10; // number of iterations of updateBadPixels
231 double deltaLim = 10.0; // acceptable value of delta, the max change due to updateBadPixels
235
236 for (int i = 0; i != niter; ++i) {
237 int const ncomp =
238 (i == 0) ? 0
239 : ((nEigenComponents == 0) ? imagePca.getEigenImages().size() : nEigenComponents);
240 double delta = imagePca.updateBadPixels(BAD | CR | INTRP, ncomp);
241 if (i > 0 && delta < deltaLim) {
242 break;
243 }
244
245 imagePca.analyze();
246 }
247
248 std::vector<std::shared_ptr<MaskedImageT>> eigenImages = imagePca.getEigenImages();
249 std::vector<double> eigenValues = imagePca.getEigenValues();
250 int const nEigen = static_cast<int>(eigenValues.size());
251
252 int const ncomp = (nEigenComponents <= 0 || nEigen < nEigenComponents) ? nEigen : nEigenComponents;
253 //
254 // Set the background level of the components to 0.0 to avoid coupling variable background
255 // levels to the form of the Psf. More precisely, we calculate the mean of an outer "annulus"
256 // of width bkg_border
257 //
258 for (int k = 0; k != ncomp; ++k) {
259 ImageT const& im = *eigenImages[k]->getImage();
260
261 int bkg_border = 2;
262 if (bkg_border > im.getWidth()) {
263 bkg_border = im.getWidth() / 2;
264 }
265 if (bkg_border > im.getHeight()) {
266 bkg_border = im.getHeight() / 2;
267 }
268
269 double sum = 0;
270 // Bottom and Top borders
271 for (int i = 0; i != bkg_border; ++i) {
272 typename ImageT::const_x_iterator ptrB = im.row_begin(i),
273 ptrT = im.row_begin(im.getHeight() - i - 1);
274 for (int j = 0; j != im.getWidth(); ++j, ++ptrB, ++ptrT) {
275 sum += *ptrB + *ptrT;
276 }
277 }
278 for (int i = bkg_border; i < im.getHeight() - bkg_border; ++i) {
279 // Left and Right borders
280 typename ImageT::const_x_iterator ptrL = im.row_begin(i),
281 ptrR = im.row_begin(i) + im.getWidth() - bkg_border;
282 for (int j = 0; j != bkg_border; ++j, ++ptrL, ++ptrR) {
283 sum += *ptrL + *ptrR;
284 }
285 }
286 sum /= 2 * (bkg_border * im.getWidth() + bkg_border * (im.getHeight() - 2 * bkg_border));
287
288 *eigenImages[k] -= sum;
289 }
290 //
291 // Now build our LinearCombinationKernel; build the lists of basis functions
292 // and spatial variation, then assemble the Kernel
293 //
294 afw::math::KernelList kernelList;
296 geom::Box2D const range = geom::Box2D(geom::Point2D(xy0), geom::Extent2D(dims));
297
298 for (int i = 0; i != ncomp; ++i) {
299 {
300 // Enforce unit sum for kernel by construction
301 // Zeroth component has unit sum
302 // Other components have zero sum by normalising and then subtracting the zeroth component
303 ImageT& image = *eigenImages[i]->getImage();
304 double sum = std::accumulate(image.begin(true), image.end(true), 0.0);
305 if (i == 0) {
306 image /= sum;
307 } else {
308 for (typename ImageT::fast_iterator ptr0 = eigenImages[0]->getImage()->begin(true),
309 ptr1 = image.begin(true), end = image.end(true);
310 ptr1 != end; ++ptr0, ++ptr1) {
311 *ptr1 = *ptr1 / sum - *ptr0;
312 }
313 }
314 }
315
316 kernelList.push_back(std::shared_ptr<afw::math::Kernel>(new afw::math::FixedKernel(
317 afw::image::Image<afw::math::Kernel::Pixel>(*eigenImages[i]->getImage(), true))));
318
319 afw::math::Kernel::SpatialFunctionPtr
320 // spatialFunction(new afw::math::PolynomialFunction2<double>(spatialOrder));
321 spatialFunction(new afw::math::Chebyshev1Function2<double>(spatialOrder, range));
322 spatialFunction->setParameter(0, 1.0); // the constant term; all others are 0
323 spatialFunctionList.push_back(spatialFunction);
324 }
325
327 new afw::math::LinearCombinationKernel(kernelList, spatialFunctionList));
328
329 return std::make_pair(psf, eigenValues);
330}
int end
T accumulate(T... args)
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
static void setWidth(int width)
Set the width of the image that getImage should return.
Definition: SpatialCell.h:136
static void setHeight(int height)
Set the height of the image that getImage should return.
Definition: SpatialCell.h:141
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
T make_pair(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::vector< std::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:462
T push_back(T... args)
T size(T... args)
Key< int > psf
Definition: Exposure.cc:65

◆ findCosmicRays()

template<typename MaskedImageT >
std::vector< std::shared_ptr< afw::detection::Footprint > > lsst::meas::algorithms::findCosmicRays ( MaskedImageT &  mimage,
afw::detection::Psf const &  psf,
double const  bkgd,
daf::base::PropertySet const &  ps,
bool const  keep 
)

Find cosmic rays in an Image, and mask and remove them.

Returns
vector of CR's Footprints

In this loop, we look for strings of CRpixels on the same row and adjoining columns; each of these becomes a Span with a unique ID.

In this loop, we look for strings of CRpixels on the same row and adjoining columns; each of these becomes a Span with a unique ID.

Parameters
mimageImage to search
psfthe Image's PSF
bkgdunsubtracted background of frame, DN
psPropertySet directing the behavior
keepif true, don't remove the CRs

Definition at line 315 of file CR.cc.

321 {
322 typedef typename MaskedImageT::Image ImageT;
323 typedef typename ImageT::Pixel ImagePixel;
324 typedef typename MaskedImageT::Mask::Pixel MaskPixel;
325
326 // Parse the PropertySet
327 double const minSigma = ps.getAsDouble("minSigma"); // min sigma over sky in pixel for CR candidate
328 double const minDn = ps.getAsDouble("min_DN"); // min number of DN in an CRs
329 double const cond3Fac = ps.getAsDouble("cond3_fac"); // fiddle factor for condition #3
330 double const cond3Fac2 = ps.getAsDouble("cond3_fac2"); // 2nd fiddle factor for condition #3
331 int const niteration = ps.getAsInt("niteration"); // Number of times to look for contaminated
332 // pixels near CRs
333 int const nCrPixelMax = ps.getAsInt("nCrPixelMax"); // maximum number of contaminated pixels
334 /*
335 * thresholds for 3rd condition
336 *
337 * Realise PSF at center of image
338 */
339 std::shared_ptr<afw::math::Kernel const> kernel = psf.getLocalKernel();
340 if (!kernel) {
341 throw LSST_EXCEPT(pexExcept::NotFoundError, "Psf is unable to return a kernel");
342 }
344 afw::detection::Psf::Image(geom::ExtentI(kernel->getWidth(), kernel->getHeight()));
345 kernel->computeImage(psfImage, true);
346
347 int const xc = kernel->getCtr().getX(); // center of PSF
348 int const yc = kernel->getCtr().getY();
349
350 double const I0 = psfImage(xc, yc);
351 double const thresH =
352 cond3Fac2 * (0.5 * (psfImage(xc - 1, yc) + psfImage(xc + 1, yc))) / I0; // horizontal
353 double const thresV = cond3Fac2 * (0.5 * (psfImage(xc, yc - 1) + psfImage(xc, yc + 1))) / I0; // vertical
354 double const thresD = cond3Fac2 *
355 (0.25 * (psfImage(xc - 1, yc - 1) + psfImage(xc + 1, yc + 1) +
356 psfImage(xc - 1, yc + 1) + psfImage(xc + 1, yc - 1))) /
357 I0; // diag
358 /*
359 * Setup desired mask planes
360 */
361 MaskPixel const badBit = mimage.getMask()->getPlaneBitMask("BAD"); // Generic bad pixels
362 MaskPixel const crBit = mimage.getMask()->getPlaneBitMask("CR"); // CR-contaminated pixels
363 MaskPixel const interpBit = mimage.getMask()->getPlaneBitMask("INTRP"); // Interpolated pixels
364 MaskPixel const saturBit = mimage.getMask()->getPlaneBitMask("SAT"); // Saturated pixels
365 MaskPixel const nodataBit = mimage.getMask()->getPlaneBitMask("NO_DATA"); // Non data pixels
366
367 MaskPixel const badMask = (badBit | interpBit | saturBit | nodataBit); // naughty pixels
368 /*
369 * Go through the frame looking at each pixel (except the edge ones which we ignore)
370 */
371 int const ncol = mimage.getWidth();
372 int const nrow = mimage.getHeight();
373
374 std::vector<CRPixel<ImagePixel>> crpixels; // storage for detected CR-contaminated pixels
375 typedef typename std::vector<CRPixel<ImagePixel>>::iterator crpixel_iter;
376 typedef typename std::vector<CRPixel<ImagePixel>>::reverse_iterator crpixel_riter;
377
378 for (int j = 1; j < nrow - 1; ++j) {
379 typename MaskedImageT::xy_locator loc = mimage.xy_at(1, j); // locator for data
380
381 for (int i = 1; i < ncol - 1; ++i, ++loc.x()) {
382 ImagePixel corr = 0;
383 if (!is_cr_pixel<MaskedImageT>(&corr, loc, minSigma, thresH, thresV, thresD, bkgd, cond3Fac)) {
384 continue;
385 }
386 /*
387 * condition #4
388 */
389 if (loc.mask() & badMask) {
390 continue;
391 }
392 if ((loc.mask(-1, 1) | loc.mask(0, 1) | loc.mask(1, 1) | loc.mask(-1, 0) | loc.mask(1, 0) |
393 loc.mask(-1, -1) | loc.mask(0, -1) | loc.mask(1, -1)) &
394 interpBit) {
395 continue;
396 }
397 /*
398 * OK, it's a CR
399 *
400 * replace CR-contaminated pixels with reasonable values as we go through
401 * image, which increases the detection rate
402 */
403 crpixels.push_back(CRPixel<ImagePixel>(i + mimage.getX0(), j + mimage.getY0(), loc.image()));
404 loc.image() = corr; /* just a preliminary estimate */
405
406 if (static_cast<int>(crpixels.size()) > nCrPixelMax) {
407 reinstateCrPixels(mimage.getImage().get(), crpixels);
408
409 throw LSST_EXCEPT(pex::exceptions::LengthError,
410 (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
411 }
412 }
413 }
414 /*
415 * We've found them on a pixel-by-pixel basis, now merge those pixels
416 * into cosmic rays
417 */
418 std::vector<int> aliases; // aliases for initially disjoint parts of CRs
419 aliases.reserve(1 + crpixels.size() / 2); // initial size of aliases
420
421 std::vector<afw::detection::IdSpan::Ptr> spans; // y:x0,x1 for objects
422 spans.reserve(aliases.capacity()); // initial size of spans
423
424 aliases.push_back(0); // 0 --> 0
425
431 int ncr = 0; // number of detected cosmic rays
432 if (!crpixels.empty()) {
433 int id; // id number for a CR
434 int x0 = -1, x1 = -1, y = -1; // the beginning and end column, and row of this span in a CR
435
436 // I am dummy
437 CRPixel<ImagePixel> dummy(0, -1, 0, -1);
438 crpixels.push_back(dummy);
439 // printf("Created dummy CR: i %i, id %i, col %i, row %i, val %g\n", dummy.get_i(), dummy.id,
440 // dummy.col, dummy.row, (double)dummy.val);
441 for (crpixel_iter crp = crpixels.begin(); crp < crpixels.end() - 1; ++crp) {
442 // printf("Looking at CR: i %i, id %i, col %i, row %i, val %g\n", crp->get_i(), crp->id, crp->col,
443 // crp->row, (double)crp->val);
444
445 if (crp->id < 0) { // not already assigned
446 crp->id = ++ncr; // a new CR
447 aliases.push_back(crp->id);
448 y = crp->row;
449 x0 = x1 = crp->col;
450 // printf(" Assigned ID %i; looking at row %i, start col %i\n", crp->id, crp->row, crp->col);
451 }
452 id = crp->id;
453 // printf(" Next CRpix has i=%i, id=%i, row %i, col %i\n", crp[1].get_i(), crp[1].id, crp[1].row,
454 // crp[1].col);
455
456 if (crp[1].row == crp[0].row && crp[1].col == crp[0].col + 1) {
457 // printf(" Adjoining! Set next CRpix id = %i; x1=%i\n", crp[1].id, x1);
458 crp[1].id = id;
459 ++x1;
460 } else {
461 assert(y >= 0 && x0 >= 0 && x1 >= 0);
462 spans.push_back(afw::detection::IdSpan::Ptr(new afw::detection::IdSpan(id, y, x0, x1)));
463 // printf(" Not adjoining; adding span id=%i, y=%i, x = [%i, %i]\n", id, y, x0, x1);
464 }
465 }
466 }
467
468 // At the end of this loop, all crpixel entries have been assigned an ID,
469 // except for the "dummy" entry at the end of the array.
470 if (crpixels.size() > 0) {
471 for (crpixel_iter cp = crpixels.begin(); cp != crpixels.end() - 1; cp++) {
472 assert(cp->id >= 0);
473 assert(cp->col >= 0);
474 assert(cp->row >= 0);
475 }
476 // dummy:
477 assert(crpixels[crpixels.size() - 1].id == -1);
478 assert(crpixels[crpixels.size() - 1].col == 0);
479 assert(crpixels[crpixels.size() - 1].row == -1);
480 }
481
482 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
483 sp++) {
484 assert((*sp)->id >= 0);
485 assert((*sp)->y >= 0);
486 assert((*sp)->x0 >= 0);
487 assert((*sp)->x1 >= (*sp)->x0);
488 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; sp2++) {
489 assert((*sp2)->y >= (*sp)->y);
490 if ((*sp2)->y == (*sp)->y) {
491 assert((*sp2)->x0 > (*sp)->x1);
492 }
493 }
494 }
495
496 /*
497 * See if spans touch each other
498 */
499 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
500 ++sp) {
501 int const y = (*sp)->y;
502 int const x0 = (*sp)->x0;
503 int const x1 = (*sp)->x1;
504
505 // this loop will probably run for only a few steps
506 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
507 if ((*sp2)->y == y) {
508 // on this row (but not adjoining columns, since it would have been merged into this span);
509 // keep looking.
510 continue;
511 } else if ((*sp2)->y != (y + 1)) {
512 // sp2 is more than one row below; can't be connected.
513 break;
514 } else if ((*sp2)->x0 > (x1 + 1)) {
515 // sp2 is more than one column away to the right; can't be connected
516 break;
517 } else if ((*sp2)->x1 >= (x0 - 1)) {
518 // touches
519 int r1 = afw::detection::resolve_alias(aliases, (*sp)->id);
520 int r2 = afw::detection::resolve_alias(aliases, (*sp2)->id);
521 aliases[r1] = r2;
522 }
523 }
524 }
525
526 /*
527 * Resolve aliases; first alias chains, then the IDs in the spans
528 */
529 for (unsigned int i = 0; i != spans.size(); ++i) {
530 spans[i]->id = afw::detection::resolve_alias(aliases, spans[i]->id);
531 }
532
533 /*
534 * Sort spans by ID, so we can sweep through them once
535 */
536 if (spans.size() > 0) {
537 std::sort(spans.begin(), spans.end(), afw::detection::IdSpanCompar());
538 }
539
540 /*
541 * Build Footprints from spans
542 */
544
545 if (spans.size() > 0) {
546 int id = spans[0]->id;
547 unsigned int i0 = 0; // initial value of i
548 for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
549 if (i == spans.size() || spans[i]->id != id) {
551
553 spanList.reserve(i - i0);
554 for (; i0 < i; ++i0) {
555 spanList.push_back(afw::geom::Span(spans[i0]->y, spans[i0]->x0, spans[i0]->x1));
556 }
557 cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(spanList)));
558 CRs.push_back(cr);
559 }
560
561 if (i < spans.size()) {
562 id = spans[i]->id;
563 }
564 }
565 }
566
567 reinstateCrPixels(mimage.getImage().get(), crpixels);
568 /*
569 * apply condition #1
570 */
571 CountsInCR<typename ImageT::Pixel> CountDN(bkgd);
572 for (std::vector<std::shared_ptr<afw::detection::Footprint>>::iterator cr = CRs.begin(), end = CRs.end();
573 cr != end; ++cr) {
574 // find the sum of pixel values within the CR
575 (*cr)->getSpans()->applyFunctor(CountDN, *mimage.getImage());
576
577 LOGL_DEBUG("TRACE4.lsst.algorithms.CR", "CR at (%d, %d) has %g DN", (*cr)->getBBox().getMinX(),
578 (*cr)->getBBox().getMinY(), CountDN.getCounts());
579 if (CountDN.getCounts() < minDn) { /* not bright enough */
580 LOGL_DEBUG("TRACE5.lsst.algorithms.CR", "Erasing CR");
581
582 cr = CRs.erase(cr);
583 --cr; // back up to previous CR (we're going to increment it)
584 --end;
585 }
586 CountDN.reset();
587 }
588 ncr = CRs.size(); /* some may have been too faint */
589 /*
590 * We've found them all, time to kill them all
591 */
592 bool const debias_values = true;
593 bool grow = false;
594 LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing initial list of CRs");
595 removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
596#if 0 // Useful to see phase 2 in display; debugging only
597 (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
598 mimage.getMask()->getPlaneBitMask("DETECTED"));
599#endif
600 /*
601 * Now that we've removed them, go through image again, examining area around
602 * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
603 *
604 * We iterate niteration times; niter==1 was sufficient for SDSS data, but megacam
605 * CCDs are different -- who knows for other devices?
606 */
607 bool too_many_crs = false; // we've seen too many CR pixels
608 int nextra = 0; // number of pixels added to list of CRs
609 for (int i = 0; i != niteration && !too_many_crs; ++i) {
610 LOGL_DEBUG("TRACE1.lsst.algorithms.CR", "Starting iteration %d", i);
612 fiter != CRs.end(); fiter++) {
614 /*
615 * Are all those `CR' pixels interpolated? If so, don't grow it
616 */
617 {
618 // this work should be taken on in DM-9538
619 // std::shared_ptr<afw::detection::Footprint> om = footprintAndMask(cr, mimage.getMask(),
620 // interpBit);
621 // Placeholder until DM-9538 then remove empty footprint
622 auto om = std::make_shared<afw::detection::Footprint>();
623 int const npix = (om) ? om->getArea() : 0;
624
625 if (static_cast<std::size_t>(npix) == cr->getArea()) {
626 continue;
627 }
628 }
629 /*
630 * No; some of the suspect pixels aren't interpolated
631 */
632 afw::detection::Footprint extra; // extra pixels added to cr
633 for (auto siter = cr->getSpans()->begin(); siter != cr->getSpans()->end(); siter++) {
634 auto const span = siter;
635
636 /*
637 * Check the lines above and below the span. We're going to check a 3x3 region around
638 * the pixels, so we need a buffer around the edge. We check the pixels just to the
639 * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
640 * column direction, but only 1 in the row direction.
641 */
642 int const y = span->getY() - mimage.getY0();
643 if (y < 2 || y >= nrow - 2) {
644 continue;
645 }
646 int x0 = span->getX0() - mimage.getX0();
647 int x1 = span->getX1() - mimage.getX0();
648 x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
649 x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;
650
651 checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
652 bkgd, 0, keep);
653 checkSpanForCRs(&extra, crpixels, y, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
654 bkgd, 0, keep);
655 checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
656 bkgd, 0, keep);
657 }
658
659 if (extra.getSpans()->size() > 0) { // we added some pixels
660 if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
661 too_many_crs = true;
662 break;
663 }
664
665 nextra += extra.getArea();
666
667 std::vector<afw::geom::Span> tmpSpanList(cr->getSpans()->begin(), cr->getSpans()->end());
668 for (auto const &spn : (*extra.getSpans())) {
669 tmpSpanList.push_back(spn);
670 }
671 cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(tmpSpanList)));
672 }
673 }
674
675 if (nextra == 0) {
676 break;
677 }
678 }
679 /*
680 * mark those pixels as CRs
681 */
682 if (!too_many_crs) {
683 for (auto const &foot : CRs) {
684 foot->getSpans()->setMask(*mimage.getMask().get(), crBit);
685 }
686 }
687 /*
688 * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
689 * first value stored (hence the uses of rbegin/rend)
690 *
691 * We have to do this if we decide _not_ to remove certain CRs,
692 * for example those which lie next to saturated pixels
693 */
694 if (keep || too_many_crs) {
695 if (crpixels.size() > 0) {
696 int const imageX0 = mimage.getX0();
697 int const imageY0 = mimage.getY0();
698
699 std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
700
701 crpixel_riter rend = crpixels.rend();
702 for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
703 if (crp->row == -1)
704 // dummy; skip it.
705 continue;
706 mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
707 }
708 }
709 } else {
710 if (true || nextra > 0) {
711 grow = true;
712 LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing final list of CRs, grow = %d", grow);
713 removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
714 }
715 /*
716 * we interpolated over all CR pixels, so set the interp bits too
717 */
718 for (auto const &foot : CRs) {
719 foot->getSpans()->setMask(*mimage.getMask().get(), static_cast<MaskPixel>(crBit | interpBit));
720 }
721 }
722
723 if (too_many_crs) { // we've cleaned up, so we can throw the exception
724 throw LSST_EXCEPT(pex::exceptions::LengthError,
725 (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
726 }
727
728 return CRs;
729}
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
int y
Definition: SpanSet.cc:48
T capacity(T... args)
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
T empty(T... args)
T end(T... args)
T erase(T... args)
T move(T... args)
int resolve_alias(const std::vector< int > &aliases, int id)
Follow a chain of aliases, returning the final resolved value.
Definition: CR.cc:96
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
lsst::afw::detection::Footprint Footprint
Definition: Source.h:57
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
T rbegin(T... args)
T rend(T... args)
T reserve(T... args)
T sort(T... args)
int id
Definition: CR.cc:143
int row
Definition: CR.cc:145
int col
Definition: CR.cc:144

◆ fitKernelParamsToImage()

template<typename Image >
std::pair< std::vector< double >, afw::math::KernelList > lsst::meas::algorithms::fitKernelParamsToImage ( afw::math::LinearCombinationKernel const &  kernel,
Image const &  image,
geom::Point2D const &  pos 
)

Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

Returns
std::pair(coefficients, std::pair(kernels, center amplitude))
Parameters
kernelthe Kernel to fit
imagethe image to be fit
posthe position of the object

Definition at line 1037 of file SpatialModelPsf.cc.

1041 {
1043
1044 afw::math::KernelList kernels = kernel.getKernelList(); // the Kernels that kernel adds together
1045 int const nKernel = kernels.size();
1046
1047 if (nKernel == 0) {
1048 throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Your kernel must have at least one component");
1049 }
1050
1051 /*
1052 * Go through all the kernels, get a copy centered at the desired sub-pixel position, and then
1053 * extract a subImage from the parent image at the same place
1054 */
1055 std::vector<std::shared_ptr<KernelT>> kernelImages = offsetKernel<KernelT>(kernel, pos[0], pos[1]);
1056 geom::BoxI bbox(kernelImages[0]->getBBox());
1057 Image const& subImage(Image(image, bbox, afw::image::PARENT, false)); // shallow copy
1058
1059 /*
1060 * Solve the linear problem subImage = sum x_i K_i + epsilon; we solve this for x_i by constructing the
1061 * normal equations, A x = b
1062 */
1063 Eigen::MatrixXd A(nKernel, nKernel);
1064 Eigen::VectorXd b(nKernel);
1065
1066 for (int i = 0; i != nKernel; ++i) {
1067 b(i) = afw::image::innerProduct(*kernelImages[i], *subImage.getImage());
1068
1069 for (int j = i; j != nKernel; ++j) {
1070 A(i, j) = A(j, i) = afw::image::innerProduct(*kernelImages[i], *kernelImages[j]);
1071 }
1072 }
1073 Eigen::VectorXd x(nKernel);
1074
1075 if (nKernel == 1) {
1076 x(0) = b(0) / A(0, 0);
1077 } else {
1078 x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
1079 }
1080
1081 // the XY0() point of the shifted Kernel basis functions
1082 geom::Point2I const xy0 = kernelImages[0]->getXY0();
1083
1084 afw::math::KernelList newKernels(nKernel);
1085 std::vector<double> params(nKernel);
1086 for (int i = 0; i != nKernel; ++i) {
1087 std::shared_ptr<afw::math::Kernel> newKernel(new afw::math::FixedKernel(*kernelImages[i]));
1088 newKernel->setCtr(xy0 + newKernel->getDimensions() / 2);
1089
1090 params[i] = x[i];
1091 newKernels[i] = newKernel;
1092 }
1093
1094 return std::make_pair(params, newKernels);
1095}
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
double x
table::Key< int > b
An integer coordinate rectangle.
Definition: Box.h:55
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition: ImagePca.cc:414

◆ fitKernelToImage()

template<typename Image >
std::pair< std::shared_ptr< afw::math::Kernel >, std::pair< double, double > > lsst::meas::algorithms::fitKernelToImage ( afw::math::LinearCombinationKernel const &  kernel,
Image const &  image,
geom::Point2D const &  pos 
)

Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

Returns
std::pair(best-fit kernel, std::pair(amp, chi^2))
Parameters
kernelthe Kernel to fit
imagethe image to be fit
posthe position of the object

Definition at line 1104 of file SpatialModelPsf.cc.

1108 {
1110 fitKernelParamsToImage(kernel, image, pos);
1111 std::vector<double> params = fit.first;
1112 afw::math::KernelList kernels = fit.second;
1113 int const nKernel = params.size();
1114 assert(kernels.size() == static_cast<unsigned int>(nKernel));
1115
1116 double amp = 0.0;
1117 for (int i = 0; i != nKernel; ++i) {
1118 std::shared_ptr<afw::math::Kernel> base = kernels[i];
1119 std::shared_ptr<afw::math::FixedKernel> k = std::static_pointer_cast<afw::math::FixedKernel>(base);
1120 amp += params[i] * k->getSum();
1121 }
1122
1124 double chisq = 0.0;
1125 outputKernel->setCtr(kernels[0]->getCtr());
1126
1127 return std::make_pair(outputKernel, std::make_pair(amp, chisq));
1128}
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:704
std::pair< std::vector< double >, afw::math::KernelList > fitKernelParamsToImage(afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

◆ fitSpatialKernelFromPsfCandidates() [1/2]

template<typename PixelT >
std::pair< bool, double > lsst::meas::algorithms::fitSpatialKernelFromPsfCandidates ( afw::math::Kernel kernel,
afw::math::SpatialCellSet const &  psfCells,
int const  nStarPerCell = -1,
double const  tolerance = 1e-5,
double const  lambda = 0.0 
)

Fit spatial kernel using full-nonlinear optimization to estimate candidate amplitudes.

Parameters
kernelthe Kernel to fit
psfCellsA SpatialCellSet containing PsfCandidates
nStarPerCellmax no. of stars per cell; <= 0 => infty
toleranceTolerance; how close chi^2 should be to true minimum
lambdafloor for variance is lambda*data

Definition at line 592 of file SpatialModelPsf.cc.

598 {
599 int const nComponents = kernel->getNKernelParameters();
600 int const nSpatialParams = kernel->getNSpatialParameters();
601 //
602 // visitor that evaluates the chi^2 of the current fit
603 //
604 evalChi2Visitor<PixelT> getChi2(*kernel, lambda);
605 //
606 // We have to unpack the Kernel coefficients into a linear array, coeffs
607 //
608 std::vector<double> coeffs; // The coefficients we want to fit
609 coeffs.assign(nComponents * nSpatialParams, 0.0);
610
611 std::vector<double> stepSize; // step sizes
612 stepSize.assign(nComponents * nSpatialParams, 100);
613 //
614 // Translate that into minuit's language
615 //
616 ROOT::Minuit2::MnUserParameters fitPar;
617 std::vector<std::string> paramNames;
618 paramNames.reserve(nComponents * nSpatialParams);
619
620 for (int i = 0, c = 0; c != nComponents; ++c) {
621 coeffs[i] = 1; // the constant part of each spatial order
622 for (int s = 0; s != nSpatialParams; ++s, ++i) {
623 paramNames.push_back((boost::format("C%d:%d") % c % s).str());
624 fitPar.Add(paramNames[i].c_str(), coeffs[i], stepSize[i]);
625 }
626 }
627 fitPar.Fix("C0:0");
628 //
629 // Create the minuit object that knows how to minimise our functor
630 //
631 MinimizeChi2<PixelT> minimizerFunc(getChi2, kernel, psfCells, nStarPerCell, nComponents, nSpatialParams);
632
633 double const errorDef = 1.0; // use +- 1sigma errors
634 minimizerFunc.setErrorDef(errorDef);
635 //
636 // tell minuit about it
637 //
638 ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
639 //
640 // And let it loose
641 //
642 int maxFnCalls = 0; // i.e. unlimited
643 ROOT::Minuit2::FunctionMinimum min =
644 migrad(maxFnCalls, tolerance / (1e-4 * errorDef)); // minuit uses 0.1*1e-3*tolerance*errorDef
645
646 float minChi2 = min.Fval();
647 bool const isValid = min.IsValid() && std::isfinite(minChi2);
648
649 if (true || isValid) { // calculate coeffs even in minuit is unhappy
650 for (int i = 0; i != nComponents * nSpatialParams; ++i) {
651 coeffs[i] = min.UserState().Value(i);
652 }
653
654 setSpatialParameters(kernel, coeffs);
655 }
656
657#if 0 // Estimate errors; we don't really need this
658 ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
659 for (int i = 0, c = 0; c != nComponents; ++c) {
660 for (int s = 0; s != nSpatialParams; ++s, ++i) {
661 char const *name = paramNames[i].c_str();
662 printf("%s %g", name, min.UserState().Value(name));
663 if (isValid && !fitPar.Parameter(fitPar.Index(name)).IsFixed()) {
664 printf(" (%g+%g)\n", minos(i).first, minos(i).second);
665 }
666 printf("\n");
667 }
668 }
669#endif
670 //
671 // One time more through the Candidates setting their chi^2 values. We'll
672 // do all the candidates this time, not just the first nStarPerCell
673 //
674 psfCells.visitAllCandidates(&getChi2, true);
675
676 return std::make_pair(isValid, minChi2);
677}
table::Key< std::string > name
Definition: Amplifier.cc:116
int min
T assign(T... args)
unsigned int getNKernelParameters() const
Return the number of kernel parameters (0 if none)
Definition: Kernel.h:246
int getNSpatialParameters() const
Return the number of spatial parameters (0 if not spatially varying)
Definition: Kernel.h:251
A class to pass around to all our PsfCandidates to evaluate the PSF fit's X^2.
bool isValid
Definition: fits.cc:399
T printf(T... args)
T isfinite(T... args)
void setSpatialParameters(afw::math::Kernel *kernel, Eigen::VectorXd const &vec)
Fit a Kernel's spatial variability from a set of stars.

◆ fitSpatialKernelFromPsfCandidates() [2/2]

template<typename PixelT >
std::pair< bool, double > lsst::meas::algorithms::fitSpatialKernelFromPsfCandidates ( afw::math::Kernel kernel,
afw::math::SpatialCellSet const &  psfCells,
bool const  doNonLinearFit,
int const  nStarPerCell = -1,
double const  tolerance = 1e-5,
double const  lambda = 0.0 
)
Parameters
kernelthe Kernel to fit
psfCellsA SpatialCellSet containing PsfCandidates
doNonLinearFitUse the full-up nonlinear fitter
nStarPerCellmax no. of stars per cell; <= 0 => infty
toleranceTolerance; how close chi^2 should be to true minimum
lambdafloor for variance is lambda*data

Definition at line 876 of file SpatialModelPsf.cc.

883 {
884 if (doNonLinearFit) {
885 return fitSpatialKernelFromPsfCandidates<PixelT>(kernel, psfCells, nStarPerCell, tolerance);
886 }
887
888 double const tau = 0; // softening for errors
889
890 afw::math::LinearCombinationKernel const* lcKernel =
891 dynamic_cast<afw::math::LinearCombinationKernel const*>(kernel);
892 if (!lcKernel) {
893 throw LSST_EXCEPT(
895 "Failed to cast Kernel to LinearCombinationKernel while building spatial PSF model");
896 }
897#if 1
898 //
899 // Set the initial amplitudes of all our candidates
900 //
901 setAmplitudeVisitor<PixelT> setAmplitude;
902 psfCells.visitAllCandidates(&setAmplitude, true);
903#endif
904 //
905 // visitor that fills out the A and b matrices (we'll solve A x = b for the coeffs, x)
906 //
907 FillABVisitor<PixelT> getAB(*lcKernel, tau);
908 //
909 // Actually visit all our candidates
910 //
911 psfCells.visitCandidates(&getAB, nStarPerCell, true);
912 //
913 // Extract A and b, and solve Ax = b
914 //
915 Eigen::MatrixXd const& A = getAB.getA();
916 Eigen::VectorXd const& b = getAB.getB();
917 Eigen::VectorXd x0(b.size()); // Solution to matrix problem
918
919 switch (b.size()) {
920 case 0: // One candidate, no spatial variability
921 break;
922 case 1: // eigen can't/won't handle 1x1 matrices
923 x0(0) = b(0) / A(0, 0);
924 break;
925 default:
926 x0 = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
927 break;
928 }
929#if 0
930 std::cout << "A " << A << std::endl;
931 std::cout << "b " << b.transpose() << std::endl;
932 std::cout << "x " << x.transpose() << std::endl;
933
934 afw::image::Image<double> img(b.size(), b.size());
935 for (int j = 0; j < b.size(); ++j) {
936 for (int i = 0; i < b.size(); ++i) {
937 img(i, j) = A(i, j);
938 }
939 }
940 img.writeFits("a.fits");
941
942 if (x.cols() >= 6) {
943 for (int i = 0; i != 6; ++i) {
944 double xcen = 25; double ycen = 35 + 35*i;
945 std::cout << "x, y " << xcen << " , " << ycen << " b "
946 << (x[3] + xcen*x[4] + ycen*x[5])/(x[0] + xcen*x[1] + ycen*x[2]) << std::endl;
947 }
948 }
949#endif
950
951 // Generate kernel parameters (including 0th component) from matrix solution
952 Eigen::VectorXd x(kernel->getNKernelParameters() * kernel->getNSpatialParameters()); // Kernel parameters
953 x(0) = 1.0;
954 std::fill(x.data() + 1, x.data() + kernel->getNSpatialParameters(), 0.0);
955 std::copy(x0.data(), x0.data() + x0.size(), x.data() + kernel->getNSpatialParameters());
956
957 setSpatialParameters(kernel, x);
958 //
959 // One time more through the Candidates setting their chi^2 values. We'll
960 // do all the candidates this time, not just the first nStarPerCell
961 //
962 // visitor that evaluates the chi^2 of the current fit
963 //
964 evalChi2Visitor<PixelT> getChi2(*kernel, lambda);
965
966 psfCells.visitAllCandidates(&getChi2, true);
967
968 return std::make_pair(true, getChi2.getValue());
969}
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
T copy(T... args)
T endl(T... args)
T fill(T... args)

◆ getOverallBBox()

geom::Box2I lsst::meas::algorithms::getOverallBBox ( std::vector< std::shared_ptr< afw::image::Image< double > > > const &  imgVector)

Definition at line 188 of file CoaddPsf.cc.

188 {
190 // Calculate the box which will contain them all
191 for (unsigned int i = 0; i < imgVector.size(); i++) {
192 std::shared_ptr<afw::image::Image<double>> componentImg = imgVector[i];
193 geom::Box2I cBBox = componentImg->getBBox();
194 bbox.include(cBBox); // JFB: this works even on empty bboxes
195 }
196 return bbox;
197}

◆ interpolateOverDefects()

template<typename MaskedImageT >
void lsst::meas::algorithms::interpolateOverDefects ( MaskedImageT &  mimage,
afw::detection::Psf const &  psf,
std::vector< Defect::Ptr > &  _badList,
double  fallbackValue = 0.0,
bool  useFallbackValueAtEdge = false 
)

Process a set of known bad pixels in an image.

Parameters
mimageImage to patch
psfthe Image's PSF
_badListList of Defects to patch
fallbackValueValue to fallback to if all else fails
useFallbackValueAtEdgeUse the fallback value at the image's edge?

Definition at line 2047 of file Interp.cc.

2052 {
2053 /*
2054 * Allow for image's origin
2055 */
2056 int const width = mimage.getWidth();
2057 int const height = mimage.getHeight();
2058
2060 badList.reserve(_badList.size());
2061 for (std::vector<Defect::Ptr>::iterator ptr = _badList.begin(), end = _badList.end(); ptr != end; ++ptr) {
2062 geom::BoxI bbox = (*ptr)->getBBox();
2063 bbox.shift(geom::ExtentI(-mimage.getX0(), -mimage.getY0())); // allow for image's origin
2064 geom::PointI min = bbox.getMin(), max = bbox.getMax();
2065 if (min.getX() >= width) {
2066 continue;
2067 } else if (min.getX() < 0) {
2068 if (max.getX() < 0) {
2069 continue;
2070 } else {
2071 min.setX(0);
2072 }
2073 }
2074
2075 if (max.getX() < 0) {
2076 continue;
2077 } else if (max.getX() >= width) {
2078 max.setX(width - 1);
2079 }
2080
2081 bbox = geom::BoxI(min, max);
2082 Defect::Ptr ndefect(new Defect(bbox));
2083 ndefect->classify((*ptr)->getPos(), (*ptr)->getType());
2084 badList.push_back(ndefect);
2085 }
2086
2087 sort(badList.begin(), badList.end(), Sort_ByX0<Defect>());
2088 /*
2089 * Go through the frame looking at each pixel (except the edge ones which we ignore)
2090 */
2091 typename MaskedImageT::Mask::Pixel const interpBit =
2092 mimage.getMask()->getPlaneBitMask("INTRP"); // interp'd pixels
2093
2094 constexpr int nUseInterp = 6; // no. of pixels to interpolate towards edge
2095 static_assert(nUseInterp < Defect::WIDE_DEFECT,
2096 "make sure that we can handle these defects using"
2097 "the full interpolation not edge code");
2098
2099 for (int y = 0; y != height; y++) {
2100 std::vector<Defect::Ptr> badList1D = classify_defects(badList, y, width);
2101
2102 do_defects(badList1D, y, *mimage.getImage(),
2104 useFallbackValueAtEdge, nUseInterp);
2105
2106 do_defects(badList1D, y, *mimage.getMask(), interpBit, useFallbackValueAtEdge, nUseInterp);
2107
2108 do_defects(badList1D, y, *mimage.getVariance(),
2110 useFallbackValueAtEdge, nUseInterp);
2111 }
2112}
int max
uint64_t * ptr
Definition: RangeSet.cc:88
Box2I BoxI
Definition: Box.h:780

◆ makeCoaddTransmissionCurve()

std::shared_ptr< afw::image::TransmissionCurve const > lsst::meas::algorithms::makeCoaddTransmissionCurve ( std::shared_ptr< afw::geom::SkyWcs const >  coaddWcs,
afw::table::ExposureCatalog const &  inputSensors 
)

Create a TransmissionCurve that represents the effective throughput on a coadd.

Parameters
[in]coaddWcsWCS that relates the coadd coordinate system to the sky.
[in]inputSensorsA catalog containing the WCSs, bounding boxes and polygons, coaddition weights (in a field called 'weight'), and TransmissionCurves of the sensor-level images that went into the coadd.
Exceptions
NotFoundErrorThrown if the 'weight' field does not exist in the schema.
InvalidParameterErrorThrown if one or more inputs do not have a TransmissionCurve or a Wcs (ValidPolygons may be null to indicate no spatial restrictions other than the bounding box).
Returns
a new TransmissionCurve object.

Definition at line 257 of file CoaddTransmissionCurve.cc.

258 {
259 return std::make_shared<CoaddTransmissionCurve>(coaddWcs, inputSensors);
260}
afw::table::Key< int > coaddWcs

◆ makeExposurePatch() [1/2]

template<typename ExposureT >
std::shared_ptr< ExposurePatch< ExposureT > > lsst::meas::algorithms::makeExposurePatch ( std::shared_ptr< ExposureT const >  exp,
afw::detection::Footprint const &  standardFoot,
geom::Point2D const &  standardCenter,
afw::geom::SkyWcs const &  standardWcs 
)

Definition at line 98 of file ExposurePatch.h.

99 {
100 return std::make_shared<ExposurePatch<ExposureT> >(exp, standardFoot, standardCenter, standardWcs);
101}

◆ makeExposurePatch() [2/2]

template<typename ExposureT >
std::shared_ptr< ExposurePatch< ExposureT > > lsst::meas::algorithms::makeExposurePatch ( std::shared_ptr< ExposureT const >  exp,
std::shared_ptr< afw::detection::Footprint const >  foot,
geom::Point2D const &  center 
)

Factory function for ExposurePatch.

Definition at line 92 of file ExposurePatch.h.

93 {
94 return std::make_shared<ExposurePatch<ExposureT> >(exp, foot, center);
95}

◆ makePsfCandidate()

template<typename PixelT >
std::shared_ptr< PsfCandidate< PixelT > > lsst::meas::algorithms::makePsfCandidate ( std::shared_ptr< afw::table::SourceRecord > const &  source,
std::shared_ptr< afw::image::Exposure< PixelT > >  image 
)

Return a PsfCandidate of the right sort.

Cf. std::make_pair

Parameters
sourceThe detected Source
imageThe image wherein lies the object

Definition at line 176 of file PsfCandidate.h.

180 {
181 return std::make_shared<PsfCandidate<PixelT>>(source, image);
182}
afw::table::Key< afw::table::Array< ImagePixelT > > image
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224

◆ PYBIND11_MODULE()

lsst::meas::algorithms::PYBIND11_MODULE ( psfCandidate  ,
mod   
)

Definition at line 77 of file psfCandidate.cc.

77 {
78 declarePsfCandidate<float>(mod, "F");
79}

◆ setSpatialParameters() [1/2]

void lsst::meas::algorithms::setSpatialParameters ( afw::math::Kernel kernel,
Eigen::VectorXd const &  vec 
)

Fit a Kernel's spatial variability from a set of stars.

Definition at line 521 of file SpatialModelPsf.cc.

521 {
522 int const nComponents = kernel->getNKernelParameters();
523 int const nSpatialParams = kernel->getNSpatialParameters();
524
525 assert(nComponents * nSpatialParams == vec.size());
526
527 std::vector<std::vector<double>> kCoeffs; // coefficients rearranged for Kernel
528 kCoeffs.reserve(nComponents);
529 for (int i = 0; i != nComponents; ++i) {
530 std::vector<double> spatialCoeffs(nSpatialParams);
531 for (int j = 0; j != nSpatialParams; ++j) {
532 spatialCoeffs[j] = vec[i * nSpatialParams + j];
533 }
534 kCoeffs.push_back(spatialCoeffs);
535 }
536
537 kernel->setSpatialParameters(kCoeffs);
538}
void setSpatialParameters(const std::vector< std::vector< double > > params)
Set the parameters of all spatial functions.
Definition: Kernel.cc:110

◆ setSpatialParameters() [2/2]

void lsst::meas::algorithms::setSpatialParameters ( afw::math::Kernel kernel,
std::vector< double > const &  coeffs 
)

Fit a Kernel's spatial variability from a set of stars.

Definition at line 500 of file SpatialModelPsf.cc.

500 {
501 int const nComponents = kernel->getNKernelParameters();
502 int const nSpatialParams = kernel->getNSpatialParameters();
503
504 assert(nComponents * nSpatialParams == static_cast<long>(coeffs.size()));
505
506 std::vector<std::vector<double>> kCoeffs; // coefficients rearranged for Kernel
507 kCoeffs.reserve(nComponents);
508 for (int i = 0; i != nComponents; ++i) {
509 kCoeffs.push_back(std::vector<double>(nSpatialParams));
510 std::copy(coeffs.begin() + i * nSpatialParams, coeffs.begin() + (i + 1) * nSpatialParams,
511 kCoeffs[i].begin());
512 }
513
514 kernel->setSpatialParameters(kCoeffs);
515}

◆ subtractPsf() [1/2]

template<typename ImageT >
double lsst::meas::algorithms::subtractPsf ( afw::detection::Psf const &  psf,
ImageT *  data,
double  x,
double  y,
double  psfFlux = std::numeric_limits< double >::quiet_NaN() 
)

◆ subtractPsf() [2/2]

template<typename MaskedImageT >
double lsst::meas::algorithms::subtractPsf ( afw::detection::Psf const &  psf,
MaskedImageT *  data,
double  x,
double  y,
double  psfFlux 
)

Subtract a PSF from an image at a given position.

Parameters
psfthe PSF to subtract
dataImage to subtract PSF from
xcolumn position
yrow position
psfFluxobject's PSF flux (if not NaN)

Definition at line 976 of file SpatialModelPsf.cc.

981 {
982 if (std::isnan(x + y)) {
984 }
985
986 //
987 // Get Psf candidate
988 //
990
991 //
992 // Now find the proper sub-Image
993 //
994 geom::BoxI bbox = kImage->getBBox();
995
997 new MaskedImageT(*data, bbox, afw::image::PARENT, false)); // shallow copy
998 //
999 // Now we've got both; find the PSF's amplitude
1000 //
1001 double lambda = 0.0; // floor for variance is lambda*data
1002 try {
1003 double chi2; // chi^2 for fit
1004 double amp; // estimate of amplitude of model at this point
1005
1006 if (std::isnan(psfFlux)) {
1007 std::pair<double, double> result = fitKernel(*kImage, *subData, lambda, true);
1008 chi2 = result.first; // chi^2 for fit
1009 amp = result.second; // estimate of amplitude of model at this point
1010 } else {
1012 amp = psfFlux / afw::math::makeStatistics(*kImage, afw::math::SUM).getValue();
1013 }
1014 //
1015 // Convert kImage to the proper type so that I can subtract it.
1016 //
1018 new typename MaskedImageT::Image(*kImage, true)); // of data's type
1019
1020 *kImageF *= amp;
1021 *subData->getImage() -= *kImageF;
1022
1023 return chi2;
1025 LSST_EXCEPT_ADD(e, (boost::format("Object at (%.2f, %.2f)") % x % y).str());
1026 throw e;
1027 }
1028}
py::object result
Definition: _schema.cc:429
char * data
Definition: BaseRecord.cc:61
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
T isnan(T... args)
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
@ SUM
find sum of pixels in the image
Definition: Statistics.h:77
T quiet_NaN(T... args)