Namespace lsst::afw::math

namespace math

Typedefs

typedef std::vector<std::shared_ptr<Kernel>> KernelList
typedef lsst::afw::image::VariancePixel WeightPixel

Enums

enum UndersampleStyle

Values:

THROW_EXCEPTION
REDUCE_INTERP_ORDER
INCREASE_NXNYSAMPLE
enum Property

control what is calculated

Values:

NOTHING = 0x0

We don’t want anything.

ERRORS = 0x1

Include errors of requested quantities.

NPOINT = 0x2

number of sample points

MEAN = 0x4

estimate sample mean

STDEV = 0x8

estimate sample standard deviation

VARIANCE = 0x10

estimate sample variance

MEDIAN = 0x20

estimate sample median

IQRANGE = 0x40

estimate sample inter-quartile range

MEANCLIP = 0x80

estimate sample N-sigma clipped mean (N set in StatisticsControl, default=3)

STDEVCLIP = 0x100

estimate sample N-sigma clipped stdev (N set in StatisticsControl, default=3)

VARIANCECLIP = 0x200

estimate sample N-sigma clipped variance (N set in StatisticsControl, default=3)

MIN = 0x400

estimate sample minimum

MAX = 0x800

estimate sample maximum

SUM = 0x1000

find sum of pixels in the image

MEANSQUARE = 0x2000

find mean value of square of pixel values

ORMASK = 0x4000

get the or-mask of all pixels used.

NCLIPPED = 0x8000

number of clipped points

NMASKED = 0x10000

number of masked points

Functions

template<typename PixelT>
std::shared_ptr<Approximate<PixelT>> makeApproximate(std::vector<double> const &x, std::vector<double> const &y, image::MaskedImage<PixelT> const &im, lsst::geom::Box2I const &bbox, ApproximateControl const &ctrl)

Construct a new Approximate object, inferring the type from the type of the given MaskedImage.

Parameters
  • x: the x-values of points

  • y: the y-values of points

  • im: The values at (x, y)

  • bbox: Range where approximation should be valid

  • ctrl: desired approximation algorithm

UndersampleStyle stringToUndersampleStyle(std::string const &style)

Conversion function to switch a string to an UndersampleStyle

template<typename ImageT>
std::shared_ptr<Background> makeBackground(ImageT const &img, BackgroundControl const &bgCtrl)

A convenience function that uses function overloading to make the correct type of Background

cf. std::make_pair()

std::shared_ptr<BoundedField> operator*(double const scale, std::shared_ptr<BoundedField const> bf)
template<typename OutImageT, typename InImageT>
void scaledPlus(OutImageT &outImage, double c1, InImageT const &inImage1, double c2, InImageT const &inImage2)

Compute the scaled sum of two images

outImage = c1 inImage1 + c2 inImage2

For example to linearly interpolate between two images set: c1 = 1.0 - fracDist c2 = fracDist where fracDist is the fractional distance of outImage from inImage1: location of outImage - location of inImage1 fracDist = —————————————- location of inImage2 - location of inImage1

Parameters
  • [out] outImage: output image

  • [in] c1: coefficient for image 1

  • [in] inImage1: input image 1

  • [in] c2: coefficient for image 2

  • [in] inImage2: input image 2

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if outImage is not same dimensions as inImage1 and inImage2.

template<typename OutImageT, typename InImageT>
OutImageT::SinglePixel convolveAtAPoint(typename InImageT::const_xy_locator inImageLocator, typename lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel>::const_xy_locator kernelLocator, int kWidth, int kHeight)

Apply convolution kernel to an image at one point

Note

This subroutine sacrifices convenience for performance. The user is expected to figure out the kernel center and adjust the supplied pixel accessors accordingly. For an example of how to do this see convolve().

Parameters
  • inImageLocator: locator for image pixel that overlaps pixel (0,0) of kernel (the origin of the kernel, not its center)

  • kernelLocator: locator for (0,0) pixel of kernel (the origin of the kernel, not its center)

  • kWidth: number of columns in kernel

  • kHeight: number of rows in kernel

template<typename OutImageT, typename InImageT>
OutImageT::SinglePixel convolveAtAPoint(typename InImageT::const_xy_locator inImageLocator, std::vector<lsst::afw::math::Kernel::Pixel> const &kernelColList, std::vector<lsst::afw::math::Kernel::Pixel> const &kernelRowList)

Apply separable convolution kernel to an image at one point

Note

This subroutine sacrifices convenience for performance. The user is expected to figure out the kernel center and adjust the supplied pixel accessors accordingly. For an example of how to do this see convolve().

Parameters
  • inImageLocator: locator for image pixel that overlaps pixel (0,0) of kernel (the origin of the kernel, not its center)

  • kernelXList: kernel column vector

  • kernelYList: kernel row vector

template<typename OutImageT, typename InImageT, typename KernelT>
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.

Various convolution kernels are available, including:

If a kernel is spatially varying, its spatial model is computed at each pixel position on the image (pixel position, not pixel index). At present (2009-09-24) this position is computed relative to the lower left corner of the sub-image, but it will almost certainly change to be the lower left corner of the parent image.

All convolution is performed in real space. This allows convolution to handle masked pixels and spatially varying kernels. Although convolution of an Image with a spatially invariant kernel could, in fact, be performed in Fourier space, the code does not do this.

Note that mask bits are smeared by convolution; all nonzero pixels in the kernel smear the mask, even pixels that have very small values. Larger kernels smear the mask more and are also slower to convolve. Use the smallest kernel that will do the job.

convolvedImage has a border of edge pixels which cannot be computed normally. Normally these pixels are set to the standard edge pixel, as returned by edgePixel(). However, if your code cannot handle nans in the image or infs in the variance, you may set doCopyEdge true, in which case the edge pixels are set to the corresponding pixels of the input image and (if there is a mask) the mask EDGE bit is set.

The border of edge pixels has size:

  • kernel.getCtrX() along the left edge

  • kernel.getCtrY() along the bottom edge

  • kernel.getWidth() - 1 - kernel.getCtrX() along the right edge

  • kernel.getHeight() - 1 - kernel.getCtrY() along the top edge You can obtain a bounding box for the good pixels in the convolved image from a bounding box for the entire image using the Kernel method shrinkBBox.

Convolution has been optimized for the various kinds of kernels, as follows (listed approximately in order of decreasing speed):

  • DeltaFunctionKernel convolution is a simple image shift.

  • SeparableKernel convolution is performed by convolving the input by one of the two functions, then the result by the other function. Thus convolution with a kernel of size nCols x nRows becomes convolution with a kernel of size nCols x 1, followed by convolution with a kernel of size 1 x nRows.

  • Convolution with spatially invariant versions of the other kernels is performed by computing the kernel image once and convolving with that. The code has been optimized for cache performance and so should be fairly efficient.

  • Convolution with a spatially varying LinearCombinationKernel is performed by convolving the image by each basis kernel and combining the result by solving the spatial model. This will be efficient provided the kernel does not contain too many or very large basis kernels.

  • Convolution with spatially varying AnalyticKernel is likely to be slow. The code simply computes the output one pixel at a time by computing the AnalyticKernel at that point and applying it to the input image. This is not favorable for cache performance (especially for large kernels) but avoids recomputing the AnalyticKernel. It is probably possible to do better.

Additional convolution functions include:

  • convolveAtAPoint(): convolve a Kernel to an Image or MaskedImage at a point.

  • basicConvolve(): convolve a Kernel with an Image or MaskedImage, but do not set the edge pixels of the output. Optimization of convolution for different types of Kernel are handled by different specializations of basicConvolve().

afw/examples offers programs that time convolution including timeConvolve and timeSpatiallyVaryingConvolve.

Parameters
  • [out] convolvedImage: convolved image; must be the same size as inImage

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if convolvedImage is not the same size as inImage

  • lsst::pex::exceptions::InvalidParameterError: if inImage is smaller than kernel in columns and/or rows.

  • std::bad_alloc: when allocation of memory fails

template<typename OutImageT, typename InImageT, typename KernelT>
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge = false)

Old, deprecated version of convolve.

Parameters
  • [out] convolvedImage: convolved image; must be the same size as inImage

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] doNormalize: if true, normalize the kernel, else use “as is”

  • [in] doCopyEdge: if false (default), set edge pixels to the standard edge pixel; if true, copy edge pixels from input and set EDGE bit of mask

template<typename ImageT>
ImageT::SinglePixel edgePixel(lsst::afw::image::detail::Image_tag)

Return an off-the-edge pixel appropriate for a given Image type

The value is quiet_NaN if that exists for the pixel type, else 0

template<typename MaskedImageT>
MaskedImageT::SinglePixel edgePixel(lsst::afw::image::detail::MaskedImage_tag)

Return an off-the-edge pixel appropriate for a given MaskedImage type

The components are:

  • image = quiet_NaN if that exists for the pixel type, else 0

  • mask = NO_DATA bit set

  • variance = infinity if that exists for the pixel type, else 0

Exceptions
  • lsst::pex::exceptions::LogicError: Thrown if the global mask plane dictionary does not have a NO_DATA bit.

template<typename OutImageT, typename InImageT>
OutImageT::SinglePixel convolveAtAPoint(typename InImageT::const_xy_locator inImageLocator, std::vector<Kernel::Pixel> const &kernelXList, std::vector<Kernel::Pixel> const &kernelYList)

Apply separable convolution kernel to an image at one point

Note

