LSSTApplications  20.0.0
LSSTDataManagementBasePackage
Namespaces | Classes | Typedefs | Functions | Variables
lsst::ip::diffim Namespace Reference

Namespaces

 dcrModel
 
 detail
 
 diaCatalogSourceSelector
 
 diaSourceAnalysis
 
 diffimLib
 
 diffimTools
 
 dipoleFitTask
 
 dipoleMeasurement
 
 getTemplate
 
 imageDecorrelation
 
 imageMapReduce
 
 imagePsfMatch
 
 kernelCandidateQa
 
 makeKernelBasisList
 
 metrics
 
 modelPsfMatch
 
 psfMatch
 
 snapPsfMatch
 
 utils
 
 version
 
 zogy
 

Classes

class  DipoleCentroidAlgorithm
 Intermediate base class for algorithms that compute a centroid. More...
 
class  DipoleCentroidControl
 
class  DipoleFluxAlgorithm
 Intermediate base class for algorithms that compute a flux. More...
 
class  DipoleFluxControl
 
class  FindSetBits
 Class to accumulate Mask bits. More...
 
class  ImageStatistics
 Class to calculate difference image statistics. More...
 
class  KernelCandidate
 Class stored in SpatialCells for spatial Kernel fitting. More...
 
class  KernelCandidateDetection
 Search through images for Footprints with no masked pixels. More...
 
class  KernelSolution
 
class  MaskedKernelSolution
 
class  MinimizeDipoleChi2
 Class to minimize PsfDipoleFlux; this is the object that Minuit minimizes. More...
 
class  NaiveDipoleCentroid
 Intermediate base class for algorithms that compute a centroid. More...
 
class  NaiveDipoleFlux
 
class  PsfDipoleFlux
 Implementation of Psf dipole flux. More...
 
class  PsfDipoleFluxControl
 C++ control object for PSF dipole fluxes. More...
 
class  RegularizedKernelSolution
 
class  SpatialKernelSolution
 
class  StaticKernelSolution
 

Typedefs

typedef float PixelT
 
typedef float InputT
 

Functions

lsst::afw::math::KernelList makeDeltaFunctionBasisList (int width, int height)
 Build a set of Delta Function basis kernels. More...
 
Eigen::MatrixXd makeRegularizationMatrix (lsst::daf::base::PropertySet const &ps)
 Build a regularization matrix for Delta function kernels. More...
 
Eigen::MatrixXd makeForwardDifferenceMatrix (int width, int height, std::vector< int > const &orders, float borderPenalty, bool fitForBackground)
 Build a forward difference regularization matrix for Delta function kernels. More...
 
Eigen::MatrixXd makeCentralDifferenceMatrix (int width, int height, int stencil, float borderPenalty, bool fitForBackground)
 Build a central difference Laplacian regularization matrix for Delta function kernels. More...
 
lsst::afw::math::KernelList renormalizeKernelList (lsst::afw::math::KernelList const &kernelListIn)
 Renormalize a list of basis kernels. More...
 
lsst::afw::math::KernelList makeAlardLuptonBasisList (int halfWidth, int nGauss, std::vector< double > const &sigGauss, std::vector< int > const &degGauss)
 Build a set of Alard/Lupton basis kernels. More...
 
template<typename PixelT , typename BackgroundT >
lsst::afw::image::MaskedImage< PixelT > convolveAndSubtract (lsst::afw::image::MaskedImage< PixelT > const &templateImage, lsst::afw::image::MaskedImage< PixelT > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, BackgroundT background, bool invert=true)
 Execute fundamental task of convolving template and subtracting it from science image. More...
 
template<typename PixelT , typename BackgroundT >
lsst::afw::image::MaskedImage< PixelT > convolveAndSubtract (lsst::afw::image::Image< PixelT > const &templateImage, lsst::afw::image::MaskedImage< PixelT > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, BackgroundT background, bool invert=true)
 Execute fundamental task of convolving template and subtracting it from science image. More...
 
template<typename PixelT >
Eigen::MatrixXd imageToEigenMatrix (lsst::afw::image::Image< PixelT > const &img)
 Turns a 2-d Image into a 2-d Eigen Matrix. More...
 
Eigen::MatrixXi maskToEigenMatrix (lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
 
template<typename PixelT >
std::shared_ptr< KernelCandidate< PixelT > > makeKernelCandidate (float const xCenter, float const yCenter, std::shared_ptr< afw::image::MaskedImage< PixelT > > const &templateMaskedImage, std::shared_ptr< afw::image::MaskedImage< PixelT > > const &scienceMaskedImage, daf::base::PropertySet const &ps)
 Return a KernelCandidate pointer of the right sort. More...
 
template<typename PixelT >
std::shared_ptr< KernelCandidate< PixelT > > makeKernelCandidate (std::shared_ptr< afw::table::SourceRecord > const &source, std::shared_ptr< afw::image::MaskedImage< PixelT > > const &templateMaskedImage, std::shared_ptr< afw::image::MaskedImage< PixelT > > const &scienceMaskedImage, daf::base::PropertySet const &ps)
 Return a KernelCandidate pointer of the right sort. More...
 
 PYBIND11_MODULE (basisLists, mod)
 
 PYBIND11_MODULE (_dipoleAlgorithms, mod)
 
 PYBIND11_MODULE (findSetBits, mod)
 
 PYBIND11_MODULE (imageStatistics, mod)
 
 PYBIND11_MODULE (imageSubtract, mod)
 
 PYBIND11_MODULE (kernelCandidate, mod)
 
 PYBIND11_MODULE (kernelCandidateDetection, mod)
 
 PYBIND11_MODULE (kernelSolution, mod)
 
Eigen::MatrixXd makeFiniteDifferenceRegularizationDeprecated (unsigned int width, unsigned int height, unsigned int order, unsigned int boundary_style, unsigned int difference_style, bool printB)
 Generate regularization matrix for delta function kernels. More...
 
int const NEGCENTXPAR (0)
 
int const NEGCENTYPAR (1)
 
int const NEGFLUXPAR (2)
 
int const POSCENTXPAR (3)
 
int const POSCENTYPAR (4)
 
int const POSFLUXPAR (5)
 
template<typename PixelT , typename BackgroundT >
afwImage::MaskedImage< PixelT > convolveAndSubtract (lsst::afw::image::MaskedImage< PixelT > const &templateImage, lsst::afw::image::MaskedImage< PixelT > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, BackgroundT background, bool invert)
 Implement fundamental difference imaging step of convolution and subtraction : D = I - (K*T + bg) where * denotes convolution. More...
 
template<typename PixelT , typename BackgroundT >
afwImage::MaskedImage< PixelT > convolveAndSubtract (lsst::afw::image::Image< PixelT > const &templateImage, lsst::afw::image::MaskedImage< PixelT > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, BackgroundT background, bool invert)
 Implement fundamental difference imaging step of convolution and subtraction : D = I - (K.x.T + bg) More...
 
template Eigen::MatrixXd imageToEigenMatrix (lsst::afw::image::Image< double > const &)
 
template lsst::afw::image::MaskedImage< float > convolveAndSubtract (lsst::afw::image::Image< float > const &templateImage, lsst::afw::image::MaskedImage< float > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, double background, bool invert)
 
template afwImage::MaskedImage< float > convolveAndSubtract (lsst::afw::image::Image< float > const &templateImage, lsst::afw::image::MaskedImage< float > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, lsst::afw::math::Function2< double > const &backgroundFunction, bool invert)
 
template lsst::afw::image::MaskedImage< float > convolveAndSubtract (lsst::afw::image::MaskedImage< float > const &templateImage, lsst::afw::image::MaskedImage< float > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, double background, bool invert)
 
template afwImage::MaskedImage< float > convolveAndSubtract (lsst::afw::image::MaskedImage< float > const &templateImage, lsst::afw::image::MaskedImage< float > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, lsst::afw::math::Function2< double > const &backgroundFunction, bool invert)
 
template lsst::afw::image::MaskedImage< double > convolveAndSubtract (lsst::afw::image::Image< double > const &templateImage, lsst::afw::image::MaskedImage< double > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, double background, bool invert)
 
template afwImage::MaskedImage< double > convolveAndSubtract (lsst::afw::image::Image< double > const &templateImage, lsst::afw::image::MaskedImage< double > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, lsst::afw::math::Function2< double > const &backgroundFunction, bool invert)
 
template lsst::afw::image::MaskedImage< double > convolveAndSubtract (lsst::afw::image::MaskedImage< double > const &templateImage, lsst::afw::image::MaskedImage< double > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, double background, bool invert)
 
template afwImage::MaskedImage< double > convolveAndSubtract (lsst::afw::image::MaskedImage< double > const &templateImage, lsst::afw::image::MaskedImage< double > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, lsst::afw::math::Function2< double > const &backgroundFunction, bool invert)
 

Variables

 NaiveDipoleCentroid
 
 Control
 
 DipoleCentroidControl
 
 executionOrder
 
 NaiveDipoleFlux
 
 DipoleFluxControl
 
 PsfDipoleFlux
 
 PsfDipoleFluxControl
 

Typedef Documentation

◆ InputT

typedef float lsst::ip::diffim::InputT

Definition at line 1572 of file KernelSolution.cc.

◆ PixelT

typedef float lsst::ip::diffim::PixelT

Definition at line 383 of file KernelCandidate.cc.

Function Documentation

◆ convolveAndSubtract() [1/12]

template lsst::afw::image::MaskedImage< double > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::Image< double > const &  templateImage,
lsst::afw::image::MaskedImage< double > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
double  background,
bool  invert 
)

