LSST Applications 24.1.5,g02d81e74bb+fa3a7a026e,g180d380827+a53a32eff8,g2079a07aa2+86d27d4dc4,g2305ad1205+c0501b3732,g295015adf3+7d3e92f0ec,g2bbee38e9b+0e5473021a,g337abbeb29+0e5473021a,g33d1c0ed96+0e5473021a,g3a166c0a6a+0e5473021a,g3ddfee87b4+5dd1654d75,g48712c4677+3bf1020dcb,g487adcacf7+065c13d9cf,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+d7ac436cfb,g5a732f18d5+53520f316c,g64a986408d+fa3a7a026e,g858d7b2824+fa3a7a026e,g8a8a8dda67+585e252eca,g99cad8db69+a5a909b84f,g9ddcbc5298+9a081db1e4,ga1e77700b3+15fc3df1f7,ga8c6da7877+4cf350ccb2,gb0e22166c9+60f28cb32d,gba4ed39666+c2a2e4ac27,gbb8dafda3b+f991a0b59f,gc120e1dc64+9ccbfdb8be,gc28159a63d+0e5473021a,gcf0d15dbbd+5dd1654d75,gd96a1ce819+42fd0ee607,gdaeeff99f8+f9a426f77a,ge6526c86ff+0d71447b4b,ge79ae78c31+0e5473021a,gee10cc3b42+585e252eca,gff1a9f87cc+fa3a7a026e
LSST Data Management Base Package
Loading...
Searching...
No Matches
Namespaces | Classes | Typedefs | Functions | Variables
lsst::ip::diffim Namespace Reference

Namespaces

namespace  computeSpatiallySampledMetrics
 
namespace  dcrModel
 
namespace  detail
 
namespace  detectAndMeasure
 
namespace  diffimLib
 
namespace  dipoleFitTask
 
namespace  dipoleMeasurement
 
namespace  getTemplate
 
namespace  imageDecorrelation
 
namespace  imageMapReduce
 
namespace  makeKernel
 
namespace  makeKernelBasisList
 
namespace  metrics
 
namespace  modelPsfMatch
 
namespace  psfMatch
 
namespace  subtractImages
 
namespace  utils
 
namespace  version
 

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  ImageStatistics
 Class to calculate difference image statistics. More...
 
class  KernelCandidate
 Class stored in SpatialCells for spatial Kernel fitting. 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.
 
Eigen::MatrixXd makeRegularizationMatrix (lsst::daf::base::PropertySet const &ps)
 Build a regularization matrix for Delta function kernels.
 
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.
 
Eigen::MatrixXd makeCentralDifferenceMatrix (int width, int height, int stencil, float borderPenalty, bool fitForBackground)
 Build a central difference Laplacian regularization matrix for Delta function kernels.
 
lsst::afw::math::KernelList renormalizeKernelList (lsst::afw::math::KernelList const &kernelListIn)
 Renormalize a list of basis kernels.
 
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.
 
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.
 
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.
 
template<typename PixelT >
Eigen::MatrixXd imageToEigenMatrix (lsst::afw::image::Image< PixelT > const &img)
 Turns a 2-d Image into a 2-d Eigen Matrix.
 
