LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
KernelSolution.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
11 #include <iterator>
12 #include <cmath>
13 #include <algorithm>
14 #include <limits>
15 
16 #include <memory>
17 #include "boost/timer.hpp"
18 
19 #include "Eigen/Core"
20 #include "Eigen/Cholesky"
21 #include "Eigen/QR"
22 #include "Eigen/LU"
23 #include "Eigen/Eigenvalues"
24 #include "Eigen/SVD"
25 
26 #include "lsst/afw/math.h"
27 #include "lsst/afw/geom.h"
28 #include "lsst/afw/image.h"
29 #include "lsst/afw/detection.h"
30 #include "lsst/geom.h"
31 #include "lsst/log/Log.h"
33 
36 
37 #include "ndarray.h"
38 #include "ndarray/eigen.h"
39 
40 #define DEBUG_MATRIX 0
41 #define DEBUG_MATRIX2 0
42 
43 namespace geom = lsst::geom;
44 namespace afwDet = lsst::afw::detection;
45 namespace afwMath = lsst::afw::math;
46 namespace afwGeom = lsst::afw::geom;
47 namespace afwImage = lsst::afw::image;
49 
50 namespace lsst {
51 namespace ip {
52 namespace diffim {
53 
54  /* Unique identifier for solution */
56 
58  Eigen::MatrixXd mMat,
59  Eigen::VectorXd bVec,
60  bool fitForBackground
61  ) :
62  _id(++_SolutionId),
63  _mMat(mMat),
64  _bVec(bVec),
65  _aVec(),
66  _solvedBy(NONE),
67  _fitForBackground(fitForBackground)
68  {};
69 
71  bool fitForBackground
72  ) :
73  _id(++_SolutionId),
74  _mMat(),
75  _bVec(),
76  _aVec(),
77  _solvedBy(NONE),
78  _fitForBackground(fitForBackground)
79  {};
80 
82  _id(++_SolutionId),
83  _mMat(),
84  _bVec(),
85  _aVec(),
86  _solvedBy(NONE),
87  _fitForBackground(true)
88  {};
89 
91  solve(_mMat, _bVec);
92  }
93 
95  return getConditionNumber(_mMat, conditionType);
96  }
97 
98  double KernelSolution::getConditionNumber(Eigen::MatrixXd const& mMat,
99  ConditionNumberType conditionType) {
100  switch (conditionType) {
101  case EIGENVALUE:
102  {
103  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
104  Eigen::VectorXd eValues = eVecValues.eigenvalues();
105  double eMax = eValues.maxCoeff();
106  double eMin = eValues.minCoeff();
107  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
108  "EIGENVALUE eMax / eMin = %.3e", eMax / eMin);
109  return (eMax / eMin);
110  break;
111  }
112  case SVD:
113  {
114  Eigen::VectorXd sValues = mMat.jacobiSvd().singularValues();
115  double sMax = sValues.maxCoeff();
116  double sMin = sValues.minCoeff();
117  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
118  "SVD eMax / eMin = %.3e", sMax / sMin);
119  return (sMax / sMin);
120  break;
121  }
122  default:
123  {
125  "Undefined ConditionNumberType : only EIGENVALUE, SVD allowed.");
126  break;
127  }
128  }
129  }
130 
131  void KernelSolution::solve(Eigen::MatrixXd const& mMat,
132  Eigen::VectorXd const& bVec) {
133 
134  if (DEBUG_MATRIX) {
135  std::cout << "M " << std::endl;
136  std::cout << mMat << std::endl;
137  std::cout << "B " << std::endl;
138  std::cout << bVec << std::endl;
139  }
140 
141  Eigen::VectorXd aVec = Eigen::VectorXd::Zero(bVec.size());
142 
143  boost::timer t;
144  t.restart();
145 
146  LOGL_DEBUG("TRACE2.ip.diffim.KernelSolution.solve",
147  "Solving for kernel");
148  _solvedBy = LU;
149  Eigen::FullPivLU<Eigen::MatrixXd> lu(mMat);
150  if (lu.isInvertible()) {
151  aVec = lu.solve(bVec);
152  } else {
153  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
154  "Unable to determine kernel via LU");
155  /* LAST RESORT */
156  try {
157 
159  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
160  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
161  Eigen::VectorXd eValues = eVecValues.eigenvalues();
162 
163  for (int i = 0; i != eValues.rows(); ++i) {
164  if (eValues(i) != 0.0) {
165  eValues(i) = 1.0/eValues(i);
166  }
167  }
168 
169  aVec = rMat * eValues.asDiagonal() * rMat.transpose() * bVec;
170  } catch (pexExcept::Exception& e) {
171 
172  _solvedBy = NONE;
173  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
174  "Unable to determine kernel via eigen-values");
175 
176  throw LSST_EXCEPT(pexExcept::Exception, "Unable to determine kernel solution");
177  }
178  }
179 
180  double time = t.elapsed();
181  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
182  "Compute time for matrix math : %.2f s", time);
183 
184  if (DEBUG_MATRIX) {
185  std::cout << "A " << std::endl;
186  std::cout << aVec << std::endl;
187  }
188 
189  _aVec = aVec;
190  }
191 
192  /*******************************************************************************************************/
193 
194  template <typename InputT>
196  lsst::afw::math::KernelList const& basisList,
197  bool fitForBackground
198  )
199  :
200  KernelSolution(fitForBackground),
201  _cMat(),
202  _iVec(),
203  _ivVec(),
204  _kernel(),
205  _background(0.0),
206  _kSum(0.0)
207  {
208  std::vector<double> kValues(basisList.size());
210  new afwMath::LinearCombinationKernel(basisList, kValues)
211  );
212  };
213 
214  template <typename InputT>
216  if (_solvedBy == KernelSolution::NONE) {
217  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
218  }
219  return _kernel;
220  }
221 
222  template <typename InputT>
224  if (_solvedBy == KernelSolution::NONE) {
225  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
226  }
228  new afwImage::Image<afwMath::Kernel::Pixel>(_kernel->getDimensions())
229  );
230  (void)_kernel->computeImage(*image, false);
231  return image;
232  }
233 
234  template <typename InputT>
236  if (_solvedBy == KernelSolution::NONE) {
237  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return background");
238  }
239  return _background;
240  }
241 
242  template <typename InputT>
244  if (_solvedBy == KernelSolution::NONE) {
245  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
246  }
247  return _kSum;
248  }
249 
250  template <typename InputT>
253  if (_solvedBy == KernelSolution::NONE) {
254  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
255  }
256 
257  return std::make_pair(_kernel, _background);
258  }
259 
260  template <typename InputT>
262  lsst::afw::image::Image<InputT> const &templateImage,
263  lsst::afw::image::Image<InputT> const &scienceImage,
265  ) {
266 
267  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
268  if (varStats.getValue(afwMath::MIN) < 0.0) {
270  "Error: variance less than 0.0");
271  }
272  if (varStats.getValue(afwMath::MIN) == 0.0) {
274  "Error: variance equals 0.0, cannot inverse variance weight");
275  }
276 
277  lsst::afw::math::KernelList basisList =
278  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
279 
280  unsigned int const nKernelParameters = basisList.size();
281  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
282  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
283 
284  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
285 
286  /* Ignore buffers around edge of convolved images :
287  *
288  * If the kernel has width 5, it has center pixel 2. The first good pixel
289  * is the (5-2)=3rd pixel, which is array index 2, and ends up being the
290  * index of the central pixel.
291  *
292  * You also have a buffer of unusable pixels on the other side, numbered
293  * width-center-1. The last good usable pixel is N-width+center+1.
294  *
295  * Example : the kernel is width = 5, center = 2
296  *
297  * ---|---|-c-|---|---|
298  *
299  * the image is width = N
300  * convolve this with the kernel, and you get
301  *
302  * |-x-|-x-|-g-|---|---| ... |---|---|-g-|-x-|-x-|
303  *
304  * g = first/last good pixel
305  * x = bad
306  *
307  * the first good pixel is the array index that has the value "center", 2
308  * the last good pixel has array index N-(5-2)+1
309  * eg. if N = 100, you want to use up to index 97
310  * 100-3+1 = 98, and the loops use i < 98, meaning the last
311  * index you address is 97.
312  */
313 
314  /* NOTE - we are accessing particular elements of Eigen arrays using
315  these coordinates, therefore they need to be in LOCAL coordinates.
316  This was written before ndarray unification.
317  */
318  geom::Box2I goodBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
319  unsigned int const startCol = goodBBox.getMinX();
320  unsigned int const startRow = goodBBox.getMinY();
321  // endCol/Row is one past the index of the last good col/row
322  unsigned int endCol = goodBBox.getMaxX() + 1;
323  unsigned int endRow = goodBBox.getMaxY() + 1;
324 
325  boost::timer t;
326  t.restart();
327 
328  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
329  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
330  startCol,
331  endRow-startRow,
332  endCol-startCol);
333  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
334  startCol,
335  endRow-startRow,
336  endCol-startCol);
337  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
338  startRow, startCol, endRow-startRow, endCol-startCol
339  ).array().inverse().matrix();
340 
341  /* Resize into 1-D for later usage */
342  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
343  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
344  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
345 
346  /* Holds image convolved with basis function */
347  afwImage::Image<PixelT> cimage(templateImage.getDimensions());
348 
349  /* Holds eigen representation of image convolved with all basis functions */
350  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
351 
352  /* Iterators over convolved image list and basis list */
353  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
354  /* Create C_i in the formalism of Alard & Lupton */
355  afwMath::ConvolutionControl convolutionControl;
356  convolutionControl.setDoNormalize(false);
357  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
358  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
359 
360  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
361  startCol,
362  endRow-startRow,
363  endCol-startCol);
364  cMat.resize(cMat.size(), 1);
365  *eiter = cMat;
366 
367  }
368 
369  double time = t.elapsed();
370  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
371  "Total compute time to do basis convolutions : %.2f s", time);
372  t.restart();
373 
374  /*
375  Load matrix with all values from convolvedEigenList : all images
376  (eigeniVariance, convolvedEigenList) must be the same size
377  */
378  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
379  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
380  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
381  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
382  cMat.col(kidxj) = eiterj->col(0);
383  }
384  /* Treat the last "image" as all 1's to do the background calculation. */
385  if (_fitForBackground)
386  cMat.col(nParameters-1).fill(1.);
387 
388  _cMat = cMat;
389  _ivVec = eigeniVariance.col(0);
390  _iVec = eigenScience.col(0);
391 
392  /* Make these outside of solve() so I can check condition number */
393  _mMat = _cMat.transpose() * (_ivVec.asDiagonal() * _cMat);
394  _bVec = _cMat.transpose() * (_ivVec.asDiagonal() * _iVec);
395  }
396 
397  template <typename InputT>
399  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.solve",
400  "mMat is %d x %d; bVec is %d; cMat is %d x %d; vVec is %d; iVec is %d",
401  _mMat.rows(), _mMat.cols(), _bVec.size(),
402  _cMat.rows(), _cMat.cols(), _ivVec.size(), _iVec.size());
403 
404  if (DEBUG_MATRIX) {
405  std::cout << "C" << std::endl;
406  std::cout << _cMat << std::endl;
407  std::cout << "iV" << std::endl;
408  std::cout << _ivVec << std::endl;
409  std::cout << "I" << std::endl;
410  std::cout << _iVec << std::endl;
411  }
412 
413  try {
415  } catch (pexExcept::Exception &e) {
416  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
417  throw e;
418  }
419  /* Turn matrices into _kernel and _background */
420  _setKernel();
421  }
422 
423  template <typename InputT>
425  if (_solvedBy == KernelSolution::NONE) {
426  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot make solution");
427  }
428 
429  unsigned int const nParameters = _aVec.size();
430  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
431  unsigned int const nKernelParameters =
432  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList().size();
433  if (nParameters != (nKernelParameters + nBackgroundParameters))
434  throw LSST_EXCEPT(pexExcept::Exception, "Mismatched sizes in kernel solution");
435 
436  /* Fill in the kernel results */
437  std::vector<double> kValues(nKernelParameters);
438  for (unsigned int idx = 0; idx < nKernelParameters; idx++) {
439  if (std::isnan(_aVec(idx))) {
441  str(boost::format("Unable to determine kernel solution %d (nan)") % idx));
442  }
443  kValues[idx] = _aVec(idx);
444  }
445  _kernel->setKernelParameters(kValues);
446 
448  new ImageT(_kernel->getDimensions())
449  );
450  _kSum = _kernel->computeImage(*image, false);
451 
452  if (_fitForBackground) {
453  if (std::isnan(_aVec(nParameters-1))) {
455  str(boost::format("Unable to determine background solution %d (nan)") %
456  (nParameters-1)));
457  }
458  _background = _aVec(nParameters-1);
459  }
460  }
461 
462 
463  template <typename InputT>
465  throw LSST_EXCEPT(pexExcept::Exception, "Uncertainty calculation not supported");
466 
467  /* Estimate of parameter uncertainties comes from the inverse of the
468  * covariance matrix (noise spectrum).
469  * N.R. 15.4.8 to 15.4.15
470  *
471  * Since this is a linear problem no need to use Fisher matrix
472  * N.R. 15.5.8
473  *
474  * Although I might be able to take advantage of the solution above.
475  * Since this now works and is not the rate limiting step, keep as-is for DC3a.
476  *
477  * Use Cholesky decomposition again.
478  * Cholkesy:
479  * Cov = L L^t
480  * Cov^(-1) = (L L^t)^(-1)
481  * = (L^T)^-1 L^(-1)
482  *
483  * Code would be:
484  *
485  * Eigen::MatrixXd cov = _mMat.transpose() * _mMat;
486  * Eigen::LLT<Eigen::MatrixXd> llt = cov.llt();
487  * Eigen::MatrixXd error2 = llt.matrixL().transpose().inverse()*llt.matrixL().inverse();
488  */
489  }
490 
491  /*******************************************************************************************************/
492 
493  template <typename InputT>
495  lsst::afw::math::KernelList const& basisList,
496  bool fitForBackground
497  ) :
498  StaticKernelSolution<InputT>(basisList, fitForBackground)
499  {};
500 
501  template <typename InputT>
503  lsst::afw::image::Image<InputT> const &templateImage,
504  lsst::afw::image::Image<InputT> const &scienceImage,
507  ) {
508 
509  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
510  if (varStats.getValue(afwMath::MIN) < 0.0) {
512  "Error: variance less than 0.0");
513  }
514  if (varStats.getValue(afwMath::MIN) == 0.0) {
516  "Error: variance equals 0.0, cannot inverse variance weight");
517  }
518 
519  /* Full footprint of all input images */
521  new afwDet::Footprint(std::make_shared<afwGeom::SpanSet>(templateImage.getBBox())));
522 
523  afwMath::KernelList basisList =
524  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
525  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
526 
527  /* Only BAD pixels marked in this mask */
528  afwImage::MaskPixel bitMask =
533 
534  /* Create a Footprint that contains all the masked pixels set above */
536  afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);
537 
538  /* And spread it by the kernel half width */
539  int growPix = (*kiter)->getCtr().getX();
540  afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
541 
542 #if 0
543  for (typename afwDet::FootprintSet::FootprintList::iterator
544  ptr = maskedFpSetGrown.getFootprints()->begin(),
545  end = maskedFpSetGrown.getFootprints()->end();
546  ptr != end;
547  ++ptr) {
548 
549  afwDet::setMaskFromFootprint(finalMask,
550  (**ptr),
552  }
553 #endif
554 
555  afwImage::Mask<afwImage::MaskPixel> finalMask(pixelMask.getDimensions());
556  for (auto const & foot : *(maskedFpSetGrown.getFootprints())) {
557  foot->getSpans()->setMask(finalMask, afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
558  }
559  pixelMask.writeFits("pixelmask.fits");
560  finalMask.writeFits("finalmask.fits");
561 
562 
563  ndarray::Array<int, 1, 1> maskArray =
564  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
565  fullFp->getSpans()->flatten(maskArray, finalMask.getArray(), templateImage.getXY0());
566  auto maskEigen = ndarray::asEigenMatrix(maskArray);
567 
568  ndarray::Array<InputT, 1, 1> arrayTemplate =
569  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
570  fullFp->getSpans()->flatten(arrayTemplate, templateImage.getArray(), templateImage.getXY0());
571  auto eigenTemplate0 = ndarray::asEigenMatrix(arrayTemplate);
572 
573  ndarray::Array<InputT, 1, 1> arrayScience =
574  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
575  fullFp->getSpans()->flatten(arrayScience, scienceImage.getArray(), scienceImage.getXY0());
576  auto eigenScience0 = ndarray::asEigenMatrix(arrayScience);
577 
578  ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance =
579  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
580  fullFp->getSpans()->flatten(arrayVariance, varianceEstimate.getArray(), varianceEstimate.getXY0());
581  auto eigenVariance0 = ndarray::asEigenMatrix(arrayVariance);
582 
583  int nGood = 0;
584  for (int i = 0; i < maskEigen.size(); i++) {
585  if (maskEigen(i) == 0.0)
586  nGood += 1;
587  }
588 
589  Eigen::VectorXd eigenTemplate(nGood);
590  Eigen::VectorXd eigenScience(nGood);
591  Eigen::VectorXd eigenVariance(nGood);
592  int nUsed = 0;
593  for (int i = 0; i < maskEigen.size(); i++) {
594  if (maskEigen(i) == 0.0) {
595  eigenTemplate(nUsed) = eigenTemplate0(i);
596  eigenScience(nUsed) = eigenScience0(i);
597  eigenVariance(nUsed) = eigenVariance0(i);
598  nUsed += 1;
599  }
600  }
601 
602 
603  boost::timer t;
604  t.restart();
605 
606  unsigned int const nKernelParameters = basisList.size();
607  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
608  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
609 
610  /* Holds image convolved with basis function */
611  afwImage::Image<InputT> cimage(templateImage.getDimensions());
612 
613  /* Holds eigen representation of image convolved with all basis functions */
614  std::vector<Eigen::VectorXd> convolvedEigenList(nKernelParameters);
615 
616  /* Iterators over convolved image list and basis list */
617  typename std::vector<Eigen::VectorXd>::iterator eiter = convolvedEigenList.begin();
618 
619  /* Create C_i in the formalism of Alard & Lupton */
620  afwMath::ConvolutionControl convolutionControl;
621  convolutionControl.setDoNormalize(false);
622  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
623  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
624 
625  ndarray::Array<InputT, 1, 1> arrayC =
626  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
627  fullFp->getSpans()->flatten(arrayC, cimage.getArray(), cimage.getXY0());
628  auto eigenC0 = ndarray::asEigenMatrix(arrayC);
629 
630  Eigen::VectorXd eigenC(nGood);
631  int nUsed = 0;
632  for (int i = 0; i < maskEigen.size(); i++) {
633  if (maskEigen(i) == 0.0) {
634  eigenC(nUsed) = eigenC0(i);
635  nUsed += 1;
636  }
637  }
638 
639  *eiter = eigenC;
640  }
641  double time = t.elapsed();
642  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.buildWithMask",
643  "Total compute time to do basis convolutions : %.2f s", time);
644  t.restart();
645 
646  /* Load matrix with all convolved images */
647  Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
648  typename std::vector<Eigen::VectorXd>::iterator eiterj = convolvedEigenList.begin();
649  typename std::vector<Eigen::VectorXd>::iterator eiterE = convolvedEigenList.end();
650  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
651  cMat.block(0, kidxj, eigenTemplate.size(), 1) =
652  Eigen::MatrixXd(*eiterj).block(0, 0, eigenTemplate.size(), 1);
653  }
654  /* Treat the last "image" as all 1's to do the background calculation. */
655  if (this->_fitForBackground)
656  cMat.col(nParameters-1).fill(1.);
657 
658  this->_cMat = cMat;
659  this->_ivVec = eigenVariance.array().inverse().matrix();
660  this->_iVec = eigenScience;
661 
662  /* Make these outside of solve() so I can check condition number */
663  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_cMat);
664  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_iVec);
665  }
666 
667 
668  template <typename InputT>
670  lsst::afw::image::Image<InputT> const &templateImage,
671  lsst::afw::image::Image<InputT> const &scienceImage,
674  ) {
675 
676  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
677  if (varStats.getValue(afwMath::MIN) < 0.0) {
679  "Error: variance less than 0.0");
680  }
681  if (varStats.getValue(afwMath::MIN) == 0.0) {
683  "Error: variance equals 0.0, cannot inverse variance weight");
684  }
685 
686  lsst::afw::math::KernelList basisList =
687  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
688  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
689 
690  /* Only BAD pixels marked in this mask */
692  afwImage::MaskPixel bitMask =
696  sMask &= bitMask;
697  /* TBD: Need to figure out a way to spread this; currently done in Python */
698 
699  unsigned int const nKernelParameters = basisList.size();
700  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
701  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
702 
703  /* NOTE - we are accessing particular elements of Eigen arrays using
704  these coordinates, therefore they need to be in LOCAL coordinates.
705  This was written before ndarray unification.
706  */
707  /* Ignore known EDGE pixels for speed */
708  geom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
709  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
710  "Limits of good pixels after convolution: %d,%d -> %d,%d (local)",
711  shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(),
712  shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());
713 
714  /* Subimages are addressed in raw pixel coordinates */
715  unsigned int startCol = shrunkLocalBBox.getMinX();
716  unsigned int startRow = shrunkLocalBBox.getMinY();
717  unsigned int endCol = shrunkLocalBBox.getMaxX();
718  unsigned int endRow = shrunkLocalBBox.getMaxY();
719  /* Needed for Eigen block slicing operation */
720  endCol += 1;
721  endRow += 1;
722 
723  boost::timer t;
724  t.restart();
725 
726  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
727  Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow,
728  startCol,
729  endRow-startRow,
730  endCol-startCol);
731 
732  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
733  startCol,
734  endRow-startRow,
735  endCol-startCol);
736  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
737  startCol,
738  endRow-startRow,
739  endCol-startCol);
740  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
741  startRow, startCol, endRow-startRow, endCol-startCol
742  ).array().inverse().matrix();
743 
744  /* Resize into 1-D for later usage */
745  eMask.resize(eMask.rows()*eMask.cols(), 1);
746  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
747  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
748  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
749 
750  /* Do the masking, slowly for now... */
751  Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
752  Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
753  Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
754  int nGood = 0;
755  for (int i = 0; i < eMask.rows(); i++) {
756  if (eMask(i, 0) == 0) {
757  maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
758  maskedEigenScience(nGood, 0) = eigenScience(i, 0);
759  maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
760  nGood += 1;
761  }
762  }
763  /* Can't resize() since all values are lost; use blocks */
764  eigenTemplate = maskedEigenTemplate.block(0, 0, nGood, 1);
765  eigenScience = maskedEigenScience.block(0, 0, nGood, 1);
766  eigeniVariance = maskedEigeniVariance.block(0, 0, nGood, 1);
767 
768 
769  /* Holds image convolved with basis function */
770  afwImage::Image<InputT> cimage(templateImage.getDimensions());
771 
772  /* Holds eigen representation of image convolved with all basis functions */
773  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
774 
775  /* Iterators over convolved image list and basis list */
776  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
777  /* Create C_i in the formalism of Alard & Lupton */
778  afwMath::ConvolutionControl convolutionControl;
779  convolutionControl.setDoNormalize(false);
780  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
781  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
782 
783  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
784  startCol,
785  endRow-startRow,
786  endCol-startCol);
787  cMat.resize(cMat.size(), 1);
788 
789  /* Do masking */
790  Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
791  int nGood = 0;
792  for (int i = 0; i < eMask.rows(); i++) {
793  if (eMask(i, 0) == 0) {
794  maskedcMat(nGood, 0) = cMat(i, 0);
795  nGood += 1;
796  }
797  }
798  cMat = maskedcMat.block(0, 0, nGood, 1);
799  *eiter = cMat;
800  }
801 
802  double time = t.elapsed();
803  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
804  "Total compute time to do basis convolutions : %.2f s", time);
805  t.restart();
806 
807  /*
808  Load matrix with all values from convolvedEigenList : all images
809  (eigeniVariance, convolvedEigenList) must be the same size
810  */
811  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
812  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
813  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
814  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
815  cMat.col(kidxj) = eiterj->col(0);
816  }
817  /* Treat the last "image" as all 1's to do the background calculation. */
818  if (this->_fitForBackground)
819  cMat.col(nParameters-1).fill(1.);
820 
821  this->_cMat = cMat;
822  this->_ivVec = eigeniVariance.col(0);
823  this->_iVec = eigenScience.col(0);
824 
825  /* Make these outside of solve() so I can check condition number */
826  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
827  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
828 
829  }
830 
831  /* NOTE - this was written before the ndarray unification. I am rewriting
832  buildSingleMask to use the ndarray stuff */
833  template <typename InputT>
835  lsst::afw::image::Image<InputT> const &templateImage,
836  lsst::afw::image::Image<InputT> const &scienceImage,
838  lsst::geom::Box2I maskBox
839  ) {
840 
841  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
842  if (varStats.getValue(afwMath::MIN) < 0.0) {
844  "Error: variance less than 0.0");
845  }
846  if (varStats.getValue(afwMath::MIN) == 0.0) {
848  "Error: variance equals 0.0, cannot inverse variance weight");
849  }
850 
851  lsst::afw::math::KernelList basisList =
852  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
853 
854  unsigned int const nKernelParameters = basisList.size();
855  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
856  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
857 
858  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
859 
860  /*
861  NOTE : If we are using these views in Afw's Image space, we need to
862  make sure and compensate for the XY0 of the image:
863 
864  geom::Box2I fullBBox = templateImage.getBBox();
865  int maskStartCol = maskBox.getMinX();
866  int maskStartRow = maskBox.getMinY();
867  int maskEndCol = maskBox.getMaxX();
868  int maskEndRow = maskBox.getMaxY();
869 
870 
871  If we are going to be doing the slicing in Eigen matrices derived from
872  the images, ignore the XY0.
873 
874  geom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);
875 
876  int maskStartCol = maskBox.getMinX() - templateImage.getX0();
877  int maskStartRow = maskBox.getMinY() - templateImage.getY0();
878  int maskEndCol = maskBox.getMaxX() - templateImage.getX0();
879  int maskEndRow = maskBox.getMaxY() - templateImage.getY0();
880 
881  */
882 
883 
884  geom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());
885 
886  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
887  "Limits of good pixels after convolution: %d,%d -> %d,%d",
888  shrunkBBox.getMinX(), shrunkBBox.getMinY(),
889  shrunkBBox.getMaxX(), shrunkBBox.getMaxY());
890 
891  unsigned int const startCol = shrunkBBox.getMinX();
892  unsigned int const startRow = shrunkBBox.getMinY();
893  unsigned int const endCol = shrunkBBox.getMaxX();
894  unsigned int const endRow = shrunkBBox.getMaxY();
895 
896  /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
897  slicing using Afw, not Eigen
898 
899  Eigen arrays have index 0,0 in the upper right, while LSST images
900  have 0,0 in the lower left. The y-axis is flipped. When translating
901  Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
902  accounted for. However, we still need to be aware of this fact if
903  addressing subregions of an Eigen matrix. This is why the slicing is
904  done in Afw, its cleaner.
905 
906  Please see examples/maskedKernel.cc for elaboration on some of the
907  tests done to make sure this indexing gets done right.
908 
909  */
910 
911 
912  /* Inner limits; of mask box */
913  int maskStartCol = maskBox.getMinX();
914  int maskStartRow = maskBox.getMinY();
915  int maskEndCol = maskBox.getMaxX();
916  int maskEndRow = maskBox.getMaxY();
917 
918  /*
919 
920  |---------------------------|
921  | Kernel Boundary |
922  | |---------------------| |
923  | | Top | |
924  | |......_________......| |
925  | | | | | |
926  | | L | Box | R | |
927  | | | | | |
928  | |......---------......| |
929  | | Bottom | |
930  | |---------------------| |
931  | |
932  |---------------------------|
933 
934  4 regions we want to extract from the pixels: top bottom left right
935 
936  */
937  geom::Box2I tBox = geom::Box2I(geom::Point2I(startCol, maskEndRow + 1),
938  geom::Point2I(endCol, endRow));
939 
940  geom::Box2I bBox = geom::Box2I(geom::Point2I(startCol, startRow),
941  geom::Point2I(endCol, maskStartRow - 1));
942 
943  geom::Box2I lBox = geom::Box2I(geom::Point2I(startCol, maskStartRow),
944  geom::Point2I(maskStartCol - 1, maskEndRow));
945 
946  geom::Box2I rBox = geom::Box2I(geom::Point2I(maskEndCol + 1, maskStartRow),
947  geom::Point2I(endCol, maskEndRow));
948 
949  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
950  "Upper good pixel region: %d,%d -> %d,%d",
951  tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
952  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
953  "Bottom good pixel region: %d,%d -> %d,%d",
954  bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
955  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
956  "Left good pixel region: %d,%d -> %d,%d",
957  lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
958  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
959  "Right good pixel region: %d,%d -> %d,%d",
960  rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());
961 
962  std::vector<geom::Box2I> boxArray;
963  boxArray.push_back(tBox);
964  boxArray.push_back(bBox);
965  boxArray.push_back(lBox);
966  boxArray.push_back(rBox);
967 
968  int totalSize = tBox.getWidth() * tBox.getHeight();
969  totalSize += bBox.getWidth() * bBox.getHeight();
970  totalSize += lBox.getWidth() * lBox.getHeight();
971  totalSize += rBox.getWidth() * rBox.getHeight();
972 
973  Eigen::MatrixXd eigenTemplate(totalSize, 1);
974  Eigen::MatrixXd eigenScience(totalSize, 1);
975  Eigen::MatrixXd eigeniVariance(totalSize, 1);
976  eigenTemplate.setZero();
977  eigenScience.setZero();
978  eigeniVariance.setZero();
979 
980  boost::timer t;
981  t.restart();
982 
983  int nTerms = 0;
984  typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
985  for (; biter != boxArray.end(); ++biter) {
986  int area = (*biter).getWidth() * (*biter).getHeight();
987 
988  afwImage::Image<InputT> siTemplate(templateImage, *biter);
989  afwImage::Image<InputT> siScience(scienceImage, *biter);
990  afwImage::Image<InputT> sVarEstimate(varianceEstimate, *biter);
991 
992  Eigen::MatrixXd eTemplate = imageToEigenMatrix(siTemplate);
993  Eigen::MatrixXd eScience = imageToEigenMatrix(siScience);
994  Eigen::MatrixXd eiVarEstimate = imageToEigenMatrix(sVarEstimate).array().inverse().matrix();
995 
996  eTemplate.resize(area, 1);
997  eScience.resize(area, 1);
998  eiVarEstimate.resize(area, 1);
999 
1000  eigenTemplate.block(nTerms, 0, area, 1) = eTemplate.block(0, 0, area, 1);
1001  eigenScience.block(nTerms, 0, area, 1) = eScience.block(0, 0, area, 1);
1002  eigeniVariance.block(nTerms, 0, area, 1) = eiVarEstimate.block(0, 0, area, 1);
1003 
1004  nTerms += area;
1005  }
1006 
1007  afwImage::Image<InputT> cimage(templateImage.getDimensions());
1008 
1009  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
1010  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
1011  /* Create C_i in the formalism of Alard & Lupton */
1012  afwMath::ConvolutionControl convolutionControl;
1013  convolutionControl.setDoNormalize(false);
1014  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
1015  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
1016  Eigen::MatrixXd cMat(totalSize, 1);
1017  cMat.setZero();
1018 
1019  int nTerms = 0;
1020  typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
1021  for (; biter != boxArray.end(); ++biter) {
1022  int area = (*biter).getWidth() * (*biter).getHeight();
1023 
1024  afwImage::Image<InputT> csubimage(cimage, *biter);
1025  Eigen::MatrixXd esubimage = imageToEigenMatrix(csubimage);
1026  esubimage.resize(area, 1);
1027  cMat.block(nTerms, 0, area, 1) = esubimage.block(0, 0, area, 1);
1028 
1029  nTerms += area;
1030  }
1031 
1032  *eiter = cMat;
1033 
1034  }
1035 
1036  double time = t.elapsed();
1037  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
1038  "Total compute time to do basis convolutions : %.2f s", time);
1039  t.restart();
1040 
1041  /*
1042  Load matrix with all values from convolvedEigenList : all images
1043  (eigeniVariance, convolvedEigenList) must be the same size
1044  */
1045  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
1046  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
1047  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
1048  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
1049  cMat.col(kidxj) = eiterj->col(0);
1050  }
1051  /* Treat the last "image" as all 1's to do the background calculation. */
1052  if (this->_fitForBackground)
1053  cMat.col(nParameters-1).fill(1.);
1054 
1055  this->_cMat = cMat;
1056  this->_ivVec = eigeniVariance.col(0);
1057  this->_iVec = eigenScience.col(0);
1058 
1059  /* Make these outside of solve() so I can check condition number */
1060  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1061  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1062  }
1063  /*******************************************************************************************************/
1064 
1065 
1066  template <typename InputT>
1068  lsst::afw::math::KernelList const& basisList,
1069  bool fitForBackground,
1070  Eigen::MatrixXd const& hMat,
1072  )
1073  :
1074  StaticKernelSolution<InputT>(basisList, fitForBackground),
1075  _hMat(hMat),
1076  _ps(ps.deepCopy())
1077  {};
1078 
1079  template <typename InputT>
1081  Eigen::MatrixXd vMat = this->_cMat.jacobiSvd().matrixV();
1082  Eigen::MatrixXd vMatvMatT = vMat * vMat.transpose();
1083 
1084  /* Find pseudo inverse of mMat, which may be ill conditioned */
1085  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(this->_mMat);
1086  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
1087  Eigen::VectorXd eValues = eVecValues.eigenvalues();
1088  double eMax = eValues.maxCoeff();
1089  for (int i = 0; i != eValues.rows(); ++i) {
1090  if (eValues(i) == 0.0) {
1091  eValues(i) = 0.0;
1092  }
1093  else if ((eMax / eValues(i)) > maxCond) {
1094  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1095  "Truncating eValue %d; %.5e / %.5e = %.5e vs. %.5e",
1096  i, eMax, eValues(i), eMax / eValues(i), maxCond);
1097  eValues(i) = 0.0;
1098  }
1099  else {
1100  eValues(i) = 1.0 / eValues(i);
1101  }
1102  }
1103  Eigen::MatrixXd mInv = rMat * eValues.asDiagonal() * rMat.transpose();
1104 
1105  std::vector<double> lambdas = _createLambdaSteps();
1106  std::vector<double> risks;
1107  for (unsigned int i = 0; i < lambdas.size(); i++) {
1108  double l = lambdas[i];
1109  Eigen::MatrixXd mLambda = this->_mMat + l * _hMat;
1110 
1111  try {
1112  KernelSolution::solve(mLambda, this->_bVec);
1113  } catch (pexExcept::Exception &e) {
1114  LSST_EXCEPT_ADD(e, "Unable to solve regularized kernel matrix");
1115  throw e;
1116  }
1117  Eigen::VectorXd term1 = (this->_aVec.transpose() * vMatvMatT * this->_aVec);
1118  if (term1.size() != 1)
1119  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1120 
1121  double term2a = (vMatvMatT * mLambda.inverse()).trace();
1122 
1123  Eigen::VectorXd term2b = (this->_aVec.transpose() * (mInv * this->_bVec));
1124  if (term2b.size() != 1)
1125  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1126 
1127  double risk = term1(0) + 2 * (term2a - term2b(0));
1128  LOGL_DEBUG("TRACE4.ip.diffim.RegularizedKernelSolution.estimateRisk",
1129  "Lambda = %.3f, Risk = %.5e",
1130  l, risk);
1131  LOGL_DEBUG("TRACE5.ip.diffim.RegularizedKernelSolution.estimateRisk",
1132  "%.5e + 2 * (%.5e - %.5e)",
1133  term1(0), term2a, term2b(0));
1134  risks.push_back(risk);
1135  }
1136  std::vector<double>::iterator it = min_element(risks.begin(), risks.end());
1137  int index = distance(risks.begin(), it);
1138  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1139  "Minimum Risk = %.3e at lambda = %.3e", risks[index], lambdas[index]);
1140 
1141  return lambdas[index];
1142 
1143  }
1144 
1145  template <typename InputT>
1146  Eigen::MatrixXd RegularizedKernelSolution<InputT>::getM(bool includeHmat) {
1147  if (includeHmat == true) {
1148  return this->_mMat + _lambda * _hMat;
1149  }
1150  else {
1151  return this->_mMat;
1152  }
1153  }
1154 
1155  template <typename InputT>
1157 
1158  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1159  "cMat is %d x %d; vVec is %d; iVec is %d; hMat is %d x %d",
1160  this->_cMat.rows(), this->_cMat.cols(), this->_ivVec.size(),
1161  this->_iVec.size(), _hMat.rows(), _hMat.cols());
1162 
1163  if (DEBUG_MATRIX2) {
1164  std::cout << "ID: " << (this->_id) << std::endl;
1165  std::cout << "C:" << std::endl;
1166  std::cout << this->_cMat << std::endl;
1167  std::cout << "Sigma^{-1}:" << std::endl;
1168  std::cout << Eigen::MatrixXd(this->_ivVec.asDiagonal()) << std::endl;
1169  std::cout << "Y:" << std::endl;
1170  std::cout << this->_iVec << std::endl;
1171  std::cout << "H:" << std::endl;
1172  std::cout << _hMat << std::endl;
1173  }
1174 
1175 
1176  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1177  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1178 
1179 
1180  /* See N.R. 18.5
1181 
1182  Matrix equation to solve is Y = C a + N
1183  Y = vectorized version of I (I = image to not convolve)
1184  C_i = K_i (x) R (R = image to convolve)
1185  a = kernel coefficients
1186  N = noise
1187 
1188  If we reweight everything by the inverse square root of the noise
1189  covariance, we get a linear model with the identity matrix for
1190  the noise. The problem can then be solved using least squares,
1191  with the normal equations
1192 
1193  C^T Y = C^T C a
1194 
1195  or
1196 
1197  b = M a
1198 
1199  with
1200 
1201  b = C^T Y
1202  M = C^T C
1203  a = (C^T C)^{-1} C^T Y
1204 
1205 
1206  We can regularize the least square problem
1207 
1208  Y = C a + lambda a^T H a (+ N, which can be ignored)
1209 
1210  or the normal equations
1211 
1212  (C^T C + lambda H) a = C^T Y
1213 
1214 
1215  The solution to the regularization of the least squares problem is
1216 
1217  a = (C^T C + lambda H)^{-1} C^T Y
1218 
1219  The approximation to Y is
1220 
1221  C a = C (C^T C + lambda H)^{-1} C^T Y
1222 
1223  with smoothing matrix
1224 
1225  S = C (C^T C + lambda H)^{-1} C^T
1226 
1227  */
1228 
1229  std::string lambdaType = _ps->getAsString("lambdaType");
1230  if (lambdaType == "absolute") {
1231  _lambda = _ps->getAsDouble("lambdaValue");
1232  }
1233  else if (lambdaType == "relative") {
1234  _lambda = this->_mMat.trace() / this->_hMat.trace();
1235  _lambda *= _ps->getAsDouble("lambdaScaling");
1236  }
1237  else if (lambdaType == "minimizeBiasedRisk") {
1238  double tol = _ps->getAsDouble("maxConditionNumber");
1239  _lambda = estimateRisk(tol);
1240  }
1241  else if (lambdaType == "minimizeUnbiasedRisk") {
1242  _lambda = estimateRisk(std::numeric_limits<double>::max());
1243  }
1244  else {
1245  throw LSST_EXCEPT(pexExcept::Exception, "lambdaType in PropertySet not recognized");
1246  }
1247 
1248  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1249  "Applying kernel regularization with lambda = %.2e", _lambda);
1250 
1251 
1252  try {
1253  KernelSolution::solve(this->_mMat + _lambda * _hMat, this->_bVec);
1254  } catch (pexExcept::Exception &e) {
1255  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
1256  throw e;
1257  }
1258  /* Turn matrices into _kernel and _background */
1260  }
1261 
1262  template <typename InputT>
1264  std::vector<double> lambdas;
1265 
1266  std::string lambdaStepType = _ps->getAsString("lambdaStepType");
1267  if (lambdaStepType == "linear") {
1268  double lambdaLinMin = _ps->getAsDouble("lambdaLinMin");
1269  double lambdaLinMax = _ps->getAsDouble("lambdaLinMax");
1270  double lambdaLinStep = _ps->getAsDouble("lambdaLinStep");
1271  for (double l = lambdaLinMin; l <= lambdaLinMax; l += lambdaLinStep) {
1272  lambdas.push_back(l);
1273  }
1274  }
1275  else if (lambdaStepType == "log") {
1276  double lambdaLogMin = _ps->getAsDouble("lambdaLogMin");
1277  double lambdaLogMax = _ps->getAsDouble("lambdaLogMax");
1278  double lambdaLogStep = _ps->getAsDouble("lambdaLogStep");
1279  for (double l = lambdaLogMin; l <= lambdaLogMax; l += lambdaLogStep) {
1280  lambdas.push_back(pow(10, l));
1281  }
1282  }
1283  else {
1284  throw LSST_EXCEPT(pexExcept::Exception, "lambdaStepType in PropertySet not recognized");
1285  }
1286  return lambdas;
1287  }
1288 
1289  /*******************************************************************************************************/
1290 
1292  lsst::afw::math::KernelList const& basisList,
1293  lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction,
1296  ) :
1297  KernelSolution(),
1298  _spatialKernelFunction(spatialKernelFunction),
1299  _constantFirstTerm(false),
1300  _kernel(),
1301  _background(background),
1302  _kSum(0.0),
1303  _ps(ps.deepCopy()),
1304  _nbases(0),
1305  _nkt(0),
1306  _nbt(0),
1307  _nt(0) {
1308 
1309  bool isAlardLupton = _ps->getAsString("kernelBasisSet") == "alard-lupton";
1310  bool usePca = _ps->getAsBool("usePcaForSpatialKernel");
1311  if (isAlardLupton || usePca) {
1312  _constantFirstTerm = true;
1313  }
1314  this->_fitForBackground = _ps->getAsBool("fitForBackground");
1315 
1316  _nbases = basisList.size();
1317  _nkt = _spatialKernelFunction->getParameters().size();
1318  _nbt = _fitForBackground ? _background->getParameters().size() : 0;
1319  _nt = 0;
1320  if (_constantFirstTerm) {
1321  _nt = (_nbases - 1) * _nkt + 1 + _nbt;
1322  } else {
1323  _nt = _nbases * _nkt + _nbt;
1324  }
1325 
1326  Eigen::MatrixXd mMat(_nt, _nt);
1327  Eigen::VectorXd bVec(_nt);
1328  mMat.setZero();
1329  bVec.setZero();
1330 
1331  this->_mMat = mMat;
1332  this->_bVec = bVec;
1333 
1335  new afwMath::LinearCombinationKernel(basisList, *_spatialKernelFunction)
1336  );
1337 
1338  LOGL_DEBUG("TRACE3.ip.diffim.SpatialKernelSolution",
1339  "Initializing with size %d %d %d and constant first term = %s",
1340  _nkt, _nbt, _nt,
1341  _constantFirstTerm ? "true" : "false");
1342 
1343  }
1344 
1345  void SpatialKernelSolution::addConstraint(float xCenter, float yCenter,
1346  Eigen::MatrixXd const& qMat,
1347  Eigen::VectorXd const& wVec) {
1348 
1349  LOGL_DEBUG("TRACE5.ip.diffim.SpatialKernelSolution.addConstraint",
1350  "Adding candidate at %f, %f", xCenter, yCenter);
1351 
1352  /* Calculate P matrices */
1353  /* Pure kernel terms */
1354  Eigen::VectorXd pK(_nkt);
1355  std::vector<double> paramsK = _spatialKernelFunction->getParameters();
1356  for (int idx = 0; idx < _nkt; idx++) { paramsK[idx] = 0.0; }
1357  for (int idx = 0; idx < _nkt; idx++) {
1358  paramsK[idx] = 1.0;
1359  _spatialKernelFunction->setParameters(paramsK);
1360  pK(idx) = (*_spatialKernelFunction)(xCenter, yCenter); /* Assume things don't vary over stamp */
1361  paramsK[idx] = 0.0;
1362  }
1363  Eigen::MatrixXd pKpKt = (pK * pK.transpose());
1364 
1365  Eigen::VectorXd pB;
1366  Eigen::MatrixXd pBpBt;
1367  Eigen::MatrixXd pKpBt;
1368  if (_fitForBackground) {
1369  pB = Eigen::VectorXd(_nbt);
1370 
1371  /* Pure background terms */
1372  std::vector<double> paramsB = _background->getParameters();
1373  for (int idx = 0; idx < _nbt; idx++) { paramsB[idx] = 0.0; }
1374  for (int idx = 0; idx < _nbt; idx++) {
1375  paramsB[idx] = 1.0;
1376  _background->setParameters(paramsB);
1377  pB(idx) = (*_background)(xCenter, yCenter); /* Assume things don't vary over stamp */
1378  paramsB[idx] = 0.0;
1379  }
1380  pBpBt = (pB * pB.transpose());
1381 
1382  /* Cross terms */
1383  pKpBt = (pK * pB.transpose());
1384  }
1385 
1386  if (DEBUG_MATRIX) {
1387  std::cout << "Spatial weights" << std::endl;
1388  std::cout << "pKpKt " << pKpKt << std::endl;
1389  if (_fitForBackground) {
1390  std::cout << "pBpBt " << pBpBt << std::endl;
1391  std::cout << "pKpBt " << pKpBt << std::endl;
1392  }
1393  }
1394 
1395  if (DEBUG_MATRIX) {
1396  std::cout << "Spatial matrix inputs" << std::endl;
1397  std::cout << "M " << qMat << std::endl;
1398  std::cout << "B " << wVec << std::endl;
1399  }
1400 
1401  /* first index to start the spatial blocks; default=0 for non-constant first term */
1402  int m0 = 0;
1403  /* how many rows/cols to adjust the matrices/vectors; default=0 for non-constant first term */
1404  int dm = 0;
1405  /* where to start the background terms; this is always true */
1406  int mb = _nt - _nbt;
1407 
1408  if (_constantFirstTerm) {
1409  m0 = 1; /* we need to manually fill in the first (non-spatial) terms below */
1410  dm = _nkt-1; /* need to shift terms due to lack of spatial variation in first term */
1411 
1412  _mMat(0, 0) += qMat(0,0);
1413  for(int m2 = 1; m2 < _nbases; m2++) {
1414  _mMat.block(0, m2*_nkt-dm, 1, _nkt) += qMat(0,m2) * pK.transpose();
1415  }
1416  _bVec(0) += wVec(0);
1417 
1418  if (_fitForBackground) {
1419  _mMat.block(0, mb, 1, _nbt) += qMat(0,_nbases) * pB.transpose();
1420  }
1421  }
1422 
1423  /* Fill in the spatial blocks */
1424  for(int m1 = m0; m1 < _nbases; m1++) {
1425  /* Diagonal kernel-kernel term; only use upper triangular part of pKpKt */
1426  _mMat.block(m1*_nkt-dm, m1*_nkt-dm, _nkt, _nkt) +=
1427  (pKpKt * qMat(m1,m1)).triangularView<Eigen::Upper>();
1428 
1429  /* Kernel-kernel terms */
1430  for(int m2 = m1+1; m2 < _nbases; m2++) {
1431  _mMat.block(m1*_nkt-dm, m2*_nkt-dm, _nkt, _nkt) += qMat(m1,m2) * pKpKt;
1432  }
1433 
1434  if (_fitForBackground) {
1435  /* Kernel cross terms with background */
1436  _mMat.block(m1*_nkt-dm, mb, _nkt, _nbt) += qMat(m1,_nbases) * pKpBt;
1437  }
1438 
1439  /* B vector */
1440  _bVec.segment(m1*_nkt-dm, _nkt) += wVec(m1) * pK;
1441  }
1442 
1443  if (_fitForBackground) {
1444  /* Background-background terms only */
1445  _mMat.block(mb, mb, _nbt, _nbt) +=
1446  (pBpBt * qMat(_nbases,_nbases)).triangularView<Eigen::Upper>();
1447  _bVec.segment(mb, _nbt) += wVec(_nbases) * pB;
1448  }
1449 
1450  if (DEBUG_MATRIX) {
1451  std::cout << "Spatial matrix outputs" << std::endl;
1452  std::cout << "mMat " << _mMat << std::endl;
1453  std::cout << "bVec " << _bVec << std::endl;
1454  }
1455 
1456  }
1457 
1460  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
1461  }
1464  );
1465  (void)_kernel->computeImage(*image, false, pos[0], pos[1]);
1466  return image;
1467  }
1468 
1470  /* Fill in the other half of mMat */
1471  for (int i = 0; i < _nt; i++) {
1472  for (int j = i+1; j < _nt; j++) {
1473  _mMat(j,i) = _mMat(i,j);
1474  }
1475  }
1476 
1477  try {
1479  } catch (pexExcept::Exception &e) {
1480  LSST_EXCEPT_ADD(e, "Unable to solve spatial kernel matrix");
1481  throw e;
1482  }
1483  /* Turn matrices into _kernel and _background */
1484  _setKernel();
1485  }
1486 
1490  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
1491  }
1492 
1493  return std::make_pair(_kernel, _background);
1494  }
1495 
1496  void SpatialKernelSolution::_setKernel() {
1497  /* Report on condition number */
1498  double cNumber = this->getConditionNumber(EIGENVALUE);
1499 
1500  if (_nkt == 1) {
1501  /* Not spatially varying; this fork is a specialization for convolution speed--up */
1502 
1503  /* Set the basis coefficients */
1504  std::vector<double> kCoeffs(_nbases);
1505  for (int i = 0; i < _nbases; i++) {
1506  if (std::isnan(_aVec(i))) {
1507  throw LSST_EXCEPT(
1509  str(boost::format(
1510  "I. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % i % cNumber));
1511  }
1512  kCoeffs[i] = _aVec(i);
1513  }
1514  lsst::afw::math::KernelList basisList =
1515  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
1516  _kernel.reset(
1517  new afwMath::LinearCombinationKernel(basisList, kCoeffs)
1518  );
1519  }
1520  else {
1521 
1522  /* Set the kernel coefficients */
1524  kCoeffs.reserve(_nbases);
1525  for (int i = 0, idx = 0; i < _nbases; i++) {
1526  kCoeffs.push_back(std::vector<double>(_nkt));
1527 
1528  /* Deal with the possibility the first term doesn't vary spatially */
1529  if ((i == 0) && (_constantFirstTerm)) {
1530  if (std::isnan(_aVec(idx))) {
1531  throw LSST_EXCEPT(
1533  str(boost::format(
1534  "II. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % idx % cNumber));
1535  }
1536  kCoeffs[i][0] = _aVec(idx++);
1537  }
1538  else {
1539  for (int j = 0; j < _nkt; j++) {
1540  if (std::isnan(_aVec(idx))) {
1541  throw LSST_EXCEPT(
1543  str(boost::format(
1544  "III. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % idx % cNumber));
1545  }
1546  kCoeffs[i][j] = _aVec(idx++);
1547  }
1548  }
1549  }
1550  _kernel->setSpatialParameters(kCoeffs);
1551  }
1552 
1553  /* Set kernel Sum */
1555  _kSum = _kernel->computeImage(*image, false);
1556 
1557  /* Set the background coefficients */
1558  std::vector<double> bgCoeffs(_fitForBackground ? _nbt : 1);
1559  if (_fitForBackground) {
1560  for (int i = 0; i < _nbt; i++) {
1561  int idx = _nt - _nbt + i;
1562  if (std::isnan(_aVec(idx))) {
1564  str(boost::format(
1565  "Unable to determine spatial background solution %d (nan)") % (idx)));
1566  }
1567  bgCoeffs[i] = _aVec(idx);
1568  }
1569  }
1570  else {
1571  bgCoeffs[0] = 0.;
1572  }
1573  _background->setParameters(bgCoeffs);
1574  }
1575 
1576 /***********************************************************************************************************/
1577 //
1578 // Explicit instantiations
1579 //
1580  typedef float InputT;
1581 
1582  template class StaticKernelSolution<InputT>;
1583  template class MaskedKernelSolution<InputT>;
1584  template class RegularizedKernelSolution<InputT>;
1585 
1586 }}} // end of namespace lsst::ip::diffim
int end
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
Image Subtraction helper functions.
Declaration of classes to store the solution for convolution kernels.
LSST DM logging module built on log4cxx.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:504
uint64_t * ptr
Definition: RangeSet.cc:88
T begin(T... args)
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
std::shared_ptr< FootprintList > getFootprints()
: Return the Footprints of detected objects
Definition: FootprintSet.h:156
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
@ BITMASK
Use (pixels & (given mask))
Definition: Threshold.h:48
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:445
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition: ImageBase.h:356
lsst::geom::Point2I getXY0() const
Return the image's origin.
Definition: ImageBase.h:323
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
void writeFits(std::string const &fileName, std::shared_ptr< lsst::daf::base::PropertySet const > metadata=std::shared_ptr< lsst::daf::base::PropertySet >(), std::string const &mode="w") const
Write a mask to a regular FITS file.
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:379
Parameters to control convolution.
Definition: ConvolveImage.h:50
void setDoNormalize(bool doNormalize)
Definition: ConvolveImage.h:66
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:213
std::shared_ptr< lsst::afw::math::Function2< double > > SpatialFunctionPtr
Definition: Kernel.h:114
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:85
void setSpatialParameters(const std::vector< std::vector< double >> params)
Set the parameters of all spatial functions.
Definition: Kernel.cc:119
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:705
A class to evaluate image statistics.
Definition: Statistics.h:215
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
Class for storing generic metadata.
Definition: PropertySet.h:67
An integer coordinate rectangle.
Definition: Box.h:55
int getMinY() const noexcept
Definition: Box.h:158
int getHeight() const noexcept
Definition: Box.h:188
int getMinX() const noexcept
Definition: Box.h:157
int getWidth() const noexcept
Definition: Box.h:187
int getMaxX() const noexcept
Definition: Box.h:161
int getMaxY() const noexcept
Definition: Box.h:162
bool _fitForBackground
Background terms included in fit.
virtual double getConditionNumber(ConditionNumberType conditionType)
Eigen::VectorXd _bVec
Derived least squares B vector.
Eigen::VectorXd _aVec
Derived least squares solution matrix.
KernelSolvedBy _solvedBy
Type of algorithm used to make solution.
Eigen::MatrixXd _mMat
Derived least squares M matrix.
static int _SolutionId
Unique identifier for solution.
Eigen::MatrixXd const & getM()
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
virtual void buildWithMask(lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &pixelMask)
virtual void buildSingleMaskOrig(lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::geom::Box2I maskBox)
MaskedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground)
virtual void buildOrig(lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > pixelMask)
RegularizedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground, Eigen::MatrixXd const &hMat, lsst::daf::base::PropertySet const &ps)
std::pair< std::shared_ptr< lsst::afw::math::LinearCombinationKernel >, lsst::afw::math::Kernel::SpatialFunctionPtr > getSolutionPair()
std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage(lsst::geom::Point2D const &pos)
SpatialKernelSolution(lsst::afw::math::KernelList const &basisList, lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction, lsst::afw::math::Kernel::SpatialFunctionPtr background, lsst::daf::base::PropertySet const &ps)
void addConstraint(float xCenter, float yCenter, Eigen::MatrixXd const &qMat, Eigen::VectorXd const &wVec)
virtual std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > getSolutionPair()
void _setKernel()
Set kernel after solution.
virtual std::shared_ptr< lsst::afw::math::Kernel > getKernel()
std::shared_ptr< lsst::afw::math::Kernel > _kernel
Derived single-object convolution kernel.
void _setKernelUncertainty()
Not implemented.
virtual std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage()
StaticKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground)
virtual void build(lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate)
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
Reports invalid arguments.
Definition: Runtime.h:66
T end(T... args)
T endl(T... args)
T isnan(T... args)
T make_pair(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
@ MIN
estimate sample minimum
Definition: Statistics.h:76
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns a 2-d Image into a 2-d Eigen Matrix.
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
T reset(T... args)
T size(T... args)
#define DEBUG_MATRIX2
#define DEBUG_MATRIX