◆ convolveAndSubtract() [2/12]

template afwImage::MaskedImage< double > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::Image< double > const &  templateImage,
lsst::afw::image::MaskedImage< double > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
lsst::afw::math::Function2< double > const &  backgroundFunction,
bool  invert 
)

◆ convolveAndSubtract() [3/12]

template lsst::afw::image::MaskedImage< float > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::Image< float > const &  templateImage,
lsst::afw::image::MaskedImage< float > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
double  background,
bool  invert 
)

◆ convolveAndSubtract() [4/12]

template afwImage::MaskedImage< float > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::Image< float > const &  templateImage,
lsst::afw::image::MaskedImage< float > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
lsst::afw::math::Function2< double > const &  backgroundFunction,
bool  invert 
)

◆ convolveAndSubtract() [5/12]

template<typename PixelT , typename BackgroundT >
afwImage::MaskedImage<PixelT> lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::Image< PixelT > const &  templateImage,
lsst::afw::image::MaskedImage< PixelT > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
BackgroundT  background,
bool  invert 
)

Implement fundamental difference imaging step of convolution and subtraction : D = I - (K.x.T + bg)

Execute fundamental task of convolving template and subtracting it from science image.

Note
The template is taken to be an Image, not a MaskedImage; it therefore has neither variance nor bad pixels
If you convolve the science image, D = (K*I + bg) - T, set invert=False
Instantiated such that background can be a double or Function2D
Returns
Difference image
Parameters
templateImageImage T to convolve with Kernel
scienceMaskedImageImage I to subtract T from
convolutionKernelPSF-matching Kernel used
backgroundDifferential background
invertInvert the output difference image

Definition at line 166 of file ImageSubtract.cc.

172  {
173 
174  boost::timer t;
175  t.restart();
176 
177  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
178  afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
179  convolutionControl.setDoNormalize(false);
180  afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
181  convolutionKernel, convolutionControl);
182 
183  /* Add in background */
184  *(convolvedMaskedImage.getImage()) += background;
185 
186  /* Do actual subtraction */
187  *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
188 
189  /* Invert */
190  if (invert) {
191  *convolvedMaskedImage.getImage() *= -1.0;
192  }
193  convolvedMaskedImage.getMask()->assign(*scienceMaskedImage.getMask());
194  convolvedMaskedImage.getVariance()->assign(*scienceMaskedImage.getVariance());
195 
196  double time = t.elapsed();
197  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
198  "Total compute time to convolve and subtract : %.2f s", time);
199 
200  return convolvedMaskedImage;
201 }

◆ convolveAndSubtract() [6/12]

template<typename PixelT , typename BackgroundT >
lsst::afw::image::MaskedImage<PixelT> lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::Image< PixelT > const &  templateImage,
lsst::afw::image::MaskedImage< PixelT > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
BackgroundT  background,
bool  invert 
)

Execute fundamental task of convolving template and subtracting it from science image.

Note
This version accepts an Image for the template, and is thus faster during convolution
Parameters
templateImageImage to apply convolutionKernel to
scienceMaskedImageMaskedImage from which convolved templateImage is subtracted
convolutionKernelKernel to apply to templateImage
backgroundBackground scalar or function to subtract after convolution
invertInvert the output difference image

Execute fundamental task of convolving template and subtracting it from science image.

Note
The template is taken to be an Image, not a MaskedImage; it therefore has neither variance nor bad pixels
If you convolve the science image, D = (K*I + bg) - T, set invert=False
Instantiated such that background can be a double or Function2D
Returns
Difference image

Definition at line 166 of file ImageSubtract.cc.

172  {
173 
174  boost::timer t;
175  t.restart();
176 
177  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
178  afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
179  convolutionControl.setDoNormalize(false);
180  afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
181  convolutionKernel, convolutionControl);
182 
183  /* Add in background */
184  *(convolvedMaskedImage.getImage()) += background;
185 
186  /* Do actual subtraction */
187  *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
188 
189  /* Invert */
190  if (invert) {
191  *convolvedMaskedImage.getImage() *= -1.0;
192  }
193  convolvedMaskedImage.getMask()->assign(*scienceMaskedImage.getMask());
194  convolvedMaskedImage.getVariance()->assign(*scienceMaskedImage.getVariance());
195 
196  double time = t.elapsed();
197  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
198  "Total compute time to convolve and subtract : %.2f s", time);
199 
200  return convolvedMaskedImage;
201 }

◆ convolveAndSubtract() [7/12]

template lsst::afw::image::MaskedImage< double > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::MaskedImage< double > const &  templateImage,
lsst::afw::image::MaskedImage< double > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
double  background,
bool  invert 
)

◆ convolveAndSubtract() [8/12]

template afwImage::MaskedImage< double > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::MaskedImage< double > const &  templateImage,
lsst::afw::image::MaskedImage< double > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
lsst::afw::math::Function2< double > const &  backgroundFunction,
bool  invert 
)

◆ convolveAndSubtract() [9/12]

template lsst::afw::image::MaskedImage< float > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::MaskedImage< float > const &  templateImage,
lsst::afw::image::MaskedImage< float > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
double  background,
bool  invert 
)

◆ convolveAndSubtract() [10/12]

template afwImage::MaskedImage< float > lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::MaskedImage< float > const &  templateImage,
lsst::afw::image::MaskedImage< float > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
lsst::afw::math::Function2< double > const &  backgroundFunction,
bool  invert 
)

◆ convolveAndSubtract() [11/12]

template<typename PixelT , typename BackgroundT >
afwImage::MaskedImage<PixelT> lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::MaskedImage< PixelT > const &  templateImage,
lsst::afw::image::MaskedImage< PixelT > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
BackgroundT  background,
bool  invert 
)

Implement fundamental difference imaging step of convolution and subtraction : D = I - (K*T + bg) where * denotes convolution.

Execute fundamental task of convolving template and subtracting it from science image.

Note
If you convolve the science image, D = (K*I + bg) - T, set invert=False
The template is taken to be an MaskedImage; this takes c 1.6 times as long as using an Image.
Instantiated such that background can be a double or Function2D
Returns
Difference image
Parameters
templateImageImage T to convolve with Kernel
scienceMaskedImageImage I to subtract T from
convolutionKernelPSF-matching Kernel used
backgroundDifferential background
invertInvert the output difference image

Definition at line 115 of file ImageSubtract.cc.

121  {
122 
123  boost::timer t;
124  t.restart();
125 
126  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
127  afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
128  convolutionControl.setDoNormalize(false);
129  afwMath::convolve(convolvedMaskedImage, templateImage,
130  convolutionKernel, convolutionControl);
131 
132  /* Add in background */
133  *(convolvedMaskedImage.getImage()) += background;
134 
135  /* Do actual subtraction */
136  convolvedMaskedImage -= scienceMaskedImage;
137 
138  /* Invert */
139  if (invert) {
140  convolvedMaskedImage *= -1.0;
141  }
142 
143  double time = t.elapsed();
144  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
145  "Total compute time to convolve and subtract : %.2f s", time);
146 
147  return convolvedMaskedImage;
148 }

