lsst::meas::astrom::ScaledPolynomialTransformFitter Class Reference

A fitter class for scaled polynomial transforms. More...

#include <ScaledPolynomialTransformFitter.h>

Classes
class	Keys

Public Member Functions
void	fit (int order=-1)
	Perform a linear least-squares fit of the polynomial coefficients.
void	updateModel ()
	Update the 'model' field in the data catalog using the current best- fit transform.
double	updateIntrinsicScatter ()
	Infer the intrinsic scatter in the offset between the data points and the best-fit model, and update the uncertainties accordingly.
double	getIntrinsicScatter () const
	Return the current estimate of the intrinsic scatter.
std::pair< double, size_t >	rejectOutliers (OutlierRejectionControl const &ctrl)
	Mark outliers in the data catalog using sigma clipping.
afw::table::BaseCatalog const &	getData () const
	Return a catalog of data points and model values for diagnostic purposes.
ScaledPolynomialTransform const &	getTransform () const
	Return the best-fit transform.
PolynomialTransform const &	getPoly () const
	Return the polynomial part of the best-fit transform.
geom::AffineTransform const &	getInputScaling () const
	Return the input scaling transform that maps input data points to [-1, 1].
geom::AffineTransform const &	getOutputScaling () const
	Return the output scaling transform that maps output data points to [-1, 1].

Static Public Member Functions
static ScaledPolynomialTransformFitter	fromMatches (int maxOrder, afw::table::ReferenceMatchVector const &matches, afw::geom::SkyWcs const &initialWcs, double intrinsicScatter)
	Initialize a fit from intermediate world coordinates to pixels using source/reference matches.
static ScaledPolynomialTransformFitter	fromGrid (int maxOrder, geom::Box2D const &bbox, int nGridX, int nGridY, ScaledPolynomialTransform const &toInvert)
	Initialize a fit that inverts an existing transform by evaluating and fitting to points on a grid.

Detailed Description

A fitter class for scaled polynomial transforms.

This class allows for iteration between actual fitting, outlier rejection, and estimation of intrinsic scatter. It also provides access to the current model for debugging via an afw::table::BaseCatalog that contains the input values, output values, and the best-fit-transformed input values.

ScaledPolynomialTransformFitter has two public construction methods:

• fromMatches fits a transform that maps intermediate world coordinates to pixel coordinates, using as inputs a match of reference coordinates to measured source positions. This sets up the fitter to do outlier rejection and intrinsic scatter estimation.
• fromGrid fits an inverse to an existing transform, using the to-be-inverted transform to populate a grid of positions. Because the to-be-inverted transform is assumed to be known exactly, fitters initialized with this method have no need for outlier rejection or scatter estimation.

In either case, the fitter creates affine transforms that map the input and output data points onto [-1, 1]. It then fits a polynomial transform that, when composed with the input scaling transform and the inverse of the output scaling transform, maps the input data points to the output data points.

The fitter can be used in an outlier-rejection loop with the following pattern (with user-defined convergence criteria):

while (true) {
    fitter.fit();
    if (converged) break;
    fitter.updateModel();
    fitter.computeIntrinsicScatter();
    fitter.rejectOutliers();
}

This pattern fits a model, uses that model to transform the input points, estimates intrinsic scatter from the differences between model-transformed points and output data points, and finally rejects outliers by sigma-clipping those differences before repeating.

ScaledPolynomialTransformFitter instances should be confined to a single thread.

Definition at line 101 of file ScaledPolynomialTransformFitter.h.

Member Function Documentation

fit()

void lsst::meas::astrom::ScaledPolynomialTransformFitter::fit

(

int

order = -1

)

Perform a linear least-squares fit of the polynomial coefficients.

Parameters

[in]

order

The maximum order of the polynomial transform. If negative (the default) the maxOrder from construction is used.

After fitting, updateModel() should be called to apply the model transform to the input points if the user wants to make use of the data catalog for diagnostics. Calling fit() alone is sufficient if the user is only interested in getting the model transform itself.

