lsst::meas::astrom Namespace Reference

Namespaces
namespace	approximateWcs
namespace	astrometry
namespace	denormalizeMatches
namespace	detail
namespace	directMatch
namespace	display
namespace	fitAffineWcs
namespace	fitSipDistortion
namespace	fitTanSipWcs
namespace	match_probabilistic_task
namespace	matcher_probabilistic
namespace	matchOptimisticB
namespace	matchOptimisticBTask
namespace	matchPessimisticB
namespace	pessimistic_pattern_matcher_b_3D
namespace	ref_match
namespace	setMatchDistance
namespace	sip
namespace	verifyWcs
namespace	version

Classes
struct	MatchOptimisticBControl
class	OutlierRejectionControl
	Control object for outlier rejection in ScaledPolynomialTransformFitter. More...
struct	PatternResult
	Result of construct_pattern_and_shift_rot_matrix(), containing the matched sources, rotation matrix, and a success flag. More...
class	PolynomialTransform
	A 2-d coordinate transform represented by a pair of standard polynomials (one for each coordinate). More...
struct	ProxyPair
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	RotationTestResult
	Result of test_rotation(), containing the rotation matrix and success flag. More...
class	ScaledPolynomialTransform
	A 2-d coordinate transform represented by a lazy composition of an AffineTransform, a PolynomialTransform, and another AffineTransform. More...
class	ScaledPolynomialTransformFitter
	A fitter class for scaled polynomial transforms. More...
struct	ShiftRotMatrixResult
	Result of create_shift_rot_matrix() containing the rotation matrix and rotation angle. 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...
class	SipTransformBase
	Base class for SIP transform objects. More...
struct	SortedArrayResult
	Result of create_sorted_arrays(), containing the sorted distances and array indexes. 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.
