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"
465 if (posRefCat.empty()) {
468 if (sourceCat.empty()) {
471 if (posRefBegInd < 0) {
474 if (static_cast<size_t>(posRefBegInd) >= posRefCat.size()) {
477 double const maxRotationRad =
geom::degToRad(control.maxRotationDeg);
482 for (
auto iter = sourceCat.begin();
iter != sourceCat.end(); ++
iter) {
485 srcCenter /= sourceCat.size();
487 sphgeom::Vector3d refCenter(0, 0, 0);
488 for (
auto iter = posRefCat.begin();
iter != posRefCat.end(); ++
iter) {
489 refCenter +=
iter->getCoord().getVector();
491 refCenter /= posRefCat.size();
501 selectPoint(sourceProxyCat, sourceCat.getSchema().find<
double>(control.sourceFluxField).
key,
502 control.numBrightStars);
507 selectPoint(posRefProxyCat, posRefCat.getSchema().find<
double>(control.refFluxField).
key,
508 sourceSubCat.
size() + 25, posRefBegInd);
515 size_t const posRefCatSubSize = posRefSubCat.
size();
516 for (
size_t i = 0; i < posRefCatSubSize - 1; i++) {
517 for (
size_t j = i + 1; j < posRefCatSubSize; j++) {
525 size_t const sourceSubCatSize = sourceSubCat.
size();
526 for (
size_t i = 0; i < sourceSubCatSize - 1; i++) {
527 for (
size_t j = i + 1; j < sourceSubCatSize; j++) {
528 sourcePairList.push_back(
ProxyPair(sourceSubCat[i], sourceSubCat[j]));
531 std::sort(sourcePairList.begin(), sourcePairList.end(), cmpPair);
537 size_t const fullShapeSize = control.numPointsForShape - 1;
538 for (
size_t ii = 0; ii < sourcePairList.size(); ii++) {
542 searchPair(posRefPairList, p, control.matchingAllowancePix, maxRotationRad);
550 for (
size_t l = 0; l < q.
size(); l++) {
553 double dpa = p.
pa - q[l].pa;
566 for (
size_t k = 0; k < sourceSubCat.
size(); k++) {
567 if (p.
first == sourceSubCat[k] || p.
second == sourceSubCat[k])
continue;
572 posRefPairList, pp, q[l], control.matchingAllowancePix, dpa, maxRotationRad);
573 if (r != posRefPairList.
end()) {
580 if (srcMatPair.
size() == fullShapeSize) {
586 bool goodMatch =
false;
587 if (srcMatPair.
size() == fullShapeSize) {
589 for (
size_t k = 1; k < catMatPair.
size(); k++) {
590 if (catMatPair[0].
first != catMatPair[k].
first) {
597 if (goodMatch && srcMatPair.
size() == fullShapeSize) {
603 for (
size_t k = 0; k < srcMatPair.
size(); k++) {
608 boost::shared_array<double>
coeff = polyfit(1, srcMat, catMat);
611 for (
size_t k = 0; k < srcMat.
size(); k++) {
612 std::cout <<
"circle(" << srcMat[k].getX() <<
"," << srcMat[k].getY()
614 std::cout <<
"circle(" << catMat[k].getX() <<
"," << catMat[k].getY()
616 std::cout <<
"line(" << srcMat[0].getX() <<
"," << srcMat[0].getY() <<
"," 617 << srcMat[k].getX() <<
"," << srcMat[k].getY()
618 <<
") # line=0 0 color=green" <<
std::endl;
619 std::cout <<
"line(" << catMat[0].getX() <<
"," << catMat[0].getY() <<
"," 620 << catMat[k].getX() <<
"," << catMat[k].getY()
621 <<
") # line=0 0 color=red" <<
std::endl;
636 std::cout <<
"Determinant (max " << control.maxDeterminant <<
"): ";
638 std::cout <<
"Angle between axes (deg; allowed 90 +/- ";
639 std::cout << control.allowedNonperpDeg <<
"): ";
642 if (
std::fabs(coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1.) > control.maxDeterminant ||
643 std::fabs(theta.asDegrees() - 90) > control.allowedNonperpDeg ||
644 std::fabs(coeff[0]) > control.maxOffsetPix ||
645 std::fabs(coeff[3]) > control.maxOffsetPix) {
651 double x0, y0, x1, y1;
655 for (
size_t i = 0; i < sourceSubCat.
size(); i++) {
656 x0 = sourceSubCat[i].getX();
657 y0 = sourceSubCat[i].getY();
660 searchNearestPoint(posRefSubCat, x1, y1, control.matchingAllowancePix);
661 if (refObjDist.first != posRefSubCat.
end()) {
666 std::cout <<
"Match: " << x0 <<
"," << y0 <<
" --> " << x1 <<
"," << y1
667 <<
" <==> " << refObjDist.first->getX() <<
"," 668 << refObjDist.first->getY() <<
std::endl;
672 if (num <= control.numPointsForShape) {
679 coeff = polyfit(1, srcMat, catMat);
682 for (
size_t i = 0; i < 6; ++i) {
688 matPair = FinalVerify(coeff, posRefProxyCat, sourceProxyCat,
689 control.matchingAllowancePix, verbose);
691 std::cout <<
"Number of matches: " << matPair.size() <<
" vs " 694 if (matPair.size() <=
static_cast<std::size_t>(control.minMatchedPairs)) {
698 if (matPair.size() > matPairSave.
size()) {
699 matPairSave = matPair;
707 if (matPairCand.
size() == 3) {
718 if (matPairCand.
size() == 0) {
721 size_t nmatch = matPairCand[0].
size();
723 for (
size_t i = 1; i < matPairCand.
size(); i++) {
724 if (matPairCand[i].size() > nmatch) {
725 nmatch = matPairCand[i].
size();
726 matPairRet = matPairCand[i];
AngleUnit constexpr radians
constant with units of radians
A class representing an angle.
table::Key< table::Array< std::uint8_t > > wcs
ProxyVector makeProxies(afw::table::SourceCatalog const &sourceCat, afw::geom::SkyWcs const &distortedWcs, afw::geom::SkyWcs const &tanWcs)
table::Key< table::Array< double > > coeff
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Point in an unspecified spherical coordinate system.
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Reports invalid arguments.
Extent< double, 2 > Extent2D