Interpolation Tools

C++ Interpolants

class Interpolant

Base class representing one-dimensional interpolant functions.

One-dimensional interpolant function. X units are in pixels and the frequency-domain u values are in cycles per pixel. This differs from the usual definition of k in radians per arcsec, hence the different letter variable.

All Interpolants are assumed symmetric so that frequency-domain values are real.

Subclassed by galsim::Cubic, galsim::Delta, galsim::Lanczos, galsim::Linear, galsim::Nearest, galsim::Quintic, galsim::SincInterpolant

Public Functions

inline Interpolant(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline Interpolant(const Interpolant &rhs)

Copy constructor: does not copy photon sampler, will need to rebuild.

inline virtual ~Interpolant()

Destructor.

virtual double xrange() const = 0

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

virtual int ixrange() const = 0

The total range as an integer. Typically xrange() == 0.5 * ixrange().

virtual double urange() const = 0

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

virtual double xval(double x) const = 0

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double xvalWrapped(double x, int N) const

Value of interpolant, wrapped at period N.

This returns sum_{j=-inf}^{inf} xval(x + jN):

Parameters:

• x – [in] Distance from sample (pixels)

• N – [in] Wrapping period (pixels)

Returns:

Value of interpolant after wrapping

void xvalMany(double *x, int N) const

Calculate xval for array of input values x.

Parameters:

• [in/out] – Each x[i] is replaces by xval[x[i]]

• How – [in] many x values to calculate

virtual double uval(double u) const = 0

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

void uvalMany(double *u, int N) const

Calculate uval for array of input values u.

Parameters:

• [in/out] – Each u[i] is replaces by uval[u[i]]

• How – [in] many u values to calculate

inline virtual bool isExactAtNodes() const

Report whether interpolation will reproduce values at samples.

This will return true if the interpolant is exact at nodes, meaning that F(0)=1 and F(n)=0 for non-zero integer n. Right now this is true for every implementation.

Returns:

True if samples are returned exactly.

inline virtual double getPositiveFlux() const

Return the integral of the positive portions of the kernel.

Should return 1 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of positive portions of kernel

inline virtual double getNegativeFlux() const

Return the (absolute value of) integral of the negative portions of the kernel.

Should return 0 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of abs value of negative portions of kernel

double getPositiveFlux2d() const

Return the positive and negative flux for a 2d application of the interpolant.

double getNegativeFlux2d() const

inline virtual void shoot(PhotonArray &photons, UniformDeviate ud) const

Return array of displacements drawn from this kernel.

Since Interpolant is 1d, will use only x array of PhotonArray. It will be assumed that photons returned are randomly ordered (no need to shuffle them). Also assumed that all photons will have nearly equal absolute value of flux. Total flux returned may not equal 1 due to shot noise in negative/positive photons, and small fluctuations in photon weights.

Parameters:

• N – [in] number of photons to shoot

• ud – [in] UniformDeviate used to generate random values

Returns:

a PhotonArray containing the vector of displacements for interpolation kernel.

virtual std::string makeStr() const = 0

class Delta : public galsim::Interpolant

Delta-function interpolant in 1d.

The interpolant for when you do not want to interpolate between samples. It is not really intended to be used for any analytic drawing because it is infinite in the x domain at the location of samples, and it extends to infinity in the u domain. But it could be useful for photon-shooting, where it is trivially implemented as no displacements. The argument in the constructor is used to make a crude box approximation to the x-space delta function and to give a large but finite urange.

Public Functions

inline Delta(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~Delta()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

inline virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

inline virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

inline virtual double getPositiveFlux() const

Return the integral of the positive portions of the kernel.

Should return 1 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of positive portions of kernel

inline virtual double getNegativeFlux() const

Return the (absolute value of) integral of the negative portions of the kernel.

Should return 0 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of abs value of negative portions of kernel

virtual void shoot(PhotonArray &photons, UniformDeviate ud) const

Return array of displacements drawn from this kernel.

Since Interpolant is 1d, will use only x array of PhotonArray. It will be assumed that photons returned are randomly ordered (no need to shuffle them). Also assumed that all photons will have nearly equal absolute value of flux. Total flux returned may not equal 1 due to shot noise in negative/positive photons, and small fluctuations in photon weights.

Parameters:

• N – [in] number of photons to shoot

• ud – [in] UniformDeviate used to generate random values

Returns:

a PhotonArray containing the vector of displacements for interpolation kernel.

virtual std::string makeStr() const

class Nearest : public galsim::Interpolant

Nearest-neighbor interpolation: boxcar.

The nearest-neighbor interpolant performs poorly as a k-space or x-space interpolant for SBInterpolatedImage. (See paper by Bernstein & Gruen, http://arxiv.org/abs/1401.2636.) The objection to its use in Fourier space does not apply when shooting photons to generate an image; in that case, the nearest-neighbor interpolant is quite efficient (but not necessarily the best choice in terms of accuracy).

Tolerance determines how far onto sinc wiggles the uval will go. Very far, by default!

Public Functions

inline Nearest(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~Nearest()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

inline virtual double getPositiveFlux() const

Return the integral of the positive portions of the kernel.

Should return 1 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of positive portions of kernel

inline virtual double getNegativeFlux() const

Return the (absolute value of) integral of the negative portions of the kernel.

Should return 0 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of abs value of negative portions of kernel

virtual void shoot(PhotonArray &photons, UniformDeviate ud) const

Return array of displacements drawn from this kernel.

Since Interpolant is 1d, will use only x array of PhotonArray. It will be assumed that photons returned are randomly ordered (no need to shuffle them). Also assumed that all photons will have nearly equal absolute value of flux. Total flux returned may not equal 1 due to shot noise in negative/positive photons, and small fluctuations in photon weights.

Parameters:

• N – [in] number of photons to shoot

• ud – [in] UniformDeviate used to generate random values

Returns:

a PhotonArray containing the vector of displacements for interpolation kernel.

virtual std::string makeStr() const

class Linear : public galsim::Interpolant

Linear interpolant.

The linear interpolant is a poor choice for FFT-based operations on SBInterpolatedImage, as it rings to high frequencies. (See paper by Bernstein & Gruen, http://arxiv.org/abs/1401.2636.) This objection does not apply when shooting photons, in which case the linear interpolant is quite efficient (but not necessarily the best choice in terms of accuracy).

Public Functions

inline Linear(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~Linear()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

inline virtual double getPositiveFlux() const

Return the integral of the positive portions of the kernel.

Should return 1 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of positive portions of kernel

inline virtual double getNegativeFlux() const

Return the (absolute value of) integral of the negative portions of the kernel.

Should return 0 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of abs value of negative portions of kernel

virtual void shoot(PhotonArray &photons, UniformDeviate ud) const

Return array of displacements drawn from this kernel.

Since Interpolant is 1d, will use only x array of PhotonArray. It will be assumed that photons returned are randomly ordered (no need to shuffle them). Also assumed that all photons will have nearly equal absolute value of flux. Total flux returned may not equal 1 due to shot noise in negative/positive photons, and small fluctuations in photon weights.

Parameters:

• N – [in] number of photons to shoot

• ud – [in] UniformDeviate used to generate random values

Returns:

a PhotonArray containing the vector of displacements for interpolation kernel.

virtual std::string makeStr() const

class Cubic : public galsim::Interpolant

Cubic interpolator exact to 3rd order Taylor expansion.

From R. G. Keys, IEEE Trans. Acoustics, Speech, & Signal Proc 29, p 1153, 1981

The cubic interpolant is a reasonable choice for a four-point interpolant for SBInterpolatedImage. (See paper by Bernstein & Gruen, http://arxiv.org/abs/1401.2636.)

Public Functions

Cubic(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~Cubic()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

inline virtual double getPositiveFlux() const

Return the integral of the positive portions of the kernel.

Should return 1 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of positive portions of kernel

inline virtual double getNegativeFlux() const

Return the (absolute value of) integral of the negative portions of the kernel.

Should return 0 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of abs value of negative portions of kernel

virtual std::string makeStr() const

class Quintic : public galsim::Interpolant

Piecewise-quintic polynomial interpolant, ideal for Fourier-space interpolation.

See paper by Bernstein & Gruen, http://arxiv.org/abs/1401.2636.

Public Functions

Quintic(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~Quintic()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

inline virtual double getPositiveFlux() const

Return the integral of the positive portions of the kernel.

Should return 1 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of positive portions of kernel

inline virtual double getNegativeFlux() const

Return the (absolute value of) integral of the negative portions of the kernel.

Should return 0 unless the kernel has negative portions. Default is to ask the numerical sampler for its stored value.

Returns:

Integral of abs value of negative portions of kernel

virtual std::string makeStr() const

class SincInterpolant : public galsim::Interpolant

Sinc interpolation: inverse of Nearest-neighbor.

The Sinc interpolant (K(x) = sin(pi x)/(pi x)) is mathematically perfect for band-limited data, introducing no spurious frequency content beyond kmax = pi/dx for input data with pixel scale dx. However, it is formally infinite in extent and, even with reasonable trunction, is still quite large. It will give exact results in SBInterpolatedImage::kValue() when it is used as a k-space interpolant, but is extremely slow. The usual compromise between sinc accuracy vs. speed is the Lanczos interpolant (see its documentation for details).

Public Functions

inline SincInterpolant(const GSParams &gsparams)

Constructor.

Parameters:

gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~SincInterpolant()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double xvalWrapped(double x, int N) const

Value of interpolant, wrapped at period N.

This returns sum_{j=-inf}^{inf} xval(x + jN):

Parameters:

• x – [in] Distance from sample (pixels)

• N – [in] Wrapping period (pixels)

Returns:

Value of interpolant after wrapping

virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

virtual void shoot(PhotonArray &photons, UniformDeviate ud) const

Return array of displacements drawn from this kernel.

Since Interpolant is 1d, will use only x array of PhotonArray. It will be assumed that photons returned are randomly ordered (no need to shuffle them). Also assumed that all photons will have nearly equal absolute value of flux. Total flux returned may not equal 1 due to shot noise in negative/positive photons, and small fluctuations in photon weights.

Parameters:

• N – [in] number of photons to shoot

• ud – [in] UniformDeviate used to generate random values

Returns:

a PhotonArray containing the vector of displacements for interpolation kernel.

virtual std::string makeStr() const

class Lanczos : public galsim::Interpolant

The Lanczos interpolation filter, nominally sinc(x)*sinc(x/n), truncated at +/-n.

The Lanczos filter is an approximation to the band-limiting sinc filter with a smooth cutoff at high x. Order n Lanczos has a range of +/- n pixels. It typically is a good compromise between kernel size and accuracy.

Note that pure Lanczos, when interpolating a set of constant-valued samples, does not return this constant. Setting conserve_dc in the constructor tweaks the function so that it approximately conserves the value of constant (DC) input data (accurate to better than 1.e-5 when used in two dimensions).

Public Functions

Lanczos(int n, bool conserve_dc, const GSParams &gsparams)

Constructor.

Parameters:

• n – [in] Filter order; must be given on input and cannot be changed.

• conserve_dc – [in] Set true to adjust filter to be more nearly correct for constant inputs.

• gsparams – [in] GSParams object storing constants that control the accuracy of operations, if different from the default.

inline ~Lanczos()

inline virtual double xrange() const

Maximum extent of interpolant from origin in x space (pixels)

Returns:

Range of non-zero values of interpolant.

inline virtual int ixrange() const

The total range as an integer. Typically xrange() == 0.5 * ixrange().

inline virtual double urange() const

Maximum extent of interpolant from origin in u space (cycles per pixel)

Returns:

Range of non-zero values of interpolant in u space

inline int getN() const

inline bool conservesDC() const

virtual double xval(double x) const

Value of interpolant in real space.

Parameters:

x – [in] Distance from sample (pixels)

Returns:

Value of interpolant

virtual double uval(double u) const

Value of interpolant in frequency space.

Parameters:

u – [in] Frequency for evaluation (cycles per pixel)

Returns:

Value of interpolant, normalized so uval(0) = 1 for flux-conserving interpolation.

virtual std::string makeStr() const

C++ Lookup Tables

class Table : public galsim::FluxDensity

A class to represent lookup tables for a function y = f(x).

Subclassed by galsim::TableBuilder

Public Types

enum interpolant

Values:

enumerator linear

enumerator floor

enumerator ceil

enumerator nearest

enumerator spline

enumerator gsinterp

Public Functions

Table(const double *args, const double *vals, int N, interpolant in)

Table from args, vals.

Table(const double *args, const double *vals, int N, const Interpolant *gsinterp)

double argMin() const

double argMax() const

size_t size() const

virtual double operator()(double a) const

interp, return double(0) if beyond bounds This is a virtual function from FluxDensity, which lets a Table be a FluxDensity.

double lookup(double a) const

interp, but exception if beyond bounds

void interpMany(const double *argvec, double *valvec, int N) const

interp many values at once

double integrate(double xmin, double xmax) const

double integrateProduct(const Table &g, double xmin, double xmax, double xfact) const

class TableBuilder : public galsim::Table

Public Functions

inline TableBuilder(interpolant in)

inline TableBuilder(const Interpolant *gsinterp)

inline bool finalized() const

inline double lookup(double a) const

inline void addEntry(double x, double f)

Insert an (x, y(x)) pair into the table.

void finalize()

class Table2D

A class to represent lookup tables for a function z = f(x, y).

Public Types

enum interpolant

Values:

enumerator linear

enumerator floor

enumerator ceil

enumerator nearest

enumerator spline

enumerator gsinterp

Public Functions

Table2D(const double *xargs, const double *yargs, const double *vals, int Nx, int Ny, interpolant in)

Table from xargs, yargs, vals.

Table2D(const double *xargs, const double *yargs, const double *vals, int Nx, int Ny, const double *dfdx, const double *dfdy, const double *d2fdxdy)

Table2D(const double *xargs, const double *yargs, const double *vals, int Nx, int Ny, const Interpolant *gsinterp)

double lookup(double x, double y) const

interp

void interpMany(const double *xvec, const double *yvec, double *valvec, int N) const

interp many values at once

void interpGrid(const double *xvec, const double *yvec, double *valvec, int Nx, int Ny) const

void gradient(double x, double y, double &dfdx, double &dfdy) const

Estimate df/dx, df/dy at a single location.

void gradientMany(const double *xvec, const double *yvec, double *dfdxvec, double *dfdyvec, int N) const

Estimate many df/dx and df/dy values.

void gradientGrid(const double *xvec, const double *yvec, double *dfdxvec, double *dfdyvec, int Nx, int Ny) const