Eigen::MatrixXi maskToEigenMatrix (lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
 
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.
 
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.
 
void wrapBasisLists (WrapperCollection &wrappers)
 
void wrapDipoleAlgorithms (WrapperCollection &wrappers)
 
void wrapImageStatistics (WrapperCollection &wrappers)
 
void wrapImageSubtract (WrapperCollection &wrappers)
 
void wrapKernelCandidate (WrapperCollection &wrappers)
 
void wrapKernelSolution (WrapperCollection &wrappers)
 
 PYBIND11_MODULE (_ipdiffimLib, 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.
 
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.
 
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)
 
template Eigen::MatrixXd imageToEigenMatrix (lsst::afw::image::Image< float > const &)
 
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
 
 deprecated
 
 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 410 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::cpu_timer t;
175
176 afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
177 afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
178 convolutionControl.setDoNormalize(false);
179 afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
180 convolutionKernel, convolutionControl);
181
182 /* Add in background */
183 *(convolvedMaskedImage.getImage()) += background;
184
185 /* Do actual subtraction */
186 *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
187
188 /* Invert */
189 if (invert) {
190 *convolvedMaskedImage.getImage() *= -1.0;
191 }
192 convolvedMaskedImage.getMask()->assign(*scienceMaskedImage.getMask());
193 convolvedMaskedImage.getVariance()->assign(*scienceMaskedImage.getVariance());
194
195 t.stop();
196 double time = 1e-9 * t.elapsed().wall;
197 LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
198 "Total compute time to convolve and subtract : %.2f s", time);
199
200 return convolvedMaskedImage;
201}
LOGL_DEBUG
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition Log.h:515
lsst::afw::image::ImageBase::getDimensions
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition ImageBase.h:356
lsst::afw::image::ImageBase::assign
void assign(ImageBase const &rhs, lsst::geom::Box2I const &bbox=lsst::geom::Box2I(), ImageOrigin origin=PARENT)
Copy pixels from another image to a specified subregion of this image.
Definition Image.cc:153
lsst::afw::image::MaskedImage
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
lsst::afw::image::MaskedImage::getVariance
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition MaskedImage.h:1052
lsst::afw::image::MaskedImage::getMask
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
Definition MaskedImage.h:1031
lsst::afw::image::MaskedImage::getImage
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition MaskedImage.h:1019
lsst::afw::math::ConvolutionControl
Parameters to control convolution.
Definition ConvolveImage.h:50
lsst::afw::math::ConvolutionControl::setDoNormalize
void setDoNormalize(bool doNormalize)
Definition ConvolveImage.h:66
lsst::afw::math::convolve
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.
Definition ConvolveImage.cc:181

◆ 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::cpu_timer t;
175
176 afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
177 afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
178 convolutionControl.setDoNormalize(false);
179 afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
180 convolutionKernel, convolutionControl);
181
182 /* Add in background */
183 *(convolvedMaskedImage.getImage()) += background;
184
185 /* Do actual subtraction */
186 *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
187
188 /* Invert */
189 if (invert) {
190 *convolvedMaskedImage.getImage() *= -1.0;
191 }
192 convolvedMaskedImage.getMask()->assign(*scienceMaskedImage.getMask());
193 convolvedMaskedImage.getVariance()->assign(*scienceMaskedImage.getVariance());
194
195 t.stop();
196 double time = 1e-9 * t.elapsed().wall;
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::cpu_timer t;
124
125 afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
126 afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
127 convolutionControl.setDoNormalize(false);
128 afwMath::convolve(convolvedMaskedImage, templateImage,
129 convolutionKernel, convolutionControl);
130
131 /* Add in background */
132 *(convolvedMaskedImage.getImage()) += background;
133
134 /* Do actual subtraction */
135 convolvedMaskedImage -= scienceMaskedImage;
136
137 /* Invert */
138 if (invert) {
139 convolvedMaskedImage *= -1.0;
140 }
141
142 t.stop();
143 double time = 1e-9 * t.elapsed().wall;
144 LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
145 "Total compute time to convolve and subtract : %.2f s", time);
146
147 return convolvedMaskedImage;
148}
lsst::afw::image::MaskedImage::getDimensions
lsst::geom::Extent2I getDimensions() const
Definition MaskedImage.h:1058

◆ 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::cpu_timer t;
124
125 afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
126 afwMath::ConvolutionControl convolutionControl = afwMath::ConvolutionControl();
127 convolutionControl.setDoNormalize(false);
128 afwMath::convolve(convolvedMaskedImage, templateImage,
129 convolutionKernel, convolutionControl);
130
131 /* Add in background */
132 *(convolvedMaskedImage.getImage()) += background;
133
134 /* Do actual subtraction */
135 convolvedMaskedImage -= scienceMaskedImage;
136
137 /* Invert */
138 if (invert) {
139 convolvedMaskedImage *= -1.0;
140 }
141
142 t.stop();
143 double time = 1e-9 * t.elapsed().wall;
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/3]

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