Definition at line 200 of file ScaledPolynomialTransformFitter.cc.

                                                   {
    int maxOrder = _transform.getPoly().getOrder();
    if (order < 0) {
        order = maxOrder;
    }
    if (order > maxOrder) {
        throw LSST_EXCEPT(
                pex::exceptions::LengthError,
                (boost::format("Order (%d) exceeded maximum order for the fitter (%d)") % order % maxOrder)
                        .str());
    }
 
    int const packedSize = detail::computePackedSize(order);
    std::size_t nGood = 0;
    if (_keys.rejected.isValid()) {
        for (auto const &record : _data) {
            if (!record.get(_keys.rejected)) {
                ++nGood;
            }
        }
    } else {
        nGood = _data.size();
    }
    // One block of the block-diagonal (2x2) unweighted design matrix M;
    // m[i,j] = u_i^{p(j)} v_i^{q(j)}.  The two nonzero blocks are the same,
    // because we're using the same polynomial basis for x and y.
    Eigen::MatrixXd m = Eigen::MatrixXd::Zero(nGood, packedSize);
    // vx, vy: (2x1) blocks of the unweighted data vector v
    Eigen::VectorXd vx = Eigen::VectorXd::Zero(nGood);
    Eigen::VectorXd vy = Eigen::VectorXd::Zero(nGood);
    // sxx, syy, sxy: (2x2) blocks of the covariance matrix S, each of which is
    // individually diagonal.
    Eigen::ArrayXd sxx(nGood);
    Eigen::ArrayXd syy(nGood);
    Eigen::ArrayXd sxy(nGood);
    Eigen::Matrix2d outS = _outputScaling.getLinear().getMatrix();
    for (std::size_t i1 = 0, i2 = 0; i1 < _data.size(); ++i1) {
        // check that the 'rejected' field (== 'not outlier-rejected') is both
        // present in the schema and not rejected.
        if (!_keys.rejected.isValid() || !_data[i1].get(_keys.rejected)) {
            geom::Point2D output = _outputScaling(_data[i1].get(_keys.output));
            vx[i2] = output.getX();
            vy[i2] = output.getY();
            m.row(i2) = _vandermonde.row(i1).head(packedSize);
            if (_keys.outputErr.isValid()) {
                Eigen::Matrix2d modelErr =
                        outS * _data[i1].get(_keys.outputErr).cast<double>() * outS.adjoint();
                sxx[i2] = modelErr(0, 0);
                sxy[i2] = modelErr(0, 1);
                syy[i2] = modelErr(1, 1);
            } else {
                sxx[i2] = 1.0;
                sxy[i2] = 0.0;
                syy[i2] = 1.0;
            }
            ++i2;
        }
    }
    // Do a blockwise inverse of S.  Note that the result F is still symmetric
    Eigen::ArrayXd fxx = 1.0 / (sxx - sxy.square() / syy);
    Eigen::ArrayXd fyy = 1.0 / (syy - sxy.square() / sxx);
    Eigen::ArrayXd fxy = -(sxy / sxx) * fyy;
#ifdef LSST_ScaledPolynomialTransformFitter_TEST_IN_PLACE
    assert((sxx * fxx + sxy * fxy).isApproxToConstant(1.0));
    assert((syy * fyy + sxy * fxy).isApproxToConstant(1.0));
    assert((sxx * fxy).isApprox(-sxy * fyy));
    assert((sxy * fxx).isApprox(-syy * fxy));
#endif
    // Now that we've got all the block quantities, we'll form the full normal equations matrix.
    // That's H = M^T F M:
    Eigen::MatrixXd h(2 * packedSize, 2 * packedSize);
    h.topLeftCorner(packedSize, packedSize) = m.adjoint() * fxx.matrix().asDiagonal() * m;
    h.topRightCorner(packedSize, packedSize) = m.adjoint() * fxy.matrix().asDiagonal() * m;
    h.bottomLeftCorner(packedSize, packedSize) = h.topRightCorner(packedSize, packedSize).adjoint();
    h.bottomRightCorner(packedSize, packedSize) = m.adjoint() * fyy.matrix().asDiagonal() * m;
    // And here's the corresponding RHS vector, g = M^T F v
    Eigen::VectorXd g(2 * packedSize);
    g.head(packedSize) = m.adjoint() * (fxx.matrix().asDiagonal() * vx + fxy.matrix().asDiagonal() * vy);
    g.tail(packedSize) = m.adjoint() * (fxy.matrix().asDiagonal() * vx + fyy.matrix().asDiagonal() * vy);
    // Solve the normal equations.
    auto lstsq = afw::math::LeastSquares::fromNormalEquations(h, g);
    auto solution = lstsq.getSolution();
    // Unpack the solution vector back into the polynomial coefficient matrices.
    for (int n = 0, j = 0; n <= order; ++n) {
        for (int p = 0, q = n; p <= n; ++p, --q, ++j) {
            _transform._poly._xCoeffs(p, q) = solution[j];
            _transform._poly._yCoeffs(p, q) = solution[j + packedSize];
        }
    }
}

