LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
SpatialModelPsf.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008, 2009, 2010 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
32#include <numeric>
33
34#if !defined(DOXYGEN)
35#include "Minuit2/FCNBase.h"
36#include "Minuit2/FunctionMinimum.h"
37#include "Minuit2/MnMigrad.h"
38#include "Minuit2/MnMinos.h"
39#include "Minuit2/MnPrint.h"
40#endif
41
42#include "Eigen/Core"
43#include "Eigen/Cholesky"
44#include "Eigen/SVD"
45
46#include "lsst/geom/Box.h"
47#include "lsst/geom/Point.h"
54
55namespace lsst {
56namespace meas {
57namespace algorithms {
58
59namespace {
60
61int const WARP_BUFFER(1); // Buffer (border) around kernel to prevent warp issues
62std::string const WARP_ALGORITHM("lanczos5"); // Name of warping algorithm to use
63
64// A class to pass around to all our PsfCandidates which builds the PcaImageSet
65template <typename PixelT>
66class SetPcaImageVisitor : public afw::math::CandidateVisitor {
67 typedef afw::image::Image<PixelT> ImageT;
68 typedef afw::image::MaskedImage<PixelT> MaskedImageT;
69 typedef afw::image::Exposure<PixelT> ExposureT;
70
71public:
72 explicit SetPcaImageVisitor(PsfImagePca<MaskedImageT>* imagePca, // Set of Images to initialise
73 unsigned int const mask = 0x0 // Ignore pixels with any of these bits set
74 )
75 : afw::math::CandidateVisitor(), _imagePca(imagePca) {
76 ;
77 }
78
79 // Called by SpatialCellSet::visitCandidates for each Candidate
80 void processCandidate(afw::math::SpatialCellCandidate* candidate) {
81 PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
82 if (imCandidate == NULL) {
84 "Failed to cast SpatialCellCandidate to PsfCandidate");
85 }
86
87 try {
88 std::shared_ptr<MaskedImageT> im = imCandidate->getOffsetImage(WARP_ALGORITHM, WARP_BUFFER);
89
90 // static int count = 0;
91 // im->writeFits(str(boost::format("cand%03d.fits") % count));
92 // count += 1;
93
94 afw::math::StatisticsControl sctrl;
95 sctrl.setNanSafe(false);
96
97 if (!std::isfinite(
98 afw::math::makeStatistics(*im->getImage(), afw::math::MAX, sctrl).getValue())) {
100 str(boost::format("Image at %d, %d contains NaN") %
101 imCandidate->getXCenter() % imCandidate->getYCenter()));
102 }
103 if (!std::isfinite(
104 afw::math::makeStatistics(*im->getVariance(), afw::math::MAX, sctrl).getValue())) {
106 str(boost::format("Variance of Image at %d, %d contains NaN") %
107 imCandidate->getXCenter() % imCandidate->getYCenter()));
108 }
109
110 _imagePca->addImage(im, imCandidate->getSource()->getPsfInstFlux());
112 return;
113 }
114 }
115
116private:
117 PsfImagePca<MaskedImageT>* _imagePca; // the ImagePca we're building
118};
119
120/************************************************************************************************************/
122template <typename PixelT>
123class countVisitor : public afw::math::CandidateVisitor {
124 typedef afw::image::MaskedImage<PixelT> MaskedImage;
125 typedef afw::image::Exposure<PixelT> Exposure;
126
127public:
128 explicit countVisitor() : afw::math::CandidateVisitor(), _n(0) {}
129
130 void reset() { _n = 0; }
131
132 // Called by SpatialCellSet::visitCandidates for each Candidate
133 void processCandidate(afw::math::SpatialCellCandidate* candidate) {
134 PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
135 if (imCandidate == NULL) {
137 "Failed to cast SpatialCellCandidate to PsfCandidate");
138 }
139
140 try {
141 imCandidate->getMaskedImage();
143 return;
144 }
145
146 ++_n;
147 }
148
149 // Return the number
150 double getN() const { return _n; }
151
152private:
153 int mutable _n; // the desired number
154};
155
160template <typename ImageT>
162 afw::math::LinearCombinationKernel const& kernel,
163 float dx, float dy
164 ) {
165 afw::math::KernelList kernels = kernel.getKernelList(); // The Kernels that kernel adds together
166 unsigned int const nKernel = kernels.size(); // Number of kernel components
167 std::vector<std::shared_ptr<ImageT>> kernelImages(nKernel); // Images of each Kernel in kernels
168 if (nKernel == 0) {
169 throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Kernel has no components");
170 }
171
172 ImageT scratch(kernel.getDimensions()); // Buffered scratch space
173 for (unsigned int i = 0; i != nKernel; ++i) {
174 kernels[i]->computeImage(scratch, false);
175 kernelImages[i] = afw::math::offsetImage(scratch, dx, dy, WARP_ALGORITHM, WARP_BUFFER);
176 }
177
178 return kernelImages;
179}
180
181} // Anonymous namespace
182
183/************************************************************************************************************/
191template <typename PixelT>
194 afw::math::SpatialCellSet const& psfCells,
195 lsst::geom::Extent2I const& dims,
196 lsst::geom::Point2I const& xy0,
197 int const nEigenComponents,
198 int const spatialOrder,
199 int const ksize,
200 int const nStarPerCell,
201 bool const constantWeight,
202 int const border
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
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
317 afw::image::Image<afw::math::Kernel::Pixel>(*eigenImages[i]->getImage(), true))));
318
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}
331
332/************************************************************************************************************/
336template <typename PixelT>
337int countPsfCandidates(afw::math::SpatialCellSet const& psfCells, int const nStarPerCell) {
338 countVisitor<PixelT> counter;
339 psfCells.visitCandidates(&counter, nStarPerCell);
340
341 return counter.getN();
342}
343
344/************************************************************************************************************/
345namespace {
351template <typename ModelImageT, typename DataImageT>
352std::pair<double, double> fitKernel(ModelImageT const& mImage, // The model image at this point
353 DataImageT const& data, // the data to fit
354 double lambda = 0.0, // floor for variance is lambda*data
355 bool detected = true, // only fit DETECTED pixels?
356 int const id = -1 // ID for this object; useful in debugging
357 ) {
358 assert(data.getDimensions() == mImage.getDimensions());
359 assert(id == id);
360 int const DETECTED = afw::image::Mask<>::getPlaneBitMask("DETECTED");
362
363 double sumMM = 0.0, sumMD = 0.0, sumDD = 0.0; // sums of model*model/variance etc.
364 int npix = 0; // number of pixels used to evaluate chi^2
365 for (int y = 0; y != data.getHeight(); ++y) {
366 typename ModelImageT::x_iterator mptr = mImage.row_begin(y);
367 for (typename DataImageT::x_iterator ptr = data.row_begin(y), end = data.row_end(y); ptr != end;
368 ++ptr, ++mptr) {
369 double const m = (*mptr)[0]; // value of model
370 double const d = ptr.image(); // value of data
371 double const var = ptr.variance() + lambda * d; // data's variance
372 if (detected && !(ptr.mask() & DETECTED)) {
373 continue;
374 }
375 if (ptr.mask() & BAD) {
376 continue;
377 }
378 if (var != 0.0) { // assume variance == 0 => infinity XXX
379 double const iVar = 1.0 / var;
380 npix++;
381 sumMM += m * m * iVar;
382 sumMD += m * d * iVar;
383 sumDD += d * d * iVar;
384 }
385 }
386 }
387
388 if (npix == 0) {
389 throw LSST_EXCEPT(lsst::pex::exceptions::RangeError, "No good pixels");
390 }
391 if (sumMM == 0.0) {
392 throw LSST_EXCEPT(lsst::pex::exceptions::RangeError, "sum(data*data)/var == 0");
393 }
394
395 double const amp = sumMD / sumMM; // estimate of amplitude of model at this point
396 double const chi2 = (sumDD - 2 * amp * sumMD + amp * amp * sumMM) / (npix - 1);
397
398#if 0
399 bool show = false; // Display the centre of the image; set from gdb
400
401 if (show) {
402 show = true; // you can jump here in gdb to set show if direct attempts fail
403 int y = data.getHeight()/2;
404 int x = data.getWidth()/2;
405 int hsize = 2;
406 printf("\ndata ");
407 for (int ii = -hsize; ii <= hsize; ++ii) {
408 for (int jj = -hsize; jj <= hsize; ++jj) {
409 printf("%7.1f ", data.at(x + jj, y - ii).image());
410 }
411 printf(" model ");
412 for (int jj = -hsize; jj <= hsize; ++jj) {
413 printf("%7.1f ", amp*(*(mImage.at(x + jj, y - ii)))[0]);
414 }
415 printf("\n ");
416 }
417 printf("%g %.1f\n", amp, chi2);
418 }
419#endif
420
421 return std::make_pair(chi2, amp);
422}
423} // namespace
424
425/************************************************************************************************************/
426/*
427 * Fit for the spatial variation of the PSF parameters over the field
428 */
430template <typename PixelT>
435
437
438public:
439 explicit evalChi2Visitor(afw::math::Kernel const& kernel, double lambda)
440 : afw::math::CandidateVisitor(),
441 _chi2(0.0),
442 _kernel(kernel),
443 _lambda(lambda),
444 _kImage(std::shared_ptr<KImage>(new KImage(kernel.getDimensions()))) {}
445
446 void reset() { _chi2 = 0.0; }
447
448 // Called by SpatialCellSet::visitCandidates for each Candidate
450 PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
451 if (imCandidate == NULL) {
453 "Failed to cast SpatialCellCandidate to PsfCandidate");
454 }
455
456 double const xcen = imCandidate->getSource()->getX();
457 double const ycen = imCandidate->getSource()->getY();
458
459 _kernel.computeImage(*_kImage, true, xcen, ycen);
461 try {
462 data = imCandidate->getOffsetImage(WARP_ALGORITHM, WARP_BUFFER);
464 return;
465 }
466
467 try {
469 fitKernel(*_kImage, *data, _lambda, false, imCandidate->getSource()->getId());
470
471 double dchi2 = result.first; // chi^2 from this object
472 double const amp = result.second; // estimate of amplitude of model at this point
473
474 imCandidate->setChi2(dchi2);
475 imCandidate->setAmplitude(amp);
476
477 _chi2 += dchi2;
482 }
483 }
484
485 // Return the computed chi^2
486 double getValue() const { return _chi2; }
487
488private:
489 double mutable _chi2; // the desired chi^2
490 afw::math::Kernel const& _kernel; // the kernel
491 double _lambda; // floor for variance is _lambda*data
492 std::shared_ptr<KImage> mutable _kImage; // The Kernel at this point; a scratch copy
493};
494
495/********************************************************************************************************/
499// Set the Kernel's spatial parameters from a vector
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}
516
520// Set the Kernel's spatial parameters from an Eigen::VectorXd
521void setSpatialParameters(afw::math::Kernel* kernel, Eigen::VectorXd const& vec) {
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}
539
540//
541// The object that minuit minimises
542//
543template <typename PixelT>
544class MinimizeChi2 : public ROOT::Minuit2::FCNBase {
545public:
547 afw::math::SpatialCellSet const& psfCells, int nStarPerCell, int nComponents,
548 int nSpatialParams)
549 : _errorDef(1.0),
550 _chi2Visitor(chi2Visitor),
551 _kernel(kernel),
552 _psfCells(psfCells),
553 _nStarPerCell(nStarPerCell),
554 _nComponents(nComponents),
555 _nSpatialParams(nSpatialParams) {}
556
563 double Up() const { return _errorDef; }
564
565 // Evaluate our cost function (in this case chi^2)
566 double operator()(const std::vector<double>& coeffs) const {
567 setSpatialParameters(_kernel, coeffs);
568
569 _psfCells.visitCandidates(&_chi2Visitor, _nStarPerCell);
570
571 return _chi2Visitor.getValue();
572 }
573
574 void setErrorDef(double def) { _errorDef = def; }
575
576private:
577 double _errorDef; // how much cost function has changed at the +- 1 error points
578
579 evalChi2Visitor<PixelT>& _chi2Visitor;
580 afw::math::Kernel* _kernel;
581 afw::math::SpatialCellSet const& _psfCells;
582 int _nStarPerCell;
583 int _nComponents;
584 int _nSpatialParams;
585};
586
587/************************************************************************************************************/
591template <typename PixelT>
593 afw::math::Kernel* kernel,
594 afw::math::SpatialCellSet const& psfCells,
595 int const nStarPerCell,
596 double const tolerance,
597 double const lambda
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}
678
679/************************************************************************************************************/
680namespace {
717template <typename PixelT>
718class FillABVisitor : public afw::math::CandidateVisitor {
720 typedef afw::image::MaskedImage<PixelT> MaskedImage;
721 typedef afw::image::Exposure<PixelT> Exposure;
722
724
725public:
726 explicit FillABVisitor(afw::math::LinearCombinationKernel const& kernel, // the Kernel we're fitting
727 double tau2 = 0.0 // floor to the per-candidate variance
728 )
729 : afw::math::CandidateVisitor(),
730 _kernel(kernel),
731 _tau2(tau2),
732 _nSpatialParams(_kernel.getNSpatialParameters()),
733 _nComponents(_kernel.getNKernelParameters()),
734 _basisImgs(),
735 _A((_nComponents - 1) * _nSpatialParams, (_nComponents - 1) * _nSpatialParams),
736 _b((_nComponents - 1) * _nSpatialParams),
737 _basisDotBasis(_nComponents, _nComponents) {
738 _basisImgs.resize(_nComponents);
739
740 _A.setZero();
741 _b.setZero();
742 //
743 // Get all the Kernel's components as Images
744 //
745 afw::math::KernelList const& kernels = _kernel.getKernelList(); // Kernel's components
746 for (int i = 0; i != _nComponents; ++i) {
747 _basisImgs[i] = std::shared_ptr<KImage>(new KImage(kernels[i]->getDimensions()));
748 kernels[i]->computeImage(*_basisImgs[i], false);
749 }
750
751 //
752 // Calculate the inner products of the Kernel components once and for all
753 //
754 for (int i = 1; i != _nComponents; ++i) { // Don't need 0th component
755 for (int j = i; j != _nComponents; ++j) {
756 _basisDotBasis(i, j) = _basisDotBasis(j, i) = afw::image::innerProduct(
757 *_basisImgs[i], *_basisImgs[j], PsfCandidate<PixelT>::getBorderWidth());
758 }
759 }
760 }
761
762 void reset() {}
763
764 // Called by SpatialCellSet::visitCandidates for each Candidate
765 void processCandidate(afw::math::SpatialCellCandidate* candidate) {
766 PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
767 if (imCandidate == NULL) {
769 "Failed to cast SpatialCellCandidate to PsfCandidate");
770 }
771
773 try {
774 data = imCandidate->getMaskedImage(_kernel.getWidth(), _kernel.getHeight());
776 return;
777 }
778 double const xcen = imCandidate->getXCenter();
779 double const ycen = imCandidate->getYCenter();
780 double const dx = afw::image::positionToIndex(xcen, true).second;
781 double const dy = afw::image::positionToIndex(ycen, true).second;
782
783#if 0
784 double const amp = imCandidate->getAmplitude();
785#else
786 /*
787 * Estimate the amplitude based on the current basis functions.
788 *
789 * N.b. you have to be a little careful here. Consider a PSF that is phi == (N0 + b*y*N1)/(1 + b*y)
790 * where the amplitude of N0 and N1 is 1.0, so a star has profile I = A*(N0 + b*y*N1)/(1 + b*y)
791 *
792 * If we set the amplitude to be A = I(0)/phi(0) (i.e. the central value of the data and best-fit phi)
793 * then the coefficient of N0 becomes 1/(1 + b*y) which makes the model non-linear in y.
794 */
796 fitKernelToImage(_kernel, *data, geom::Point2D(xcen, ycen));
797 double const amp = ret.second.first;
798#endif
799
800 double const var = imCandidate->getVar();
801 double const ivar = 1 / (var + _tau2); // Allow for floor on variance
802
803 // Spatial params of all the components
804 std::vector<std::vector<double>> params(_nComponents);
805 for (int ic = 1; ic != _nComponents; ++ic) { // Don't need params[0]
806 params[ic] = _kernel.getSpatialFunction(ic)->getDFuncDParameters(xcen, ycen);
807 }
808
809 std::vector<std::shared_ptr<KImage>> basisImages = offsetKernel<KImage>(_kernel, dx, dy);
810
811 // Prepare values for basis dot data
812 // Scale data and subtract 0th component as part of unit kernel sum construction
813 std::shared_ptr<Image> dataImage(new Image(*data->getImage(), true));
814 typename KImage::fast_iterator bPtr = basisImages[0]->begin(true);
815 for (typename Image::fast_iterator dPtr = dataImage->begin(true), end = dataImage->end(true);
816 dPtr != end; ++dPtr, ++bPtr) {
817 *dPtr = *dPtr / amp - *bPtr;
818 }
819
820 for (int i = 0, ic = 1; ic != _nComponents; ++ic) { // Don't need 0th component now
821 double const basisDotData = afw::image::innerProduct(*basisImages[ic], *dataImage,
823 for (int is = 0; is != _nSpatialParams; ++is, ++i) {
824 _b(i) += ivar * params[ic][is] * basisDotData;
825
826 for (int j = i, jc = ic; jc != _nComponents; ++jc) {
827 for (int js = (i == j) ? is : 0; js != _nSpatialParams; ++js, ++j) {
828 _A(i, j) += ivar * params[ic][is] * params[jc][js] * _basisDotBasis(ic, jc);
829 _A(j, i) = _A(i, j); // could do this after _A is fully calculated
830 }
831 }
832 }
833 }
834 }
835
836 Eigen::MatrixXd const& getA() const { return _A; }
837 Eigen::VectorXd const& getB() const { return _b; }
838
839private:
840 afw::math::LinearCombinationKernel const& _kernel; // the kernel
841 double _tau2; // variance floor added in quadrature to true candidate variance
842 int const _nSpatialParams; // number of spatial parameters
843 int const _nComponents; // number of basis functions
844 std::vector<std::shared_ptr<KImage>> _basisImgs; // basis function images from _kernel
845 Eigen::MatrixXd _A; // We'll solve the matrix equation A x = b for the Kernel's coefficients
846 Eigen::VectorXd _b;
847 Eigen::MatrixXd _basisDotBasis; // the inner products of the Kernel components
848};
849
851template <typename PixelT>
852class setAmplitudeVisitor : public afw::math::CandidateVisitor {
853 typedef afw::image::MaskedImage<PixelT> MaskedImage;
854 typedef afw::image::Exposure<PixelT> Exposure;
855
856public:
857 // Called by SpatialCellSet::visitCandidates for each Candidate
858 void processCandidate(afw::math::SpatialCellCandidate* candidate) {
859 PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
860 if (imCandidate == NULL) {
862 "Failed to cast SpatialCellCandidate to PsfCandidate");
863 }
864 imCandidate->setAmplitude(
865 afw::math::makeStatistics(*imCandidate->getMaskedImage()->getImage(), afw::math::MAX)
866 .getValue());
867 }
868};
869
870} // namespace
871
875template <typename PixelT>
877 afw::math::Kernel* kernel,
878 afw::math::SpatialCellSet const& psfCells,
879 bool const doNonLinearFit,
880 int const nStarPerCell,
881 double const tolerance,
882 double const lambda
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}
970
971/************************************************************************************************************/
975template <typename MaskedImageT>
977 MaskedImageT* data,
978 double x,
979 double y,
980 double psfFlux
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}
1029
1030/************************************************************************************************************/
1036template <typename Image>
1039 Image const& image,
1040 geom::Point2D const& pos
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}
1096
1097/************************************************************************************************************/
1103template <typename Image>
1106 Image const& image,
1107 geom::Point2D const& pos
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}
1129
1130/************************************************************************************************************/
1131//
1132// Explicit instantiations
1133//
1135typedef float Pixel;
1136
1138createKernelFromPsfCandidates<Pixel>(afw::math::SpatialCellSet const&, geom::Extent2I const&,
1139 geom::Point2I const&, int const, int const, int const, int const,
1140 bool const, int const);
1141template int countPsfCandidates<Pixel>(afw::math::SpatialCellSet const&, int const);
1142
1143template std::pair<bool, double> fitSpatialKernelFromPsfCandidates<Pixel>(afw::math::Kernel*,
1145 int const, double const,
1146 double const);
1147template std::pair<bool, double> fitSpatialKernelFromPsfCandidates<Pixel>(afw::math::Kernel*,
1149 bool const, int const, double const,
1150 double const);
1151
1152template double subtractPsf(afw::detection::Psf const&, afw::image::MaskedImage<float>*, double, double,
1153 double);
1154
1157 geom::Point2D const&);
1158
1161 geom::Point2D const&);
1163
1164} // namespace algorithms
1165} // namespace meas
1166} // namespace lsst
py::object result
Definition _schema.cc:429
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
char * data
Definition BaseRecord.cc:61
int min
int end
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition Exception.h:54
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< MaskPixelT > > mask
Class used by SpatialCell for spatial PSF fittig.
std::uint64_t * ptr
Definition RangeSet.cc:95
int y
Definition SpanSet.cc:48
int m
Definition SpanSet.cc:48
Class used by SpatialCell for spatial PSF fittig.
table::Key< int > b
T accumulate(T... args)
T assign(T... args)
T begin(T... args)
A polymorphic base class for representing an image's Point Spread Function.
Definition Psf.h:82
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
x_iterator fast_iterator
A fast STL compliant iterator for contiguous images N.b.
Definition ImageBase.h:137
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
void writeFits(std::string const &fileName, daf::base::PropertySet const *metadata=nullptr, std::string const &mode="w") const
Write an image to a regular FITS file.
virtual double updateBadPixels(unsigned long mask, int const ncomp)
Update the bad pixels (i.e.
Definition ImagePca.cc:385
ImageList const & getEigenImages() const
Return Eigen images.
Definition ImagePca.h:100
std::vector< double > const & getEigenValues() const
Return Eigen values.
Definition ImagePca.h:98
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:82
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
A kernel created from an Image.
Definition Kernel.h:471
void setParameter(unsigned int ind, double newValue)
Set one function parameter without range checking.
Definition Function.h:143
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition Kernel.h:110
unsigned int getNKernelParameters() const
Return the number of kernel parameters (0 if none)
Definition Kernel.h:246
void setSpatialParameters(const std::vector< std::vector< double > > params)
Set the parameters of all spatial functions.
Definition Kernel.cc:110
int getNSpatialParameters() const
Return the number of spatial parameters (0 if not spatially varying)
Definition Kernel.h:251
double computeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize, double x=0.0, double y=0.0) const
Compute an image (pixellized representation of the kernel) in place.
Definition Kernel.cc:76
A kernel that is a linear combination of fixed basis kernels.
Definition Kernel.h:704
virtual KernelList const & getKernelList() const
Get the fixed basis kernels.
Base class for candidate objects in a SpatialCell.
Definition SpatialCell.h:70
void setStatus(Status status)
Set the candidate's status.
static void setWidth(int width)
Set the width of the image that getImage should return.
static void setHeight(int height)
Set the height of the image that getImage should return.
void setChi2(double chi2)
Set the candidate's chi^2.
A collection of SpatialCells covering an entire image.
void visitAllCandidates(CandidateVisitor *visitor, bool const ignoreExceptions=false)
Call the visitor's processCandidate method for every Candidate in the SpatialCellSet.
void visitCandidates(CandidateVisitor *visitor, int const nMaxPerCell=-1, bool const ignoreExceptions=false)
Call the visitor's processCandidate method for each Candidate in the SpatialCellSet.
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
A floating-point coordinate rectangle geometry.
Definition Box.h:413
An integer coordinate rectangle.
Definition Box.h:55
double operator()(const std::vector< double > &coeffs) const
MinimizeChi2(evalChi2Visitor< PixelT > &chi2Visitor, afw::math::Kernel *kernel, afw::math::SpatialCellSet const &psfCells, int nStarPerCell, int nComponents, int nSpatialParams)
double Up() const
Error definition of the function.
Class stored in SpatialCells for spatial Psf fitting.
static int getBorderWidth()
Return the number of pixels being ignored around the candidate image's edge.
void setAmplitude(double amplitude)
Set the best-fit amplitude.
std::shared_ptr< afw::table::SourceRecord > getSource() const
Return the original Source.
std::shared_ptr< afw::image::MaskedImage< PixelT > > getOffsetImage(std::string const algorithm, unsigned int buffer) const
Return an offset version of the image of the source.
virtual void analyze()
Generate eigenimages that are normalised and background-subtracted.
Definition ImagePca.cc:41
A class to pass around to all our PsfCandidates to evaluate the PSF fit's X^2.
evalChi2Visitor(afw::math::Kernel const &kernel, double lambda)
void processCandidate(afw::math::SpatialCellCandidate *candidate)
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports when the result of an operation cannot be represented by the destination type.
Definition Runtime.h:115
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T copy(T... args)
T endl(T... args)
T fill(T... args)
bool isValid
Definition fits.cc:404
T printf(T... args)
T isfinite(T... args)
T isnan(T... args)
T make_pair(T... args)
Class for doing PCA on PSF stars.
Definition __init__.py:1
show(frame=None)
Definition ds9.py:88
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition ImagePca.cc:413
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition ImageUtils.h:69
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:361
@ MAX
estimate sample maximum
Definition Statistics.h:67
@ SUM
find sum of pixels in the image
Definition Statistics.h:68
std::shared_ptr< ImageT > offsetImage(ImageT const &image, float dx, float dy, std::string const &algorithmName="lanczos5", unsigned int buffer=0)
Return an image offset by (dx, dy) using the specified algorithm.
std::vector< std::shared_ptr< Kernel > > KernelList
Definition Kernel.h:462
std::pair< std::shared_ptr< afw::math::LinearCombinationKernel >, std::vector< double > > createKernelFromPsfCandidates(afw::math::SpatialCellSet const &psfCells, geom::Extent2I const &dims, geom::Point2I const &xy0, int const nEigenComponents, int const spatialOrder, int const ksize, int const nStarPerCell=-1, bool const constantWeight=true, int const border=3)
Return a Kernel pointer and a list of eigenvalues resulting from analysing the provided SpatialCellSe...
int countPsfCandidates(afw::math::SpatialCellSet const &psfCells, int const nStarPerCell=-1)
Count the number of candidates in use.
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.
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.
std::pair< bool, double > 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.
void setSpatialParameters(afw::math::Kernel *kernel, std::vector< double > const &coeffs)
Fit a Kernel's spatial variability from a set of stars.
double subtractPsf(afw::detection::Psf const &psf, ImageT *data, double x, double y, double psfFlux=std::numeric_limits< double >::quiet_NaN())
STL namespace.
T push_back(T... args)
T quiet_NaN(T... args)
T reserve(T... args)
T size(T... args)
Key< int > psf
Definition Exposure.cc:65
g2d::python::Image< double > Image
Definition test_image.cc:14