◆ convolveAndSubtract() [12/12]

template<typename PixelT , typename BackgroundT >
lsst::afw::image::MaskedImage<PixelT> lsst::ip::diffim::convolveAndSubtract ( lsst::afw::image::MaskedImage< PixelT > const &  templateImage,
lsst::afw::image::MaskedImage< PixelT > const &  scienceMaskedImage,
lsst::afw::math::Kernel const &  convolutionKernel,
BackgroundT  background,
bool  invert 
)

Execute fundamental task of convolving template and subtracting it from science image.

Note
This version accepts a MaskedImage for the template
Parameters
templateImageMaskedImage to apply convolutionKernel to
scienceMaskedImageMaskedImage from which convolved templateImage is subtracted
convolutionKernelKernel to apply to templateImage
backgroundBackground scalar or function to subtract after convolution
invertInvert the output difference image

Execute fundamental task of convolving template and subtracting it from science image.

Note
If you convolve the science image, D = (K*I + bg) - T, set invert=False
The template is taken to be an MaskedImage; this takes c 1.6 times as long as using an Image.
Instantiated such that background can be a double or Function2D
Returns
Difference image

Definition at line 115 of file ImageSubtract.cc.

121  {
122 
123  boost::timer t;
124  t.restart();
125 
126  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
127  afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
128  convolutionControl.setDoNormalize(false);
129  afwMath::convolve(convolvedMaskedImage, templateImage,
130  convolutionKernel, convolutionControl);
131 
132  /* Add in background */
133  *(convolvedMaskedImage.getImage()) += background;
134 
135  /* Do actual subtraction */
136  convolvedMaskedImage -= scienceMaskedImage;
137 
138  /* Invert */
139  if (invert) {
140  convolvedMaskedImage *= -1.0;
141  }
142 
143  double time = t.elapsed();
144  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
145  "Total compute time to convolve and subtract : %.2f s", time);
146 
147  return convolvedMaskedImage;
148 }

◆ imageToEigenMatrix() [1/2]

template Eigen::MatrixXd lsst::ip::diffim::imageToEigenMatrix ( lsst::afw::image::Image< double > const &  )

◆ imageToEigenMatrix() [2/2]

template<typename PixelT >
template Eigen::MatrixXd lsst::ip::diffim::imageToEigenMatrix ( lsst::afw::image::Image< PixelT > const &  img)

Turns a 2-d Image into a 2-d Eigen Matrix.

Turns Image into a 2-D Eigen Matrix.

Parameters
imgImage whose pixel values are read into an Eigen::MatrixXd

Definition at line 62 of file ImageSubtract.cc.

64  {
65  unsigned int rows = img.getHeight();
66  unsigned int cols = img.getWidth();
67  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(rows, cols);
68  for (int y = 0; y != img.getHeight(); ++y) {
69  int x = 0;
70  for (typename afwImage::Image<PixelT>::x_iterator ptr = img.row_begin(y);
71  ptr != img.row_end(y); ++ptr, ++x) {
72  // M is addressed row, col. Need to invert y-axis.
73  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
74  M(rows-y-1,x) = *ptr;
75  }
76  }
77  return M;
78 }

◆ makeAlardLuptonBasisList()

lsst::afw::math::KernelList lsst::ip::diffim::makeAlardLuptonBasisList ( int  halfWidth,
int  nGauss,
std::vector< double > const &  sigGauss,
std::vector< int > const &  degGauss 
)

Build a set of Alard/Lupton basis kernels.

Generate an Alard-Lupton basis set of Kernels.

Note
Should consider implementing as SeparableKernels for additional speed, but this will make the normalization a bit more complicated
Parameters
halfWidthsize is 2*N + 1
nGaussnumber of gaussians
sigGaussWidths of the Gaussian Kernels
degGaussLocal spatial variation of bases
Note
Should consider implementing as SeparableKernels for additional speed, but this will make the normalization a bit more complicated
Returns
Vector of Alard-Lupton Kernels.
Parameters
halfWidthsize is 2*N + 1
nGaussnumber of gaussians
sigGausswidth of the gaussians
degGausslocal spatial variation of gaussians

Definition at line 78 of file BasisLists.cc.

83  {
84  typedef afwMath::Kernel::Pixel Pixel;
85  typedef afwImage::Image<Pixel> Image;
86 
87  if (halfWidth < 1) {
88  throw LSST_EXCEPT(pexExcept::Exception, "halfWidth must be positive");
89  }
90  if (nGauss != static_cast<int>(sigGauss.size())) {
91  throw LSST_EXCEPT(pexExcept::Exception, "sigGauss does not have enough entries");
92  }
93  if (nGauss != static_cast<int>(degGauss.size())) {
94  throw LSST_EXCEPT(pexExcept::Exception, "degGauss does not have enough entries");
95  }
96  int fullWidth = 2 * halfWidth + 1;
97  Image image(geom::Extent2I(fullWidth, fullWidth));
98 
99  afwMath::KernelList kernelBasisList;
100  for (int i = 0; i < nGauss; i++) {
101  /*
102  sigma = FWHM / ( 2 * sqrt(2 * ln(2)) )
103  */
104  double sig = sigGauss[i];
105  int deg = degGauss[i];
106 
107  LOGL_DEBUG("TRACE1.ip.diffim.BasisLists.makeAlardLuptonBasisList",
108  "Gaussian %d : sigma %.2f degree %d", i, sig, deg);
109 
110  afwMath::GaussianFunction2<Pixel> gaussian(sig, sig);
111  afwMath::AnalyticKernel kernel(fullWidth, fullWidth, gaussian);
112  afwMath::PolynomialFunction2<Pixel> polynomial(deg);
113 
114  for (int j = 0, n = 0; j <= deg; j++) {
115  for (int k = 0; k <= (deg - j); k++, n++) {
116  /* for 0th order term, skip polynomial */
117  (void)kernel.computeImage(image, true);
118  if (n == 0) {
119  std::shared_ptr<afwMath::Kernel>
120  kernelPtr(new afwMath::FixedKernel(image));
121  kernelBasisList.push_back(kernelPtr);
122  continue;
123  }
124 
125  /* gaussian to be modified by this term in the polynomial */
126  polynomial.setParameter(n, 1.);
127  (void)kernel.computeImage(image, true);
128  for (int y = 0, v = -halfWidth; y < image.getHeight(); y++, v++) {
129  int u = -halfWidth;
130  for (Image::xy_locator ptr = image.xy_at(0, y),
131  end = image.xy_at(image.getWidth(), y);
132  ptr != end; ++ptr.x(), u++) {
133  /* Evaluate from -1 to 1 */
134  *ptr = *ptr * polynomial(u/static_cast<double>(halfWidth),
135  v/static_cast<double>(halfWidth));
136  }
137  }
138  std::shared_ptr<afwMath::Kernel>
139  kernelPtr(new afwMath::FixedKernel(image));
140  kernelBasisList.push_back(kernelPtr);
141  polynomial.setParameter(n, 0.);
142  }
143  }
144  }
145  return renormalizeKernelList(kernelBasisList);
146  }

◆ makeCentralDifferenceMatrix()

Eigen::MatrixXd lsst::ip::diffim::makeCentralDifferenceMatrix ( int  width,
int  height,
int  stencil,
float  borderPenalty,
bool  fitForBackground 
)

Build a central difference Laplacian regularization matrix for Delta function kernels.

Generate regularization matrix for delta function kernels.

Parameters
widthWidth of basis set you want to regularize
heightHeight of basis set you want to regularize
stencilWhich type of Laplacian approximation to use
borderPenaltyAmount of penalty (if any) to apply to border pixels; > 0
fitForBackgroundFit for differential background?

Definition at line 209 of file BasisLists.cc.