fromGrid()

ScaledPolynomialTransformFitter lsst::meas::astrom::ScaledPolynomialTransformFitter::fromGrid ( int maxOrder, geom::Box2D const & bbox, int nGridX, int nGridY, ScaledPolynomialTransform const & toInvert )

static

Initialize a fit that inverts an existing transform by evaluating and fitting to points on a grid.

Parameters

[in]	maxOrder	Maximum polynomial order for the fit.
[in]	bbox	Bounding box for the grid in the input coordinates of the toInvert transform (or equivalently the output coordinates of the transform to be fit).
[in]	nGridX	Number of grid points in the X direction.
[in]	nGridY	Number of grid points in the Y direction.
[in]	toInvert	Transform to invert

This initializes the data catalog with the following fields:

• output: grid positions passed as input to the toInvert transform, treated as output data points when fitting its inverse.
• input: grid positions generated as the output of the toInvert transform, treated as input data points when fitting its inverse.
• model: best-fit-transformed input points.

The updateIntrinsicScatter and rejectOutliers methods cannot be used on fitters initialized using this method.

Definition at line 150 of file ScaledPolynomialTransformFitter.cc.

                                                   {
    Keys const &keys = Keys::forGrid();
    afw::table::BaseCatalog catalog(keys.schema);
    catalog.reserve(nGridX * nGridY);
    geom::Extent2D dx(bbox.getWidth() / nGridX, 0.0);
    geom::Extent2D dy(0.0, bbox.getHeight() / nGridY);
    for (int iy = 0; iy < nGridY; ++iy) {
        for (int ix = 0; ix < nGridX; ++ix) {
            geom::Point2D point = bbox.getMin() + dx * ix + dy * iy;
            auto record = catalog.addNew();
            record->set(keys.output, point);
            record->set(keys.input, toInvert(point));
        }
    }
    return ScaledPolynomialTransformFitter(catalog, keys, maxOrder, 0.0, computeScaling(catalog, keys.input),
                                           computeScaling(catalog, keys.output));
}

fromMatches()

ScaledPolynomialTransformFitter lsst::meas::astrom::ScaledPolynomialTransformFitter::fromMatches ( int maxOrder, afw::table::ReferenceMatchVector const & matches, afw::geom::SkyWcs const & initialWcs, double intrinsicScatter )

static

Initialize a fit from intermediate world coordinates to pixels using source/reference matches.

Parameters

[in]	maxOrder	Maximum polynomial order for the fit.
[in]	matches	Vector of source-reference matches. The centroids and centroid errors from the sources are used, along with the coords from the reference objects.
[in]	initialWcs	Initial WCS, used to transform reference positions to intermediate world coordinates and populate the "ref" field in the data catalog.
[in]	intrinsicScatter	Initial intrinsic scatter to be added in quadrature with source measurement errors in setting the uncertainty on each match.

This initializes the data catalog with the following fields:

• ref_id: reference object ID
• src_id: source ID from the match
• src_{x,y}: source position in pixels
• src_{xErr,yErr,x_y_Cov}: uncertainty on source positions, including measurement errors and the inferred intrinsic scatter.
• intermediate_{x,y}: reference positions in intermediate world coordinates.
• initial_{x,y}: reference positions transformed by initialWcs.
• model_{x,y}: reference positions transformed by current best-fit distortion.
• rejected Flag that is set for objects that have been rejected as outliers.

