LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Namespaces | Classes | Typedefs | Functions
lsst::meas::astrom Namespace Reference

Namespaces

 anetAstrometry
 
 anetBasicAstrometry
 
 approximateWcs
 
 astrometry
 
 astrometryNetDataConfig
 
 catalogStarSelector
 
 detail
 
 display
 
 fitTanSipWcs
 
 loadAstrometryNetObjects
 
 matchOptimisticB
 
 multiindex
 
 setMatchDistance
 
 sip
 
 verifyWcs
 
 version
 

Classes

struct  RecordProxy
 
struct  ProxyPair
 
struct  MatchOptimisticBControl
 

Typedefs

typedef std::vector< RecordProxy > ProxyVector
 

Functions

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)
 

Typedef Documentation

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

Definition at line 53 of file matchOptimisticB.h.

Function Documentation

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

Definition at line 492 of file matchOptimisticB.cc.

496  {
497  ProxyVector r;
498  r.reserve(sourceCat.size());
499  for (afwTable::SourceCatalog::const_iterator sourcePtr = sourceCat.begin();
500  sourcePtr != sourceCat.end(); ++sourcePtr) {
501  r.push_back(RecordProxy(sourcePtr,
502  tanWcs.skyToPixel(*distortedWcs.pixelToSky(sourcePtr->getCentroid()))));
503  }
504  return r;
505  }
lsst::meas::astrom::ProxyVector
std::vector< RecordProxy > ProxyVector
Definition: matchOptimisticB.h:53
lsst.afw.table::SortedCatalogT::const_iterator
Base::const_iterator const_iterator
Definition: SortedCatalog.h:49
ProxyVector lsst::meas::astrom::makeProxies ( lsst::afw::table::SimpleCatalog const &  posRefCat,
afw::image::Wcs const &  tanWcs 
)

Definition at line 507 of file matchOptimisticB.cc.

510  {
511  auto coordKey = afwTable::CoordKey(posRefCat.getSchema()["coord"]);
512  ProxyVector r;
513  r.reserve(posRefCat.size());
514  for (afwTable::SimpleCatalog::const_iterator posRefPtr = posRefCat.begin();
515  posRefPtr != posRefCat.end(); ++posRefPtr) {
516  r.push_back(RecordProxy(posRefPtr, tanWcs.skyToPixel(posRefPtr->get(coordKey))));
517  }
518  return r;
519  }
lsst::meas::astrom::ProxyVector
std::vector< RecordProxy > ProxyVector
Definition: matchOptimisticB.h:53
lsst.afw.table::SortedCatalogT::const_iterator
Base::const_iterator const_iterator
Definition: SortedCatalog.h:49
lsst.afw.table::CoordKey
A FunctorKey used to get or set celestial coordiantes from a pair of Angle keys.
Definition: aggregates.h:119
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]posRefCatcatalog of position reference stars; fields read:
  • "coord"
  • control.refFluxField
[in]sourceCatcatalog of detected sources; fields read:
  • "Centroid_x"
  • "Centroid_y"
  • control.refFluxField
[in]wcsestimated WCS
[in]controlcontrol object
[in]posRefBegIndindex of first start to use in posRefCat
[in]verbosetrue 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 521 of file matchOptimisticB.cc.