215  {
216 
217  /* 5- or 9-point stencil to approximate the Laplacian; i.e. this is a second
218  * order central finite difference.
219  *
220  * 5-point approximation ignores off-diagonal elements, and is
221  *
222  * f(x-1,y) + f(x+1,y) + f(x,y-1) + f(x,y+1) - 4 f(x,y)
223  *
224  * 0 1 0
225  * 1 -4 1
226  * 0 1 0
227  *
228  *
229  * 9-point approximation includes diagonals and is
230  *
231  * 1 4 1
232  * 4 -20 4 / 6
233  * 1 4 1
234  *
235  */
236 
237  if (borderPenalty < 0)
238  throw LSST_EXCEPT(pexExcept::Exception, "Only border penalty of >= 0 allowed");
239 
240  std::vector<std::vector<float> >
241  coeffs(3, std::vector<float>(3,0));
242 
243  if (stencil == 5) {
244  /* http://www.physics.arizona.edu/~restrepo/475B/Notes/source/node51.html
245  *
246  * This is equivalent to a second order central difference summed along
247  * the x-axis, and then summed along the y-axis.
248  *
249  * http://www.holoborodko.com/pavel/?page_id=239
250  *
251  */
252  coeffs[0][1] = 1.;
253  coeffs[1][0] = 1.;
254  coeffs[1][1] = -4.;
255  coeffs[1][2] = 1.;
256  coeffs[2][1] = 1.;
257  }
258  else if (stencil == 9) {
259  /* http://www.physics.arizona.edu/~restrepo/475B/Notes/source/node52.html */
260  coeffs[0][0] = 1. / 6.;
261  coeffs[0][1] = 4. / 6.;
262  coeffs[0][2] = 1. / 6.;
263  coeffs[1][0] = 4. / 6.;
264  coeffs[1][1] = -20. / 6.;
265  coeffs[1][2] = 4. / 6.;
266  coeffs[2][0] = 1. / 6.;
267  coeffs[2][1] = 4. / 6.;
268  coeffs[2][2] = 1. / 6.;
269  }
270  else {
271  throw LSST_EXCEPT(pexExcept::Exception, "Only 5- or 9-point Laplacian stencils allowed");
272  }
273 
274  int nBgTerms = fitForBackground ? 1 : 0;
275  Eigen::MatrixXd bMat = Eigen::MatrixXd::Zero(width * height + nBgTerms, width * height + nBgTerms);
276 
277  for (int i = 0; i < width*height; i++) {
278  int const x0 = i % width; // the x coord in the kernel image
279  int const y0 = i / width; // the y coord in the kernel image
280  int const distX = width - x0 - 1; // distance from edge of image
281  int const distY = height - y0 - 1; // distance from edge of image
282 
283  if ( (x0 > 0) && (y0 > 0) && (distX > 0) && (distY > 0) ) {
284  for (int dx = -1; dx < 2; dx += 1) {
285  for (int dy = -1; dy < 2; dy += 1) {
286  bMat(i, i + dx + dy * width) += coeffs[dx+1][dy+1];
287  }
288  }
289  }
290  else {
291  bMat(i, i) = borderPenalty;
292  }
293  }
294 
295  if (fitForBackground) {
296  /* Last row / col should have no regularization since its the background term */
297  if (bMat.col(width*height).sum() != 0.) {
298  throw LSST_EXCEPT(pexExcept::Exception, "Error 1 in regularization matrix");
299  }
300  if (bMat.row(width*height).sum() != 0.) {
301  throw LSST_EXCEPT(pexExcept::Exception, "Error 2 in regularization matrix");
302  }
303  }
304 
305  return bMat;
306  }

◆ makeDeltaFunctionBasisList()

lsst::afw::math::KernelList lsst::ip::diffim::makeDeltaFunctionBasisList ( int  width,
int  height 
)

Build a set of Delta Function basis kernels.

Generate a basis set of delta function Kernels.

Note
Total number of basis functions is width*height
Parameters
widthWidth of basis set (cols)
heightHeight of basis set (rows)

Generates a vector of Kernels sized nCols * nRows, where each Kernel has a unique pixel set to value 1.0 with the other pixels valued 0.0. This is the "delta function" basis set.

Returns
Vector of orthonormal delta function Kernels.
Exceptions
lsst::pex::exceptions::DomainErrorif nRows or nCols not positive
Parameters
widthnumber of columns in the set
heightnumber of rows in the set

Definition at line 47 of file BasisLists.cc.

50  {
51  if ((width < 1) || (height < 1)) {
52  throw LSST_EXCEPT(pexExcept::Exception, "nRows and nCols must be positive");
53  }
54  const int signedWidth = static_cast<int>(width);
55  const int signedHeight = static_cast<int>(height);
56  afwMath::KernelList kernelBasisList;
57  for (int row = 0; row < signedHeight; ++row) {
58  for (int col = 0; col < signedWidth; ++col) {
59  std::shared_ptr<afwMath::Kernel>
60  kernelPtr(new afwMath::DeltaFunctionKernel(width, height, geom::Point2I(col,row)));
61  kernelBasisList.push_back(kernelPtr);
62  }
63  }
64  return kernelBasisList;
65  }

◆ makeFiniteDifferenceRegularizationDeprecated()

Eigen::MatrixXd lsst::ip::diffim::makeFiniteDifferenceRegularizationDeprecated ( unsigned int  width,
unsigned int  height,
unsigned int  order,
unsigned int  boundary_style,
unsigned int  difference_style,
bool  printB 
)

Generate regularization matrix for delta function kernels.

Definition at line 509 of file BasisLists.cc.