template<typename MatchT >
afw::math::Statistics	makeMatchStatisticsInPixels (afw::geom::SkyWcs 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.
template<typename MatchT >
afw::math::Statistics	makeMatchStatisticsInRadians (afw::geom::SkyWcs 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.
ProxyVector	makeProxies (afw::table::SourceCatalog const &sourceCat, afw::geom::SkyWcs const &distortedWcs, afw::geom::SkyWcs const &tanWcs)
ProxyVector	makeProxies (afw::table::SimpleCatalog const &posRefCat, afw::geom::SkyWcs const &tanWcs)
afw::table::ReferenceMatchVector	matchOptimisticB (afw::table::SimpleCatalog const &posRefCat, afw::table::SourceCatalog const &sourceCat, MatchOptimisticBControl const &control, afw::geom::SkyWcs const &wcs, int posRefBegInd=0, bool verbose=false)
	Match sources to stars in a position reference catalog using optimistic pattern matching B.
PatternResult	construct_pattern_and_shift_rot_matrix (ndarray::Array< double, 2, 1 > src_pattern_array, ndarray::Array< double, 2, 1 > src_delta_array, ndarray::Array< double, 1, 1 > src_dist_array, ndarray::Array< float, 1, 1 > dist_array, ndarray::Array< uint16_t, 2, 1 > id_array, ndarray::Array< double, 2, 1 > reference_array, size_t n_match, double max_cos_theta_shift, double max_cos_rot_sq, double max_dist_rad)
	Test an input source pattern against the reference catalog.
RotationTestResult	test_rotation (ndarray::Array< double, 1, 1 > const &src_center, ndarray::Array< double, 1, 1 > const &ref_center, ndarray::Array< double, 1, 1 > const &src_delta, ndarray::Array< double, 1, 1 > const &ref_delta, double cos_shift, double max_cos_rot_sq)
	Test if the rotation implied between the source pattern and reference pattern is within tolerance.
SortedArrayResult	create_sorted_arrays (ndarray::Array< double, 1, 1 > const &ref_center, ndarray::Array< double, 2, 1 > const &reference_array)
	Create arrays sorted on the distance between the center of this candidate reference object and all the reference objects.
std::vector< std::pair< size_t, size_t > >	create_pattern_spokes (ndarray::Array< double, 1, 1 > const &src_ctr, ndarray::Array< double, 2, 1 > const &src_delta_array, ndarray::Array< double, 1, 1 > const &src_dist_array, ndarray::Array< double, 1, 1 > const &ref_ctr, Eigen::Vector3d const &proj_ref_ctr_delta, std::vector< float > const &ref_dist_array, std::vector< uint16_t > const &ref_id_array, ndarray::Array< double, 2, 1 > const &reference_array, double max_dist_rad, size_t n_match)
	Create the individual spokes that make up the pattern now that the shift and rotation are within tolerance.
ShiftRotMatrixResult	create_shift_rot_matrix (double cos_rot_sq, Eigen::Matrix3d const &shift_matrix, ndarray::Array< double, 1, 1 > const &src_delta, ndarray::Array< double, 1, 1 > const &ref_ctr, ndarray::Array< double, 1, 1 > const &ref_delta)
	Create the complete spherical rotation matrix.
Eigen::Matrix3d	create_spherical_rotation_matrix (Eigen::Vector3d const &rot_axis, double cos_rotation, double sin_rotion)
	Construct a generalized 3D rotation matrix about a given axis.
bool	intermediate_verify_comparison (std::vector< Eigen::Vector3d > const &src_pattern, std::vector< Eigen::Vector3d > const &ref_pattern, Eigen::Matrix3d const &shift_rot_matrix, double max_dist_rad)
	Test the input rotation matrix by comparing the rotated src pattern to the ref pattern.
int	check_spoke (double cos_theta_src, double sin_theta_src, ndarray::Array< double, 1, 1 > const &ref_ctr, Eigen::Vector3d const &proj_ref_ctr_delta, double proj_ref_ctr_dist_sq, std::pair< size_t, size_t > const &candidate_range, std::vector< uint16_t > const &ref_id_array, ndarray::Array< double, 2, 1 > const &reference_array, double src_sin_tol)
	Check the opening angle between the first spoke of our pattern for the source object against the reference object.
PolynomialTransform	compose (geom::AffineTransform const &t1, PolynomialTransform const &t2)
	Return a PolynomialTransform that is equivalent to the composition t1(t2())
PolynomialTransform	compose (PolynomialTransform const &t1, geom::AffineTransform const &t2)
	Return a PolynomialTransform that is equivalent to the composition t1(t2())
std::shared_ptr< afw::geom::SkyWcs >	makeWcs (SipForwardTransform const &sipForward, SipReverseTransform const &sipReverse, geom::SpherePoint const &skyOrigin)
	Create a new TAN SIP Wcs from a pair of SIP transforms and the sky origin.
std::shared_ptr< afw::geom::SkyWcs >	transformWcsPixels (afw::geom::SkyWcs const &wcs, geom::AffineTransform const &s)
	Create a new SkyWcs whose pixel coordinate system has been transformed via an affine transform.
std::shared_ptr< afw::geom::SkyWcs >	rotateWcsPixelsBy90 (afw::geom::SkyWcs const &wcs, int nQuarter, geom::Extent2I const &dimensions)
	Return a new SkyWcs that represents a rotation of the image it corresponds to about the image's center.
void	wrapPolynomialTransform (WrapperCollection &wrappers)
void	wrapScaledPolynomialTransformFitter (WrapperCollection &wrappers)
void	wrapSipTransform (WrapperCollection &wrappers)
void	wrapMatchOptimisticB (WrapperCollection &wrappers)
void	wrapMakeMatchStatistics (WrapperCollection &wrappers)
void	wrapPessimisticPatternMatcherUtils (WrapperCollection &wrappers)
	PYBIND11_MODULE (_measAstromLib, mod)
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::geom::SkyWcs 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::geom::SkyWcs 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::geom::SkyWcs 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::geom::SkyWcs const &wcs, std::vector< afw::table::SourceMatch > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl)
std::pair< size_t, size_t >	find_candidate_reference_pair_range (float src_dist, ndarray::Array< float, 1, 1 > const &ref_dist_array, double max_dist_rad)
	Find the range of reference pairs within the distance tolerance of our source pair spoke.
std::pair< size_t, size_t >	find_candidate_reference_pair_range (float src_dist, std::vector< float > const &ref_dist_array, double max_dist_rad)
	Find the range of reference pairs within the distance tolerance of our source pair spoke.

Typedef Documentation

ProxyVector

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

Definition at line 52 of file matchOptimisticB.h.

Function Documentation

check_spoke()

int lsst::meas::astrom::check_spoke	(	double	cos_theta_src,
		double	sin_theta_src,
		ndarray::Array< double, 1, 1 > const &	ref_ctr,
		Eigen::Vector3d const &	proj_ref_ctr_delta,
		double	proj_ref_ctr_dist_sq,
		std::pair< size_t, size_t > const &	candidate_range,
		std::vector< uint16_t > const &	ref_id_array,
		ndarray::Array< double, 2, 1 > const &	reference_array,
		double	src_sin_tol )

Check the opening angle between the first spoke of our pattern for the source object against the reference object.

This method makes heavy use of the small angle approximation to perform the comparison.

Parameters

[in]	cos_theta_src	Cosine of the angle between the current candidate source spoke and the first spoke.
[in]	sin_theta_src	Sine of the angle between the current candidate source spoke and the first spoke.
[in]	ref_ctr	3 vector of the candidate reference center
[in]	proj_ref_ctr_delta	Plane projected first spoke in the reference pattern using the pattern center as normal.
[in]	proj_ref_ctr_dist_sq	Squared length of the projected vector.
[in]	candidate_range	Min and max index locations in ref_id_array that have pair lengths within the tolerance range.
[in]	ref_id_array	Array of id lookups into the master reference array that our center id object is paired with. uint16 type is to limit the amount of memory used and is set in the pessimistic_pattern_matcher_3d python class with reference catalogs trimmed by the matchPessimisticB runner class.
[in]	reference_array	Array of three vectors representing the locations of all reference objects.
[in]	src_sin_tol	Sine of tolerance allowed between source and reference spoke opening angles.

Returns

ID of the candidate reference object successfully matched or -1 if no match is found.

Definition at line 368 of file pessimisticPatternMatcherUtils.cc.

                                                                                       {
    // Loop over our candidate reference objects. candidate_range is the min
    // and max of for pair candidates and are view into ref_id_array. Here we
    // start from the midpoint of min and max values and step outward.
    size_t midpoint = (candidate_range.first + candidate_range.second) / 2;
    for (size_t idx = 0; idx < candidate_range.second - candidate_range.first; idx++) {
        // TODO DM-33514: cleanup this loop to use an iterator that handles the "inside-out" iteration.
        if (idx % 2 == 0) {
            midpoint = midpoint + idx;
        } else {
            midpoint = midpoint - idx;
        }
        // Compute the delta vector from the pattern center.
        uint16_t ref_id = ref_id_array[midpoint];
        ndarray::Array<double, 1, 1> ref_delta = copy(reference_array[ref_id] - ref_ctr);
 
        double ref_dot = ndarray::asEigenMatrix(ref_delta).dot(ndarray::asEigenMatrix(ref_ctr));
        ndarray::Array<double, 1, 1> proj_ref_delta = copy(ref_delta - ref_dot * ref_ctr);
        // Compute the cos between our "center" reference vector and the
        // current reference candidate.
        auto proj_ref_delta_eigen = ndarray::asEigenMatrix(proj_ref_delta);
        double proj_delta_dist_sq = proj_ref_delta_eigen.dot(proj_ref_delta_eigen);
        double geom_dist_ref = sqrt(proj_ref_ctr_dist_sq * proj_delta_dist_sq);
        double cos_theta_ref = proj_ref_delta_eigen.dot(proj_ref_ctr_delta) / geom_dist_ref;
 
        // Make sure we can safely make the comparison in case
        // our "center" and candidate vectors are mostly aligned.
        double cos_sq_comparison;
        if (cos_theta_ref * cos_theta_ref < (1 - src_sin_tol * src_sin_tol)) {
            cos_sq_comparison = (cos_theta_src - cos_theta_ref) * (cos_theta_src - cos_theta_ref) /
                                (1 - cos_theta_ref * cos_theta_ref);
        } else {
            cos_sq_comparison = (cos_theta_src - cos_theta_ref) * (cos_theta_src - cos_theta_ref) /
                                (src_sin_tol * src_sin_tol);
        }
        // Test the difference of the cosine of the reference angle against
        // the source angle. Assumes that the delta between the two is
        // small.
        if (cos_sq_comparison > src_sin_tol * src_sin_tol) {
            continue;
        }
        // The cosine tests the magnitude of the angle but not
        // its direction. To do that we need to know the sine as well.
        // This cross product calculation does that.
        Eigen::Vector3d cross_ref = proj_ref_delta_eigen.head<3>().cross(proj_ref_ctr_delta) / geom_dist_ref;
        double sin_theta_ref = cross_ref.dot(ndarray::asEigenMatrix(ref_ctr));
        // Check the value of the cos again to make sure that it is not
        // near zero.
        double sin_comparison;
        if (abs(cos_theta_src) < src_sin_tol) {
            sin_comparison = (sin_theta_src - sin_theta_ref) / src_sin_tol;
        } else {
            sin_comparison = (sin_theta_src - sin_theta_ref) / cos_theta_ref;
        }
 
        // Return the correct id of the candidate we found.
        if (abs(sin_comparison) < src_sin_tol) {
            return ref_id;
        }
    }
    return -1;
}

compose() [1/2]

PolynomialTransform lsst::meas::astrom::compose	(	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 214 of file PolynomialTransform.cc.

                                                                                          {
    typedef geom::AffineTransform AT;
    PolynomialTransform result(t2.getOrder());
 
    result._xCoeffs.deep() = t2._xCoeffs * t1[AT::XX] + t2._yCoeffs * t1[AT::XY];
    result._yCoeffs.deep() = 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;
}

compose() [2/2]

PolynomialTransform lsst::meas::astrom::compose	(	PolynomialTransform const &	t1,
		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 225 of file PolynomialTransform.cc.

                                                                                          {
    typedef 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;
}

construct_pattern_and_shift_rot_matrix()

PatternResult lsst::meas::astrom::construct_pattern_and_shift_rot_matrix	(	ndarray::Array< double, 2, 1 >	src_pattern_array,
		ndarray::Array< double, 2, 1 >	src_delta_array,
		ndarray::Array< double, 1, 1 >	src_dist_array,
		ndarray::Array< float, 1, 1 >	dist_array,
		ndarray::Array< uint16_t, 2, 1 >	id_array,
		ndarray::Array< double, 2, 1 >	reference_array,
		size_t	n_match,
		double	max_cos_theta_shift,
		double	max_cos_rot_sq,
		double	max_dist_rad )

Test an input source pattern against the reference catalog.

Parameters

[in]	src_pattern_array	Sub selection of source 3 vectors to create a pattern from.
[in]	src_delta_array	Deltas of pairs from the source pattern center.
[in]	src_dist_array	Distances of the pairs in src_delta_array.
[in]	dist_array	Set of all distances between pairs of reference objects.
[in]	id_array	Pairs of indices creating the reference pairs. Sorted the same as dist_array.
[in]	reference_array	Set of all reference points.
[in]	n_match	Number of points to attempt to create a pattern from. Must be <= len(src_pattern_array)
[in]	max_cos_theta_shift	Maximum shift allowed between two patterns' centers. Equivalent to the shift on the sky.
[in]	max_cos_rot_sq	Maximum rotation between two patterns that have been shifted to have their centers on top of each other.
[in]	max_dist_rad	Maximum difference allowed between the source and reference pair distances to consider the reference pair a candidate for the source pair. Also sets the tolerance between the opening angles of the spokes when compared to the reference.

Returns

The result of the pattern matching.

Definition at line 40 of file pessimisticPatternMatcherUtils.cc.

                                                                                {
    // Our first test. We search the reference dataset for pairs that have the same length as our first source
    // pairs to within plus/minus the max_dist tolerance.
    std::pair<size_t, size_t> candidate_range =
            find_candidate_reference_pair_range(src_dist_array[0], dist_array, max_dist_rad);
    size_t ref_dist_idx = (candidate_range.first + candidate_range.second) / 2;
 
    // Start our loop over the candidate reference objects. Looping from the inside (minimum difference to our
    // source dist) to the outside.
    for (size_t idx = 0; idx < candidate_range.second - candidate_range.first; idx++) {
        // TODO DM-33514: cleanup this loop to use an iterator that handles the "inside-out" iteration.
        if (idx % 2 == 0) {
            ref_dist_idx = ref_dist_idx + idx;
        } else {
            ref_dist_idx = ref_dist_idx - idx;
        }
        // We have two candidates for which reference object corresponds with the source at the center of our
        // pattern. As such we loop over and test both possibilities.
        ndarray::Array<uint16_t, 1, 1> tmp_ref_pair_list = id_array[ref_dist_idx];
        for (uint16_t ref_pair_idx = 0; ref_pair_idx < 2; ref_pair_idx++) {
            uint16_t ref_id = id_array[ref_dist_idx][ref_pair_idx];
            std::vector<std::pair<uint16_t, uint16_t>> candidate_pairs;
            // Test the angle between our candidate ref center and the source center of our pattern. This
            // angular distance also defines the shift we will later use.
            ndarray::Array<double, 1, 1> ref_center = reference_array[ref_id];
            double cos_shift =
                    ndarray::asEigenMatrix(src_pattern_array[0]).dot(ndarray::asEigenMatrix(ref_center));
            if (cos_shift < max_cos_theta_shift) {
                continue;
            }
            // We can now append this one as a candidate.
            candidate_pairs.push_back(std::make_pair(ref_id, 0));
            ndarray::Array<double, 1, 1> ref_delta;
            // Test to see which reference object to use in the pair.
            if (ref_id == *tmp_ref_pair_list.begin()) {
                candidate_pairs.push_back(std::make_pair(tmp_ref_pair_list[1], 1));
                ref_delta = copy(reference_array[tmp_ref_pair_list[1]] - ref_center);
            } else {
                candidate_pairs.push_back(std::make_pair(tmp_ref_pair_list[0], 1));
                ref_delta = copy(reference_array[tmp_ref_pair_list[0]] - ref_center);
            }
            // For dense fields it will be faster to compute the absolute rotation this pair suggests first
            // rather than saving it after all the spokes are found. We then compute the cos^2 of the rotation
            // and first part of the rotation matrix from source to reference frame.
            RotationTestResult test_rot_result =
                    test_rotation(src_pattern_array[0], ref_center, src_delta_array[0], ref_delta, cos_shift,
                                  max_cos_rot_sq);
            if (!test_rot_result.success) {
                continue;
            }
            // Now that we have a candidate first spoke and reference pattern center, we mask our future
            // search to only those pairs that contain our candidate reference center.
            SortedArrayResult sorted_array_struct = create_sorted_arrays(ref_center, reference_array);
            // Now we feed this sub data to match the spokes of our pattern.
            std::vector<std::pair<size_t, size_t>> pattern_spokes =
                    create_pattern_spokes(src_pattern_array[0], src_delta_array, src_dist_array, ref_center,
                                          test_rot_result.proj_ref_ctr_delta, sorted_array_struct.dists,
                                          sorted_array_struct.ids, reference_array, max_dist_rad, n_match);
            // If we don't find enough candidates we can continue to the next reference center pair.
            if (pattern_spokes.size() < n_match - 2) {
                continue;
            }
            // If we have the right number of matched ids we store these.
            candidate_pairs.reserve(candidate_pairs.size() + pattern_spokes.size());
            candidate_pairs.insert(candidate_pairs.end(), pattern_spokes.begin(), pattern_spokes.end());
            // We can now create our full matrix for both the shift and rotation. The shift aligns
            // the pattern centers, while the rotation rotates the spokes on top of each other.
            ShiftRotMatrixResult shift_rot_result =
                    create_shift_rot_matrix(test_rot_result.cos_rot_sq, test_rot_result.shift_matrix,
                                            src_delta_array[0], ref_center, ref_delta);
 
            // concatenate our final patterns.
            Eigen::Matrix3d shift_rot_matrix = shift_rot_result.shift_rot_matrix;
            std::vector<Eigen::Vector3d> ref_pattern;
            std::vector<Eigen::Vector3d> src_pattern;
            ref_pattern.reserve(n_match);
            src_pattern.reserve(n_match);
            for (size_t idx = 0; idx < n_match; idx++) {
                ref_pattern[idx] = ndarray::asEigenMatrix(reference_array[candidate_pairs[idx].first]);
                src_pattern[idx] = ndarray::asEigenMatrix(src_pattern_array[candidate_pairs[idx].second]);
            }
            // TODO: DM-32985 Implement least squares linear fitting here. Look at python example linked in
            // ticket for how to implement.
            // Test that the all points in each pattern are within tolerance after shift/rotation.
            bool passed =
                    intermediate_verify_comparison(src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad);
            if (passed) {
                return PatternResult(candidate_pairs, shift_rot_matrix, cos_shift, shift_rot_result.sin_rot);
            }
        }
    }
    // failed fit
    return PatternResult();
}

create_pattern_spokes()

std::vector< std::pair< size_t, size_t > > lsst::meas::astrom::create_pattern_spokes	(	ndarray::Array< double, 1, 1 > const &	src_ctr,
		ndarray::Array< double, 2, 1 > const &	src_delta_array,
		ndarray::Array< double, 1, 1 > const &	src_dist_array,
		ndarray::Array< double, 1, 1 > const &	ref_ctr,
		Eigen::Vector3d const &	proj_ref_ctr_delta,
		std::vector< float > const &	ref_dist_array,
		std::vector< uint16_t > const &	ref_id_array,
		ndarray::Array< double, 2, 1 > const &	reference_array,
		double	max_dist_rad,
		size_t	n_match )

Create the individual spokes that make up the pattern now that the shift and rotation are within tolerance.

If we can't create a valid pattern we exit early with a partial result.

Parameters

[in]	src_ctr	3 vector of the source pinwheel center
[in]	src_delta_array	Array of 3 vector deltas between the source center and the pairs that make up the remaining spokes of the pinwheel.
[in]	src_dist_array	Array of the distances of each src_delta in the pinwheel.
[in]	ref_ctr	3 vector of the candidate reference center.
[in]	proj_ref_ctr_delta	Plane projected 3 vector formed from the center point of the candidate pin-wheel and the second point in the pattern to create the first spoke pair. This is the candidate pair that was matched in the main _construct_pattern_and_shift_rot_matrix loop.
[in]	ref_dist_array	Array of vector distances for each of the reference pairs
[in]	ref_id_array	Array of id lookups into the master reference array that our center id object is paired with.
[in]	reference_array	Full 3 vector data for the reference catalog.
[in]	max_dist_rad	Maximum search radius for distances.
[in]	n_match	Number of source deltas that must be matched into the reference deltas in order to consider this a successful pattern match.

Returns

Return pairs of reference ids and their matched src ids.

Append the successful indices to our list. The src_idx needs an extra iteration to skip the first and second source objects.

Append the successful indices to our list. The src_idx needs an extra iteration to skip the first and second source objects.

Definition at line 278 of file pessimisticPatternMatcherUtils.cc.

                                             {
    // Struct where we will be putting our results.
    std::vector<std::pair<size_t, size_t>> output_spokes;
 
    // Counter for number of spokes we failed to find a reference
    // candidate for. We break the loop if we haven't found enough.
    size_t n_fail = 0;
 
    // Plane project the center/first spoke of the source pattern using
    // the center vector of the pattern as normal.
    auto src_ctr_eigen = ndarray::asEigenMatrix(src_ctr);
    double src_delta_ctr_dot = ndarray::asEigenMatrix(src_delta_array[0]).dot(src_ctr_eigen);
 
    ndarray::Array<double, 1, 1> proj_src_ctr_delta = copy(src_delta_array[0] - src_delta_ctr_dot * src_ctr);
    auto proj_src_ctr_delta_eigen = ndarray::asEigenMatrix(proj_src_ctr_delta);
    double proj_src_ctr_dist_sq = proj_src_ctr_delta_eigen.dot(proj_src_ctr_delta_eigen);
 
    // Pre - compute the squared length of the projected reference vector.
    double proj_ref_ctr_dist_sq = proj_ref_ctr_delta.dot(proj_ref_ctr_delta);
 
    // Value of sin where sin(theta) ~= theta to within 0.1%. Used to make
    // sure we are still in small angle approximation.
    double max_sin_tol = 0.0447;
    // Loop over the source pairs.
    for (size_t src_idx = 1; src_idx < src_dist_array.size(); src_idx++) {
        if (n_fail > src_dist_array.size() - (n_match - 1)) {
            break;
        }
        // Find the reference pairs that include our candidate pattern center
        // and sort them in increasing delta. Check this first so we don't
        // compute anything else if no candidates exist.
        std::pair<size_t, size_t> candidate_range =
                find_candidate_reference_pair_range(src_dist_array[src_idx], ref_dist_array, max_dist_rad);
        if (candidate_range.first == candidate_range.second) {
            n_fail++;
            continue;
        }
 
        // Given our length tolerance we can use it to compute a tolerance
        // on the angle between our spoke.
        double src_sin_tol = max_dist_rad / (src_dist_array[src_idx] + max_dist_rad);
 
        // Test if the small angle approximation will still hold. This is
        // defined as when sin(theta) ~= theta to within 0.1% of each
        // other. If the implied opening angle is too large we set it to
        // the 0.1% threshold.
        if (src_sin_tol > max_sin_tol) {
            src_sin_tol = max_sin_tol;
        }
 
        // Plane project the candidate source spoke and compute the cosine
        // and sine of the opening angle.
        double proj_src_delta_dot = ndarray::asEigenMatrix(src_delta_array[src_idx]).dot(src_ctr_eigen);
 
        ndarray::Array<double, 1, 1> proj_src_delta =
                copy(src_delta_array[src_idx] - proj_src_delta_dot * src_ctr);
        auto proj_src_delta_eigen = ndarray::asEigenMatrix(proj_src_delta);
        double geom_dist_src = sqrt(proj_src_delta_eigen.dot(proj_src_delta_eigen) * proj_src_ctr_dist_sq);
 
        // Compute cosine and sine of the delta vector opening angle.
        double cos_theta_src = proj_src_delta_eigen.dot(proj_src_ctr_delta_eigen) / geom_dist_src;
        Eigen::Vector3d cross_src =
                proj_src_delta_eigen.head<3>().cross(proj_src_ctr_delta_eigen.head<3>()) / geom_dist_src;
        double sin_theta_src = cross_src.dot(src_ctr_eigen);
 
        // Test the spokes and return the id of the reference object.
        // Return -1 if no match is found.
        int ref_id =
                check_spoke(cos_theta_src, sin_theta_src, ref_ctr, proj_ref_ctr_delta, proj_ref_ctr_dist_sq,
                            candidate_range, ref_id_array, reference_array, src_sin_tol);
        if (ref_id < 0) {
            n_fail++;
            continue;
        }
 
        output_spokes.push_back(std::make_pair(static_cast<size_t>(ref_id), src_idx + 1));
        if (output_spokes.size() >= n_match - 2) {
            break;
        }
    }
    return output_spokes;
}

create_shift_rot_matrix()

ShiftRotMatrixResult lsst::meas::astrom::create_shift_rot_matrix	(	double	cos_rot_sq,
		Eigen::Matrix3d const &	shift_matrix,
		ndarray::Array< double, 1, 1 > const &	src_delta,
		ndarray::Array< double, 1, 1 > const &	ref_ctr,
		ndarray::Array< double, 1, 1 > const &	ref_delta )

Create the complete spherical rotation matrix.

Parameters

[in]	cos_rot_sq	cosine of the rotation needed to align our source and reference candidate patterns.
[in]	shift_matrix	3x3 rotation matrix for shifting the source pattern center on top of the candidate reference pattern center.
[in]	src_delta	3 vector delta representing the first spoke of the source pattern
[in]	ref_ctr	3 vector on the unit-sphere representing the center of our reference pattern.
[in]	ref_delta	3 vector delta made by the first pair of the reference pattern.

Returns

Struct containing constructed matrix and implied rotation.

Definition at line 237 of file pessimisticPatternMatcherUtils.cc.

                                                                                          {
    double cos_rot = sqrt(cos_rot_sq);
    Eigen::Vector3d src_delta_eigen = ndarray::asEigenMatrix(src_delta);
    Eigen::Vector3d rot_src_delta = shift_matrix * src_delta_eigen;
    Eigen::Vector3d tmp_cross = rot_src_delta.cross(ndarray::asEigenMatrix(ref_delta).head<3>());
    double delta_dot_cross = tmp_cross.dot(ndarray::asEigenMatrix(ref_ctr));
 
    double sin_rot = sgn(delta_dot_cross) * sqrt(1 - cos_rot_sq);
    Eigen::Matrix3d rot_matrix =
            create_spherical_rotation_matrix(ndarray::asEigenMatrix(ref_ctr), cos_rot, sin_rot);
 
    Eigen::Matrix3d shift_rot_matrix = rot_matrix * shift_matrix;
 
    ShiftRotMatrixResult result;
    result.sin_rot = sin_rot;
    result.shift_rot_matrix = shift_rot_matrix;
    return result;
}

create_sorted_arrays()

SortedArrayResult lsst::meas::astrom::create_sorted_arrays	(	ndarray::Array< double, 1, 1 > const &	ref_center,
		ndarray::Array< double, 2, 1 > const &	reference_array )

Create arrays sorted on the distance between the center of this candidate reference object and all the reference objects.

Parameters

[in]	ref_center	Center point of the candidate pattern.
[in]	reference_array	Array of all reference object points.

Returns

Sorted distances and indexes of reference objects.

Definition at line 139 of file pessimisticPatternMatcherUtils.cc.

                                                                                          {
    SortedArrayResult result;
    // NOTE: this algorithm is quadratic in the length of reference_array. It might be worth using std::sort
    // instead of this approach, if reference_array is long, but the length at which that matters should
    // be tested because it would require caching of distances.
    for (uint16_t idx = 0; idx < reference_array.getShape()[0]; idx++) {
        Eigen::Vector3d diff = ndarray::asEigenMatrix(copy(reference_array[idx] - ref_center));
        double dist = sqrt(diff.dot(diff));
        auto dists_itr = std::lower_bound(result.dists.begin(), result.dists.end(), dist);
        auto ids_itr = result.ids.begin() + (dists_itr - result.dists.begin());
        result.ids.insert(ids_itr, idx);
        result.dists.insert(dists_itr, dist);
    }
    return result;
}

create_spherical_rotation_matrix()

Eigen::Matrix3d lsst::meas::astrom::create_spherical_rotation_matrix	(	Eigen::Vector3d const &	rot_axis,
		double	cos_rotation,
		double	sin_rotion )

Construct a generalized 3D rotation matrix about a given axis.

param[in] rot_axis 3 vector defining the axis to rotate about. param[in] cos_rotation cosine of the rotation angle. param[in] sin_rotion sine of the rotation angle.

Returns

3x3 spherical, rotation matrix.

Definition at line 223 of file pessimisticPatternMatcherUtils.cc.

                                                                      {
    Eigen::Matrix3d rot_cross_matrix;
    // clang-format off
    rot_cross_matrix << 0., -rot_axis[2], rot_axis[1],
                        rot_axis[2], 0., -rot_axis[0],
                        -rot_axis[1], rot_axis[0], 0.;
    // clang-format on
    Eigen::Matrix3d shift_matrix = cos_rotation * Eigen::Matrix3d::Identity() +
                                   sin_rotation * rot_cross_matrix +
                                   (1 - cos_rotation) * rot_axis * rot_axis.transpose();
    return shift_matrix;
}

find_candidate_reference_pair_range() [1/2]

std::pair< size_t, size_t > lsst::meas::astrom::find_candidate_reference_pair_range	(	float	src_dist,
		ndarray::Array< float, 1, 1 > const &	ref_dist_array,
		double	max_dist_rad )

Find the range of reference pairs within the distance tolerance of our source pair spoke.

Parameters

[in]	src_dist	Value of the distance we would like to search for in the reference array in radians.
[in]	ref_dist_array	sorted array of distances in radians.
[in]	max_dist_rad	maximum plus/minus search to find in the reference array in radians.

Returns

pair of indices for the min and max range of indices spanning src_dist +/- max_dist_rad to search.

Definition at line 156 of file pessimisticPatternMatcherUtils.cc.

                                                                                              {
    auto itr =
            std::lower_bound(ref_dist_array.begin(), ref_dist_array.end(), src_dist - max_dist_rad - 1e-16);
    auto itrEnd =
            std::upper_bound(ref_dist_array.begin(), ref_dist_array.end(), src_dist + max_dist_rad + 1e-16);
 
    size_t startIdx = itr - ref_dist_array.begin();
    size_t endIdx = itrEnd - ref_dist_array.begin();
    return std::make_pair(startIdx, endIdx);
}

find_candidate_reference_pair_range() [2/2]

std::pair< size_t, size_t > lsst::meas::astrom::find_candidate_reference_pair_range	(	float	src_dist,
		std::vector< float > const &	ref_dist_array,
		double	max_dist_rad )

Find the range of reference pairs within the distance tolerance of our source pair spoke.

Parameters

[in]	src_dist	Value of the distance we would like to search for in the reference array in radians.
[in]	ref_dist_array	sorted array of distances in radians.
[in]	max_dist_rad	maximum plus/minus search to find in the reference array in radians.

Returns

pair of indices for the min and max range of indices spanning src_dist +/- max_dist_rad to search.

Definition at line 168 of file pessimisticPatternMatcherUtils.cc.

                                                                                   {
    auto itr =
            std::lower_bound(ref_dist_array.begin(), ref_dist_array.end(), src_dist - max_dist_rad - 1e-16);
    auto itrEnd =
            std::upper_bound(ref_dist_array.begin(), ref_dist_array.end(), src_dist + max_dist_rad + 1e-16);
 
    size_t startIdx = itr - ref_dist_array.begin();
    size_t endIdx = itrEnd - ref_dist_array.begin();
    return std::make_pair(startIdx, endIdx);
}

intermediate_verify_comparison()

bool lsst::meas::astrom::intermediate_verify_comparison	(	std::vector< Eigen::Vector3d > const &	src_pattern,
		std::vector< Eigen::Vector3d > const &	ref_pattern,
		Eigen::Matrix3d const &	shift_rot_matrix,
		double	max_dist_rad )

Test the input rotation matrix by comparing the rotated src pattern to the ref pattern.

If every point in the pattern after rotation is within a distance of max_dist_rad to its candidate point in the other pattern, we return True.

Parameters

[in]	src_pattern	Vector of 3 vectors representing the points that make up our source pinwheel pattern.
[in]	ref_pattern	Vector of 3 vectors representing our candidate reference pinwheel pattern.
[in]	shift_rot_matrix	3x3 rotation matrix that takes the source objects and rotates them onto the frame of the reference objects.
[in]	max_dist_rad	Maximum distance allowed to consider two objects the same.

Returns

Whether or not the rotation matrix gives results within tolerances.

Definition at line 259 of file pessimisticPatternMatcherUtils.cc.

                                                                                                {
    double max_dist_sq = max_dist_rad * max_dist_rad;
    bool passed = true;
    auto iSrc = src_pattern.begin();
    auto iRef = ref_pattern.begin();
    for (auto iSrc = src_pattern.begin(), iRef = ref_pattern.begin();
         iSrc != src_pattern.end() && iRef != ref_pattern.end(); iSrc++, iRef++) {
        Eigen::Vector3d rot_src_vect = shift_rot_matrix * *iSrc;
        Eigen::Vector3d diff_vect = rot_src_vect - *iRef;
        if (max_dist_sq < diff_vect.dot(diff_vect)) {
            passed = false;
            break;
        }
    }
    return passed;
}

makeMatchStatistics()

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);
}

makeMatchStatistics< afw::table::ReferenceMatch >()

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 )

