LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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 
55 namespace lsst {
56 namespace meas {
57 namespace algorithms {
58 
59 namespace {
60 
61 int const WARP_BUFFER(1); // Buffer (border) around kernel to prevent warp issues
62 std::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
65 template <typename PixelT>
66 class 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 
71 public:
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 
116 private:
117  PsfImagePca<MaskedImageT>* _imagePca; // the ImagePca we're building
118 };
119 
120 /************************************************************************************************************/
122 template <typename PixelT>
123 class countVisitor : public afw::math::CandidateVisitor {
124  typedef afw::image::MaskedImage<PixelT> MaskedImage;
125  typedef afw::image::Exposure<PixelT> Exposure;
126 
127 public:
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 
152 private:
153  int mutable _n; // the desired number
154 };
155 
160 template <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 /************************************************************************************************************/
191 template <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 /************************************************************************************************************/
336 template <typename PixelT>
337 int 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 /************************************************************************************************************/
345 namespace {
351 template <typename ModelImageT, typename DataImageT>
352 std::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  */
430 template <typename PixelT>
435 
437 
438 public:
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;
478  } catch (lsst::pex::exceptions::RangeError& e) {
482  }
483  }
484 
485  // Return the computed chi^2
486  double getValue() const { return _chi2; }
487 
488 private:
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
521 void 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 //
543 template <typename PixelT>
544 class MinimizeChi2 : public ROOT::Minuit2::FCNBase {
545 public:
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 
576 private:
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 /************************************************************************************************************/
591 template <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 /************************************************************************************************************/
683 namespace {
720 template <typename PixelT>
721 class FillABVisitor : public afw::math::CandidateVisitor {
722  typedef afw::image::Image<PixelT> Image;
723  typedef afw::image::MaskedImage<PixelT> MaskedImage;
724  typedef afw::image::Exposure<PixelT> Exposure;
725 
727 
728 public:
729  explicit FillABVisitor(afw::math::LinearCombinationKernel const& kernel, // the Kernel we're fitting
730  double tau2 = 0.0 // floor to the per-candidate variance
731  )
732  : afw::math::CandidateVisitor(),
733  _kernel(kernel),
734  _tau2(tau2),
735  _nSpatialParams(_kernel.getNSpatialParameters()),
736  _nComponents(_kernel.getNKernelParameters()),
737  _basisImgs(),
738  _A((_nComponents - 1) * _nSpatialParams, (_nComponents - 1) * _nSpatialParams),
739  _b((_nComponents - 1) * _nSpatialParams),
740  _basisDotBasis(_nComponents, _nComponents) {
741  _basisImgs.resize(_nComponents);
742 
743  _A.setZero();
744  _b.setZero();
745  //
746  // Get all the Kernel's components as Images
747  //
748  afw::math::KernelList const& kernels = _kernel.getKernelList(); // Kernel's components
749  for (int i = 0; i != _nComponents; ++i) {
750  _basisImgs[i] = std::shared_ptr<KImage>(new KImage(kernels[i]->getDimensions()));
751  kernels[i]->computeImage(*_basisImgs[i], false);
752  }
753 
754  //
755  // Calculate the inner products of the Kernel components once and for all
756  //
757  for (int i = 1; i != _nComponents; ++i) { // Don't need 0th component
758  for (int j = i; j != _nComponents; ++j) {
759  _basisDotBasis(i, j) = _basisDotBasis(j, i) = afw::image::innerProduct(
760  *_basisImgs[i], *_basisImgs[j], PsfCandidate<PixelT>::getBorderWidth());
761  }
762  }
763  }
764 
765  void reset() {}
766 
767  // Called by SpatialCellSet::visitCandidates for each Candidate
768  void processCandidate(afw::math::SpatialCellCandidate* candidate) {
769  PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
770  if (imCandidate == NULL) {
772  "Failed to cast SpatialCellCandidate to PsfCandidate");
773  }
774 
776  try {
777  data = imCandidate->getMaskedImage(_kernel.getWidth(), _kernel.getHeight());
779  return;
780  }
781  double const xcen = imCandidate->getXCenter();
782  double const ycen = imCandidate->getYCenter();
783  double const dx = afw::image::positionToIndex(xcen, true).second;
784  double const dy = afw::image::positionToIndex(ycen, true).second;
785 
786 #if 0
787  double const amp = imCandidate->getAmplitude();
788 #else
789  /*
790  * Estimate the amplitude based on the current basis functions.
791  *
792  * N.b. you have to be a little careful here. Consider a PSF that is phi == (N0 + b*y*N1)/(1 + b*y)
793  * where the amplitude of N0 and N1 is 1.0, so a star has profile I = A*(N0 + b*y*N1)/(1 + b*y)
794  *
795  * If we set the amplitude to be A = I(0)/phi(0) (i.e. the central value of the data and best-fit phi)
796  * then the coefficient of N0 becomes 1/(1 + b*y) which makes the model non-linear in y.
797  */
799  fitKernelToImage(_kernel, *data, geom::Point2D(xcen, ycen));
800  double const amp = ret.second.first;
801 #endif
802 
803  double const var = imCandidate->getVar();
804  double const ivar = 1 / (var + _tau2); // Allow for floor on variance
805 
806  // Spatial params of all the components
807  std::vector<std::vector<double>> params(_nComponents);
808  for (int ic = 1; ic != _nComponents; ++ic) { // Don't need params[0]
809  params[ic] = _kernel.getSpatialFunction(ic)->getDFuncDParameters(xcen, ycen);
810  }
811 
812  std::vector<std::shared_ptr<KImage>> basisImages = offsetKernel<KImage>(_kernel, dx, dy);
813 
814  // Prepare values for basis dot data
815  // Scale data and subtract 0th component as part of unit kernel sum construction
816  std::shared_ptr<Image> dataImage(new Image(*data->getImage(), true));
817  typename KImage::fast_iterator bPtr = basisImages[0]->begin(true);
818  for (typename Image::fast_iterator dPtr = dataImage->begin(true), end = dataImage->end(true);
819  dPtr != end; ++dPtr, ++bPtr) {
820  *dPtr = *dPtr / amp - *bPtr;
821  }
822 
823  for (int i = 0, ic = 1; ic != _nComponents; ++ic) { // Don't need 0th component now
824  double const basisDotData = afw::image::innerProduct(*basisImages[ic], *dataImage,
826  for (int is = 0; is != _nSpatialParams; ++is, ++i) {
827  _b(i) += ivar * params[ic][is] * basisDotData;
828 
829  for (int j = i, jc = ic; jc != _nComponents; ++jc) {
830  for (int js = (i == j) ? is : 0; js != _nSpatialParams; ++js, ++j) {
831  _A(i, j) += ivar * params[ic][is] * params[jc][js] * _basisDotBasis(ic, jc);
832  _A(j, i) = _A(i, j); // could do this after _A is fully calculated
833  }
834  }
835  }
836  }
837  }
838 
839  Eigen::MatrixXd const& getA() const { return _A; }
840  Eigen::VectorXd const& getB() const { return _b; }
841 
842 private:
843  afw::math::LinearCombinationKernel const& _kernel; // the kernel
844  double _tau2; // variance floor added in quadrature to true candidate variance
845  int const _nSpatialParams; // number of spatial parameters
846  int const _nComponents; // number of basis functions
847  std::vector<std::shared_ptr<KImage>> _basisImgs; // basis function images from _kernel
848  Eigen::MatrixXd _A; // We'll solve the matrix equation A x = b for the Kernel's coefficients
849  Eigen::VectorXd _b;
850  Eigen::MatrixXd _basisDotBasis; // the inner products of the Kernel components
851 };
852 
854 template <typename PixelT>
855 class setAmplitudeVisitor : public afw::math::CandidateVisitor {
856  typedef afw::image::MaskedImage<PixelT> MaskedImage;
857  typedef afw::image::Exposure<PixelT> Exposure;
858 
859 public:
860  // Called by SpatialCellSet::visitCandidates for each Candidate
861  void processCandidate(afw::math::SpatialCellCandidate* candidate) {
862  PsfCandidate<PixelT>* imCandidate = dynamic_cast<PsfCandidate<PixelT>*>(candidate);
863  if (imCandidate == NULL) {
865  "Failed to cast SpatialCellCandidate to PsfCandidate");
866  }
867  imCandidate->setAmplitude(
868  afw::math::makeStatistics(*imCandidate->getMaskedImage()->getImage(), afw::math::MAX)
869  .getValue());
870  }
871 };
872 
873 } // namespace
874 
875 template <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 /************************************************************************************************************/
975 template <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;
1024  } catch (lsst::pex::exceptions::RangeError& e) {
1025  LSST_EXCEPT_ADD(e, (boost::format("Object at (%.2f, %.2f)") % x % y).str());
1026  throw e;
1027  }
1028 }
1029 
1030 /************************************************************************************************************/
1036 template <typename Image>
1038  afw::math::LinearCombinationKernel const& kernel,
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 /************************************************************************************************************/
1103 template <typename Image>
1105  afw::math::LinearCombinationKernel const& kernel,
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 //
1135 typedef float Pixel;
1136 
1138 createKernelFromPsfCandidates<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);
1141 template int countPsfCandidates<Pixel>(afw::math::SpatialCellSet const&, int const);
1142 
1143 template std::pair<bool, double> fitSpatialKernelFromPsfCandidates<Pixel>(afw::math::Kernel*,
1145  int const, double const,
1146  double const);
1147 template std::pair<bool, double> fitSpatialKernelFromPsfCandidates<Pixel>(afw::math::Kernel*,
1149  bool const, int const, double const,
1150  double const);
1151 
1152 template 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
table::Key< std::string > name
Definition: Amplifier.cc:116
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
char * data
Definition: BaseRecord.cc:61
int min
int end
double x
#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.
uint64_t * ptr
Definition: RangeSet.cc:88
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:76
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, std::shared_ptr< lsst::daf::base::PropertySet const > metadata=std::shared_ptr< lsst::daf::base::PropertySet const >(), 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:386
std::vector< double > const & getEigenValues() const
Return Eigen values.
Definition: ImagePca.h:98
ImageList const & getEigenImages() const
Return Eigen images.
Definition: ImagePca.h:100
static MaskPixelT getPlaneBitMask(const std::vector< std::string > &names)
Return the bitmask corresponding to a vector of plane names OR'd together.
Definition: Mask.cc:372
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
A kernel created from an Image.
Definition: Kernel.h:471
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
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
void setSpatialParameters(const std::vector< std::vector< double >> params)
Set the parameters of all spatial functions.
Definition: Kernel.cc:110
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.
Definition: SpatialCell.cc:53
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
void setChi2(double chi2)
Set the candidate's chi^2.
Definition: SpatialCell.h:148
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:383
void visitAllCandidates(CandidateVisitor *visitor, bool const ignoreExceptions=false)
Call the visitor's processCandidate method for every Candidate in the SpatialCellSet.
Definition: SpatialCell.cc:398
void visitCandidates(CandidateVisitor *visitor, int const nMaxPerCell=-1, bool const ignoreExceptions=false)
Call the visitor's processCandidate method for each Candidate in the SpatialCellSet.
Definition: SpatialCell.cc:379
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
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.
Definition: PsfCandidate.h:55
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.
Definition: PsfCandidate.h:114
std::shared_ptr< afw::table::SourceRecord > getSource() const
Return the original Source.
Definition: PsfCandidate.h:108
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:399
T printf(T... args)
T isfinite(T... args)
T isnan(T... args)
T make_pair(T... args)
Class for doing PCA on PSF stars.
def show(frame=None)
Definition: ds9.py:88
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition: ImagePca.cc:414
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:359
@ MAX
estimate sample maximum
Definition: Statistics.h:76
@ SUM
find sum of pixels in the image
Definition: Statistics.h:77
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.
Definition: offsetImage.cc:41
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())
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
A base class for image defects.
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