516  {
517 
518  if (order > 2) throw LSST_EXCEPT(pexExcept::Exception, "Only orders 0..2 allowed");
519 
520  if (boundary_style > 2) {
521  throw LSST_EXCEPT(pexExcept::Exception, "Boundary styles 0..2 defined");
522  }
523  if (difference_style > 1) {
524  throw LSST_EXCEPT(pexExcept::Exception, "Only forward(0), and central(1) difference types defined");
525  }
526 
527  /* what works, and what doesn't */
528  // == good job ==
529  // - order 0, wrapped, forward
530  // - order 1, wrapped or tapered, central or forward
531  // - order 2, wrapped or tapered, forward
532  // == bad job (usually diagongal stripes) ==
533  // - all others
534 
535 
536  /*
537  Instead of Taylor expanding the forward difference approximation of
538  derivatives (N.R. section 18.5) lets just hard code in the expansion of
539  the 1st through 3rd derivatives, which will try and enforce smoothness of
540  0 through 2nd derivatives.
541 
542  A property of the basic "finite difference regularization" is that their
543  rows (column multipliers) sum to 0.
544 
545  Another consideration is to use *multiple* finite difference operators as
546  a constraint.
547 
548  */
549 
550 
551  // ===========================================================================
552  // Get the coeffs for the finite differencing
553  // note: The coeffs are stored 2D although they are essentially 1D entities.
554  // The 2d design was chosen to allow cross-terms to be included,
555  // though they are not yet implemented.
556  //
557  std::vector<std::vector<std::vector<float> > >
558  coeffs(3, std::vector<std::vector<float> >(5, std::vector<float>(5,0)));
559  unsigned int x_cen = 0, y_cen = 0; // center of reqested order coeffs
560  unsigned int x_cen1 = 0, y_cen1 = 0; // center of order 1 coeffs
561  unsigned int x_cen2 = 0, y_cen2 = 0; // center of order 2 coeffs
562  unsigned int x_size = 0, y_size = 0;
563 
564  // forward difference coefficients
565  if (difference_style == 0) {
566 
567  y_cen = x_cen = 0;
568  x_cen1 = y_cen1 = 0;
569  x_cen2 = y_cen2 = 0;
570 
571  x_size = y_size = order + 2;
572 
573  // default forward difference suggested in NR chap 18
574  // 0th order
575  coeffs[0][0][0] = -2;
576  coeffs[0][0][1] = 1;
577  coeffs[0][1][0] = 1;
578  coeffs[0][1][1] = 0;
579 
580  // 1st 2
581  coeffs[1][0][0] = -2;
582  coeffs[1][0][1] = 2;
583  coeffs[1][0][2] = -1;
584  coeffs[1][1][0] = 2;
585  coeffs[1][1][1] = 0;
586  coeffs[1][1][2] = 0;
587  coeffs[1][2][0] = -1;
588  coeffs[1][2][1] = 0;
589  coeffs[1][2][2] = 0;
590 
591  // 2nd 2
592  coeffs[2][0][0] = -2;
593  coeffs[2][0][1] = 3;
594  coeffs[2][0][2] = -3;
595  coeffs[2][0][3] = 1;
596  coeffs[2][1][0] = 3;
597  coeffs[2][1][1] = 0;
598  coeffs[2][1][2] = 0;
599  coeffs[2][1][3] = 0;
600  coeffs[2][2][0] = -3;
601  coeffs[2][2][1] = 0;
602  coeffs[2][2][2] = 0;
603  coeffs[2][2][3] = 0;
604  coeffs[2][3][0] = 1;
605  coeffs[2][3][1] = 0;
606  coeffs[2][3][2] = 0;
607  coeffs[2][3][3] = 0;
608  }
609 
610  // central difference coefficients
611  if (difference_style == 1) {
612 
613  // this is asymmetric and produces diagonal banding in the kernel
614  // from: http://www.holoborodko.com/pavel/?page_id=239
615  if (order == 0) {
616  y_cen = x_cen = 1;
617  x_size = y_size = 3;
618  }
619  coeffs[0][0][0] = 0;
620  coeffs[0][0][1] = -1;
621  coeffs[0][0][2] = 0;
622 
623  coeffs[0][1][0] = -1;
624  coeffs[0][1][1] = 0;
625  coeffs[0][1][2] = 1;
626 
627  coeffs[0][2][0] = 0;
628  coeffs[0][2][1] = 1;
629  coeffs[0][2][2] = 0;
630 
631 
632  // this works well and is largely the same as order=1 forward-diff.
633  // from: http://www.holoborodko.com/pavel/?page_id=239
634  if (order == 1) {
635  y_cen = x_cen = 1;
636  x_size = y_size = 3;
637  }
638  y_cen1 = x_cen1 = 1;
639  coeffs[1][0][0] = 0;
640  coeffs[1][0][1] = 1;
641  coeffs[1][0][2] = 0;
642 
643  coeffs[1][1][0] = 1;
644  coeffs[1][1][1] = -4;
645  coeffs[1][1][2] = 1;
646 
647  coeffs[1][2][0] = 0;
648  coeffs[1][2][1] = 1;
649  coeffs[1][2][2] = 0;
650 
651 
652  // asymmetric and produces diagonal banding in the kernel
653  // from http://www.holoborodko.com/pavel/?page_id=239
654  if (order == 2) {
655  y_cen = x_cen = 2;
656  x_size = y_size = 5;
657  }
658  y_cen2 = x_cen2 = 2;
659  coeffs[2][0][0] = 0;
660  coeffs[2][0][1] = 0;
661  coeffs[2][0][2] = -1;
662  coeffs[2][0][3] = 0;
663  coeffs[2][0][4] = 0;
664 
665  coeffs[2][1][0] = 0;
666  coeffs[2][1][1] = 0;
667  coeffs[2][1][2] = 2;
668  coeffs[2][1][3] = 0;
669  coeffs[2][1][4] = 0;
670 
671  coeffs[2][2][0] = -1;
672  coeffs[2][2][1] = 2;
673  coeffs[2][2][2] = 0;
674  coeffs[2][2][3] = -2;
675  coeffs[2][2][4] = 1;
676 
677  coeffs[2][3][0] = 0;
678  coeffs[2][3][1] = 0;
679  coeffs[2][3][2] = -2;
680  coeffs[2][3][3] = 0;
681  coeffs[2][3][4] = 0;
682 
683  coeffs[2][4][0] = 0;
684  coeffs[2][4][1] = 0;
685  coeffs[2][4][2] = 1;
686  coeffs[2][4][3] = 0;
687  coeffs[2][4][4] = 0;
688  }
689 
690 
691  /* Note we have to add 1 extra (empty) term here because of the differential
692  * background fitting */
693  Eigen::MatrixXd bMat = Eigen::MatrixXd::Zero(width*height+1, width*height+1);
694 
695  /* Forward difference approximation */
696  for (unsigned int i = 0; i < width*height; i++) {
697 
698  unsigned int const x0 = i % width; // the x coord in the kernel image
699  unsigned int const y0 = i / width; // the y coord in the kernel image
700 
701  unsigned int x_edge_distance = (x0 > (width - x0 - 1)) ? width - x0 - 1 : x0;
702  unsigned int y_edge_distance = (y0 > (height - y0 - 1)) ? height - y0 - 1 : y0;
703  unsigned int edge_distance = (x_edge_distance < y_edge_distance) ? x_edge_distance :
704  y_edge_distance;
705 
706  for (unsigned int dx = 0; dx < x_size; dx++) {
707  for (unsigned int dy = 0; dy < y_size; dy++) {
708 
709  // determine where to put this coeff
710 
711  // handle the boundary condition
712  // note: adding width and height in the sum prevents negatives
713  unsigned int x = 0;
714  unsigned int y = 0;
715  double this_coeff = 0;
716 
717  // no-wrapping at edges
718  if (boundary_style == 0) {
719  x = x0 + dx - x_cen;
720  y = y0 + dy - y_cen;
721  if (y > height - 1 || x > width - 1) {
722  continue;
723  }
724  this_coeff = coeffs[order][dx][dy];
725 
726  // wrapping at edges
727  } else if (boundary_style == 1) {
728  x = (width + x0 + dx - x_cen) % width;
729  y = (height + y0 + dy - y_cen) % height;
730  this_coeff = coeffs[order][dx][dy];
731 
732  // order tapering to the edge (just clone wrapping for now)
733  // - use the lowest order possible
734  } else if (boundary_style == 2) {
735 
736  // edge rows and columns ... set to constant
737  if (edge_distance == 0) {
738  x = x0;
739  y = y0;
740  this_coeff = 1;
741  }
742  // in one from edge, use 1st order
743  else if (edge_distance == 1 && order > 0) {
744  x = (width + x0 + dx - x_cen1) % width;
745  y = (height + y0 + dy - y_cen1) % height;
746  if ((dx < 3) && (dy < 3)) { this_coeff = coeffs[1][dx][dy]; }
747  }
748  // in two from edge, use 2st order if order > 1
749  else if (edge_distance == 2 && order > 1){
750  x = (width + x0 + dx - x_cen2) % width;
751  y = (height + y0 + dy - y_cen2) % height;
752  if ((dx < 5) && (dy < 5)) { this_coeff = coeffs[2][dx][dy]; }
753  }
754  // if we're somewhere in the middle
755  else if (edge_distance > order) {
756  x = (width + x0 + dx - x_cen) % width;
757  y = (height + y0 + dy - y_cen) % height;
758  this_coeff = coeffs[order][dx][dy];
759  }
760 
761  }
762 
763  bMat(i, y*width + x) = this_coeff;
764 
765  }
766 
767  }
768 
769  }
770 
771  if (printB) {
772  std::cout << bMat << std::endl;
773  }
774 
775  Eigen::MatrixXd hMat = bMat.transpose() * bMat;
776  return hMat;
777  }

◆ makeForwardDifferenceMatrix()

Eigen::MatrixXd lsst::ip::diffim::makeForwardDifferenceMatrix ( int  width,
int  height,
std::vector< int > const &  orders,
float  borderPenalty,
bool  fitForBackground 
)

Build a forward difference regularization matrix for Delta function kernels.

Generate regularization matrix for delta function kernels.

Parameters
widthWidth of basis set you want to regularize
heightHeight of basis set you want to regularize
ordersWhich derivatives to penalize (1,2,3)
borderPenaltyAmount of penalty (if any) to apply to border pixels; > 0
fitForBackgroundFit for differential background?

Definition at line 311 of file BasisLists.cc.

