lsst.meas.astrom Namespace Reference

Namespaces
	anetAstrometry
	anetBasicAstrometry
	approximateWcs
	astrometry
	astrometryNetDataConfig
	catalogStarSelector
	createMatchMetadata
	detail
	directMatch
	display
	fitSipDistortion
	fitTanSipWcs
	loadAstrometryNetObjects
	matchOptimisticB
	multiindex
	ref_match
	setMatchDistance
	sip
	verifyWcs
	version

Classes
struct	RecordProxy
	A wrapper around a SimpleRecord or SourceRecord that allows us to record a pixel position in a way that is independent of the record type. More...
struct	ProxyPair
struct	MatchOptimisticBControl
class	PolynomialTransform
	A 2-d coordinate transform represented by a pair of standard polynomials (one for each coordinate). More...
class	ScaledPolynomialTransform
	A 2-d coordinate transform represented by a lazy composition of an AffineTransform, a PolynomialTransform, and another AffineTransform. More...
class	OutlierRejectionControl
	Control object for outlier rejection in ScaledPolynomialTransformFitter. More...
class	ScaledPolynomialTransformFitter
	A fitter class for scaled polynomial transforms. More...
class	SipTransformBase
	Base class for SIP transform objects. More...
class	SipForwardTransform
	A transform that maps pixel coordinates to intermediate world coordinates according to the SIP convention. More...
class	SipReverseTransform
	A transform that maps intermediate world coordinates to pixel coordinates according to the SIP convention. More...

Typedefs
typedef std::vector< RecordProxy >	ProxyVector