makeMatchStatistics< afw::table::SourceMatch >()

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 )

makeMatchStatisticsInPixels()

template<typename MatchT >

afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInPixels	(	afw::geom::SkyWcs 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 48 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);
}

makeMatchStatisticsInPixels< afw::table::ReferenceMatch >()

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

makeMatchStatisticsInPixels< afw::table::SourceMatch >()

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

makeMatchStatisticsInRadians()

template<typename MatchT >

afw::math::Statistics lsst::meas::astrom::makeMatchStatisticsInRadians	(	afw::geom::SkyWcs 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 71 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.separation(srcCoord);
        val.push_back(angSep.asRadians());
    }
    return afw::math::makeStatistics(val, flags, sctrl);
}

makeMatchStatisticsInRadians< afw::table::ReferenceMatch >()

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

makeMatchStatisticsInRadians< afw::table::SourceMatch >()

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

makeProxies() [1/2]

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

Definition at line 448 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;
}

makeProxies() [2/2]

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

Definition at line 436 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;
}

makeWcs()

std::shared_ptr< afw::geom::SkyWcs > lsst::meas::astrom::makeWcs	(	SipForwardTransform const &	sipForward,
		SipReverseTransform const &	sipReverse,
		geom::SpherePoint 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	ICRS 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 148 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(ndarray::asEigenMatrix(sipForward.getPoly().getXCoeffs()));
    Eigen::MatrixXd sipB(ndarray::asEigenMatrix(sipForward.getPoly().getYCoeffs()));
    Eigen::MatrixXd sipAP(ndarray::asEigenMatrix(sipReverse.getPoly().getXCoeffs()));
    Eigen::MatrixXd sipBP(ndarray::asEigenMatrix(sipReverse.getPoly().getYCoeffs()));
 
    return afw::geom::makeTanSipWcs(sipForward.getPixelOrigin(), skyOrigin,
                                    sipForward.getCdMatrix().getMatrix(), sipA, sipB, sipAP, sipBP);
}