528  {
529  control.validate();
530  if (posRefBegInd < 0) {
531  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd < 0");
532  }
533  if (posRefBegInd >= posRefCat.size()) {
534  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd too big");
535  }
536  double const maxRotationRad = afw::geom::degToRad(control.maxRotationDeg);
537 
538  CONST_PTR(afw::image::Wcs) tanWcs; // Undistorted Wcs, providing tangent plane
539  if (wcs.hasDistortion()) {
540  // Create an undistorted Wcs to project everything with
541  // We'll anchor it at the center of the area.
542  afw::geom::Extent2D srcCenter(0, 0);
543  for (auto iter = sourceCat.begin(); iter != sourceCat.end(); ++iter) {
544  srcCenter += afw::geom::Extent2D(iter->getCentroid());
545  }
546  srcCenter /= sourceCat.size();
547 
548  afw::geom::Extent3D refCenter(0, 0, 0);
549  for (auto iter = posRefCat.begin(); iter != posRefCat.end(); ++iter) {
550  refCenter += afw::geom::Extent3D(iter->getCoord().toIcrs().getVector());
551  }
552  refCenter /= posRefCat.size();
553 
554  tanWcs = boost::make_shared<afw::image::Wcs>(
555  afw::coord::IcrsCoord(afw::geom::Point3D(refCenter)).getPosition(),
556  afw::geom::Point2D(srcCenter),
557  wcs.getCDMatrix()
558  );
559  }
560 
561  ProxyVector posRefProxyCat = makeProxies(posRefCat, tanWcs ? *tanWcs : wcs);
562  ProxyVector sourceProxyCat = makeProxies(sourceCat, wcs, tanWcs ? *tanWcs : wcs);
563 
564  // sourceSubCat contains at most the numBrightStars brightest sources, sorted by decreasing flux
565  ProxyVector sourceSubCat = selectPoint(
566  sourceProxyCat,
567  sourceCat.getSchema().find<double>(control.sourceFluxField).key,
568  control.numBrightStars);
569 
570  // posRefSubCat skips the initial posRefBegInd brightest reference objects and contains
571  // at most the next len(sourceSubCat) + 25 brightest reference objects, sorted by decreasing flux
572  ProxyVector posRefSubCat = selectPoint(
573  posRefProxyCat,
574  posRefCat.getSchema().find<double>(control.refFluxField).key,
575  sourceSubCat.size()+25,
576  posRefBegInd);
577  if (verbose) {
578  std::cout << "Catalog sizes: " << sourceSubCat.size() << " " << posRefSubCat.size() << std::endl;
579  }
580 
581  // Construct a list of pairs of position reference stars sorted by increasing separation
582  std::vector<ProxyPair> posRefPairList;
583  size_t const posRefCatSubSize = posRefSubCat.size();
584  for (size_t i = 0; i < posRefCatSubSize-1; i++) {
585  for (size_t j = i+1; j < posRefCatSubSize; j++) {
586  posRefPairList.push_back(ProxyPair(posRefSubCat[i], posRefSubCat[j]));
587  }
588  }
589  std::sort(posRefPairList.begin(), posRefPairList.end(), cmpPair);
590 
591  // Construct a list of pairs of sources sorted by increasing separation
592  std::vector<ProxyPair> sourcePairList;
593  size_t const sourceSubCatSize = sourceSubCat.size();
594  for (size_t i = 0; i < sourceSubCatSize-1; i++) {
595  for (size_t j = i+1; j < sourceSubCatSize; j++) {
596  sourcePairList.push_back(ProxyPair(sourceSubCat[i], sourceSubCat[j]));
597  }
598  }
599  std::sort(sourcePairList.begin(), sourcePairList.end(), cmpPair);
600 
601  afwTable::ReferenceMatchVector matPair;
602  afwTable::ReferenceMatchVector matPairSave;
603  std::vector<afwTable::ReferenceMatchVector> matPairCand;
604 
605  size_t const fullShapeSize = control.numPointsForShape - 1; // Max size of shape array
606  for (size_t ii = 0; ii < sourcePairList.size(); ii++) {
607  ProxyPair p = sourcePairList[ii];
608 
609  std::vector<ProxyPair> q = searchPair(posRefPairList, p,
610  control.matchingAllowancePix, maxRotationRad);
611 
612  // If candidate pairs are found
613  if (q.size() != 0) {
614 
615  std::vector<ProxyPair> srcMatPair;
616  std::vector<ProxyPair> catMatPair;
617 
618  // Go through candidate pairs
619  for (size_t l = 0; l < q.size(); l++) {
620  // sign matters, so don't use deltaAng; no need to wrap because
621  // the result is used with deltaAng later
622  double dpa = p.pa - q[l].pa;
623 
624  srcMatPair.clear();
625  catMatPair.clear();
626 
627  srcMatPair.push_back(p);
628  catMatPair.push_back(q[l]);
629 
630  if (verbose) {
631  std::cout << "p dist: " << p.distance << " pa: " << p.pa << std::endl;
632  std::cout << "q dist: " << q[l].distance << " pa: " << q[l].pa << std::endl;
633  }
634 
635  for (size_t k = 0; k < sourceSubCat.size(); k++) {
636  if (p.first == sourceSubCat[k] || p.second == sourceSubCat[k]) continue;
637 
638  ProxyPair pp(p.first, sourceSubCat[k]);
639 
640  std::vector<ProxyPair>::iterator r = searchPair3(posRefPairList, pp, q[l],
641  control.matchingAllowancePix, dpa, maxRotationRad);
642  if (r != posRefPairList.end()) {
643  srcMatPair.push_back(pp);
644  catMatPair.push_back(*r);
645  if (verbose) {
646  std::cout << " p dist: " << pp.distance << " pa: " << pp.pa << std::endl;
647  std::cout << " r dist: " << (*r).distance << " pa: " << (*r).pa << std::endl;
648  }
649  if (srcMatPair.size() == fullShapeSize) {
650  break;
651  }
652  }
653  }
654 
655  bool goodMatch = false;
656  if (srcMatPair.size() == fullShapeSize) {
657  goodMatch = true;
658  for (size_t k = 1; k < catMatPair.size(); k++) {
659  if (catMatPair[0].first != catMatPair[k].first) {
660  goodMatch = false;
661  break;
662  }
663  }
664  }
665 
666  if (goodMatch && srcMatPair.size() == fullShapeSize) {
667 
668  ProxyVector srcMat;
669  ProxyVector catMat;
670 
671  srcMat.push_back(srcMatPair[0].first);
672  catMat.push_back(catMatPair[0].first);
673  for (size_t k = 0; k < srcMatPair.size(); k++) {
674  srcMat.push_back(srcMatPair[k].second);
675  catMat.push_back(catMatPair[k].second);
676  }
677 
678  boost::shared_array<double> coeff = polyfit(1, srcMat, catMat);
679 
680  if (verbose) {
681  for (size_t k = 0; k < srcMat.size(); k++) {
682  std::cout << "circle(" << srcMat[k].getX() << ","
683  << srcMat[k].getY() << ",10) # color=green" << std::endl;
684  std::cout << "circle(" << catMat[k].getX() << ","
685  << catMat[k].getY() << ",10) # color=red" << std::endl;
686  std::cout << "line(" << srcMat[0].getX() << "," << srcMat[0].getY() << ","
687  << srcMat[k].getX() << "," << srcMat[k].getY()
688  << ") # line=0 0 color=green" << std::endl;
689  std::cout << "line(" << catMat[0].getX() << "," << catMat[0].getY() << ","
690  << catMat[k].getX() << "," << catMat[k].getY()
691  << ") # line=0 0 color=red" << std::endl;
692  }
693  }
694 
695  double a = coeff[1];
696  double b = coeff[2];
697  double c = coeff[4];
698  double d = coeff[5];
699  afw::geom::Angle const theta(std::acos((a*b+c*d)/
700  (std::sqrt(a*a+c*c)*std::sqrt(b*b+d*d))),
701  afw::geom::radians);
702  if (verbose) {
703  std::cout << "Linear fit from match:" << std::endl;
704  std::cout << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
705  std::cout << coeff[3] << " " << coeff[4] << " " << coeff[5] << std::endl;
706  std::cout << "Determinant (max " << control.maxDeterminant << "): ";
707  std::cout << coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1. << std::endl;
708  std::cout << "Angle between axes (deg; allowed 90 +/- ";
709  std::cout << control.allowedNonperpDeg << "): ";
710  std::cout << theta.asDegrees() << std::endl;
711  }
712  if (std::fabs(coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1.)
713  > control.maxDeterminant ||
714  std::fabs(theta.asDegrees() - 90) > control.allowedNonperpDeg ||
715  std::fabs(coeff[0]) > control.maxOffsetPix ||
716  std::fabs(coeff[3]) > control.maxOffsetPix) {
717  if (verbose) {
718  std::cout << "Bad; continuing" << std::endl;
719  }
720  continue;
721  } else {
722  double x0, y0, x1, y1;
723  int num = 0;
724  srcMat.clear();
725  catMat.clear();
726  for (size_t i = 0; i < sourceSubCat.size(); i++) {
727  x0 = sourceSubCat[i].getX();
728  y0 = sourceSubCat[i].getY();
729  transform(1, coeff, x0, y0, &x1, &y1);
730  auto refObjDist = searchNearestPoint(posRefSubCat, x1, y1,
731  control.matchingAllowancePix);
732  if (refObjDist.first != posRefSubCat.end()) {
733  num++;
734  srcMat.push_back(sourceSubCat[i]);
735  catMat.push_back(*refObjDist.first);
736  if (verbose) {
737  std::cout << "Match: " << x0 << "," << y0 << " --> "
738  << x1 << "," << y1 << " <==> "
739  << refObjDist.first->getX() << "," << refObjDist.first->getY()
740  << std::endl;
741  }
742  }
743  }
744  if (num <= control.numPointsForShape) {
745  // Can get matrix = 0,0,0,0; everything matches a single catalog object
746  if (verbose) {
747  std::cout << "Insufficient initial matches; continuing" << std::endl;
748  }
749  continue;
750  }
751  coeff = polyfit(1, srcMat, catMat);
752  if (verbose) {
753  std::cout << "Coefficients from initial matching:" << std::endl;
754  for (size_t i = 0; i < 6; ++i) {
755  std::cout << coeff[i] << " ";
756  }
757  std::cout << std::endl;
758  }
759 
760  matPair = FinalVerify(coeff, posRefProxyCat, sourceProxyCat,
761  control.matchingAllowancePix, verbose);
762  if (verbose) {
763  std::cout << "Number of matches: " << matPair.size() << " vs " <<
764  control.minMatchedPairs << std::endl;
765  }
766  if (matPair.size() <= control.minMatchedPairs) {
767  if (verbose) {
768  std::cout << "Insufficient final matches; continuing" << std::endl;
769  }
770  if (matPair.size() > matPairSave.size()) {
771  matPairSave = matPair;
772  }
773  continue;
774  } else {
775  if (verbose) {
776  std::cout << "Finish" << std::endl;
777  }
778  matPairCand.push_back(matPair);
779  if (matPairCand.size() == 3) {
780  goto END;
781  }
782  }
783  }
784  }
785  }
786  }
787  }
788 
789  END:
790  if (matPairCand.size() == 0) {
791  return matPairSave;
792  } else {
793  size_t nmatch = matPairCand[0].size();
794  afwTable::ReferenceMatchVector matPairRet = matPairCand[0];
795  for (size_t i = 1; i < matPairCand.size(); i++) {
796  if (matPairCand[i].size() > nmatch) {
797  nmatch = matPairCand[i].size();
798  matPairRet = matPairCand[i];
799  }
800  }
801  return matPairRet;
802  }
803  }
iter
int iter
Definition: GaussianCentroid.cc:73
lsst.afw.geom::radians
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
lsst.afw.geom::Extent2D
Extent< double, 2 > Extent2D
Definition: Extent.h:358
CONST_PTR
#define CONST_PTR(...)
Definition: base.h:47
wcs
tbl::Key< int > wcs
Definition: CoaddBoundedField.cc:131
lsst.afw.table::ReferenceMatchVector
std::vector< ReferenceMatch > ReferenceMatchVector
Definition: fwd.h:97
lsst.afw.table::IcrsCoord
lsst::afw::coord::IcrsCoord IcrsCoord
Definition: misc.h:38
lsst.afw.display::x0
int const x0
Definition: saturated.cc:45
lsst::meas::astrom::ProxyVector
std::vector< RecordProxy > ProxyVector
Definition: matchOptimisticB.h:53
lsst.afw.geom::Extent3D
Extent< double, 3 > Extent3D
Definition: Extent.h:359
lsst.afw.display::if
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
lsst.afw.geom::Point3D
Point< double, 3 > Point3D
Definition: Point.h:287
Wcs
lsst::afw::image::Wcs Wcs
Definition: Wcs.cc:60
degToRad
double degToRad(long double angleInDegrees)
Definition: RaDecStr.cc:67
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
lsst.afw.table::Angle
lsst::afw::geom::Angle Angle
Definition: misc.h:39
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
lsst.afw.geom::Point2D
Point< double, 2 > Point2D
Definition: Point.h:286
lsst::meas::astrom::makeProxies
ProxyVector makeProxies(lsst::afw::table::SourceCatalog const &sourceCat, afw::image::Wcs const &distortedWcs, afw::image::Wcs const &tanWcs)
Definition: matchOptimisticB.cc:492
lsst.afw.display::y0
int const y0
Definition: saturated.cc:45
lsst.afw.coord::Coord::getPosition
lsst::afw::geom::Point2D getPosition(lsst::afw::geom::AngleUnit unit=lsst::afw::geom::degrees) const
Return our contents in a Point2D object.
Definition: Coord.cc:477
Generated on Wed Sep 16 2015 13:35:51 for LSSTApplications by   doxygen 1.8.5