317  {
318 
319  /*
320  Instead of Taylor expanding the forward difference approximation of
321  derivatives (N.R. section 18.5) lets just hard code in the expansion of
322  the 1st through 3rd derivatives, which will try and enforce smoothness of
323  0 through 2nd derivatives.
324 
325  A property of the basic "finite difference regularization" is that their
326  rows (column multipliers) sum to 0.
327 
328  We taper the order of the constraint as you get close to the boundary.
329  The actual perimeter values are left unconstrained but we might want to
330  consider giving these the same penalties as the other border pixels,
331  which is +1 (since the value gets squared).
332 
333  */
334 
335  if (borderPenalty < 0)
336  throw LSST_EXCEPT(pexExcept::Exception, "Only border penalty of >= 0 allowed");
337 
338  std::vector<std::vector<float> >
339  coeffs(4, std::vector<float>(4,0));
340 
341  // penalize border? this is along the diagonal and gets squared, so sign does not matter
342  coeffs[0][0] = borderPenalty;
343 
344  // penalize zeroth derivative
345  coeffs[1][0] = -1.;
346  coeffs[1][1] = +1.;
347 
348  // penalize first derivative
349  coeffs[2][0] = -1.;
350  coeffs[2][1] = +2.;
351  coeffs[2][2] = -1.;
352 
353  // penalize second derivative
354  coeffs[3][0] = -1.;
355  coeffs[3][1] = +3.;
356  coeffs[3][2] = -3.;
357  coeffs[3][3] = +1.;
358 
359  int nBgTerms = fitForBackground ? 1 : 0;
360  Eigen::MatrixXd bTot = Eigen::MatrixXd::Zero(width * height + nBgTerms, width * height + nBgTerms);
361 
362  std::vector<int>::const_iterator order;
363  for (order = orders.begin(); order != orders.end(); order++) {
364  if ((*order < 1) || (*order > 3))
365  throw LSST_EXCEPT(pexExcept::Exception, "Only orders 1..3 allowed");
366 
367  Eigen::MatrixXd bMatX = Eigen::MatrixXd::Zero(width * height + nBgTerms,
368  width * height + nBgTerms);
369  Eigen::MatrixXd bMatY = Eigen::MatrixXd::Zero(width * height + nBgTerms,
370  width * height + nBgTerms);
371 
372  for (int i = 0; i < width*height; i++) {
373  int const x0 = i % width; // the x coord in the kernel image
374  int const y0 = i / width; // the y coord in the kernel image
375 
376  int distX = width - x0 - 1; // distance from edge of image
377  int orderToUseX = std::min(distX, *order);
378  for (int j = 0; j < orderToUseX+1; j++) {
379  bMatX(i, i + j) = coeffs[orderToUseX][j];
380  }
381 
382  int distY = height - y0 - 1; // distance from edge of image
383  int orderToUseY = std::min(distY, *order);
384  for (int j = 0; j < orderToUseY+1; j++) {
385  bMatY(i, i + j * width) = coeffs[orderToUseY][j];
386  }
387  }
388  bTot += bMatX;
389  bTot += bMatY;
390  }
391 
392  if (fitForBackground) {
393  /* Last row / col should have no regularization since its the background term */
394  if (bTot.col(width*height).sum() != 0.) {
395  throw LSST_EXCEPT(pexExcept::Exception, "Error in regularization matrix");
396  }
397  if (bTot.row(width*height).sum() != 0.) {
398  throw LSST_EXCEPT(pexExcept::Exception, "Error in regularization matrix");
399  }
400  }
401 
402  return bTot;
403  }

◆ makeKernelCandidate() [1/2]

template<typename PixelT >
std::shared_ptr<KernelCandidate<PixelT> > lsst::ip::diffim::makeKernelCandidate ( float const  xCenter,
float const  yCenter,
std::shared_ptr< afw::image::MaskedImage< PixelT > > const &  templateMaskedImage,
std::shared_ptr< afw::image::MaskedImage< PixelT > > const &  scienceMaskedImage,
daf::base::PropertySet const &  ps 
)

Return a KernelCandidate pointer of the right sort.

Parameters
xCenterX-center of candidate
yCenterY-center of candidate
templateMaskedImageTemplate subimage
scienceMaskedImageScience image subimage
psPropertySet for creation of rating

Definition at line 220 of file KernelCandidate.h.

224  {
225 
226  return std::shared_ptr<KernelCandidate<PixelT>>(new KernelCandidate<PixelT>(xCenter, yCenter,
227  templateMaskedImage,
228  scienceMaskedImage,
229  ps));
230  }

◆ makeKernelCandidate() [2/2]

template<typename PixelT >
std::shared_ptr<KernelCandidate<PixelT> > lsst::ip::diffim::makeKernelCandidate ( std::shared_ptr< afw::table::SourceRecord > const &  source,
std::shared_ptr< afw::image::MaskedImage< PixelT > > const &  templateMaskedImage,
std::shared_ptr< afw::image::MaskedImage< PixelT > > const &  scienceMaskedImage,
daf::base::PropertySet const &  ps 
)

Return a KernelCandidate pointer of the right sort.

Parameters
sourceafw::table::SourceRecord used to construct the KernelCandidate
templateMaskedImageTemplate subimage
scienceMaskedImageScience image subimage
psPropertySet for creation of rating

Definition at line 245 of file KernelCandidate.h.

248  {
249 
250  return std::shared_ptr<KernelCandidate<PixelT>>(new KernelCandidate<PixelT>(source,
251  templateMaskedImage,
252  scienceMaskedImage,
253  ps));
254  }

◆ makeRegularizationMatrix()

Eigen::MatrixXd lsst::ip::diffim::makeRegularizationMatrix ( lsst::daf::base::PropertySet const &  ps)

Build a regularization matrix for Delta function kernels.

Parameters
psPropertySet dictating which type of matrix to make
Note
Calls either makeForwardDifferenceMatrix or makeCentralDifferenceMatrix based on the PropertySet config.

Definition at line 149 of file BasisLists.cc.

151  {
152 
153  /* NOTES
154  *
155  * The 3-point first derivative central difference (Laplacian) yields
156  * diagonal stripes and should not be used.
157 
158  coeffs[0][1] = -1. / 2.;
159  coeffs[1][0] = -1. / 2.;
160  coeffs[1][2] = 1. / 2.;
161  coeffs[2][1] = 1. / 2.;
162 
163  * The 5-point second derivative central difference (Laplacian) looks great.
164  * For smaller lambdas you need a higher border penalty. In general, if you
165  * decrease the regularization strength you should increase the border
166  * penalty to avoid noise in the border pixels. This has been used to value
167  * 10 and things still look OK.
168 
169  * The 9-point second derivative central difference (Laplacian with off
170  * diagonal terms) also looks great. Seems to have slightly higher-valued
171  * border pixels so make boundary larger if using this. E.g. 1.5.
172 
173  * The forward finite difference, first derivative term works good
174  *
175  * The forward finite difference, second derivative term is slightly banded
176  * in LLC
177  *
178  * The forward finite difference, third derivative term is highly banded in
179  * LLC and should not be used.
180  *
181  */
182 
183  std::string regularizationType = ps.getAsString("regularizationType");
184  int width = ps.getAsInt("kernelSize");
185  int height = ps.getAsInt("kernelSize");
186  float borderPenalty = ps.getAsDouble("regularizationBorderPenalty");
187  bool fitForBackground = ps.getAsBool("fitForBackground");
188 
189  Eigen::MatrixXd bMat;
190  if (regularizationType == "centralDifference") {
191  int stencil = ps.getAsInt("centralRegularizationStencil");
192  bMat = makeCentralDifferenceMatrix(width, height, stencil, borderPenalty, fitForBackground);
193  }
194  else if (regularizationType == "forwardDifference") {
195  std::vector<int> orders = ps.getArray<int>("forwardRegularizationOrders");
196  bMat = makeForwardDifferenceMatrix(width, height, orders, borderPenalty, fitForBackground);
197  }
198  else {
199  throw LSST_EXCEPT(pexExcept::Exception, "regularizationType not recognized");
200  }
201 
202  Eigen::MatrixXd hMat = bMat.transpose() * bMat;
203  return hMat;
204  }

◆ maskToEigenMatrix()

Eigen::MatrixXi lsst::ip::diffim::maskToEigenMatrix ( lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &  mask)

Definition at line 80 of file ImageSubtract.cc.