This subroutine sacrifices convenience for performance. The user is expected to figure out the kernel center and adjust the supplied pixel accessors accordingly. For an example of how to do this see convolve().

Parameters
  • inImageLocator: locator for image pixel that overlaps pixel (0,0) of kernel (the origin of the kernel, not its center)

  • kernelXList: kernel column vector

  • kernelYList: kernel row vector

template<class UF>
UF::result_type int1d(UF const &func, IntRegion<typename UF::result_type> &reg, typename UF::result_type const &abserr = DEFABSERR, typename UF::result_type const &relerr = DEFRELERR)

Front end for the 1d integrator

template<class BF, class YREG>
BF::result_type int2d(BF const &func, IntRegion<typename BF::result_type> &reg, YREG const &yreg, typename BF::result_type const &abserr = DEFABSERR, typename BF::result_type const &relerr = DEFRELERR)

Front end for the 2d integrator

template<class TF, class YREG, class ZREG>
TF::result_type int3d(TF const &func, IntRegion<typename TF::result_type> &reg, YREG const &yreg, ZREG const &zreg, typename TF::result_type const &abserr = DEFABSERR, typename TF::result_type const &relerr = DEFRELERR)

Front end for the 3d integrator

template<class BF>
BF::result_type int2d(BF const &func, IntRegion<typename BF::result_type> &reg, IntRegion<typename BF::result_type> &yreg, typename BF::result_type const &abserr = DEFABSERR, typename BF::result_type const &relerr = DEFRELERR)

Front end for the 2d integrator

template<class TF>
TF::result_type int3d(TF const &func, IntRegion<typename TF::result_type> &reg, IntRegion<typename TF::result_type> &yreg, IntRegion<typename TF::result_type> &zreg, typename TF::result_type const &abserr = DEFABSERR, typename TF::result_type const &relerr = DEFRELERR)

Front end for the 3d integrator

template<typename UnaryFunctionT>
UnaryFunctionT::result_type integrate(UnaryFunctionT func, typename UnaryFunctionT::argument_type const a, typename UnaryFunctionT::argument_type const b, double eps = 1.0e-6)

The 1D integrator

Note

This simply wraps the int1d function above and handles the instantiation of the intRegion.

template<typename BinaryFunctionT>
BinaryFunctionT::result_type integrate2d(BinaryFunctionT func, typename BinaryFunctionT::first_argument_type const x1, typename BinaryFunctionT::first_argument_type const x2, typename BinaryFunctionT::second_argument_type const y1, typename BinaryFunctionT::second_argument_type const y2, double eps = 1.0e-6)

The 2D integrator

Note

Adapted from RHL’s SDSS code

std::shared_ptr<Interpolate> makeInterpolate(std::vector<double> const &x, std::vector<double> const &y, Interpolate::Style const style = Interpolate::AKIMA_SPLINE)

A factory function to make Interpolate objects

Parameters
  • x: the x-values of points

  • y: the values at x[]

  • style: desired interpolator

std::shared_ptr<Interpolate> makeInterpolate(ndarray::Array<double const, 1> const &x, ndarray::Array<double const, 1> const &y, Interpolate::Style const style = Interpolate::AKIMA_SPLINE)
Interpolate::Style stringToInterpStyle(std::string const &style)

Conversion function to switch a string to an Interpolate::Style.

Parameters
  • style: desired type of interpolation

Interpolate::Style lookupMaxInterpStyle(int const n)

Get the highest order Interpolation::Style available for ‘n’ points.

Parameters
  • n: Number of points

int lookupMinInterpPoints(Interpolate::Style const style)

Get the minimum number of points needed to use the requested interpolation style

Parameters
  • style: The style in question

void printKernel(lsst::afw::math::Kernel const &kernel, bool doNormalize, double x = 0, double y = 0, std::string pixelFmt = "%7.3f")

Print the pixel values of a Kernel to std::cout

Rows increase upward and columns to the right; thus the lower left pixel is (0,0).

Parameters
  • kernel: the kernel

  • doNormalize: if true, normalize kernel

  • x: x at which to evaluate kernel

  • y: y at which to evaluate kernel

  • pixelFmt: format for pixel values

template<typename ReturnT>
FitResults minimize(lsst::afw::math::Function1<ReturnT> const &function, std::vector<double> const &initialParameterList, std::vector<double> const &stepSizeList, std::vector<double> const &measurementList, std::vector<double> const &varianceList, std::vector<double> const &xPositionList, double errorDef)

Find the minimum of a function(x)

Uses the Minuit fitting package with a standard definition of chiSq (see MinimizerFunctionBase1).

Return

true if minimum is valid, false otherwise

Parameters
  • function: function(x) to be minimized

  • initialParameterList: initial guess for parameters

  • stepSizeList: step size for each parameter; use 0.0 to fix a parameter

  • measurementList: measured values

  • varianceList: variance for each measurement

  • xPositionList: x position of each measurement

  • errorDef: what is this?

To do:

  • Document stepSizeList better

  • Document errorDef

  • Compute stepSize automatically? (if so, find a different way to fix parameters)

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if any input vector is the wrong length

template<typename ReturnT>
FitResults minimize(lsst::afw::math::Function2<ReturnT> const &function, std::vector<double> const &initialParameterList, std::vector<double> const &stepSizeList, std::vector<double> const &measurementList, std::vector<double> const &varianceList, std::vector<double> const &xPositionList, std::vector<double> const &yPositionList, double errorDef)

Find the minimum of a function(x, y)

Uses the Minuit fitting package with a standard definition of chiSq. (see MinimizerFunctionBase2).

To do:

  • Document stepSizeList better

  • Document errorDef

  • Compute stepSize automatically? (if so, find a different way to fix parameters)

Return

true if minimum is valid, false otherwise

Parameters
  • function: function(x,y) to be minimized

  • initialParameterList: initial guess for parameters

  • stepSizeList: step size for each parameter; use 0.0 to fix a parameter

  • measurementList: measured values

  • varianceList: variance for each measurement

  • xPositionList: x position of each measurement

  • yPositionList: y position of each measurement

  • errorDef: what is this?

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if any input vector is the wrong length

template<typename ImageT>
std::shared_ptr<ImageT> offsetImage(ImageT const &image, float dx, float dy, std::string const &algorithmName = "lanczos5", unsigned int buffer = 0)

Return an image offset by (dx, dy) using the specified algorithm

Note

The image pixels are always offset by a fraction of a pixel and the image origin (XY0) picks is modified to handle the integer portion of the offset. In the special case that the offset in both x and y lies in the range (-1, 1) the origin is not changed. Otherwise the pixels are shifted by (-0.5, 0.5] pixels and the origin shifted accordingly.

Parameters
  • image: The image to offset

  • dx: move the image this far in the column direction

  • dy: move the image this far in the row direction

  • algorithmName: Type of resampling Kernel to use

  • buffer: Width of buffer (border) around kernel image to allow for warping edge effects (pixels). Values < 0 are treated as 0. This is only used during computation; the final image has the same dimensions as the kernel.

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if the algorithm is invalid

template<typename ImageT>
std::shared_ptr<ImageT> rotateImageBy90(ImageT const &image, int nQuarter)

Rotate an image by an integral number of quarter turns

Parameters
  • image: The image to rotate

  • nQuarter: the desired number of quarter turns

template<typename ImageT>
std::shared_ptr<ImageT> flipImage(ImageT const &inImage, bool flipLR, bool flipTB)

Flip an image leftright and/or topbottom

Parameters
  • inImage: The image to flip

  • flipLR: Flip left <> right?

  • flipTB: Flip top <> bottom?

template<typename ImageT>
std::shared_ptr<ImageT> binImage(ImageT const &inImage, int const binX, int const binY, lsst::afw::math::Property const flags = lsst::afw::math::MEAN)

Parameters
  • inImage: The image to bin

  • binX: Output pixels are binX*binY input pixels

  • binY: Output pixels are binX*binY input pixels

  • flags: how to generate super-pixels

template<typename ImageT>
std::shared_ptr<ImageT> binImage(ImageT const &inImage, int const binsize, lsst::afw::math::Property const flags = lsst::afw::math::MEAN)

Parameters
  • inImage: The image to bin

  • binsize: Output pixels are binsize*binsize input pixels

  • flags: how to generate super-pixels

template<typename ImageT>
void randomUniformImage(ImageT *image, Random &rand)

Set image to random numbers uniformly distributed in the range [0, 1)

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

template<typename ImageT>
void randomUniformPosImage(ImageT *image, Random &rand)

Set image to random numbers uniformly distributed in the range (0, 1)

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

template<typename ImageT>
void randomUniformIntImage(ImageT *image, Random &rand, unsigned long n)

Set image to random integers uniformly distributed in the range 0 … n - 1

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

  • [in] n: (exclusive) upper limit for random variates

template<typename ImageT>
void randomFlatImage(ImageT *image, Random &rand, double const a, double const b)

Set image to random numbers uniformly distributed in the range [a, b)

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

  • [in] a: (inclusive) lower limit for random variates

  • [in] b: (exclusive) upper limit for random variates

template<typename ImageT>
void randomGaussianImage(ImageT *image, Random &rand)

Set image to random numbers with a gaussian N(0, 1) distribution

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

template<typename ImageT>
void randomChisqImage(ImageT *image, Random &rand, double const nu)

Set image to random numbers with a chi^2_{nu} distribution

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

  • [in] nu: number of degrees of freedom