matchOptimisticB()

afwTable::ReferenceMatchVector lsst::meas::astrom::matchOptimisticB	(	afw::table::SimpleCatalog const &	posRefCat,
		afw::table::SourceCatalog const &	sourceCat,
		MatchOptimisticBControl const &	control,
		afw::geom::SkyWcs 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 459 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 = geom::degToRad(control.maxRotationDeg);
 
    // Create an undistorted Wcs to project everything with
    // We'll anchor it at the center of the area.
    geom::Extent2D srcCenter(0, 0);
    for (auto iter = sourceCat.begin(); iter != sourceCat.end(); ++iter) {
        srcCenter += geom::Extent2D(iter->getCentroid());
    }
    srcCenter /= sourceCat.size();
 
    sphgeom::Vector3d refCenter(0, 0, 0);
    for (auto iter = posRefCat.begin(); iter != posRefCat.end(); ++iter) {
        refCenter += iter->getCoord().getVector();
    }
    refCenter /= posRefCat.size();
 
    auto tanWcs =
            afw::geom::makeSkyWcs(geom::Point2D(srcCenter), geom::SpherePoint(refCenter), wcs.getCdMatrix());
 
    ProxyVector posRefProxyCat = makeProxies(posRefCat, *tanWcs);
    ProxyVector sourceProxyCat = makeProxies(sourceCat, wcs, *tanWcs);
 
    // 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];
                    geom::Angle const theta(std::acos((a * b + c * d) /
                                                      (std::sqrt(a * a + c * c) * std::sqrt(b * b + d * d))),
                                            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;
    }
}