Functions
template<typename MatchT >
afw::math::Statistics	makeMatchStatistics (std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
	Compute statistics of the distance field of a match list. More...
template<typename MatchT >
afw::math::Statistics	makeMatchStatisticsInPixels (afw::image::Wcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
	Compute statistics of on-detector radial separation for a match list, in pixels. More...
template<typename MatchT >
afw::math::Statistics	makeMatchStatisticsInRadians (afw::image::Wcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
	Compute statistics of on-sky radial separation for a match list, in radians. More...
ProxyVector	makeProxies (lsst::afw::table::SourceCatalog const &sourceCat, afw::image::Wcs const &distortedWcs, afw::image::Wcs const &tanWcs)
ProxyVector	makeProxies (lsst::afw::table::SimpleCatalog const &posRefCat, afw::image::Wcs const &tanWcs)
lsst::afw::table::ReferenceMatchVector	matchOptimisticB (lsst::afw::table::SimpleCatalog const &posRefCat, lsst::afw::table::SourceCatalog const &sourceCat, MatchOptimisticBControl const &control, afw::image::Wcs const &wcs, int posRefBegInd=0, bool verbose=false)
	Match sources to stars in a position reference catalog using optimistic pattern matching B. More...
PolynomialTransform	compose (afw::geom::AffineTransform const &t1, PolynomialTransform const &t2)
	Return a PolynomialTransform that is equivalent to the composition t1(t2()) More...
PolynomialTransform	compose (PolynomialTransform const &t1, afw::geom::AffineTransform const &t2)
	Return a PolynomialTransform that is equivalent to the composition t1(t2()) More...
std::shared_ptr < afw::image::TanWcs >	makeWcs (SipForwardTransform const &sipForward, SipReverseTransform const &sipReverse, afw::coord::Coord const &skyOrigin)
	Create a new TAN SIP Wcs from a pair of SIP transforms and the sky origin. More...
template afw::math::Statistics	makeMatchStatistics< afw::table::ReferenceMatch > (std::vector< afw::table::ReferenceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)
template afw::math::Statistics	makeMatchStatisticsInPixels< afw::table::ReferenceMatch > (afw::image::Wcs const &wcs, std::vector< afw::table::ReferenceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)
template afw::math::Statistics	makeMatchStatisticsInRadians< afw::table::ReferenceMatch > (afw::image::Wcs const &wcs, std::vector< afw::table::ReferenceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)
template afw::math::Statistics	makeMatchStatistics< afw::table::SourceMatch > (std::vector< afw::table::SourceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)
template afw::math::Statistics	makeMatchStatisticsInPixels< afw::table::SourceMatch > (afw::image::Wcs const &wcs, std::vector< afw::table::SourceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)
template afw::math::Statistics	makeMatchStatisticsInRadians< afw::table::SourceMatch > (afw::image::Wcs const &wcs, std::vector< afw::table::SourceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)

Typedef Documentation

typedef std::vector<RecordProxy> lsst::meas::astrom.ProxyVector

Definition at line 53 of file matchOptimisticB.h.

Function Documentation

PolynomialTransform lsst::meas::astrom::compose	(	afw::geom::AffineTransform const &	t1,
		PolynomialTransform const &	t2
	)

Return a PolynomialTransform that is equivalent to the composition t1(t2())

The returned composition would be exact in ideal arithmetic, but may suffer from significant round-off error for high-order polynomials.

Definition at line 254 of file PolynomialTransform.cc.

                                                                                               {
    typedef afw::geom::AffineTransform AT;
    PolynomialTransform result(t2.getOrder());
    result._xCoeffs = t2._xCoeffs*t1[AT::XX] + t2._yCoeffs*t1[AT::XY];
    result._yCoeffs = t2._xCoeffs*t1[AT::YX] + t2._yCoeffs*t1[AT::YY];
    result._xCoeffs(0, 0) += t1[AT::X];
    result._yCoeffs(0, 0) += t1[AT::Y];
    return result;
}

PolynomialTransform lsst::meas::astrom::compose	(	PolynomialTransform const &	t1,
		afw::geom::AffineTransform const &	t2
	)

Return a PolynomialTransform that is equivalent to the composition t1(t2())

The returned composition would be exact in ideal arithmetic, but may suffer from significant round-off error for high-order polynomials.

Definition at line 264 of file PolynomialTransform.cc.

                                                                                               {
    typedef afw::geom::AffineTransform AT;
    int const order = t1.getOrder();
    if (order < 1) {
        PolynomialTransform t1a(1);
        t1a._xCoeffs(0, 0) = t1._xCoeffs(0, 0);
        t1a._yCoeffs(0, 0) = t1._yCoeffs(0, 0);
        return compose(t1a, t2);
    }
    detail::BinomialMatrix binomial(order);
    // For each of these, (e.g.) a[n] == pow(a, n)
    auto const t2u = detail::computePowers(t2[AT::X], order);
    auto const t2v = detail::computePowers(t2[AT::Y], order);
    auto const t2uu = detail::computePowers(t2[AT::XX], order);
    auto const t2uv = detail::computePowers(t2[AT::XY], order);
    auto const t2vu = detail::computePowers(t2[AT::YX], order);
    auto const t2vv = detail::computePowers(t2[AT::YY], order);
    PolynomialTransform result(order);
    for (int p = 0; p <= order; ++p) {
        for (int m = 0; m <= p; ++m) {
            for (int j = 0; j <= m; ++j) {
                for (int q = 0; p + q <= order; ++q) {
                    for (int n = 0; n <= q; ++n) {
                        for (int k = 0; k <= n; ++k) {
                            double z = binomial(p,m) * t2u[p-m] * binomial(m,j) * t2uu[j] * t2uv[m-j] *
                                       binomial(q,n) * t2v[q-n] * binomial(n,k) * t2vu[k] * t2vv[n-k];
                            result._xCoeffs(j + k, m + n - j - k) += t1._xCoeffs(p, q) * z;
                            result._yCoeffs(j + k, m + n - j - k) += t1._yCoeffs(p, q) * z;
                        } // k
                    } // n
                } // q
            } // j
        } // m
    } // p
    return result;
}

template<typename MatchT >

afw::math::Statistics lsst::meas::astrom::makeMatchStatistics	(	std::vector< MatchT > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl = afw::math::StatisticsControl()
	)

Compute statistics of the distance field of a match list.

Parameters

[in]	matchList	list of matchList between reference objects and sources; fields read: distance: distance between source and reference object, in arbitrary units; the resulting statistics have the same units as distance
[in]	flags	what to calculate; OR constants such as lsst::afw::math::MEAN, MEANCLIP, STDDEV, MEDIAN, defined in lsst/afw/math/Statitics.h's Property enum
[in]	sctrl	statistics configuration

Definition at line 33 of file makeMatchStatistics.cc.

  {
    if (matchList.empty()) {
        throw LSST_EXCEPT(pexExcept::RuntimeError, "matchList is empty");
    }
    std::vector<double> val;
    val.reserve(matchList.size());

    for (auto const & match: matchList) {
        val.push_back(match.distance);
    }
    return afw::math::makeStatistics(val, flags, sctrl);
}

template afw::math::Statistics lsst::meas::astrom::makeMatchStatistics< afw::table::ReferenceMatch >	(	std::vector< afw::table::ReferenceMatch > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl
	)

template afw::math::Statistics lsst::meas::astrom::makeMatchStatistics< afw::table::SourceMatch >	(	std::vector< afw::table::SourceMatch > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl
	)

template<typename MatchT >

afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInPixels	(	afw::image::Wcs const &	wcs,
		std::vector< MatchT > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl = afw::math::StatisticsControl()
	)

Compute statistics of on-detector radial separation for a match list, in pixels.

Parameters

[in]	wcs	WCS describing pixel to sky transformation
[in]	matchList	list of matchList between reference objects and sources; fields read: first: reference object; only the coord is read second: source; only the centroid is read
[in]	flags	what to calculate; OR constants such as lsst::afw::math::MEAN, MEANCLIP, STDDEV, MEDIAN, defined in lsst/afw/math/Statitics.h's Property enum
[in]	sctrl	statistics configuration

Definition at line 51 of file makeMatchStatistics.cc.

  {
    if (matchList.empty()) {
        throw LSST_EXCEPT(pexExcept::RuntimeError, "matchList is empty");
    }
    std::vector<double> val;
    val.reserve(matchList.size());

    for (auto const & match: matchList) {
        auto refPtr = match.first;
        auto srcPtr = match.second;
        auto srcX = srcPtr->getX();
        auto srcY = srcPtr->getY();
        auto refPos = wcs.skyToPixel(refPtr->getCoord());
        auto refX = refPos[0];
        auto refY = refPos[1];
        val.push_back(::hypot(srcX - refX, srcY - refY));
    }
    return afw::math::makeStatistics(val, flags, sctrl);
}

template afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInPixels< afw::table::ReferenceMatch >	(	afw::image::Wcs const &	wcs,
		std::vector< afw::table::ReferenceMatch > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl
	)

template afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInPixels< afw::table::SourceMatch >	(	afw::image::Wcs const &	wcs,
		std::vector< afw::table::SourceMatch > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl
	)

template<typename MatchT >

afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInRadians	(	afw::image::Wcs const &	wcs,
		std::vector< MatchT > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl = afw::math::StatisticsControl()
	)

Compute statistics of on-sky radial separation for a match list, in radians.

Parameters

[in]	wcs	WCS describing pixel to sky transformation
[in]	matchList	list of matchList between reference objects and sources; fields read: first: reference object; only the coord is read second: source; only the centroid is read
[in]	flags	what to calculate; OR constants such as lsst::afw::math::MEAN, MEANCLIP, STDDEV, MEDIAN, defined in lsst/afw/math/Statitics.h's Property enum
[in]	sctrl	statistics configuration

Definition at line 77 of file makeMatchStatistics.cc.

  {
    if (matchList.empty()) {
        throw LSST_EXCEPT(pexExcept::RuntimeError, "matchList is empty");
    }
    std::vector<double> val;
    val.reserve(matchList.size());

    for (auto const & match: matchList) {
        auto refPtr = match.first;
        auto srcPtr = match.second;
        auto refCoord = refPtr->getCoord();
        auto srcCoord = wcs.pixelToSky(srcPtr->getCentroid());
        auto angSep = refCoord.angularSeparation(*srcCoord);
        val.push_back(angSep.asRadians());
    }
    return afw::math::makeStatistics(val, flags, sctrl);
}

template afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInRadians< afw::table::ReferenceMatch >	(	afw::image::Wcs const &	wcs,
		std::vector< afw::table::ReferenceMatch > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl
	)

template afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInRadians< afw::table::SourceMatch >	(	afw::image::Wcs const &	wcs,
		std::vector< afw::table::SourceMatch > const &	matchList,
		int const	flags,
		afw::math::StatisticsControl const &	sctrl
	)

ProxyVector lsst::meas::astrom::makeProxies	(	lsst::afw::table::SourceCatalog const &	sourceCat,
		afw::image::Wcs const &	distortedWcs,
		afw::image::Wcs const &	tanWcs
	)

Definition at line 491 of file matchOptimisticB.cc.

    {
        ProxyVector r;
        r.reserve(sourceCat.size());
        for (afwTable::SourceCatalog::const_iterator sourcePtr = sourceCat.begin();
            sourcePtr != sourceCat.end(); ++sourcePtr) {
            r.push_back(RecordProxy(sourcePtr,
                                    tanWcs.skyToPixel(*distortedWcs.pixelToSky(sourcePtr->getCentroid()))));
        }
        return r;
    }

ProxyVector lsst::meas::astrom::makeProxies	(	lsst::afw::table::SimpleCatalog const &	posRefCat,
		afw::image::Wcs const &	tanWcs
	)

Definition at line 506 of file matchOptimisticB.cc.

    {
        auto coordKey = afwTable::CoordKey(posRefCat.getSchema()["coord"]);
        ProxyVector r;
        r.reserve(posRefCat.size());
        for (afwTable::SimpleCatalog::const_iterator posRefPtr = posRefCat.begin();
            posRefPtr != posRefCat.end(); ++posRefPtr) {
            r.push_back(RecordProxy(posRefPtr, tanWcs.skyToPixel(posRefPtr->get(coordKey))));
        }
        return r;
    }

boost::shared_ptr< afw::image::TanWcs > lsst::meas::astrom::makeWcs	(	SipForwardTransform const &	sipForward,
		SipReverseTransform const &	sipReverse,
		afw::coord::Coord const &	skyOrigin
	)

Create a new TAN SIP Wcs from a pair of SIP transforms and the sky origin.

Parameters

[in]	sipForward	Mapping from pixel coordinates to intermediate world coordinates.
[in]	sipReverse	Mapping from intermediate world coordinates to pixel coordinates.
[in]	skyOrigin	Position of the gnomonic projection that maps sky coordinates to intermediate world coordinates (CRVAL).

Exceptions

pex::exceptions::InvalidParameterError

if the forward and reverse SIP transforms have different CRPIX values or CD matrices.

Definition at line 149 of file SipTransform.cc.

  {
    if (!sipForward.getPixelOrigin().asEigen().isApprox(sipReverse.getPixelOrigin().asEigen())) {
        std::ostringstream oss;
        oss << "SIP forward and reverse transforms have inconsistent CRPIX: "
            << sipForward.getPixelOrigin() << " != " << sipReverse.getPixelOrigin();
        throw LSST_EXCEPT(
            pex::exceptions::InvalidParameterError,
            oss.str()
        );
    }
    if (!sipForward.getCDMatrix().getMatrix().isApprox(sipReverse.getCDMatrix().getMatrix())) {
        std::ostringstream oss;
        oss << "SIP forward and reverse transforms have inconsistent CD matrix: "
            << sipForward.getCDMatrix() << "\n!=\n" << sipReverse.getCDMatrix();
        throw LSST_EXCEPT(
            pex::exceptions::InvalidParameterError,
            oss.str()
        );
    }
    Eigen::MatrixXd sipA(sipForward.getPoly().getXCoeffs().asEigen());
    Eigen::MatrixXd sipB(sipForward.getPoly().getYCoeffs().asEigen());
    Eigen::MatrixXd sipAP(sipReverse.getPoly().getXCoeffs().asEigen());
    Eigen::MatrixXd sipBP(sipReverse.getPoly().getYCoeffs().asEigen());
    // TanWcs uses strings for coordinate systems, while Coord uses an enum.
    // Frustratingly, there's no way to convert from the enum to the string.
    std::string coordSys;
    switch (skyOrigin.getCoordSystem()) {
    case afw::coord::ICRS:
        coordSys = "ICRS";
        break;
    case afw::coord::FK5:
        coordSys = "FK5";
        break;
    default:
        throw LSST_EXCEPT(
            pex::exceptions::InvalidParameterError,
            "Coordinate system not supported"
        );
    }
    return std::make_shared<afw::image::TanWcs>(
        skyOrigin.getPosition(afw::geom::degrees),
        sipForward.getPixelOrigin(),
        sipForward.getCDMatrix().getMatrix(),
        sipA, sipB, sipAP, sipBP,
        skyOrigin.getEpoch(),
        coordSys
    );
}

afwTable::ReferenceMatchVector lsst::meas::astrom::matchOptimisticB	(	lsst::afw::table::SimpleCatalog const &	posRefCat,
		lsst::afw::table::SourceCatalog const &	sourceCat,
		MatchOptimisticBControl const &	control,
		afw::image::Wcs const &	wcs,
		int	posRefBegInd = 0,
		bool	verbose = false
	)

Match sources to stars in a position reference catalog using optimistic pattern matching B.

Optimistic Pattern Matching is described in V. Tabur 2007, PASA, 24, 189-198 "Fast Algorithms for Matching CCD Images to a Stellar Catalogue"

Parameters

[in]	posRefCat	catalog of position reference stars; fields read: "coord" control.refFluxField
[in]	sourceCat	catalog of detected sources; fields read: "Centroid_x" "Centroid_y" control.refFluxField
[in]	wcs	estimated WCS
[in]	control	control object
[in]	posRefBegInd	index of first start to use in posRefCat
[in]	verbose	true to print diagnostic information to std::cout

Returns

a list of matches; the d field may be set, but you should not rely on it

Definition at line 520 of file matchOptimisticB.cc.

      {
        control.validate();
        if (posRefCat.empty()) {
            throw LSST_EXCEPT(pexExcept::InvalidParameterError, "no entries in posRefCat");
        }
        if (sourceCat.empty()) {
            throw LSST_EXCEPT(pexExcept::InvalidParameterError, "no entries in sourceCat");
        }
        if (posRefBegInd < 0) {
            throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd < 0");
        }
        if (static_cast<size_t>(posRefBegInd) >= posRefCat.size()) {
            throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd too big");
        }
        double const maxRotationRad = afw::geom::degToRad(control.maxRotationDeg);

        CONST_PTR(afw::image::Wcs) tanWcs; // Undistorted Wcs, providing tangent plane
        if (wcs.hasDistortion()) {
            // Create an undistorted Wcs to project everything with
            // We'll anchor it at the center of the area.
            afw::geom::Extent2D srcCenter(0, 0);
            for (auto iter = sourceCat.begin(); iter != sourceCat.end(); ++iter) {
                srcCenter += afw::geom::Extent2D(iter->getCentroid());
            }
            srcCenter /= sourceCat.size();

            afw::geom::Extent3D refCenter(0, 0, 0);
            for (auto iter = posRefCat.begin(); iter != posRefCat.end(); ++iter) {
                refCenter += afw::geom::Extent3D(iter->getCoord().toIcrs().getVector());
            }
            refCenter /= posRefCat.size();

            tanWcs = std::make_shared<afw::image::Wcs>(
                afw::coord::IcrsCoord(afw::geom::Point3D(refCenter)).getPosition(),
                afw::geom::Point2D(srcCenter),
                wcs.getCDMatrix()
                );
        }

        ProxyVector posRefProxyCat = makeProxies(posRefCat, tanWcs ? *tanWcs : wcs);
        ProxyVector sourceProxyCat = makeProxies(sourceCat, wcs, tanWcs ? *tanWcs : wcs);

        // sourceSubCat contains at most the numBrightStars brightest sources, sorted by decreasing flux
        ProxyVector sourceSubCat = selectPoint(
            sourceProxyCat,
            sourceCat.getSchema().find<double>(control.sourceFluxField).key,
            control.numBrightStars);

        // posRefSubCat skips the initial posRefBegInd brightest reference objects and contains
        // at most the next len(sourceSubCat) + 25 brightest reference objects, sorted by decreasing flux
        ProxyVector posRefSubCat = selectPoint(
            posRefProxyCat,
            posRefCat.getSchema().find<double>(control.refFluxField).key,
            sourceSubCat.size()+25,
            posRefBegInd);
        if (verbose) {
            std::cout << "Catalog sizes: " << sourceSubCat.size() << " " << posRefSubCat.size() << std::endl;
        }

        // Construct a list of pairs of position reference stars sorted by increasing separation
        std::vector<ProxyPair> posRefPairList;
        size_t const posRefCatSubSize = posRefSubCat.size();
        for (size_t i = 0; i < posRefCatSubSize-1; i++) {
            for (size_t j = i+1; j < posRefCatSubSize; j++) {
                posRefPairList.push_back(ProxyPair(posRefSubCat[i], posRefSubCat[j]));
            }
        }
        std::sort(posRefPairList.begin(), posRefPairList.end(), cmpPair);

        // Construct a list of pairs of sources sorted by increasing separation
        std::vector<ProxyPair> sourcePairList;
        size_t const sourceSubCatSize = sourceSubCat.size();
        for (size_t i = 0; i < sourceSubCatSize-1; i++) {
            for (size_t j = i+1; j < sourceSubCatSize; j++) {
                sourcePairList.push_back(ProxyPair(sourceSubCat[i], sourceSubCat[j]));
            }
        }
        std::sort(sourcePairList.begin(), sourcePairList.end(), cmpPair);

        afwTable::ReferenceMatchVector matPair;
        afwTable::ReferenceMatchVector matPairSave;
        std::vector<afwTable::ReferenceMatchVector> matPairCand;

        size_t const fullShapeSize = control.numPointsForShape - 1;   // Max size of shape array
        for (size_t ii = 0; ii < sourcePairList.size(); ii++) {
            ProxyPair p = sourcePairList[ii];

            std::vector<ProxyPair> q = searchPair(posRefPairList, p,
                control.matchingAllowancePix, maxRotationRad);

            // If candidate pairs are found
            if (q.size() != 0) {

                std::vector<ProxyPair> srcMatPair;
                std::vector<ProxyPair> catMatPair;

                // Go through candidate pairs
                for (size_t l = 0; l < q.size(); l++) {
                    // sign matters, so don't use deltaAng; no need to wrap because
                    // the result is used with deltaAng later
                    double dpa = p.pa - q[l].pa;

                    srcMatPair.clear();
                    catMatPair.clear();

                    srcMatPair.push_back(p);
                    catMatPair.push_back(q[l]);

                    if (verbose) {
                        std::cout << "p dist: " << p.distance << " pa: " << p.pa << std::endl;
                        std::cout << "q dist: " << q[l].distance << " pa: " << q[l].pa << std::endl;
                    }

                    for (size_t k = 0; k < sourceSubCat.size(); k++) {
                        if (p.first == sourceSubCat[k] || p.second == sourceSubCat[k]) continue;

                        ProxyPair pp(p.first, sourceSubCat[k]);

                        std::vector<ProxyPair>::iterator r = searchPair3(posRefPairList, pp, q[l],
                            control.matchingAllowancePix, dpa, maxRotationRad);
                        if (r != posRefPairList.end()) {
                            srcMatPair.push_back(pp);
                            catMatPair.push_back(*r);
                            if (verbose) {
                                std::cout << "  p dist: " << pp.distance << " pa: " << pp.pa << std::endl;
                                std::cout << "  r dist: " << (*r).distance << " pa: " << (*r).pa << std::endl;
                            }
                            if (srcMatPair.size() == fullShapeSize) {
                                break;
                            }
                        }
                    }

                    bool goodMatch = false;
                    if (srcMatPair.size() == fullShapeSize) {
                        goodMatch = true;
                        for (size_t k = 1; k < catMatPair.size(); k++) {
                            if (catMatPair[0].first != catMatPair[k].first) {
                                goodMatch = false;
                                break;
                            }
                        }
                    }

                    if (goodMatch && srcMatPair.size() == fullShapeSize) {

                        ProxyVector srcMat;
                        ProxyVector catMat;

                        srcMat.push_back(srcMatPair[0].first);
                        catMat.push_back(catMatPair[0].first);
                        for (size_t k = 0; k < srcMatPair.size(); k++) {
                            srcMat.push_back(srcMatPair[k].second);
                            catMat.push_back(catMatPair[k].second);
                        }

                        boost::shared_array<double> coeff = polyfit(1, srcMat, catMat);

                        if (verbose) {
                            for (size_t k = 0; k < srcMat.size(); k++) {
                                std::cout << "circle(" << srcMat[k].getX() << ","
                                          << srcMat[k].getY() << ",10) # color=green" << std::endl;
                                std::cout << "circle(" << catMat[k].getX() << ","
                                          << catMat[k].getY() << ",10) # color=red" << std::endl;
                                std::cout << "line(" << srcMat[0].getX() << "," << srcMat[0].getY() << ","
                                          << srcMat[k].getX() << "," << srcMat[k].getY()
                                          << ") # line=0 0 color=green" << std::endl;
                                std::cout << "line(" << catMat[0].getX() << "," << catMat[0].getY() << ","
                                          << catMat[k].getX() << "," << catMat[k].getY()
                                          << ") # line=0 0 color=red" << std::endl;
                            }
                        }

                        double a = coeff[1];
                        double b = coeff[2];
                        double c = coeff[4];
                        double d = coeff[5];
                        afw::geom::Angle const theta(std::acos((a*b+c*d)/
                                                               (std::sqrt(a*a+c*c)*std::sqrt(b*b+d*d))),
                                                     afw::geom::radians);
                        if (verbose) {
                            std::cout << "Linear fit from match:" << std::endl;
                            std::cout << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
                            std::cout << coeff[3] << " " << coeff[4] << " " << coeff[5] << std::endl;
                            std::cout << "Determinant (max " << control.maxDeterminant << "): ";
                            std::cout << coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1. << std::endl;
                            std::cout << "Angle between axes (deg; allowed 90 +/- ";
                            std::cout << control.allowedNonperpDeg << "): ";
                            std::cout << theta.asDegrees() << std::endl;
                        }
                        if (std::fabs(coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1.)
                                > control.maxDeterminant ||
                            std::fabs(theta.asDegrees() - 90) > control.allowedNonperpDeg ||
                            std::fabs(coeff[0]) > control.maxOffsetPix ||
                            std::fabs(coeff[3]) > control.maxOffsetPix) {
                            if (verbose) {
                                std::cout << "Bad; continuing" << std::endl;
                            }
                            continue;
                        } else {
                            double x0, y0, x1, y1;
                            int num = 0;
                            srcMat.clear();
                            catMat.clear();
                            for (size_t i = 0; i < sourceSubCat.size(); i++) {
                                x0 = sourceSubCat[i].getX();
                                y0 = sourceSubCat[i].getY();
                                transform(1, coeff, x0, y0, &x1, &y1);
                                auto refObjDist = searchNearestPoint(posRefSubCat, x1, y1,
                                    control.matchingAllowancePix);
                                if (refObjDist.first != posRefSubCat.end()) {
                                    num++;
                                    srcMat.push_back(sourceSubCat[i]);
                                    catMat.push_back(*refObjDist.first);
                                    if (verbose) {
                                        std::cout << "Match: " << x0 << "," << y0 << " --> "
                                            << x1 << "," << y1 << " <==> "
                                            << refObjDist.first->getX() << "," << refObjDist.first->getY()
                                            << std::endl;
                                    }
                                }
                            }
                            if (num <= control.numPointsForShape) {
                                // Can get matrix = 0,0,0,0; everything matches a single catalog object
                                if (verbose) {
                                    std::cout << "Insufficient initial matches; continuing" << std::endl;
                                }
                                continue;
                            }
                            coeff = polyfit(1, srcMat, catMat);
                            if (verbose) {
                                std::cout << "Coefficients from initial matching:" << std::endl;
                                for (size_t i = 0; i < 6; ++i) {
                                    std::cout << coeff[i] << " ";
                                }
                                std::cout << std::endl;
                            }

                            matPair = FinalVerify(coeff, posRefProxyCat, sourceProxyCat,
                                control.matchingAllowancePix, verbose);
                            if (verbose) {
                                std::cout << "Number of matches: " << matPair.size() << " vs " <<
                                    control.minMatchedPairs << std::endl;
                            }
                            if (matPair.size() <= static_cast<std::size_t>(control.minMatchedPairs)) {
                                if (verbose) {
                                    std::cout << "Insufficient final matches; continuing" << std::endl;
                                }
                                if (matPair.size() > matPairSave.size()) {
                                    matPairSave = matPair;
                                }
                                continue;
                            } else {
                                if (verbose) {
                                    std::cout << "Finish" << std::endl;
                                }
                                matPairCand.push_back(matPair);
                                if (matPairCand.size() == 3) {
                                    goto END;
                                }
                            }
                        }
                    }
                }
            }
        }

     END:
        if (matPairCand.size() == 0) {
            return matPairSave;
        } else {
            size_t nmatch = matPairCand[0].size();
            afwTable::ReferenceMatchVector matPairRet = matPairCand[0];
            for (size_t i = 1; i < matPairCand.size(); i++) {
                if (matPairCand[i].size() > nmatch) {
                    nmatch = matPairCand[i].size();
                    matPairRet = matPairCand[i];
                }
            }
            return matPairRet;
        }
    }