template<typename ImageT>
void randomPoissonImage(ImageT *image, Random &rand, double const mu)

Set image to random numbers with a Poisson distribution with mean mu (n.b. not per-pixel)

Parameters
  • [out] image: The image to set

  • [inout] rand: definition of random number algorithm, seed, etc.

  • [in] mu: mean of distribution

template<typename PixelT>
std::shared_ptr<lsst::afw::image::Image<PixelT>> statisticsStack(std::vector<std::shared_ptr<lsst::afw::image::Image<PixelT>>> &images, Property flags, StatisticsControl const &sctrl = StatisticsControl(), std::vector<lsst::afw::image::VariancePixel> const &wvector = std::vector<lsst::afw::image::VariancePixel>(0))

Parameters
  • images: Images to process

  • flags: statistics requested

  • sctrl: Control structure

  • wvector: vector containing weights

A function to compute some statistics of a stack of Images

template<typename PixelT>
void statisticsStack(lsst::afw::image::Image<PixelT> &out, std::vector<std::shared_ptr<lsst::afw::image::Image<PixelT>>> &images, Property flags, StatisticsControl const &sctrl = StatisticsControl(), std::vector<lsst::afw::image::VariancePixel> const &wvector = std::vector<lsst::afw::image::VariancePixel>(0))

Parameters
  • out: Output image

  • images: Images to process

  • flags: statistics requested

  • sctrl: Control structure

  • wvector: vector containing weights

@ brief compute statistical stack of Image. Write to output image in-situ

template<typename PixelT>
std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>> statisticsStack(std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>> &images, Property flags, StatisticsControl const &sctrl, std::vector<lsst::afw::image::VariancePixel> const &wvector, image::MaskPixel clipped, std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const &maskMap)

A function to compute some statistics of a stack of Masked Images

If none of the input images are valid for some pixel, the afwMath::StatisticsControl::getNoGoodPixelsMask() bit(s) are set.

Parameters
  • [in] images: MaskedImages to process.

  • [in] flags: Statistics requested.

  • [in] sctrl: Control structure.

  • [in] wvector: Vector of weights.

  • [in] clipped: Mask to set for pixels that were clipped (NOT rejected due to masks).

  • [in] maskMap: Vector of pairs of mask pixel values; any pixel on an input with any of the bits in .first will result in all of the bits in .second being set on the corresponding pixel on the output.

All the work is done in the function computeMaskedImageStack.

template<typename PixelT>
std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>> statisticsStack(std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>> &images, Property flags, StatisticsControl const &sctrl = StatisticsControl(), std::vector<lsst::afw::image::VariancePixel> const &wvector = std::vector<lsst::afw::image::VariancePixel>(0), image::MaskPixel clipped = 0, image::MaskPixel excuse = 0)

Parameters
  • images: MaskedImages to process

  • flags: statistics requested

  • sctrl: control structure

  • wvector: vector containing weights

  • clipped: bitmask to set if any input was clipped or masked

  • excuse: bitmask to excuse from marking as clipped

A function to compute some statistics of a stack of Masked Images

If none of the input images are valid for some pixel, the afwMath::StatisticsControl::getNoGoodPixelsMask() bit(s) are set.

Delegates to the more general version of statisticsStack taking a maskMap.

template<typename PixelT>
void statisticsStack(lsst::afw::image::MaskedImage<PixelT> &out, std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>> &images, Property flags, StatisticsControl const &sctrl, std::vector<lsst::afw::image::VariancePixel> const &wvector, image::MaskPixel clipped, std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const &maskMap)

A function to compute some statistics of a stack of Masked Images

If none of the input images are valid for some pixel, the afwMath::StatisticsControl::getNoGoodPixelsMask() bit(s) are set.

Parameters
  • [out] out: Output MaskedImage.

  • [in] images: MaskedImages to process.

  • [in] flags: Statistics requested.

  • [in] sctrl: Control structure.

  • [in] wvector: Vector of weights.

  • [in] clipped: Mask to set for pixels that were clipped (NOT rejected due to masks).

  • [in] maskMap: Vector of pairs of mask pixel values; any pixel on an input with any of the bits in .first will result in all of the bits in .second being set on the corresponding pixel on the output.

All the work is done in the function computeMaskedImageStack.

template<typename PixelT>
void statisticsStack(lsst::afw::image::MaskedImage<PixelT> &out, std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>> &images, Property flags, StatisticsControl const &sctrl = StatisticsControl(), std::vector<lsst::afw::image::VariancePixel> const &wvector = std::vector<lsst::afw::image::VariancePixel>(0), image::MaskPixel clipped = 0, image::MaskPixel excuse = 0)

Parameters
  • out: Output image

  • images: MaskedImages to process

  • flags: statistics requested

  • sctrl: control structure

  • wvector: vector containing weights

  • clipped: bitmask to set if any input was clipped or masked

  • excuse: bitmask to excuse from marking as clipped

@ brief compute statistical stack of MaskedImage. Write to output image in-situ

template<typename PixelT>
std::vector<PixelT> statisticsStack(std::vector<std::vector<PixelT>> &vectors, Property flags, StatisticsControl const &sctrl = StatisticsControl(), std::vector<lsst::afw::image::VariancePixel> const &wvector = std::vector<lsst::afw::image::VariancePixel>(0))

Parameters
  • vectors: Vectors to process

  • flags: statistics requested

  • sctrl: control structure

  • wvector: vector containing weights

A function to compute some statistics of a stack of std::vectors

template<typename PixelT>
std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>> statisticsStack(lsst::afw::image::Image<PixelT> const &image, Property flags, char dimension, StatisticsControl const &sctrl = StatisticsControl())

A function to compute statistics on the rows or columns of an image

template<typename PixelT>
std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>> statisticsStack(lsst::afw::image::MaskedImage<PixelT> const &image, Property flags, char dimension, StatisticsControl const &sctrl = StatisticsControl())

A function to compute statistics on the rows or columns of an image

Property stringToStatisticsProperty(std::string const property)

Conversion function to switch a string to a Property (see Statistics.h)