PYBIND11_MODULE()

lsst::meas::astrom::PYBIND11_MODULE	(	_measAstromLib	,
		mod	)

Definition at line 25 of file _measAstromLib.cc.

                                     {
    WrapperCollection wrappers(mod, "lsst.meas.astrom");
 
    wrappers.addInheritanceDependency("lsst.afw.math");
 
    sip::wrapLeastSqFitter1d(wrappers);
    sip::wrapCreateWcsWithSip(wrappers);
    sip::wrapMatchSrcToCatalogue(wrappers);
    sip::wrapLeastSqFitter2d(wrappers);
 
    wrapPolynomialTransform(wrappers);
    wrapScaledPolynomialTransformFitter(wrappers);
    wrapSipTransform(wrappers);
    wrapMatchOptimisticB(wrappers);
    wrapMakeMatchStatistics(wrappers);
    wrapPessimisticPatternMatcherUtils(wrappers);
    wrappers.finish();
};

rotateWcsPixelsBy90()

std::shared_ptr< afw::geom::SkyWcs > lsst::meas::astrom::rotateWcsPixelsBy90	(	afw::geom::SkyWcs const &	wcs,
		int	nQuarter,
		geom::Extent2I const &	dimensions )

Return a new SkyWcs that represents a rotation of the image it corresponds to about the image's center.