◆ imageToEigenMatrix() [2/3]

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

◆ imageToEigenMatrix() [3/3]

template<typename PixelT >
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}
ptr
uint64_t * ptr
Definition RangeSet.cc:95
y
int y
Definition SpanSet.cc:48
lsst::afw::image::ImageBase::getWidth
int getWidth() const
Return the number of columns in the image.
Definition ImageBase.h:294
lsst::afw::image::ImageBase::getHeight
int getHeight() const
Return the number of rows in the image.
Definition ImageBase.h:296
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:401
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:404
lsst::afw::image::Image
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 76 of file BasisLists.cc.

81 {
82 typedef afwMath::Kernel::Pixel Pixel;
83 typedef afwImage::Image<Pixel> Image;
84
85 if (halfWidth < 1) {
86 throw LSST_EXCEPT(pexExcept::Exception, "halfWidth must be positive");
87 }
88 if (nGauss != static_cast<int>(sigGauss.size())) {
89 throw LSST_EXCEPT(pexExcept::Exception, "sigGauss does not have enough entries");
90 }
91 if (nGauss != static_cast<int>(degGauss.size())) {
92 throw LSST_EXCEPT(pexExcept::Exception, "degGauss does not have enough entries");
93 }
94 int fullWidth = 2 * halfWidth + 1;
95 Image image(geom::Extent2I(fullWidth, fullWidth));
96
97 afwMath::KernelList kernelBasisList;
98 for (int i = 0; i < nGauss; i++) {
99 /*
100 sigma = FWHM / ( 2 * sqrt(2 * ln(2)) )
101 */
102 double sig = sigGauss[i];
103 int deg = degGauss[i];
104
105 LOGL_DEBUG("TRACE1.ip.diffim.BasisLists.makeAlardLuptonBasisList",
106 "Gaussian %d : sigma %.2f degree %d", i, sig, deg);
107
108 afwMath::GaussianFunction2<Pixel> gaussian(sig, sig);
109 afwMath::AnalyticKernel kernel(fullWidth, fullWidth, gaussian);
110 afwMath::PolynomialFunction2<Pixel> polynomial(deg);
111
112 for (int j = 0, n = 0; j <= deg; j++) {
113 for (int k = 0; k <= (deg - j); k++, n++) {
114 /* for 0th order term, skip polynomial */
115 (void)kernel.computeImage(image, true);
116 if (n == 0) {
117 std::shared_ptr<afwMath::Kernel>
118 kernelPtr(new afwMath::FixedKernel(image));
119 kernelBasisList.push_back(kernelPtr);
120 continue;
121 }
122
123 /* gaussian to be modified by this term in the polynomial */
124 polynomial.setParameter(n, 1.);
125 (void)kernel.computeImage(image, true);
126 for (int y = 0, v = -halfWidth; y < image.getHeight(); y++, v++) {
127 int u = -halfWidth;
128 for (auto ptr = image.row_begin(y), end = image.row_end(y); ptr != end; ++ptr, ++u) {
129 /* Evaluate from -1 to 1 */
130 *ptr = *ptr * polynomial(u/static_cast<double>(halfWidth),
131 v/static_cast<double>(halfWidth));
132 }
133 }
134 std::shared_ptr<afwMath::Kernel>
135 kernelPtr(new afwMath::FixedKernel(image));
136 kernelBasisList.push_back(kernelPtr);
137 polynomial.setParameter(n, 0.);
138 }
139 }
140 }
141 return renormalizeKernelList(kernelBasisList);
142 }
end
int end
Definition BoundedField.cc:105
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
image
afw::table::Key< afw::table::Array< ImagePixelT > > image
Definition HeavyFootprint.cc:212
lsst::afw::math::AnalyticKernel
A kernel described by a function.
Definition Kernel.h:535
lsst::afw::math::FixedKernel
A kernel created from an Image.
Definition Kernel.h:471
lsst::afw::math::PolynomialFunction2
2-dimensional polynomial function with cross terms
Definition FunctionLibrary.h:447
lsst.pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
lsst::ip::diffim::renormalizeKernelList
lsst::afw::math::KernelList renormalizeKernelList(lsst::afw::math::KernelList const &kernelListIn)
Renormalize a list of basis kernels.
Definition BasisLists.cc:410
std::vector::push_back
T push_back(T... args)
std::vector::size
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 205 of file BasisLists.cc.

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

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