Definition at line 127 of file ScaledPolynomialTransformFitter.cc.

                                 {
    Keys const &keys = Keys::forMatches();
    afw::table::BaseCatalog catalog(keys.schema);
    catalog.reserve(matches.size());
    float var2 = intrinsicScatter * intrinsicScatter;
    auto initialIwcToSky = getIntermediateWorldCoordsToSky(initialWcs);
    for (auto const &match : matches) {
        auto record = catalog.addNew();
        record->set(keys.refId, match.first->getId());
        record->set(keys.srcId, match.second->getId());
        record->set(keys.input, initialIwcToSky->applyInverse(match.first->getCoord()));
        record->set(keys.initial, initialWcs.skyToPixel(match.first->getCoord()));
        record->set(keys.output, match.second->getCentroid());
        record->set(keys.outputErr, match.second->getCentroidErr() + var2 * Eigen::Matrix2f::Identity());
        record->set(keys.rejected, false);
    }
    return ScaledPolynomialTransformFitter(catalog, keys, maxOrder, intrinsicScatter,
                                           computeScaling(catalog, keys.input),
                                           computeScaling(catalog, keys.output));
}

getData()

afw::table::BaseCatalog const & lsst::meas::astrom::ScaledPolynomialTransformFitter::getData ( ) const

inline

Return a catalog of data points and model values for diagnostic purposes.

The values in the returned catalog should not be modified by the user.

For information about the schema, either introspect it programmatically or see fromMatches and fromGrid.

Definition at line 249 of file ScaledPolynomialTransformFitter.h.

{ return _data; }

getInputScaling()

geom::AffineTransform const & lsst::meas::astrom::ScaledPolynomialTransformFitter::getInputScaling ( ) const

inline

Return the input scaling transform that maps input data points to [-1, 1].

Definition at line 270 of file ScaledPolynomialTransformFitter.h.

{ return _transform.getInputScaling(); }

getIntrinsicScatter()

double lsst::meas::astrom::ScaledPolynomialTransformFitter::getIntrinsicScatter ( ) const

inline

Return the current estimate of the intrinsic scatter.

Because the intrinsic scatter is included in the uncertainty used in both the fit and outlier rejection, this may be called either before updateIntrinsicScatter() to obtain the scatter used in the last such operation or after updateIntrinsicScatter() to obtain the scatter to be used in the next operation.

Definition at line 219 of file ScaledPolynomialTransformFitter.h.

{ return _intrinsicScatter; }

getOutputScaling()

geom::AffineTransform const & lsst::meas::astrom::ScaledPolynomialTransformFitter::getOutputScaling ( ) const

inline

Return the output scaling transform that maps output data points to [-1, 1].

Definition at line 275 of file ScaledPolynomialTransformFitter.h.

{ return _outputScaling; }

getPoly()

PolynomialTransform const & lsst::meas::astrom::ScaledPolynomialTransformFitter::getPoly ( ) const

inline

Return the polynomial part of the best-fit transform.

Definition at line 265 of file ScaledPolynomialTransformFitter.h.

{ return _transform.getPoly(); }

getTransform()

ScaledPolynomialTransform const & lsst::meas::astrom::ScaledPolynomialTransformFitter::getTransform ( ) const

inline

Return the best-fit transform.

If f is an instance of this class, then for an arbitrary point p:

f.getTransform()(p) == f.getOutputScaling().inverted()(f.getPoly()(f.getInputScaling()(p)))

Definition at line 260 of file ScaledPolynomialTransformFitter.h.

{ return _transform; }

rejectOutliers()

std::pair< double, std::size_t > lsst::meas::astrom::ScaledPolynomialTransformFitter::rejectOutliers

(

OutlierRejectionControl const &

ctrl

)

Mark outliers in the data catalog using sigma clipping.

A data point is declare an outlier if:

• the offset between the model prediction and the data position, weighted by the uncertainty (including inferred intrinsic scatter) of that point exceeds ctrl.nSigma
• the number of points already with larger weighted offsets is less than ctrl.nClipMax. In addition, the worst (in weighted offset) ctrl.nClipMin data points are always rejected.

Returns

A pair of the smallest rejected weighted offset (units of sigma) and the number of point rejected.

Should be called after updateIntrinsicScatter() and before fit() if the fitter is being used in an outlier rejection loop.