Parameters

[in]	wcs	Original SkyWcs to be rotated.
[in]	nQuarter	Number of 90 degree rotations (positive is counterclockwise).
[in]	dimensions	Width and height of the image.

Definition at line 179 of file SipTransform.cc.

                                                                                       {
    geom::Extent2D offset;
    switch (nQuarter % 4) {
        case 0:
            offset = geom::Extent2D(0, 0);
            break;
        case 1:
            offset = geom::Extent2D(dimensions.getY() - 1, 0);
            break;
        case 2:
            offset = geom::Extent2D(dimensions - geom::Extent2I(1, 1));
            break;
        case 3:
            offset = geom::Extent2D(0, dimensions.getX() - 1);
            break;
    }
    auto rot = geom::LinearTransform::makeRotation(nQuarter * 90.0 * geom::degrees);
    return transformWcsPixels(wcs, geom::AffineTransform(rot, offset));
}

test_rotation()

RotationTestResult lsst::meas::astrom::test_rotation	(	ndarray::Array< double, 1, 1 > const &	src_center,
		ndarray::Array< double, 1, 1 > const &	ref_center,
		ndarray::Array< double, 1, 1 > const &	src_delta,
		ndarray::Array< double, 1, 1 > const &	ref_delta,
		double	cos_shift,
		double	max_cos_rot_sq )