82  {
83  unsigned int rows = mask.getHeight();
84  unsigned int cols = mask.getWidth();
85  Eigen::MatrixXi M = Eigen::MatrixXi::Zero(rows, cols);
86  for (int y = 0; y != mask.getHeight(); ++y) {
87  int x = 0;
88  for (afwImage::Mask<lsst::afw::image::MaskPixel>::x_iterator ptr = mask.row_begin(y);
89  ptr != mask.row_end(y); ++ptr, ++x) {
90  // M is addressed row, col. Need to invert y-axis.
91  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
92  M(rows-y-1,x) = *ptr;
93  }
94  }
95  return M;
96 }

◆ NEGCENTXPAR()

int const lsst::ip::diffim::NEGCENTXPAR ( 0  )

◆ NEGCENTYPAR()

int const lsst::ip::diffim::NEGCENTYPAR ( 1  )

◆ NEGFLUXPAR()

int const lsst::ip::diffim::NEGFLUXPAR ( 2  )

◆ POSCENTXPAR()

int const lsst::ip::diffim::POSCENTXPAR ( 3  )

◆ POSCENTYPAR()

int const lsst::ip::diffim::POSCENTYPAR ( 4  )

◆ POSFLUXPAR()

int const lsst::ip::diffim::POSFLUXPAR ( 5  )

◆ PYBIND11_MODULE() [1/8]

lsst::ip::diffim::PYBIND11_MODULE ( _dipoleAlgorithms  ,
mod   
)

Definition at line 140 of file dipoleAlgorithms.cc.

140  {
141  py::module::import("lsst.afw.table");
142  py::module::import("lsst.meas.base");
143  py::module::import("lsst.pex.config");
144 
145  declareDipoleCentroidControl(mod);
146  declareDipoleFluxControl(mod);
147  declareDipolePsfFluxControl(mod);
148  declareDipoleCentroidAlgorithm(mod);
149  declareDipoleFluxAlgorithm(mod);
150  declareNaiveDipoleFlux(mod);
151  declareNaiveDipoleCentroid(mod);
152  declarePsfDipoleFlux(mod);
153 }

◆ PYBIND11_MODULE() [2/8]

lsst::ip::diffim::PYBIND11_MODULE ( basisLists  ,
mod   
)

Definition at line 39 of file basisLists.cc.

39  {
40  py::module::import("lsst.afw.math");
41  py::module::import("lsst.daf.base");
42 
43  mod.def("makeDeltaFunctionBasisList", &makeDeltaFunctionBasisList, "width"_a, "height"_a);
44  mod.def("makeRegularizationMatrix", &makeRegularizationMatrix, "ps"_a);
45  mod.def("makeForwardDifferenceMatrix", &makeForwardDifferenceMatrix, "width"_a, "height"_a, "orders"_a,
46  "borderPenalty"_a, "fitForBackground"_a);
47  mod.def("makeCentralDifferenceMatrix", &makeCentralDifferenceMatrix, "width"_a, "height"_a, "stencil"_a,
48  "borderPenalty"_a, "fitForBackground"_a);
49  mod.def("renormalizeKernelList", &renormalizeKernelList, "kernelListIn"_a);
50  mod.def("makeAlardLuptonBasisList", &makeAlardLuptonBasisList, "halfWidth"_a, "nGauss"_a, "sigGauss"_a,
51  "degGauss"_a);
52 }

◆ PYBIND11_MODULE() [3/8]

lsst::ip::diffim::PYBIND11_MODULE ( findSetBits  ,
mod   
)

Definition at line 58 of file findSetBits.cc.

58  {
59  py::module::import("lsst.afw.image");
60 
61  declareFindSetBits<afw::image::Mask<afw::image::MaskPixel>>(mod, "U");
62 }

◆ PYBIND11_MODULE() [4/8]

lsst::ip::diffim::PYBIND11_MODULE ( imageStatistics  ,
mod   
)

Definition at line 73 of file imageStatistics.cc.

73  {
74  py::module::import("lsst.daf.base");
75  py::module::import("lsst.afw.image");
76 
77  declareImageStatistics<int>(mod, "I");
78  declareImageStatistics<float>(mod, "F");
79  declareImageStatistics<double>(mod, "D");
80 }

◆ PYBIND11_MODULE() [5/8]

lsst::ip::diffim::PYBIND11_MODULE ( imageSubtract  ,
mod   
)

Definition at line 72 of file imageSubtract.cc.

72  {
73  py::module::import("lsst.afw.image");
74  py::module::import("lsst.afw.math");
75 
76  declareConvolveAndSubtract<float, double>(mod);
77  declareConvolveAndSubtract<float, afw::math::Function2<double> const &>(mod);
78 }

◆ PYBIND11_MODULE() [6/8]

lsst::ip::diffim::PYBIND11_MODULE ( kernelCandidate  ,
mod   
)

Definition at line 116 of file kernelCandidate.cc.

116  {
117  py::module::import("lsst.afw.image");
118  py::module::import("lsst.afw.math");
119  py::module::import("lsst.afw.table");
120  py::module::import("lsst.daf.base");
121 
122  declareKernelCandidate<float>(mod, "F");
123 }

◆ PYBIND11_MODULE() [7/8]

lsst::ip::diffim::PYBIND11_MODULE ( kernelCandidateDetection  ,
mod   
)

Definition at line 62 of file kernelCandidateDetection.cc.

62  {
63  py::module::import("lsst.afw.image");
64  py::module::import("lsst.afw.detection");
65  py::module::import("lsst.daf.base");
66 
67  declareKernelCandidateDetection<float>(mod, "F");
68 }

◆ PYBIND11_MODULE() [8/8]

lsst::ip::diffim::PYBIND11_MODULE ( kernelSolution  ,
mod   
)

Definition at line 176 of file kernelSolution.cc.

176  {
177  py::module::import("lsst.afw.geom");
178  py::module::import("lsst.afw.image");
179  py::module::import("lsst.afw.math");
180  py::module::import("lsst.daf.base");
181 
182  declareKernelSolution(mod);
183  declareStaticKernelSolution<float>(mod, "F");
184  declareMaskedKernelSolution<float>(mod, "F");
185  declareRegularizedKernelSolution<float>(mod, "F");
186  declareSpatialKernelSolution(mod);
187 }

◆ renormalizeKernelList()

lsst::afw::math::KernelList lsst::ip::diffim::renormalizeKernelList ( lsst::afw::math::KernelList const &  kernelListIn)

Renormalize a list of basis kernels.

Rescale an input set of kernels.

Note
Renormalization means make Ksum_0 = 1.0, Ksum_i = 0.0, K_i.dot.K_i = 1.0
Output list of shared pointers to FixedKernels
Parameters
kernelListIninput list of basis kernels
Note
Images are checked for their current kernel sum. If it is larger than std::numeric_limits<double>::epsilon(), the kernel is first divided by the kernel sum, giving it a kSum of 1.0, and then the first (normalized) component is subtracted from it, giving it a kSum of 0.0.
Returns
Vector of renormalized kernels
Parameters
kernelListInInput list to be renormalized

Definition at line 414 of file BasisLists.cc.