template<typename Pixel>
Statistics makeStatistics(lsst::afw::image::Image<Pixel> const &img, lsst::afw::image::Mask<image::MaskPixel> const &msk, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Handle a watered-down front-end to the constructor (no variance)

template<typename ImageT, typename MaskT, typename VarianceT>
Statistics makeStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Handle a straight front-end to the constructor

template<typename Pixel>
Statistics makeStatistics(lsst::afw::image::MaskedImage<Pixel> const &mimg, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Handle MaskedImages, just pass the getImage() and getMask() values right on through.

template<typename Pixel>
Statistics makeStatistics(lsst::afw::image::MaskedImage<Pixel> const &mimg, lsst::afw::image::Image<WeightPixel> const &weights, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Handle MaskedImages, just pass the getImage() and getMask() values right on through.

Statistics makeStatistics(lsst::afw::image::Mask<lsst::afw::image::MaskPixel> const &msk, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Specialization to handle Masks

Parameters
  • msk: Image (or MaskedImage) whose properties we want

  • flags: Describe what we want to calculate

  • sctrl: Control how things are calculated

template<typename Pixel>
Statistics makeStatistics(lsst::afw::image::Image<Pixel> const &img, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Parameters
  • img: Image (or Image) whose properties we want

  • flags: Describe what we want to calculate

  • sctrl: Control calculation

The makeStatistics() overload to handle regular (non-masked) Images

template<typename EntryT>
Statistics makeStatistics(std::vector<EntryT> const &v, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Parameters
  • v: Image (or MaskedImage) whose properties we want

  • flags: Describe what we want to calculate

  • sctrl: Control calculation

The makeStatistics() overload to handle std::vector<>

template<typename EntryT>
Statistics makeStatistics(std::vector<EntryT> const &v, std::vector<WeightPixel> const &vweights, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Parameters
  • v: Image (or MaskedImage) whose properties we want

  • vweights: Weights

  • flags: Describe what we want to calculate

  • sctrl: Control calculation

The makeStatistics() overload to handle std::vector<>

template<typename EntryT>
Statistics makeStatistics(lsst::afw::math::MaskedVector<EntryT> const &mv, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Parameters
  • mv: MaskedVector

  • flags: Describe what we want to calculate

  • sctrl: Control calculation

The makeStatistics() overload to handle lsst::afw::math::MaskedVector<>

template<typename EntryT>
Statistics makeStatistics(lsst::afw::math::MaskedVector<EntryT> const &mv, std::vector<WeightPixel> const &vweights, int const flags, StatisticsControl const &sctrl = StatisticsControl())

Parameters
  • mv: MaskedVector

  • vweights: weights

  • flags: Describe what we want to calculate

  • sctrl: Control calculation

The makeStatistics() overload to handle lsst::afw::math::MaskedVector<>

std::shared_ptr<SeparableKernel> makeWarpingKernel(std::string name)

Return a warping kernel given its name.

Intended for use with warpImage() and warpExposure().

Allowed names are:

A warping kernel is a subclass of SeparableKernel with the following properties (though for the sake of speed few, if any, of these are enforced):

  • Width and height are even. This is unusual for a kernel, but it is more efficient because if the extra pixel was included it would always have value 0.

  • The center pixels should be adjacent to the kernel center. Again, this avoids extra pixels that are sure to have value 0.

  • It has two parameters: fractional x and fractional row position on the source image. The fractional position is the offset of the pixel position on the source from the center of a nearby source pixel:

    • The pixel whose center is just below or to the left of the source position: 0 <= fractional x and y < 0 and the kernel center is the default (size-1)/2.

    • The pixel whose center is just above or to the right of the source position: -1.0 < fractional x and y <= 0 and the kernel center must be set to (size+1)/2.

template<typename DestExposureT, typename SrcExposureT>
int warpExposure(DestExposureT &destExposure, SrcExposureT const &srcExposure, WarpingControl const &control, typename DestExposureT::MaskedImageT::SinglePixel padValue = lsst::afw::math::edgePixel<typename DestExposureT::MaskedImageT>(typename lsst::afw::image::detail::image_traits<typename DestExposureT::MaskedImageT>::image_category()))

Parameters
  • destExposure: Remapped exposure. Wcs and xy0 are read, MaskedImage is set, and PhotoCalib, Filter and VisitInfo are copied from srcExposure. All other attributes are left alone (including Detector and Psf)

  • srcExposure: Source exposure

  • control: control parameters

  • padValue: use this value for undefined (edge) pixels

Warp (remap) one exposure to another.

This is a convenience wrapper around warpImage().

template<typename DestImageT, typename SrcImageT>
int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage, geom::SkyWcs const &srcWcs, WarpingControl const &control, typename DestImageT::SinglePixel padValue = lsst::afw::math::edgePixel<DestImageT>(typename lsst::afw::image::detail::image_traits<DestImageT>::image_category()))

Warp an Image or MaskedImage to a new Wcs. See also convenience function warpExposure() to warp an Exposure.

Parameters
  • destImage: remapped image

  • destWcs: WCS of remapped image

  • srcImage: source image

  • srcWcs: WCS of source image

  • control: control parameters

  • padValue: use this value for undefined (edge) pixels

Edge pixels are set to padValue; these are pixels that cannot be computed because they are too near the edge of srcImage or miss srcImage entirely.

Algorithm Without Interpolation:

Return

the number of valid pixels in destImage (those that are not edge pixels).

For each integer pixel position in the remapped Exposure:

  • The associated pixel position on srcImage is determined using the destination and source WCS

  • The warping kernel’s parameters are set based on the fractional part of the pixel position on srcImage

  • The warping kernel is applied to srcImage at the integer portion of the pixel position to compute the remapped pixel value

  • A flux-conservation factor is determined from the source and destination WCS and is applied to the remapped pixel

The scaling of intensity for relative area of source and destination uses two minor approximations:

  • The area of the sky marked out by a pixel on the destination image corresponds to a parallellogram on the source image.

  • The area varies slowly enough across the image that we can get away with computing the source area shifted by half a pixel up and to the left of the true area.

Algorithm With Interpolation:

Interpolation simply reduces the number of times WCS is used to map between destination and source pixel position. This computation is only made at a grid of points on the destination image, separated by interpLen pixels along rows and columns. All other source pixel positions are determined by linear interpolation between those grid points. Everything else remains the same.

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if destImage overlaps srcImage

  • std::bad_alloc: when allocation of CPU memory fails

Warning

The code that tests for image overlap is not guranteed to work correctly, based on the C++ standard. It is, in theory, possible for the code to report a “false positive”, meaning that it may claim that images overlap when they do not. We don’t believe that any of our current compilers have this problem. If, in the future, this becomes a problem then we will probably have to remove the test and rely on users being careful.

template<typename DestImageT, typename SrcImageT>
int warpImage(DestImageT &destImage, SrcImageT const &srcImage, geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, typename DestImageT::SinglePixel padValue = lsst::afw::math::edgePixel<DestImageT>(typename lsst::afw::image::detail::image_traits<DestImageT>::image_category()))

A variant of warpImage that uses a Transform<Point2Endpoint, Point2Endpoint> instead of a pair of WCS to describe the transformation.

Return

the number of good pixels

Parameters
  • [inout] destImage: Destination image; all pixels are set

  • [in] srcImage: Source image

  • [in] srcToDest: Transformation from source to destination pixels, in parent coordinates; the inverse must be defined (and is the only direction used). makeWcsPairTransform(srcWcs, destWcs) is one way to compute this transform.

  • [in] control: Warning control parameters

  • [in] padValue: Value used for pixels in the destination image that are outside the region of pixels that can be computed from the source image

template<typename DestImageT, typename SrcImageT>
int warpCenteredImage(DestImageT &destImage, SrcImageT const &srcImage, lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, WarpingControl const &control, typename DestImageT::SinglePixel padValue = lsst::afw::math::edgePixel<DestImageT>(typename lsst::afw::image::detail::image_traits<DestImageT>::image_category()))

Parameters
  • destImage: remapped image

  • srcImage: source image

  • linearTransform: linear transformation to apply

  • centerPosition: pixel position for location of linearTransform

  • control: control parameters

  • padValue: use this value for undefined (edge) pixels

Warp an image with a LinearTranform about a specified point.

This enables warping an image of e.g. a PSF without translating the centroid.

Variables

template<typename T>
bool constexpr IS_NOTHROW_INIT = noexcept(static_cast<T>(1.0))

Test that a Function’s return value is nothrow-castable to T

std::complex is an example of a numeric type that does not satisfy this requirement.

double const MOCK_INF = 1.e10
double const DEFABSERR = 1.e-15
double const DEFRELERR = 1.e-6
generic_kernel_tag generic_kernel_tag_
deltafunction_kernel_tag deltafunction_kernel_tag_
class AnalyticKernel : public lsst::afw::table::io::PersistableFacade<AnalyticKernel>, public lsst::afw::math::Kernel
#include <Kernel.h>

A kernel described by a function.

The function’s x, y arguments are as follows:

  • -getCtr() for the lower left corner pixel

  • 0, 0 for the center pixel

  • (getDimensions() - 1) - getCtr() for the upper right pixel

Note: each pixel is set to the value of the kernel function at the center of the pixel (rather than averaging the function over the area of the pixel).

template<typename PixelT>
class Approximate
#include <Approximate.h>

Approximate values for a MaskedImage

class ApproximateControl
#include <Approximate.h>

Control how to make an approximation

Note

the x- and y-order must be the same, due to a limitation of Chebyshev1Function2

class Background
#include <Background.h>

A virtual base class to evaluate image background levels

Subclassed by lsst::afw::math::BackgroundMI

class BackgroundControl
#include <Background.h>

Pass parameters to a Background object

class BackgroundMI : public lsst::afw::math::Background
#include <Background.h>

A class to evaluate image background levels

Break an image up into nx*ny sub-images and use a statistical to estimate the background levels in each square. Then use a user-specified or algorithm to estimate background at a given pixel coordinate.

Methods are available to return the background at a point (inefficiently), or an entire background image. BackgroundControl contains a public StatisticsControl member to allow user control of how the backgrounds are computed.

math::BackgroundControl bctrl(7, 7);  // number of sub-image squares in {x,y}-dimensions
bctrl.sctrl.setNumSigmaClip(5.0);     // use 5-sigma clipping for the sub-image means
std::shared_ptr<math::Background> backobj = math::makeBackground(img, bctrl);
// get a whole background image
Image<PixelT> back = backobj->getImage<PixelT>(math::Interpolate::NATURAL_SPLINE);

// get the background at a pixel at i_x,i_y
double someValue = backobj.getPixel(math::Interpolate::LINEAR, i_x, i_y);

template<typename ReturnT>
class BasePolynomialFunction2 : public lsst::afw::math::Function2<ReturnT>
#include <Function.h>

Base class for 2-dimensional polynomials of the form:

f(x,y) =   c0 f0(x) f0(y)                                        (0th order)
         + c1 f1(x) f0(x) + c2 f0(x) f1(y)                       (1st order)
         + c3 f2(x) f0(y) + c4 f1(x) f1(y) + c5 f0(x) f2(y)      (2nd order)
         + ...

and typically f0(x) = 1

Subclassed by lsst::afw::math::Chebyshev1Function2< ReturnT >, lsst::afw::math::PolynomialFunction2< ReturnT >

class BilinearWarpingKernel : public lsst::afw::math::SeparableKernel
#include <warpExposure.h>

Bilinear warping: fast; good for undersampled data.

The kernel size is 2 x 2.

For more information about warping kernels see makeWarpingKernel

class BoundedField : public lsst::afw::table::io::PersistableFacade<BoundedField>, public lsst::afw::table::io::Persistable
#include <BoundedField.h>

An abstract base class for 2-d functions defined on an integer bounding boxes

Integer bounding boxes (lsst::geom::Box2I) are inclusive of the end pixels (integer positions correspond to the centers of the pixels and include the entirety of those pixels). Thus a BoundedField defined on the box [x0, x1] x [y0, y1] actually covers the range [x0 - 0.5, x1 + 0.5] x [y0 - 0.5, y1 + 0.5].

BoundedField provides a number of ways of accessing the function, all delegating to a single evaluate-at-a-point implementation. The base class does not mandate anything about how the field is constructed, so it’s appropriate for use with e.g. model-fitting results, interpolation results points, or functions known a priori.

Usually, BoundedField will be used to represent functions that correspond to images, for quantities such as aperture corrections, photometric scaling, PSF model parameters, or backgrounds, and its bounding box will be set to match the PARENT bounding box of the image.

Subclassed by lsst::afw::math::ChebyshevBoundedField, lsst::afw::math::PixelAreaBoundedField, lsst::afw::math::ProductBoundedField, lsst::afw::math::TransformBoundedField, lsst::meas::algorithms::CoaddBoundedField

class CandidateVisitor

Subclassed by lsst::ip::diffim::detail::AssessSpatialKernelVisitor< PixelT >, lsst::ip::diffim::detail::BuildSingleKernelVisitor< PixelT >, lsst::ip::diffim::detail::BuildSpatialKernelVisitor< PixelT >, lsst::ip::diffim::detail::KernelPcaVisitor< PixelT >, lsst::ip::diffim::detail::KernelSumVisitor< PixelT >

template<typename ReturnT>
class Chebyshev1Function1 : public lsst::afw::math::Function1<ReturnT>
#include <FunctionLibrary.h>

1-dimensional weighted sum of Chebyshev polynomials of the first kind.

f(x) = c0 T0(x’) + c1 T1(x’) + c2 T2(x’) + … = c0 + c1 T1(x’) + c2 T2(x’) + … where:

  • Tn(x) is the nth Chebyshev function of the first kind: T0(x) = 1 T1(x) = 2 Tn+1(x) = 2xTn(x) + Tn-1(x)

  • x’ is x offset and scaled to range [-1, 1] as x ranges over [minX, maxX]

The function argument must be in the range [minX, maxX].

template<typename ReturnT>
class Chebyshev1Function2 : public lsst::afw::math::BasePolynomialFunction2<ReturnT>
#include <FunctionLibrary.h>

2-dimensional weighted sum of Chebyshev polynomials of the first kind.

f(x,y) = c0 T0(x’) T0(y’) # order 0

  • c1 T1(x’) T0(y’) + c2 T0(x’) T1(y’) # order 1

  • c3 T2(x’) T0(y’) + c4 T1(x’) T1(y’) + c5 T0(x’) T2(y’) # order 2

= c0 # order 0

  • c1 T1(x’) + c2 T1(y’) # order 1

  • c3 T2(x’) + c4 T1(x’) T1(y’) + c5 T2(y’) # order 2

where:

  • Tn(x) is the nth Chebyshev function of the first kind: T0(x) = 1 T1(x) = x Tn+1(x) = 2xTn(x) + Tn-1(x)

  • x’ is x offset and scaled to range [-1, 1] as x ranges over [minX, maxX]

  • y’ is y offset and scaled to range [-1, 1] as y ranges over [minY, maxY]

Return value is incorrect if function arguments are not in the range [minX, maxX], [minY, maxY].

class ChebyshevBoundedField : public lsst::afw::table::io::PersistableFacade<ChebyshevBoundedField>, public lsst::afw::math::BoundedField
#include <ChebyshevBoundedField.h>

A BoundedField based on 2-d Chebyshev polynomials of the first kind.

The 2-d Chebyshev polynomial used here is defined as:

\[ f(x,y) = \sum_i \sum_j a_{i,j} T_i(x) T_j(y) \]

where \(T_n(x)\) is the n-th order Chebyshev polynomial of \(x\) and \(a_{i,j}\) is the corresponding coefficient of the (i,j) polynomial term.

ChebyshevBoundedField supports fitting to gridded and non-gridded data, as well coefficient matrices with different x- and y-order.

There is currently quite a bit of duplication of functionality between ChebyshevBoundedField, ApproximateChebyshev, and Chebyshev1Function2; the intent is that ChebyshevBoundedField will ultimately replace ApproximateChebyshev and should be preferred over Chebyshev1Function2 when the parametrization interface that is part of the Function2 class is not needed.

class ChebyshevBoundedFieldControl
#include <ChebyshevBoundedField.h>

A control object used when fitting ChebyshevBoundedField to data (see ChebyshevBoundedField::fit)

class ConvolutionControl
#include <ConvolveImage.h>

Parameters to control convolution

template<typename T>
class Covariogram
#include <GaussianProcess.h>

The parent class of covariogram functions for use in Gaussian Process interpolation

Each instantiation of a Covariogram will store its own hyper parameters

Subclassed by lsst::afw::math::NeuralNetCovariogram< T >, lsst::afw::math::SquaredExpCovariogram< T >

struct deltafunction_kernel_tag : public lsst::afw::math::generic_kernel_tag
#include <traits.h>

Kernel has only one non-zero pixel.

class DeltaFunctionKernel : public lsst::afw::table::io::PersistableFacade<DeltaFunctionKernel>, public lsst::afw::math::Kernel
#include <Kernel.h>

A kernel that has only one non-zero pixel (of value 1)

It has no adjustable parameters and so cannot be spatially varying.

template<typename ReturnT>
class DoubleGaussianFunction2 : public lsst::afw::math::Function2<ReturnT>
#include <FunctionLibrary.h>

double Guassian (sum of two Gaussians)

Intended for use as a PSF model: the main Gaussian represents the core and the second Gaussian represents the wings.

f(x,y) = A (e^(-r^2 / 2 sigma1^2) + ampl2 e^(-r^2 / 2 sigma2^2)) where:

  • A = 1 / (2 pi (sigma1^2 + ampl2 sigma2^2))

  • r^2 = x^2 + y^2 coefficients c[0] = sigma1, c[1] = sigma2, c[2] = ampl2

struct FitResults
#include <minimize.h>

Results from minimizing a function

class FixedKernel : public lsst::afw::table::io::PersistableFacade<FixedKernel>, public lsst::afw::math::Kernel
#include <Kernel.h>

A kernel created from an Image

It has no adjustable parameters and so cannot be spatially varying.

template<typename ReturnT>
class Function : public lsst::afw::table::io::PersistableFacade<Function<ReturnT>>, public lsst::afw::table::io::Persistable
#include <Function.h>

Basic Function class.

Function objects are functions whose parameters may be read and changed using getParameters and setParameters. They were designed for use with the Kernel class.

These are simple functors with the restrictions that:

  • Function arguments and parameters are double precision

  • The return type is templated

To create a function for a particular equation, subclass Function or (much more likely) Function1 or Function2. Your subclass must:

  • Have one or more constructors, all of which must initialize _params

  • Define operator() with code to compute the function using this->_params or this->getParams() to reference the parameters

  • If the function is a linear combination of parameters then override the function isLinearCombination.

If you wish to cache any information you may use the _isCacheValid flag; this is automatically set false whenever parameters are changed.

Design Notes: The reason these functions exist (rather than using a pre-existing function class, such as Functor in VisualWorkbench) is because the Kernel class requires function objects with a standard interface for setting and getting function parameters.

The reason isLinearCombination exists is to support refactoring LinearCombinationKernels.

Subclassed by lsst::afw::math::Function1< ReturnT >, lsst::afw::math::Function2< ReturnT >

template<typename ReturnT>
class Function1 : public lsst::afw::table::io::PersistableFacade<Function1<ReturnT>>, public lsst::afw::math::Function<ReturnT>
#include <Function.h>

A Function taking one argument.

Subclass and override operator() to do useful work.

Subclassed by lsst::afw::math::Chebyshev1Function1< ReturnT >, lsst::afw::math::GaussianFunction1< ReturnT >, lsst::afw::math::IntegerDeltaFunction1< ReturnT >, lsst::afw::math::LanczosFunction1< ReturnT >, lsst::afw::math::NullFunction1< ReturnT >, lsst::afw::math::PolynomialFunction1< ReturnT >

template<typename ReturnT>
class Function2 : public lsst::afw::table::io::PersistableFacade<Function2<ReturnT>>, public lsst::afw::math::Function<ReturnT>
#include <Function.h>

A Function taking two arguments.

Subclass and override operator() to do useful work.

Subclassed by lsst::afw::math::BasePolynomialFunction2< ReturnT >, lsst::afw::math::DoubleGaussianFunction2< ReturnT >, lsst::afw::math::GaussianFunction2< ReturnT >, lsst::afw::math::IntegerDeltaFunction2< ReturnT >, lsst::afw::math::LanczosFunction2< ReturnT >, lsst::afw::math::NullFunction2< ReturnT >

template<typename ReturnT>
class GaussianFunction1 : public lsst::afw::math::Function1<ReturnT>
#include <FunctionLibrary.h>

1-dimensional Gaussian

f(x) = A e^(-x^2 / 2 sigma^2) where:

  • A = 1 / (sqrt(2 pi) xSigma) coefficient c0 = sigma

template<typename ReturnT>
class GaussianFunction2 : public lsst::afw::math::Function2<ReturnT>
#include <FunctionLibrary.h>

2-dimensional Gaussian

f(x,y) = A e^((-pos1^2 / 2 sigma1^2) - (pos2^2 / 2 sigma2^2)) where:

  • A = 1 / (2 pi sigma1 sigma2)

  • pos1 = cos(angle) x + sin(angle) y

  • pos2 = -sin(angle) x + cos(angle) y coefficients c0 = sigma1, c1 = sigma2, c2 = angle

Note

if sigma1 > sigma2 then angle is the angle of the major axis

template<typename T>
class GaussianProcess
#include <GaussianProcess.h>

Stores values of a function sampled on an image and allows you to interpolate the function to unsampled points

The data will be stored in a KD Tree for easy nearest neighbor searching when interpolating.

The array _function[] will contain the values of the function being interpolated. You can provide a two dimensional array _function[][] if you wish to interpolate a vector of functions. In this case _function[i][j] is the jth function associated with the ith data point. Note: presently, the covariance matrices do not relate elements of _function[i][] to each other, so the variances returned will be identical for all functions evaluated at the same point in parameter space.

_data[i][j] will be the jth component of the ith data point.

_max and _min contain the maximum and minimum values of each dimension in parameter space (if applicable) so that data points can be normalized by _max-_min to keep distances between points reasonable. This is an option specified by calling the relevant constructor.

class GaussianProcessTimer
#include <GaussianProcess.h>

This is a structure for keeping track of how long the interpolation methods spend on different parts of the interpolation

_eigenTime keeps track of how much time is spent using Eigen’s linear algebra packages

_iterationTime keeps track of how much time is spent iterating over matrix indices (this is also a catch-all for time that does not obviously fit in the other categories)

_searchTime keeps track of how much time is spent on nearest neighbor searches (when applicable)

_varianceTime keeps track of how much time is spent calculating the variance of our interpolated function value (note: time spent using Eigen packages for this purpose is tallied here, not in _eigenTime)

_totalTime keeps track of how much time total is spent on interpolations

_interpolationCount keeps track of how many points have been interpolated

struct generic_kernel_tag
#include <traits.h>

Tags carrying information about Kernels Kernel with no special properties.

Subclassed by lsst::afw::math::deltafunction_kernel_tag

template<typename ValueT>
class ImageImposter
#include <Statistics.h>

A vector wrapper to provide a vector with the necessary methods and typedefs to be processed by Statistics as though it were an Image.

template<typename ValueT>
class infinite_iterator : public boost::iterator_adaptor<infinite_iterator<ValueT>, const ValueT *, const ValueT, boost::forward_traversal_tag>
#include <Statistics.h>

This iterator will never increment. It is returned by row_begin() in the MaskImposter class (below) to allow phony mask pixels to be iterated over for non-mask images within Statistics.

template<typename ReturnT>
class IntegerDeltaFunction1 : public lsst::afw::math::Function1<ReturnT>
#include <FunctionLibrary.h>

1-dimensional integer delta function.

f(x) = 1 if x == xo, 0 otherwise.

For use as a kernel function be sure to handle the offset for row and column center; see examples/deltaFunctionKernel for an example.

template<typename ReturnT>
class IntegerDeltaFunction2 : public lsst::afw::math::Function2<ReturnT>
#include <FunctionLibrary.h>

2-dimensional integer delta function.

f(x) = 1 if x == xo and y == yo, 0 otherwise.

For use as a kernel function be sure to handle the offset for row and column center; see examples/deltaFunctionKernel for an example.

template<typename KernelT>
struct is_analyticKernel
#include <traits.h>

traits class to determine if a Kernel is represented as an analytic function

Subclassed by lsst::afw::math::is_analyticKernel< KernelT * >, lsst::afw::math::is_analyticKernel< std::shared_ptr< KernelT > >

template<typename T>
class KdTree
#include <GaussianProcess.h>

The data for GaussianProcess is stored in a KD tree to facilitate nearest-neighbor searches

Note: I have removed the ability to arbitrarily specify a distance function. The KD Tree nearest neighbor search algorithm only makes sense in the case of Euclidean distances, so I have forced KdTree to use Euclidean distances.

class Kernel : public lsst::afw::table::io::PersistableFacade<Kernel>, public lsst::afw::table::io::Persistable
#include <Kernel.h>

Kernels are used for convolution with MaskedImages and (eventually) Images

Kernel is a virtual base class; it cannot be instantiated. The following notes apply to Kernel and to its subclasses.

The template type should usually be float or double; integer kernels should be used with caution because they do not normalize well.

The center pixel of a Kernel is at index: (width-1)/2, (height-1)/2. Thus it is centered along columns/rows if the kernel has an odd number of columns/rows and shifted 1/2 pixel towards 0 otherwise. A kernel should have an odd number of columns and rows unless it is intended to shift an image.

Spatially Varying Kernels

Kernels may optionally vary spatially (so long as they have any kernel parameters). To make a spatially varying kernel, specify a spatial function at construction (you cannot change your mind once the kernel is constructed). You must also specify a set of spatial parameters, and you may do this at construction and/or later by calling setSpatialParameters. The spatial parameters are a vector (one per kernel function parameter) of spatial function parameters. In other words the spatial parameters are a vector of vectors indexed as [kernel parameter][spatial parameter]. The one spatial function is used to compute the kernel parameters at a given spatial position by computing each kernel parameter using its associated vector of spatial function parameters.

The convolve function computes the spatial function at the pixel position (not index) of the image. See the convolve function for details.

Note that if a kernel is spatially varying then you may not set the kernel parameters directly; that is the job of the spatial function! However, you may change the spatial parameters at any time.

Design Notes

The basic design is to use the same kernel class for both spatially varying and spatially invariant kernels. The user either does or does not supply a function describing the spatial variation at creation time. In addition, analytic kernels are described by a user-supplied function of the same basic type as the spatial variation function.

Several other designs were considered, including: A) Use different classes for spatially varying and spatially invariant versions of each kernel. Thus instead of three basic kernel classes (FixedKernel, AnalyticKernel and LinearCombinationKernel) we would have five (since FixedKernel cannot be spatially varying). Robert Lupton argued that was a needless expansion of the class hiearchy and I agreed. B) Construct analytic kernels by defining a subclass of AnalyticKernel that is specific to the desired functional (e.g. GaussianAnalyticKernel). If spatial models are handled the same way then this creates a serious proliferation of kernel classes (even if we only have two different spatial models, e.g. polynomial and Chebyshev polynomial). I felt it made more sense to define the spatial model by some kind of function class (often called a “functor”), and since we needed such a class, I chose to use it for the analytic kernel as well.

However, having a separate function class does introduce some potential inefficiencies. If a function is part of the class it can potentially be evaluated more quickly than calling a function for each pixel or spatial position.

A possible variant on the current design is to define the spatial model and analytic kernel by specifying the functions as template parameters. This has the potential to regain some efficiency in evaluating the functions. However, it would be difficult or impossible to pre-instantiate the desired template classes, a requirement of the LSST coding standards.

Subclassed by lsst::afw::math::AnalyticKernel, lsst::afw::math::DeltaFunctionKernel, lsst::afw::math::FixedKernel, lsst::afw::math::LinearCombinationKernel, lsst::afw::math::SeparableKernel

template<typename KernelT>
struct kernel_traits
#include <traits.h>

template trait class with information about Kernels

template<typename ReturnT>
class LanczosFunction1 : public lsst::afw::math::Function1<ReturnT>
#include <FunctionLibrary.h>

1-dimensional Lanczos function

f(x) = sinc(pi x’) sinc(pi x’ / n) where x’ = x - xOffset and coefficient c0 = xOffset

Warning: the Lanczos function is sometimes forced to 0 if |x’| > n but this implementation does not perform that truncation so as to improve Lanczos kernels.

template<typename ReturnT>
class LanczosFunction2 : public lsst::afw::math::Function2<ReturnT>
#include <FunctionLibrary.h>

2-dimensional separable Lanczos function

f(x, y) = sinc(pi x’) sinc(pi x’ / n) sinc(pi y’) sinc(pi y’ / n) where x’ = x - xOffset and y’ = y - yOffset and coefficients c0 = xOffset, c1 = yOffset

Warning

the Lanczos function is sometimes forced to 0 if |x’| > n or |y’| > n but this implementation does not perform that truncation so as to improve Lanczos kernels.

class LanczosWarpingKernel : public lsst::afw::math::SeparableKernel
#include <warpExposure.h>

Lanczos warping: accurate but slow and can introduce ringing artifacts.

This kernel is the product of two 1-dimensional Lanczos functions. The number of minima and maxima in the 1-dimensional Lanczos function is 2*order + 1. The kernel has one pixel per function minimum or maximum; but as applied to warping, the first or last pixel is always zero and can be omitted. Thus the kernel size is 2*order x 2*order.

For more information about warping kernels see makeWarpingKernel

class LeastSquares
#include <LeastSquares.h>

Solver for linear least-squares problems.

Linear least-squares problems are defined as finding the vector \(x\) that minimizes \(\left|A x - b\right|_2\), with the number of rows of \(A\) generally greater than the number of columns. We call \(A\) the design matrix, \(b\) the data vector, and \(x\) the solution vector. When the rank of \(A\) is equal to the number of columns, we can obtain using the solution using the normal equations:

\[ A^T A x = A^T b \]
If \(A\) is not full-rank, the problem is underconstrained, and we usually wish to solve the minimum-norm least-squares problem, which also minimizes \(|x|_2\). This can be done by computing a pseudo-inverse of \(A\) using an SVD or complete orthogonal factorization of \(A\), or by performing an Eigen decomposition of \(A^T A\).

This class can be constructed from the design matrix and data vector, or from the two terms in the normal equations (below, we call the matrix \(A^TA\) the Fisher matrix, and the vector \(A^T b\) simply the “right-hand side” (RHS) vector). If initialized with the design matrix and data vector, we can still use the normal equations to solve it. The solution via the normal equations is more susceptible to round-off error, but it is also faster, and if the normal equation terms can be computed directly it can be significantly less expensive in terms of memory. The Fisher matrix is a symmetric matrix, and it should be exactly symmetric when provided as input, because which triangle will be used is an implementation detail that is subject to change and may depend on the storage order.

The solver always operates in double precision, and returns all results in double precision. However, it can be initialized from single precision inputs. It isn’t usually a good idea to construct the normal equations in single precision, however, even when the data are single precision.

class LinearCombinationKernel : public lsst::afw::table::io::PersistableFacade<LinearCombinationKernel>, public lsst::afw::math::Kernel
#include <Kernel.h>

A kernel that is a linear combination of fixed basis kernels.

Convolution may be performed by first convolving the image with each fixed kernel, then adding the resulting images using the (possibly spatially varying) kernel coefficients.

The basis kernels are cloned (deep copied) so you may safely modify your own copies.

Warnings:

  • This class does not normalize the individual basis kernels; they are used “as is”.

template<typename ValueT>
class MaskImposter
#include <Statistics.h>

A Mask wrapper to provide an infinite_iterator for Mask::row_begin(). This allows a fake Mask to be passed in to Statistics with a regular (non-masked) Image.

class NearestWarpingKernel : public lsst::afw::math::SeparableKernel
#include <warpExposure.h>

Nearest neighbor warping: fast; good for undersampled data.

The kernel size is 2 x 2.

For more information about warping kernels see makeWarpingKernel

template<typename T>
class NeuralNetCovariogram : public lsst::afw::math::Covariogram<T>
#include <GaussianProcess.h>

a Covariogram that recreates a neural network with one hidden layer and infinite units in that layer

Contains two hyper parameters (_sigma0 and _sigma1) that characterize the expected variance of the function being interpolated

see Rasmussen and Williams (2006) http://www.gaussianprocess.org/gpml/ equation 4.29

template<typename ReturnT>
class NullFunction1 : public lsst::afw::math::Function1<ReturnT>
#include <Function.h>

a class used in function calls to indicate that no Function1 is being provided

template<typename ReturnT>
class NullFunction2 : public lsst::afw::math::Function2<ReturnT>
#include <Function.h>

a class used in function calls to indicate that no Function2 is being provided

class PixelAreaBoundedField : public lsst::afw::math::BoundedField
#include <PixelAreaBoundedField.h>

A BoundedField that evaluate the pixel area of a SkyWcs in angular units.

Typically used to move an image or source flux between surface brightness and fluence.

template<typename ReturnT>
class PolynomialFunction1 : public lsst::afw::math::Function1<ReturnT>
#include <FunctionLibrary.h>

1-dimensional polynomial function.

f(x) = c0 + c1 x + c2 x^2 + … cn-1 x^(n-1)

template<typename ReturnT>
class PolynomialFunction2 : public lsst::afw::math::BasePolynomialFunction2<ReturnT>
#include <FunctionLibrary.h>

2-dimensional polynomial function with cross terms

f(x,y) = c0 (0th order)

  • c1 x + c2 y (1st order)

  • c3 x^2 + c4 x y + c5 y^2 (2nd order)

  • c6 x^3 + c7 x^2 y + c8 x y^2 + c9 y^3 (3rd order)

Intermediate products for the most recent y are cached, so when computing for a set of x, y it is more efficient to change x before you change y.

class ProductBoundedField : public lsst::afw::table::io::PersistableFacade<ProductBoundedField>, public lsst::afw::math::BoundedField
#include <ProductBoundedField.h>

A BoundedField that lazily multiplies a sequence of other BoundedFields.

class Random
#include <Random.h>

A class that can be used to generate sequences of random numbers according to a number of different algorithms. Support for generating random variates from the uniform, Gaussian, Poisson, and chi-squared distributions is provided.

This class is a thin wrapper for the random number generation facilities of GSL, which supports many additional distributions that can easily be added to this class as the need arises.

See

Random number generation in GSL

See

Random number distributions in GSL

class SeparableKernel : public lsst::afw::table::io::PersistableFacade<SeparableKernel>, public lsst::afw::math::Kernel
#include <Kernel.h>

A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)

The function’s x, y arguments are as follows:

  • -getCtr() for the lower left corner pixel

  • 0, 0 for the center pixel

  • (getDimensions() - 1) - getCtr() for the upper right pixel

Note: each pixel is set to the value of the kernel function at the center of the pixel (rather than averaging the function over the area of the pixel).

Subclassed by lsst::afw::math::BilinearWarpingKernel, lsst::afw::math::LanczosWarpingKernel, lsst::afw::math::NearestWarpingKernel

class SpatialCell
#include <SpatialCell.h>

Class to ensure constraints for spatial modeling

A given image is divided up into cells, with each cell represented by an instance of this class. Each cell itself contains a list of instances of classes derived from SpatialCellCandidate. One class member from each cell will be chosen to fit to a spatial model. In case of a poor fit, the next class instance in the list will be fit for. If all instances in a list are rejected from the spatial model, the best one will be used.

See

The SpatialCellSet example in the module documentation.

class SpatialCellCandidate
#include <SpatialCell.h>

Base class for candidate objects in a SpatialCell

Subclassed by lsst::afw::math::SpatialCellImageCandidate

class SpatialCellCandidateIterator
#include <SpatialCell.h>

An iterator that only returns usable members of the SpatialCell

class SpatialCellImageCandidate : public lsst::afw::math::SpatialCellCandidate
#include <SpatialCell.h>

Base class for candidate objects in a SpatialCell that are able to return an Image of some sort (e.g. a PSF or a DIA kernel)

Subclassed by lsst::ip::diffim::KernelCandidate< _PixelT >, lsst::meas::algorithms::PsfCandidate< PixelT >

class SpatialCellSet
#include <SpatialCell.h>

A collection of SpatialCells covering an entire image

template<typename T>
class SquaredExpCovariogram : public lsst::afw::math::Covariogram<T>
#include <GaussianProcess.h>

SquaredExpCovariogram

a Covariogram that falls off as the negative exponent of the square of the distance between the points

Contains one hyper parameter (_ellSquared) encoding the square of the characteristic length scale of the covariogram

class Statistics
#include <Statistics.h>

A class to evaluate image statistics

The basic strategy is to construct a Statistics object from an Image and a statement of what we want to know. The desired results can then be returned using Statistics methods. A StatisticsControl object is used to pass parameters. The statistics currently implemented are listed in the enum Properties in Statistics.h.

 // sets NumSigclip (3.0), and NumIter (3) for clipping
 lsst::afw::math::StatisticsControl sctrl(3.0, 3);

 sctrl.setNumSigmaClip(4.0);            // reset number of standard deviations for N-sigma clipping
 sctrl.setNumIter(5);                   // reset number of iterations for N-sigma clipping
 sctrl.setAndMask(0x1);                 // ignore pixels with these mask bits set
 sctrl.setNanSafe(true);                // check for NaNs & Infs, a bit slower (default=true)

 lsst::afw::math::Statistics statobj =
     lsst::afw::math::makeStatistics(*img, afwMath::NPOINT |
                                           afwMath::MEAN | afwMath::MEANCLIP, sctrl);
 double const n = statobj.getValue(lsst::afw::math::NPOINT);
 std::pair<double, double> const mean =
                                  statobj.getResult(lsst::afw::math::MEAN); // Returns (value, error)
 double const meanError = statobj.getError(lsst::afw::math::MEAN);                // just the error

Note

Factory function: We used a helper function, makeStatistics, rather that the constructor directly so that the compiler could deduce the types cf. std::make_pair()

Note

Inputs: The class Statistics is templated, and makeStatistics() can take either: (1) an image, (2) a maskedImage, or (3) a std::vector<> Overloaded makeStatistics() functions then wrap what they were passed in Image/Mask-like classes and call the Statistics constructor.

Note

Clipping: The clipping is done iteratively with numSigmaClip and numIter specified in the StatisticsControl object. The first clip (ie. the first iteration) is performed at: median +/- numSigmaClip*IQ_TO_STDEV*IQR, where IQ_TO_STDEV=~0.74 is the conversion factor between the IQR and sigma for a Gaussian distribution. All subsequent iterations perform clips at mean +/- numSigmaClip*stdev.

class StatisticsControl
#include <Statistics.h>

Pass parameters to a Statistics object

A class to pass parameters which control how the stats are calculated.

class TransformBoundedField : public lsst::afw::table::io::PersistableFacade<TransformBoundedField>, public lsst::afw::math::BoundedField
#include <TransformBoundedField.h>

A BoundedField based on geom::Transform<Poin2Endpoint, GenericEndpoint<1>>.

TransformBoundedField supports arbitrary transforms.

class WarpingControl
#include <warpExposure.h>

Parameters to control convolution

Note

padValue is not member of this class to avoid making this a templated class.

namespace detail

Functions

template<typename OutImageT, typename InImageT>
void basicConvolve(OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::Kernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)

Low-level convolution function that does not set edge pixels.

convolvedImage must be the same size as inImage. convolvedImage has a border in which the output pixels are not set. This border has size:

  • kernel.getCtrX() along the left edge

  • kernel.getCtrY() along the bottom edge

  • kernel.getWidth() - 1 - kernel.getCtrX() along the right edge

  • kernel.getHeight() - 1 - kernel.getCtrY() along the top edge

Parameters
  • [out] convolvedImage: convolved image

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if convolvedImage dimensions != inImage dimensions

  • lsst::pex::exceptions::InvalidParameterError: if inImage smaller than kernel in width or height

  • lsst::pex::exceptions::InvalidParameterError: if kernel width or height < 1

  • std::bad_alloc: when allocation of CPU memory fails

template<typename OutImageT, typename InImageT>
void basicConvolve(OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::DeltaFunctionKernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)

A version of basicConvolve that should be used when convolving delta function kernels

Parameters
  • [out] convolvedImage: convolved image

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

template<typename OutImageT, typename InImageT>
void basicConvolve(OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::LinearCombinationKernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)

A version of basicConvolve that should be used when convolving a LinearCombinationKernel

The Algorithm:

  • If the kernel is spatially varying and contains only DeltaFunctionKernels then convolves the input Image by each basis kernel in turn, solves the spatial model for that component and adds in the appropriate amount of the convolved image.

  • In all other cases uses normal convolution

Parameters
  • [out] convolvedImage: convolved image

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if convolvedImage dimensions != inImage dimensions

  • lsst::pex::exceptions::InvalidParameterError: if inImage smaller than kernel in width or height

  • lsst::pex::exceptions::InvalidParameterError: if kernel width or height < 1

  • std::bad_alloc: when allocation of CPU memory fails

template<typename OutImageT, typename InImageT>
void basicConvolve(OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::SeparableKernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)

A version of basicConvolve that should be used when convolving separable kernels

Parameters
  • [out] convolvedImage: convolved image

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

template<typename OutImageT, typename InImageT>
void convolveWithBruteForce(OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::Kernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)

Convolve an Image or MaskedImage with a Kernel by computing the kernel image at every point. (If the kernel is not spatially varying then only compute it once).

convolvedImage must be the same size as inImage. convolvedImage has a border in which the output pixels are not set. This border has size:

  • kernel.getCtrX() along the left edge

  • kernel.getCtrY() along the bottom edge

  • kernel.getWidth() - 1 - kernel.getCtrX() along the right edge

  • kernel.getHeight() - 1 - kernel.getCtrY() along the top edge

Warning

Low-level convolution function that does not set edge pixels.

Parameters
  • [out] convolvedImage: convolved image

  • [in] inImage: image to convolve

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if convolvedImage dimensions != inImage dimensions

  • lsst::pex::exceptions::InvalidParameterError: if inImage smaller than kernel in width or height

  • lsst::pex::exceptions::InvalidParameterError: if kernel width or height < 1

  • std::bad_alloc: when allocation of CPU memory fails

template<typename OutImageT, typename InImageT>
void convolveWithInterpolation(OutImageT &outImage, InImageT const &inImage, lsst::afw::math::Kernel const &kernel, ConvolutionControl const &convolutionControl)

Convolve an Image or MaskedImage with a spatially varying Kernel using linear interpolation.

This is a low-level convolution function that does not set edge pixels.

The algorithm is as follows:

  • divide the image into regions whose size is no larger than maxInterpolationDistance

  • for each region:

    • convolve it using convolveRegionWithInterpolation (which see)

Note that this routine will also work with spatially invariant kernels, but not efficiently.

Parameters
  • [out] outImage: convolved image = inImage convolved with kernel

  • [in] inImage: input image

  • [in] kernel: convolution kernel

  • [in] convolutionControl: convolution control parameters

Exceptions
  • lsst::pex::exceptions::InvalidParameterError: if outImage is not the same size as inImage

template<typename OutImageT, typename InImageT>
void convolveRegionWithInterpolation(OutImageT &outImage, InImageT const &inImage, KernelImagesForRegion const &region, ConvolveWithInterpolationWorkingImages &workingImages)

Convolve a region of an Image or MaskedImage with a spatially varying Kernel using interpolation.

This is a low-level convolution function that does not set edge pixels.

Warning

: this is a low-level routine that performs no bounds checking.

Parameters
  • [out] outImage: convolved image = inImage convolved with kernel

  • [in] inImage: input image

  • [in] region: kernel image region over which to convolve

  • [in] workingImages: working kernel images

struct ConvolveWithInterpolationWorkingImages
#include <Convolve.h>

kernel images used by convolveRegionWithInterpolation

class KernelImagesForRegion
#include <Convolve.h>

A collection of Kernel images for special locations on a rectangular region of an image

See the Location enum for a list of those special locations.

This is a low-level helper class for recursive convolving with interpolation. Many of these objects may be created during a convolution, and many will share kernel images. It uses shared pointers to kernels and kernel images for increased speed and decreased memory usage (at the expense of safety). Note that null pointers are NOT acceptable for the constructors!

Warning

The kernel images along the top and right edges are computed one row or column past the bounding box. This allows abutting KernelImagesForRegion to share corner and edge kernel images, which is useful when dividing a KernelImagesForRegion into subregions.

Warning

The bounding box for the region applies to the parent image.

Also note that it uses lazy evaluation: images are computed when they are wanted.

class RowOfKernelImagesForRegion
#include <Convolve.h>

A row of KernelImagesForRegion

Intended for iterating over subregions of a KernelImagesForRegion using computeNextRow.

class Spline

Subclassed by lsst::afw::math::detail::SmoothedSpline, lsst::afw::math::detail::TautSpline

struct TrapezoidalPacker
#include <TrapezoidalPacker.h>

A helper class ChebyshevBoundedField, for mapping trapezoidal matrices to 1-d arrays.

This class is not Swigged, and should not be included by any other .h files (including lsst/afw/math/detail.h); it’s for internal use by ChebyshevBoundedField only, and it’s only in a header file instead of that .cc file only so it can be unit tested.

We characterize the matrices by their number of columns (nx) and rows (ny), and the number of complete rows minus one (m).

This splits up the matrix into a rectangular part, in which the number of columns is the same for each row, and a wide trapezoidal or triangular part, in which the number of columns decreases by one for each row.

Here are some examples of how this class handles different kinds of matrices:

A wide trapezoidal matrix with orderX=4, orderY=3: nx=5, ny=4, m=0

0 1 2 3 4 5 6 7 8 9 10 11 12 13

A tall trapezoidal matrix with orderX=2, orderY=4 nx=3, ny=5, m=2

0 1 2 3 4 5 6 7 8 9 10 11

A triangular matrix with orderX=3, orderY=3 nx=4, ny=5, m=0

0 1 2 3 4 5 6 7 8 9

A wide rectangular matrix with orderX=3, orderY=2 nx=4, ny=3, m=3

0 1 2 3 4 5 6 7 8 9 10 11

A tall rectangular matrix with orderX=2, orderY=3 nx=3, ny=4, m=4

0 1 2 3 4 5 6 7 8 9 10 11

template<typename DestImageT, typename SrcImageT>
class WarpAtOnePoint
#include <WarpAtOnePoint.h>

A functor that computes one warped pixel

namespace details

Functions

template<class T>
T norm(const T &x)
template<class T>
T real(const T &x)
template<class T>
T Epsilon()
template<class T>
T MinRep()
template<class T>
T rescale_error(T err, T const &resabs, T const &resasc)
template<class UF>
bool intGKPNA(UF const &func, IntRegion<typename UF::result_type> &reg, typename UF::result_type const epsabs, typename UF::result_type const epsrel, std::map<typename UF::result_type, typename UF::result_type> *fxmap = 0)

Non-adaptive integration of the function f over the region ‘reg’.

Note

The algorithm computes first a Gaussian quadrature value then successive Kronrod/Patterson extensions to this result. The functions terminates when the difference between successive approximations (rescaled according to rescale_error) is less than either epsabs or epsrel * I, where I is the latest estimate of the integral. The order of the Gauss/Kronron/Patterson scheme is determined by which file is included above. Currently schemes starting with order 1 and order 10 are calculated. There seems to be little practical difference in the integration times using the two schemes, so I haven’t bothered to calculate any more.

template<class UF>
void intGKP(UF const &func, IntRegion<typename UF::result_type> &reg, typename UF::result_type const epsabs, typename UF::result_type const epsrel, std::map<typename UF::result_type, typename UF::result_type> *fxmap = 0)

An adaptive integration algorithm which computes the integral of f over the region reg.

The area and estimated error are returned as reg.Area() and reg.Err() If desired, *retx and *retf return std::vectors of x,f(x) respectively They only include the evaluations in the non-adaptive pass, so they do not give an accurate estimate of the number of function evaluations.

Note

First the non-adaptive GKP algorithm is tried. If that is not accurate enough (according to the absolute and relative accuracies, epsabs and epsrel), the region is split in half, and each new region is integrated. The routine continues by successively splitting the subregion which gave the largest absolute error until the integral converges.

template<class UF>
AuxFunc1<UF> Aux1(UF uf)

Auxiliary function 1

template<class UF>
AuxFunc2<UF> Aux2(UF uf)

Auxiliary function 2

template<class BF, class Tp>
binder2_1<BF> bind21(const BF &oper, const Tp &x)
template<class TF, class Tp>
binder3_1<TF> bind31(const TF &oper, const Tp &x)
int gkp_n(int level)
template<class T>
const std::vector<T> &gkp_x(int level)
template<class T>
const std::vector<T> &gkp_wa(int level)
template<class T>
const std::vector<T> &gkp_wb(int level)

Variables

const int NGKPLEVELS = 5
template<class UF>
struct AuxFunc1 : public std::unary_function<UF::argument_type, UF::result_type>
#include <Integrate.h>

Auxiliary struct 1

template<class T>
struct ConstantReg1 : public std::unary_function<T, IntRegion<T>>
#include <Integrate.h>

Helpers for constant regions for int2d, int3d:

template<typename BinaryFunctionT>
class FunctionWrapper : public std::unary_function<BinaryFunctionT::second_argument_type, BinaryFunctionT::result_type>
#include <Integrate.h>

Wrap an integrand in a call to a 1D integrator: romberg()

When romberg2D() is called, it wraps the integrand it was given in a FunctionWrapper functor. This wrapper calls romberg() on the integrand to get a 1D (along the x-coord, for constant y) result . romberg2D() then calls romberg() with the FunctionWrapper functor as an integrand.