LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
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>
433  typedef afw::image::MaskedImage<PixelT> MaskedImage;
434  typedef afw::image::Exposure<PixelT> Exposure;
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  )
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 
775  CONST_PTR(MaskedImage) data;
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  int const x0 = kernelImages[0]->getX0(), y0 = kernelImages[0]->getY0();
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->setCtrX(x0 + static_cast<int>(newKernel->getWidth() / 2));
1089  newKernel->setCtrY(y0 + static_cast<int>(newKernel->getHeight() / 2));
1090 
1091  params[i] = x[i];
1092  newKernels[i] = newKernel;
1093  }
1094 
1095  return std::make_pair(params, newKernels);
1096 }
1097 
1098 /************************************************************************************************************/
1104 template <typename Image>
1106  afw::math::LinearCombinationKernel const& kernel,
1107  Image const& image,
1108  geom::Point2D const& pos
1109  ) {
1111  fitKernelParamsToImage(kernel, image, pos);
1112  std::vector<double> params = fit.first;
1113  afw::math::KernelList kernels = fit.second;
1114  int const nKernel = params.size();
1115  assert(kernels.size() == static_cast<unsigned int>(nKernel));
1116 
1117  double amp = 0.0;
1118  for (int i = 0; i != nKernel; ++i) {
1121  amp += params[i] * k->getSum();
1122  }
1123 
1125  double chisq = 0.0;
1126  outputKernel->setCtrX(kernels[0]->getCtrX());
1127  outputKernel->setCtrY(kernels[0]->getCtrY());
1128 
1129  return std::make_pair(outputKernel, std::make_pair(amp, chisq));
1130 }
1131 
1132 /************************************************************************************************************/
1133 //
1134 // Explicit instantiations
1135 //
1137 typedef float Pixel;
1138 
1140 createKernelFromPsfCandidates<Pixel>(afw::math::SpatialCellSet const&, geom::Extent2I const&,
1141  geom::Point2I const&, int const, int const, int const, int const,
1142  bool const, int const);
1143 template int countPsfCandidates<Pixel>(afw::math::SpatialCellSet const&, int const);
1144 
1145 template std::pair<bool, double> fitSpatialKernelFromPsfCandidates<Pixel>(afw::math::Kernel*,
1146  afw::math::SpatialCellSet const&,
1147  int const, double const,
1148  double const);
1149 template std::pair<bool, double> fitSpatialKernelFromPsfCandidates<Pixel>(afw::math::Kernel*,
1150  afw::math::SpatialCellSet const&,
1151  bool const, int const, double const,
1152  double const);
1153 
1154 template double subtractPsf(afw::detection::Psf const&, afw::image::MaskedImage<float>*, double, double,
1155  double);
1156 
1159  geom::Point2D const&);
1160 
1163  geom::Point2D const&);
1165 
1166 } // namespace algorithms
1167 } // namespace meas
1168 } // namespace lsst
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
double subtractPsf(afw::detection::Psf const &psf, ImageT *data, double x, double y, double psfFlux=std::numeric_limits< double >::quiet_NaN())
Class used by SpatialCell for spatial PSF fittig.
void setStatus(Status status)
Set the candidate&#39;s status.
Definition: SpatialCell.cc:54
uint64_t * ptr
Definition: RangeSet.cc:88
T copy(T... args)
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...
double operator()(const std::vector< double > &coeffs) const
void setAmplitude(double amplitude)
Set the best-fit amplitude.
Definition: PsfCandidate.h:116
x_iterator fast_iterator
A fast STL compliant iterator for contiguous images N.b.
Definition: ImageBase.h:138
static void setHeight(int height)
Set the height of the image that getImage should return.
Definition: SpatialCell.h:145
int countPsfCandidates(afw::math::SpatialCellSet const &psfCells, int const nStarPerCell=-1)
Count the number of candidates in use.
A floating-point coordinate rectangle geometry.
Definition: Box.h:305
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition: ImageUtils.h:69
std::vector< double > const & getEigenValues() const
Return Eigen values.
Definition: ImagePca.h:99
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
A coordinate class intended to represent absolute positions.
A class to contain the data, WCS, and other information needed to describe an image of the sky...
Definition: Exposure.h:72
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition: ImagePca.cc:415
table::Key< int > b
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
_const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition: ImageBase.h:142
T endl(T... args)
float getXCenter() const
Return the object&#39;s column-centre.
Definition: SpatialCell.h:90
py::object result
Definition: schema.cc:418
int y
Definition: SpanSet.cc:49
STL namespace.
virtual void analyze()
Generate eigenimages that are normalised and background-subtracted.
Definition: ImagePca.cc:41
estimate sample maximum
Definition: Statistics.h:77
std::shared_ptr< Image > computeImage(lsst::geom::Point2D position=makeNullPoint(), image::Color color=image::Color(), ImageOwnerEnum owner=COPY) const
Return an Image of the PSF, in a form that can be compared directly with star images.
Definition: Psf.cc:103
int min
boost::shared_ptr< afw::table::SourceRecord > getSource() const
Return the original Source.
Definition: PsfCandidate.h:110
STL class.
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.
T push_back(T... args)
double getVar() const
Return the variance in use when fitting this object.
Definition: PsfCandidate.h:119
void setSpatialParameters(afw::math::Kernel *kernel, std::vector< double > const &coeffs)
Fit a Kernel&#39;s spatial variability from a set of stars.
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:387
def show(frame=None)
Definition: ds9.py:89
double getAmplitude() const
Return the best-fit amplitude.
Definition: PsfCandidate.h:113
A class to pass around to all our PsfCandidates to evaluate the PSF fit&#39;s X^2.
A base class for image defects.
Class used by SpatialCell for spatial PSF fittig.
char * data
Definition: BaseRecord.cc:62
Class for doing PCA on PSF stars.
bool isValid
Definition: fits.cc:380
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:354
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:741
boost::shared_ptr< afw::image::MaskedImage< PixelT > const > getMaskedImage() const
Return the image at the position of the Source, without any sub-pixel shifts to put the centre of the...
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74
T make_pair(T... args)
T isfinite(T... args)
T static_pointer_cast(T... args)
static MaskPixelT getPlaneBitMask(const std::vector< std::string > &names)
Return the bitmask corresponding to a vector of plane names OR&#39;d together.
Definition: Mask.cc:383
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
table::Box2IKey bbox
Definition: Detector.cc:169
A coordinate class intended to represent offsets and dimensions.
solver_t * s
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
double x
unsigned int getNKernelParameters() const
Return the number of kernel parameters (0 if none)
Definition: Kernel.h:264
void setSpatialParameters(const std::vector< std::vector< double >> params)
Set the parameters of all spatial functions.
Definition: Kernel.cc:128
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< MaskPixelT > > mask
T assign(T... args)
STL class.
static void setWidth(int width)
Set the width of the image that getImage should return.
Definition: SpatialCell.h:140
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
T begin(T... args)
float getYCenter() const
Return the object&#39;s row-centre.
Definition: SpatialCell.h:93
std::vector< std::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:500
void processCandidate(afw::math::SpatialCellCandidate *candidate)
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...
boost::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.
double Up() const
Error definition of the function.
int getNSpatialParameters() const
Return the number of spatial parameters (0 if not spatially varying)
Definition: Kernel.h:269
MinimizeChi2(evalChi2Visitor< PixelT > &chi2Visitor, afw::math::Kernel *kernel, afw::math::SpatialCellSet const &psfCells, int nStarPerCell, int nComponents, int nSpatialParams)
T isnan(T... args)
virtual KernelList const & getKernelList() const
Get the fixed basis kernels.
int m
Definition: SpanSet.cc:49
ImageList const & getEigenImages() const
Return Eigen images.
Definition: ImagePca.h:101
void visitAllCandidates(CandidateVisitor *visitor, bool const ignoreExceptions=false)
Call the visitor&#39;s processCandidate method for every Candidate in the SpatialCellSet.
Definition: SpatialCell.cc:399
T quiet_NaN(T... args)
Base class for candidate objects in a SpatialCell.
Definition: SpatialCell.h:70
lsst::afw::image::Image< ImagePixelT > Image
Definition: MaskedImage.h:86
find sum of pixels in the image
Definition: Statistics.h:78
T fill(T... args)
T accumulate(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:76
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it...
Definition: Exception.h:54
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:112
Key< int > psf
Definition: Exposure.cc:65
An integer coordinate rectangle.
Definition: Box.h:54
void visitCandidates(CandidateVisitor *visitor, int const nMaxPerCell=-1, bool const ignoreExceptions=false)
Call the visitor&#39;s processCandidate method for each Candidate in the SpatialCellSet.
Definition: SpatialCell.cc:380
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:59
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
int end
Class stored in SpatialCells for spatial Psf fitting.
Definition: PsfCandidate.h:57
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...
void setChi2(double chi2)
Set the candidate&#39;s chi^2.
Definition: SpatialCell.h:152
virtual double updateBadPixels(unsigned long mask, int const ncomp)
Update the bad pixels (i.e.
Definition: ImagePca.cc:387
A kernel created from an Image.
Definition: Kernel.h:509
T reserve(T... args)
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
evalChi2Visitor(afw::math::Kernel const &kernel, double lambda)