Definition at line 367 of file ScaledPolynomialTransformFitter.cc.

                                             {
    // If the 'rejected' field isn't present in the schema (because the fitter
    // was constructed with fromGrid), we can't do outlier rejection.
    if (!_keys.rejected.isValid()) {
        throw LSST_EXCEPT(pex::exceptions::LogicError,
                          "Cannot reject outliers on fitter initialized with fromGrid.");
    }
    if (static_cast<std::size_t>(ctrl.nClipMin) >= _data.size()) {
        throw LSST_EXCEPT(
                pex::exceptions::LogicError,
                (boost::format("Not enough values (%d) to clip %d.") % _data.size() % ctrl.nClipMin).str());
    }
    std::map<double, afw::table::BaseRecord *> rankings;
    for (auto &record : _data) {
        Eigen::Matrix2d cov = record.get(_keys.outputErr).cast<double>();
        Eigen::Vector2d d = (record.get(_keys.output) - record.get(_keys.model)).asEigen();
        double r2 = d.dot(cov.inverse() * d);
        rankings.insert(std::make_pair(r2, &record));
    }
    auto cutoff = rankings.upper_bound(ctrl.nSigma * ctrl.nSigma);
    int nClip = 0, nGood = 0;
    for (auto iter = rankings.begin(); iter != cutoff; ++iter) {
        iter->second->set(_keys.rejected, false);
        ++nGood;
    }
    for (auto iter = cutoff; iter != rankings.end(); ++iter) {
        iter->second->set(_keys.rejected, true);
        ++nClip;
    }
    assert(static_cast<std::size_t>(nGood + nClip) == _data.size());
    while (nClip < ctrl.nClipMin) {
        --cutoff;
        cutoff->second->set(_keys.rejected, true);
        ++nClip;
    }
    while (nClip > ctrl.nClipMax && cutoff != rankings.end()) {
        cutoff->second->set(_keys.rejected, false);
        ++cutoff;
        --nClip;
    }
    std::pair<double, std::size_t> result(ctrl.nSigma, nClip);
    if (cutoff != rankings.end()) {
        result.first = std::sqrt(cutoff->first);
    }
    return result;
}

updateIntrinsicScatter()

double lsst::meas::astrom::ScaledPolynomialTransformFitter::updateIntrinsicScatter

(

)

Infer the intrinsic scatter in the offset between the data points and the best-fit model, and update the uncertainties accordingly.

The scatter is calculated as the RMS scatter that maximizes the likelihood of the current positional offsets (assumed fixed) assuming the per-data-point uncertainties are the quadrature sum of the intrinsic scatter and that source's centroid measurement uncertainty.

Should be called after updateModel() and before rejectOutliers() if the fitter is being used in an outlier-rejection loop.

Definition at line 297 of file ScaledPolynomialTransformFitter.cc.

                                                               {
    if (!_keys.rejected.isValid()) {
        throw LSST_EXCEPT(pex::exceptions::LogicError,
                          "Cannot compute intrinsic scatter on fitter initialized with fromGrid.");
    }
    double newIntrinsicScatter = computeIntrinsicScatter();
    float varDiff = newIntrinsicScatter * newIntrinsicScatter - _intrinsicScatter * _intrinsicScatter;
    for (auto &record : _data) {
        record.set(_keys.outputErr, record.get(_keys.outputErr) + varDiff * Eigen::Matrix2f::Identity());
    }
    _intrinsicScatter = newIntrinsicScatter;
    return _intrinsicScatter;
}

updateModel()

void lsst::meas::astrom::ScaledPolynomialTransformFitter::updateModel

(

)

Update the 'model' field in the data catalog using the current best- fit transform.

Should be called after fit() and before updateIntrinsicScatter() if the fitter is being used in an outlier-rejection loop.

Definition at line 291 of file ScaledPolynomialTransformFitter.cc.

                                                  {
    for (auto &record : _data) {
        record.set(_keys.model, _transform(record.get(_keys.input)));
    }
}

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-8.0.0/Linux64/meas_astrom/gff1a9f87cc+86cf3d8bc9/include/lsst/meas/astrom/ScaledPolynomialTransformFitter.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-8.0.0/Linux64/meas_astrom/gff1a9f87cc+86cf3d8bc9/src/ScaledPolynomialTransformFitter.cc