416  {
417  typedef afwMath::Kernel::Pixel Pixel;
418  typedef afwImage::Image<Pixel> Image;
419  double kSum;
420 
421  /*
422 
423  We want all the bases except for the first to sum to 0.0. This allows
424  us to achieve kernel flux conservation (Ksum) across the image since all
425  the power will be in the first term, which will not vary spatially.
426 
427  K(x,y) = Ksum * B_0 + Sum_i : a(x,y) * B_i
428 
429  To do this, normalize all Kernels to sum = 1. and subtract B_0 from all
430  subsequent kenrels.
431 
432  To get an idea of the relative contribution of each of these basis
433  functions later on down the line, lets also normalize them such that
434 
435  Sum(B_i) == 0.0 *and*
436  B_i * B_i == 1.0
437 
438  For completeness
439 
440  Sum(B_0) == 1.0
441  B_0 * B_0 != 1.0
442 
443  */
444  afwMath::KernelList kernelListOut;
445  if (kernelListIn.size() == 0) {
446  return kernelListOut;
447  }
448 
449  Image image0(kernelListIn[0]->getDimensions());
450  for (unsigned int i = 0; i < kernelListIn.size(); i++) {
451  if (i == 0) {
452  /* Make sure that it is normalized to kSum 1. */
453  (void)kernelListIn[i]->computeImage(image0, true);
454  std::shared_ptr<afwMath::Kernel>
455  kernelPtr(new afwMath::FixedKernel(image0));
456  kernelListOut.push_back(kernelPtr);
457 
458  continue;
459  }
460 
461  /* Don't normalize here */
462  Image image(kernelListIn[i]->getDimensions());
463  (void)kernelListIn[i]->computeImage(image, false);
464  /* image.writeFits(str(boost::format("in_k%d.fits") % i)); */
465 
466  /* Check the kernel sum; if its close to zero don't do anything */
467  kSum = 0.;
468  for (int y = 0; y < image.getHeight(); y++) {
469  for (Image::xy_locator ptr = image.xy_at(0, y), end = image.xy_at(image.getWidth(), y);
470  ptr != end; ++ptr.x()) {
471  kSum += *ptr;
472  }
473  }
474 
475  /* std::numeric_limits<float>::epsilon() ~ e-7
476  std::numeric_limits<double>::epsilon() ~ e-16
477 
478  If we end up with 2e-16 kernel sum, this still blows up the kernel values.
479  Even though the kernels are double, use the float limits instead
480  */
481  if (fabs(kSum) > std::numeric_limits<float>::epsilon()) {
482  image /= kSum;
483  image -= image0;
484  }
485 
486  /* Finally, rescale such that the inner product is 1 */
487  kSum = 0.;
488  for (int y = 0; y < image.getHeight(); y++) {
489  for (Image::xy_locator ptr = image.xy_at(0, y), end = image.xy_at(image.getWidth(), y);
490  ptr != end; ++ptr.x()) {
491  kSum += *ptr * *ptr;
492  }
493  }
494  image /= std::sqrt(kSum);
495  /* image.writeFits(str(boost::format("out_k%d.fits") % i)); */
496 
497  std::shared_ptr<afwMath::Kernel>
498  kernelPtr(new afwMath::FixedKernel(image));
499  kernelListOut.push_back(kernelPtr);
500  }
501  return kernelListOut;
502  }

Variable Documentation

◆ Control

lsst.ip.diffim.Control

Definition at line 49 of file __init__.py.

◆ DipoleCentroidControl

lsst.ip.diffim.DipoleCentroidControl

Definition at line 49 of file __init__.py.

◆ DipoleFluxControl

lsst.ip.diffim.DipoleFluxControl

Definition at line 50 of file __init__.py.

◆ executionOrder

lsst.ip.diffim.executionOrder

Definition at line 49 of file __init__.py.

◆ NaiveDipoleCentroid

lsst.ip.diffim.NaiveDipoleCentroid

Definition at line 49 of file __init__.py.

◆ NaiveDipoleFlux

lsst.ip.diffim.NaiveDipoleFlux

Definition at line 50 of file __init__.py.

◆ PsfDipoleFlux

lsst.ip.diffim.PsfDipoleFlux

Definition at line 51 of file __init__.py.

◆ PsfDipoleFluxControl

lsst.ip.diffim.PsfDipoleFluxControl

Definition at line 51 of file __init__.py.

y
int y
Definition: SpanSet.cc:49
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::afw::math::ConvolutionControl::setDoNormalize
void setDoNormalize(bool doNormalize)
Definition: ConvolveImage.h:66
std::string
STL class.
std::shared_ptr
STL class.
lsst::afw::image::Mask
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
std::fabs
T fabs(T... args)
lsst::ip::diffim::makeRegularizationMatrix
Eigen::MatrixXd makeRegularizationMatrix(lsst::daf::base::PropertySet const &ps)
Build a regularization matrix for Delta function kernels.
Definition: BasisLists.cc:149
lsst::afw::math::DeltaFunctionKernel
A kernel that has only one non-zero pixel (of value 1)
Definition: Kernel.h:690
lsst::meas::modelfit::Pixel
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
std::vector< std::shared_ptr< Kernel > >
std::vector::size
T size(T... args)
lsst::afw::math::FixedKernel
A kernel created from an Image.
Definition: Kernel.h:518
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition: ImageBase.h:438
lsst::afw::math::GaussianFunction2
2-dimensional Gaussian
Definition: FunctionLibrary.h:201
mask
afw::table::Key< afw::table::Array< MaskPixelT > > mask
Definition: HeavyFootprint.cc:217
lsst::afw::image::ImageBase::getHeight
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:335
std::sqrt
T sqrt(T... args)
end
int end
Definition: BoundedField.cc:105
std::vector::push_back
T push_back(T... args)
lsst::afw::image::MaskedImage< PixelT >
ast::detail::source
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224
lsst.pipe.drivers.visualizeVisit.background
background
Definition: visualizeVisit.py:37
image
afw::table::Key< afw::table::Array< ImagePixelT > > image
Definition: HeavyFootprint.cc:216
lsst::ip::diffim::makeDeltaFunctionBasisList
lsst::afw::math::KernelList makeDeltaFunctionBasisList(int width, int height)
Build a set of Delta Function basis kernels.
Definition: BasisLists.cc:47
std::cout
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y'th row.
Definition: ImageBase.h:441
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::afw::math.convolveImage.convolveImageContinued.convolve
convolve
Definition: convolveImageContinued.py:28
ptr
uint64_t * ptr
Definition: RangeSet.cc:88
std::numeric_limits::epsilon
T epsilon(T... args)
lsst::afw::image::MaskedImage::getVariance
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1090
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::min
T min(T... args)
lsst::ip::diffim::makeCentralDifferenceMatrix
Eigen::MatrixXd makeCentralDifferenceMatrix(int width, int height, int stencil, float borderPenalty, bool fitForBackground)
Build a central difference Laplacian regularization matrix for Delta function kernels.
Definition: BasisLists.cc:209
lsst::afw::math::AnalyticKernel
A kernel described by a function.
Definition: Kernel.h:582
lsst::afw::image::ImageBase::getDimensions
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition: ImageBase.h:393
std::endl
T endl(T... args)
lsst::afw::math::ConvolutionControl
Parameters to control convolution.
Definition: ConvolveImage.h:50
lsst::ip::diffim::makeAlardLuptonBasisList
lsst::afw::math::KernelList makeAlardLuptonBasisList(int halfWidth, int nGauss, std::vector< double > const &sigGauss, std::vector< int > const &degGauss)
Build a set of Alard/Lupton basis kernels.
Definition: BasisLists.cc:78
row
int row
Definition: CR.cc:145
LOGL_DEBUG
#define LOGL_DEBUG(logger, message...)
Definition: Log.h:504
std::vector::begin
T begin(T... args)
col
int col
Definition: CR.cc:144
lsst::geom::Point< int, 2 >
lsst::afw::image::MaskedImage::getImage
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1057
std::time
T time(T... args)
lsst::pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
std::vector::end
T end(T... args)
lsst::afw::math::PolynomialFunction2
2-dimensional polynomial function with cross terms
Definition: FunctionLibrary.h:449
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
lsst::ip::diffim::renormalizeKernelList
lsst::afw::math::KernelList renormalizeKernelList(lsst::afw::math::KernelList const &kernelListIn)
Renormalize a list of basis kernels.
Definition: BasisLists.cc:414
lsst::afw::image::MaskedImage::getMask
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
Definition: MaskedImage.h:1069
lsst::geom::Extent< int, 2 >
lsst::afw::image::MaskedImage::getDimensions
lsst::geom::Extent2I getDimensions() const
Definition: MaskedImage.h:1096
lsst::afw::image::ImageBase::getWidth
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:333
lsst::sphgeom::invert
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
Definition: Relationship.h:55
lsst::ip::diffim::makeForwardDifferenceMatrix
Eigen::MatrixXd makeForwardDifferenceMatrix(int width, int height, std::vector< int > const &orders, float borderPenalty, bool fitForBackground)
Build a forward difference regularization matrix for Delta function kernels.
Definition: BasisLists.cc:311
lsst::afw::math::Kernel::Pixel
double Pixel
Definition: Kernel.h:113
Generated on Wed Jun 24 2020 18:10:34 for LSSTApplications by   doxygen 1.8.18