48 {
49 if ((width < 1) || (height < 1)) {
50 throw LSST_EXCEPT(pexExcept::Exception, "nRows and nCols must be positive");
51 }
52 const int signedWidth = static_cast<int>(width);
53 const int signedHeight = static_cast<int>(height);
54 afwMath::KernelList kernelBasisList;
55 for (int row = 0; row < signedHeight; ++row) {
56 for (int col = 0; col < signedWidth; ++col) {
57 std::shared_ptr<afwMath::Kernel>
58 kernelPtr(new afwMath::DeltaFunctionKernel(width, height, geom::Point2I(col,row)));
59 kernelBasisList.push_back(kernelPtr);
60 }
61 }
62 return kernelBasisList;
63 }
lsst::afw::math::DeltaFunctionKernel
A kernel that has only one non-zero pixel (of value 1)
Definition Kernel.h:643
row
int row
Definition CR.cc:145
col
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 503 of file BasisLists.cc.

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

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

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

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

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

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

◆ 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()

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

Definition at line 47 of file _ipdiffimLib.cc.

47 {
48 lsst::utils::python::WrapperCollection wrappers(mod, "lsst.ip.diffim");
49 wrappers.addInheritanceDependency("lsst.meas.base");
50 wrappers.addInheritanceDependency("lsst.afw.math");
51 wrappers.addSignatureDependency("lsst.afw.table");
52 wrappers.addSignatureDependency("lsst.afw.image");
53 wrappers.addSignatureDependency("lsst.afw.geom");
54 wrappers.addSignatureDependency("lsst.daf.base");
55 wrappers.addSignatureDependency("lsst.afw.detection");
56
57 wrapBasisLists(wrappers);
58 wrapDipoleAlgorithms(wrappers);
59 wrapImageStatistics(wrappers);
60 wrapImageSubtract(wrappers);
61 wrapKernelCandidate(wrappers);
62 wrapKernelSolution(wrappers);
63 detail::wrapAssessSpatialKernelVisitor(wrappers);
64 detail::wrapBuildSingleKernelVisitor(wrappers);
65 detail::wrapBuildSpatialKernelVisitor(wrappers);
66 detail::wrapKernelPca(wrappers);
67 detail::wrapKernelSumVisitor(wrappers);
68 wrappers.finish();
69}
lsst::ip::diffim::wrapImageStatistics
void wrapImageStatistics(WrapperCollection &wrappers)
Definition imageStatistics.cc:76
lsst::ip::diffim::wrapBasisLists
void wrapBasisLists(WrapperCollection &wrappers)
Definition basisLists.cc:40
lsst::ip::diffim::wrapKernelCandidate
void wrapKernelCandidate(WrapperCollection &wrappers)
Definition kernelCandidate.cc:120
lsst::ip::diffim::wrapDipoleAlgorithms
void wrapDipoleAlgorithms(WrapperCollection &wrappers)
Definition dipoleAlgorithms.cc:154
lsst::ip::diffim::wrapImageSubtract
void wrapImageSubtract(WrapperCollection &wrappers)
Definition imageSubtract.cc:74
lsst::ip::diffim::wrapKernelSolution
void wrapKernelSolution(WrapperCollection &wrappers)
Definition kernelSolution.cc:186

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

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

◆ wrapBasisLists()

void lsst::ip::diffim::wrapBasisLists ( WrapperCollection & wrappers)

Definition at line 40 of file basisLists.cc.

