LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Namespaces | Classes | Typedefs | Functions | Variables
lsst::ip::diffim Namespace Reference

Namespaces

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

Classes

class  DipoleCentroidControl
 
class  DipoleFluxControl
 
class  PsfDipoleFluxControl
 C++ control object for PSF dipole fluxes. More...
 
class  DipoleCentroidAlgorithm
 Intermediate base class for algorithms that compute a centroid. More...
 
class  DipoleFluxAlgorithm
 Intermediate base class for algorithms that compute a flux. More...
 
class  NaiveDipoleFlux
 
class  NaiveDipoleCentroid
 Intermediate base class for algorithms that compute a centroid. More...
 
class  PsfDipoleFlux
 Implementation of Psf dipole flux. More...
 
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  StaticKernelSolution
 
class  MaskedKernelSolution
 
class  RegularizedKernelSolution
 
class  SpatialKernelSolution
 
class  MinimizeDipoleChi2
 Class to minimize PsfDipoleFlux; this is the object that Minuit minimizes. More...
 

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< PixelTconvolveAndSubtract (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< PixelTconvolveAndSubtract (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< PixelTconvolveAndSubtract (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< PixelTconvolveAndSubtract (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 1580 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());
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 }
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition: ImageBase.h:356
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1051
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
Definition: MaskedImage.h:1030
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1018
Parameters to control convolution.
Definition: ConvolveImage.h:50
void setDoNormalize(bool doNormalize)
Definition: ConvolveImage.h:66
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.
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
Definition: Relationship.h:55
T time(T... args)

◆ 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
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());
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());
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 }
lsst::geom::Extent2I getDimensions() const
Definition: MaskedImage.h:1057

◆ 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
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());
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;
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 }
double x
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:48
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:294
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:296
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition: ImageBase.h:401
x_iterator row_end(int y) const
Return an x_iterator to the end of the y'th row.
Definition: ImageBase.h:404
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51

◆ 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  {
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) {
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 (auto ptr = image.row_begin(y), end = image.row_end(y); ptr != end; ++ptr, ++u) {
131  /* Evaluate from -1 to 1 */
132  *ptr = *ptr * polynomial(u/static_cast<double>(halfWidth),
133  v/static_cast<double>(halfWidth));
134  }
135  }
137  kernelPtr(new afwMath::FixedKernel(image));
138  kernelBasisList.push_back(kernelPtr);
139  polynomial.setParameter(n, 0.);
140  }
141  }
142  }
143  return renormalizeKernelList(kernelBasisList);
144  }
int end
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
A kernel described by a function.
Definition: Kernel.h:535
A kernel created from an Image.
Definition: Kernel.h:471
2-dimensional polynomial function with cross terms
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
lsst::afw::math::KernelList renormalizeKernelList(lsst::afw::math::KernelList const &kernelListIn)
Renormalize a list of basis kernels.
Definition: BasisLists.cc:412
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
T push_back(T... args)
T size(T... args)

◆ 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 207 of file BasisLists.cc.

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

◆ 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) {
60  kernelPtr(new afwMath::DeltaFunctionKernel(width, height, geom::Point2I(col,row)));
61  kernelBasisList.push_back(kernelPtr);
62  }
63  }
64  return kernelBasisList;
65  }
A kernel that has only one non-zero pixel (of value 1)
Definition: Kernel.h:643
int row
Definition: CR.cc:145
int col
Definition: CR.cc:144

◆ 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 505 of file BasisLists.cc.

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

◆ 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 309 of file BasisLists.cc.

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

◆ 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  }
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224

◆ 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 147 of file BasisLists.cc.

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

◆ 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;
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 }
afw::table::Key< afw::table::Array< MaskPixelT > > mask
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77

◆ NEGCENTXPAR()

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

◆ NEGCENTYPAR()

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

◆ NEGFLUXPAR()

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

◆ POSCENTXPAR()

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

◆ POSCENTYPAR()

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

◆ POSFLUXPAR()

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

◆ 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 }
lsst::afw::math::KernelList makeDeltaFunctionBasisList(int width, int height)
Build a set of Delta Function basis kernels.
Definition: BasisLists.cc:47
Eigen::MatrixXd makeRegularizationMatrix(lsst::daf::base::PropertySet const &ps)
Build a regularization matrix for Delta function kernels.
Definition: BasisLists.cc:147
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

◆ 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 412 of file BasisLists.cc.

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

Variable Documentation

◆ Control

lsst.ip.diffim.Control

Definition at line 48 of file __init__.py.

◆ DipoleCentroidControl

Definition at line 48 of file __init__.py.

◆ DipoleFluxControl

Definition at line 49 of file __init__.py.

◆ executionOrder

lsst.ip.diffim.executionOrder

Definition at line 48 of file __init__.py.

◆ NaiveDipoleCentroid

Definition at line 48 of file __init__.py.

◆ NaiveDipoleFlux

Definition at line 49 of file __init__.py.

◆ PsfDipoleFlux

Definition at line 50 of file __init__.py.

◆ PsfDipoleFluxControl

Definition at line 50 of file __init__.py.