Test if the rotation implied between the source pattern and reference pattern is within tolerance.

To test this we need to create the first part of our spherical rotation matrix which we also return for use later.

Parameters

[in]	src_center	3 vector defining the center of the candidate source pinwheel pattern.
[in]	ref_center	3 vector defining the center of the candidate reference pinwheel pattern.
[in]	src_delta	3 vector delta between the source pattern center and the end of the pinwheel spoke.
[in]	ref_delta	3 vector delta of the candidate matched reference pair
[in]	cos_shift	Cosine of the angle between the source and reference candidate centers.
[in]	max_cos_rot_sq	Maximum allowed rotation of the pinwheel pattern.

Returns

Result of the test rotation.

Definition at line 181 of file pessimisticPatternMatcherUtils.cc.

                                                        {
    // Make sure the sine is a real number.
    if (cos_shift > 1.0) {
        cos_shift = 1;
    } else if (cos_shift < -1.0) {
        cos_shift = -1;
    }
    double sin_shift = sqrt(1 - cos_shift * cos_shift);
 
    // If the sine of our shift is zero we only need to use the identity matrix for the shift. Else we
    // construct the rotation matrix for shift.
    auto ref_center_eigen = ndarray::asEigenMatrix(ref_center).head<3>();
    Eigen::Matrix3d shift_matrix;
    if (sin_shift > 0) {
        Eigen::Vector3d rot_axis = ndarray::asEigenMatrix(src_center).head<3>().cross(ref_center_eigen);
        rot_axis /= sin_shift;
        shift_matrix = create_spherical_rotation_matrix(rot_axis, cos_shift, sin_shift);
    } else {
        shift_matrix = Eigen::Matrix3d::Identity();
    }
 
    // Now that we have our shift we apply it to the src delta vector and check the rotation.
    Eigen::Vector3d rot_src_delta = shift_matrix * ndarray::asEigenMatrix(src_delta);
    Eigen::Vector3d proj_src_delta = rot_src_delta - rot_src_delta.dot(ref_center_eigen) * ref_center_eigen;
    Eigen::Vector3d ref_delta_eigen = ndarray::asEigenMatrix(ref_delta);
 
    Eigen::Vector3d proj_ref_delta =
            ref_delta_eigen - ref_delta_eigen.dot(ref_center_eigen) * ref_center_eigen;
    double proj_src_delta_sq = proj_src_delta.dot(proj_ref_delta);
    double cos_rot_sq = proj_src_delta_sq * proj_src_delta_sq /
                        (proj_src_delta.dot(proj_src_delta) * proj_ref_delta.dot(proj_ref_delta));
    if (cos_rot_sq < max_cos_rot_sq) {
        // Return failure if the rotation isn't within tolerance.
        return RotationTestResult();
    }
    return RotationTestResult(cos_rot_sq, proj_ref_delta, shift_matrix);
}