40 {
41 auto &mod = wrappers.module;
42 mod.def("makeDeltaFunctionBasisList", &makeDeltaFunctionBasisList, "width"_a, "height"_a);
43 mod.def("makeRegularizationMatrix", &makeRegularizationMatrix, "ps"_a);
44 mod.def("makeForwardDifferenceMatrix", &makeForwardDifferenceMatrix, "width"_a, "height"_a, "orders"_a,
45 "borderPenalty"_a, "fitForBackground"_a);
46 mod.def("makeCentralDifferenceMatrix", &makeCentralDifferenceMatrix, "width"_a, "height"_a, "stencil"_a,
47 "borderPenalty"_a, "fitForBackground"_a);
48 mod.def("renormalizeKernelList", &renormalizeKernelList, "kernelListIn"_a);
49 mod.def("makeAlardLuptonBasisList", &makeAlardLuptonBasisList, "halfWidth"_a, "nGauss"_a, "sigGauss"_a,
50 "degGauss"_a);
51}

◆ wrapDipoleAlgorithms()

void lsst::ip::diffim::wrapDipoleAlgorithms ( WrapperCollection & wrappers)

Definition at line 154 of file dipoleAlgorithms.cc.

154 {
155 declareDipoleCentroidControl(wrappers);
156 declareDipoleFluxControl(wrappers);
157 declareDipolePsfFluxControl(wrappers);
158 declareDipoleCentroidAlgorithm(wrappers);
159 declareDipoleFluxAlgorithm(wrappers);
160 declareNaiveDipoleFlux(wrappers);
161 declareNaiveDipoleCentroid(wrappers);
162 declarePsfDipoleFlux(wrappers);
163}

◆ wrapImageStatistics()

void lsst::ip::diffim::wrapImageStatistics ( WrapperCollection & wrappers)
Examples
imageStatistics.cc.

Definition at line 76 of file imageStatistics.cc.

76 {
77 declareImageStatistics<int>(wrappers, "I");
78 declareImageStatistics<float>(wrappers, "F");
79 declareImageStatistics<double>(wrappers, "D");
80}

◆ wrapImageSubtract()

void lsst::ip::diffim::wrapImageSubtract ( WrapperCollection & wrappers)

Definition at line 74 of file imageSubtract.cc.

74 {
75 declareConvolveAndSubtract<float, double>(wrappers);
76 declareConvolveAndSubtract<float, afw::math::Function2<double> const &>(wrappers);
77}

◆ wrapKernelCandidate()

void lsst::ip::diffim::wrapKernelCandidate ( WrapperCollection & wrappers)

Definition at line 120 of file kernelCandidate.cc.

120 {
121 declareKernelCandidate<float>(wrappers, "F");
122}

◆ wrapKernelSolution()

void lsst::ip::diffim::wrapKernelSolution ( WrapperCollection & wrappers)

Definition at line 186 of file kernelSolution.cc.

186 {
187 declareKernelSolution(wrappers);
188 declareStaticKernelSolution<float>(wrappers, "F");
189 declareMaskedKernelSolution<float>(wrappers, "F");
190 declareRegularizedKernelSolution<float>(wrappers, "F");
191 declareSpatialKernelSolution(wrappers);
192}

Variable Documentation

◆ Control

lsst.ip.diffim.Control

Definition at line 46 of file __init__.py.

◆ deprecated

lsst.ip.diffim.deprecated

Definition at line 47 of file __init__.py.

◆ DipoleCentroidControl

lsst.ip.diffim.DipoleCentroidControl

Definition at line 46 of file __init__.py.

◆ DipoleFluxControl

lsst.ip.diffim.DipoleFluxControl

Definition at line 48 of file __init__.py.

◆ executionOrder

lsst.ip.diffim.executionOrder

Definition at line 46 of file __init__.py.

◆ NaiveDipoleCentroid

lsst.ip.diffim.NaiveDipoleCentroid

Definition at line 46 of file __init__.py.

◆ NaiveDipoleFlux

lsst.ip.diffim.NaiveDipoleFlux

Definition at line 48 of file __init__.py.

◆ PsfDipoleFlux

lsst.ip.diffim.PsfDipoleFlux

Definition at line 50 of file __init__.py.

◆ PsfDipoleFluxControl

lsst.ip.diffim.PsfDipoleFluxControl

Definition at line 50 of file __init__.py.

Generated on Thu Apr 25 2024 16:27:27 for LSST Applications by doxygen 1.10.0