transformWcsPixels()

std::shared_ptr< afw::geom::SkyWcs > lsst::meas::astrom::transformWcsPixels	(	afw::geom::SkyWcs const &	wcs,
		geom::AffineTransform const &	s )

Create a new SkyWcs whose pixel coordinate system has been transformed via an affine transform.

Parameters

[in]	wcs	Original SkyWcs object.
[in]	s	AffineTransform to apply to the pixel coordinate system.

Returns

a new Wcs that satisfies the following:

newWcs = transformWcsPixels(wcs, s);
assert(newWcs.skyToPixel(sky), s(wcs.skyToPixel(sky)));
assert(newWcs.pixelToSky(pixel), wcs.pixelToSky(s.inverted()(pixel)));

for all sky coordinates sky and pixel coordinates pixel.

Definition at line 173 of file SipTransform.cc.

                                                                                    {
    auto affineTransform22 = afw::geom::makeTransform(s);
    return afw::geom::makeModifiedWcs(*affineTransform22->inverted(), wcs, true);
}

wrapMakeMatchStatistics()

void lsst::meas::astrom::wrapMakeMatchStatistics

(

WrapperCollection &

wrappers

)

Definition at line 48 of file makeMatchStatistics.cc.

                                                                           {
    auto &mod = wrappers.module;
 
    declareMakeMatchStatistics<afw::table::ReferenceMatch>(mod);
    declareMakeMatchStatistics<afw::table::SourceMatch>(mod);
}

wrapMatchOptimisticB()

void lsst::meas::astrom::wrapMatchOptimisticB

(

WrapperCollection &

wrappers

)

Definition at line 95 of file matchOptimisticB.cc.

                                                                        {
    declareRecordProxy(wrappers);
    declareProxyPair(wrappers);
    declareMatchOptimisticBControl(wrappers);
 
    wrappers.module.def("makeProxies",
            (ProxyVector(*)(afw::table::SourceCatalog const &, afw::geom::SkyWcs const &,
                            afw::geom::SkyWcs const &)) &
                    makeProxies,
            "sourceCat"_a, "distortedWcs"_a, "tanWcs"_a);
    wrappers.module.def("makeProxies",
            (ProxyVector(*)(afw::table::SimpleCatalog const &, afw::geom::SkyWcs const &)) & makeProxies,
            "posRefCat"_a, "tanWcs"_a);
 
    wrappers.module.def("matchOptimisticB", &matchOptimisticB, "posRefCat"_a, "sourceCat"_a, "control"_a, "wcs"_a,
            "posRefBegInd"_a = 0, "verbose"_a = false);
}

wrapPessimisticPatternMatcherUtils()

void lsst::meas::astrom::wrapPessimisticPatternMatcherUtils

(

WrapperCollection &

wrappers

)

Definition at line 41 of file pessimisticPatternMatcherUtils.cc.

                                                                                      {
    wrappers.wrapType(py::class_<PatternResult>(wrappers.module, "PatternResult"), [](auto &mod, auto &cls) {
        cls.def_readonly("candidate_pairs", &PatternResult::candidate_pairs);
        cls.def_readonly("shift_rot_matrix", &PatternResult::shift_rot_matrix);
        cls.def_readonly("cos_shift", &PatternResult::cos_shift);
        cls.def_readonly("sin_rot", &PatternResult::sin_rot);
        cls.def_readonly("success", &PatternResult::success);
        mod.def("construct_pattern_and_shift_rot_matrix", &construct_pattern_and_shift_rot_matrix,
                "src_pattern_array"_a, "src_delta_array"_a, "src_dist_array"_a, "dist_array"_a, "id_array"_a,
                "reference_array"_a, "n_match"_a, "max_cos_theta_shift"_a, "max_cos_rot_sq"_a, "max_dist_rad"_a);
    });
}

wrapPolynomialTransform()

void lsst::meas::astrom::wrapPolynomialTransform

(

WrapperCollection &

wrappers

)

Definition at line 106 of file polynomialTransform.cc.

                                                                           {
    declarePolynomialTransform(wrappers);
    declareScaledPolynomialTransform(wrappers);
 
    wrappers.module.def("compose",
            (PolynomialTransform(*)(geom::AffineTransform const &, PolynomialTransform const &)) & compose,
            "t1"_a, "t2"_a);
    wrappers.module.def("compose",
            (PolynomialTransform(*)(PolynomialTransform const &, geom::AffineTransform const &)) & compose,
            "t1"_a, "t2"_a);
}

wrapScaledPolynomialTransformFitter()

void lsst::meas::astrom::wrapScaledPolynomialTransformFitter

(

WrapperCollection &

wrappers

)

Definition at line 75 of file scaledPolynomialTransformFitter.cc.

                                                                                       {
    declareOutlierRejectionControl(wrappers);
    declareScaledPolynomialTransformFitter(wrappers);
}

wrapSipTransform()

void lsst::meas::astrom::wrapSipTransform

(

WrapperCollection &

wrappers

)

Definition at line 111 of file sipTransform.cc.

                                                                    {
    declareSipTransformBase(wrappers);
    declareSipForwardTransform(wrappers);
    declareSipReverseTransform(wrappers);
 
    wrappers.module.def("makeWcs", makeWcs, "sipForward"_a, "sipReverse"_a, "skyOrigin"_a);
    wrappers.module.def("transformWcsPixels", transformWcsPixels, "wcs"_a, "s"_a);
    wrappers.module.def("rotateWcsPixelsBy90", rotateWcsPixelsBy90, "wcs"_a, "nQuarter"_a, "dimensions"_a);
}