LSSTApplications  20.0.0
LSSTDataManagementBasePackage
Namespaces | Classes | Typedefs | Functions
lsst::meas::astrom Namespace Reference

Namespaces

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

Classes

struct  MatchOptimisticBControl
 
class  OutlierRejectionControl
 Control object for outlier rejection in ScaledPolynomialTransformFitter. 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...
 
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...
 
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...
 

Typedefs

typedef std::vector< RecordProxyProxyVector
 

Functions

template<typename MatchT >
afw::math::Statistics makeMatchStatistics (std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
 Compute statistics of the distance field of a match list. More...
 
template<typename MatchT >
afw::math::Statistics makeMatchStatisticsInPixels (afw::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. More...
 
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. More...
 
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. More...
 
PolynomialTransform compose (geom::AffineTransform const &t1, PolynomialTransform const &t2)
 Return a PolynomialTransform that is equivalent to the composition t1(t2()) More...
 
PolynomialTransform compose (PolynomialTransform const &t1, geom::AffineTransform const &t2)
 Return a PolynomialTransform that is equivalent to the composition t1(t2()) More...
 
std::shared_ptr< afw::geom::SkyWcsmakeWcs (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. More...
 
std::shared_ptr< afw::geom::SkyWcstransformWcsPixels (afw::geom::SkyWcs const &wcs, geom::AffineTransform const &s)
 Create a new SkyWcs whose pixel coordinate system has been transformed via an affine transform. More...
 
std::shared_ptr< afw::geom::SkyWcsrotateWcsPixelsBy90 (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. More...
 
 PYBIND11_MODULE (makeMatchStatistics, mod)
 
 PYBIND11_MODULE (matchOptimisticB, mod)
 
 PYBIND11_MODULE (polynomialTransform, mod)
 
 PYBIND11_MODULE (scaledPolynomialTransformFitter, mod)
 
 PYBIND11_MODULE (sipTransform, 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)
 

Typedef Documentation

◆ ProxyVector

Definition at line 52 of file matchOptimisticB.h.

Function Documentation

◆ 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.

214  {
215  typedef geom::AffineTransform AT;
216  PolynomialTransform result(t2.getOrder());
217 
218  result._xCoeffs.deep() = t2._xCoeffs * t1[AT::XX] + t2._yCoeffs * t1[AT::XY];
219  result._yCoeffs.deep() = t2._xCoeffs * t1[AT::YX] + t2._yCoeffs * t1[AT::YY];
220  result._xCoeffs(0, 0) += t1[AT::X];
221  result._yCoeffs(0, 0) += t1[AT::Y];
222  return result;
223 }

◆ 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.

225  {
226  typedef geom::AffineTransform AT;
227  int const order = t1.getOrder();
228  if (order < 1) {
229  PolynomialTransform t1a(1);
230  t1a._xCoeffs(0, 0) = t1._xCoeffs(0, 0);
231  t1a._yCoeffs(0, 0) = t1._yCoeffs(0, 0);
232  return compose(t1a, t2);
233  }
234  detail::BinomialMatrix binomial(order);
235  // For each of these, (e.g.) a[n] == pow(a, n)
236  auto const t2u = detail::computePowers(t2[AT::X], order);
237  auto const t2v = detail::computePowers(t2[AT::Y], order);
238  auto const t2uu = detail::computePowers(t2[AT::XX], order);
239  auto const t2uv = detail::computePowers(t2[AT::XY], order);
240  auto const t2vu = detail::computePowers(t2[AT::YX], order);
241  auto const t2vv = detail::computePowers(t2[AT::YY], order);
242  PolynomialTransform result(order);
243  for (int p = 0; p <= order; ++p) {
244  for (int m = 0; m <= p; ++m) {
245  for (int j = 0; j <= m; ++j) {
246  for (int q = 0; p + q <= order; ++q) {
247  for (int n = 0; n <= q; ++n) {
248  for (int k = 0; k <= n; ++k) {
249  double z = binomial(p, m) * t2u[p - m] * binomial(m, j) * t2uu[j] * t2uv[m - j] *
250  binomial(q, n) * t2v[q - n] * binomial(n, k) * t2vu[k] * t2vv[n - k];
251  result._xCoeffs(j + k, m + n - j - k) += t1._xCoeffs(p, q) * z;
252  result._yCoeffs(j + k, m + n - j - k) += t1._yCoeffs(p, q) * z;
253  } // k
254  } // n
255  } // q
256  } // j
257  } // m
258  } // p
259  return result;
260 }

◆ 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]matchListlist 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]flagswhat to calculate; OR constants such as lsst::afw::math::MEAN, MEANCLIP, STDDEV, MEDIAN, defined in lsst/afw/math/Statitics.h's Property enum
[in]sctrlstatistics configuration

Definition at line 33 of file makeMatchStatistics.cc.

34  {
35  if (matchList.empty()) {
36  throw LSST_EXCEPT(pexExcept::RuntimeError, "matchList is empty");
37  }
39  val.reserve(matchList.size());
40 
41  for (auto const& match : matchList) {
42  val.push_back(match.distance);
43  }
44  return afw::math::makeStatistics(val, flags, sctrl);
45 }

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

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

◆ 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]wcsWCS describing pixel to sky transformation
[in]matchListlist 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]flagswhat to calculate; OR constants such as lsst::afw::math::MEAN, MEANCLIP, STDDEV, MEDIAN, defined in lsst/afw/math/Statitics.h's Property enum
[in]sctrlstatistics configuration

Definition at line 48 of file makeMatchStatistics.cc.

50  {
51  if (matchList.empty()) {
52  throw LSST_EXCEPT(pexExcept::RuntimeError, "matchList is empty");
53  }
55  val.reserve(matchList.size());
56 
57  for (auto const& match : matchList) {
58  auto refPtr = match.first;
59  auto srcPtr = match.second;
60  auto srcX = srcPtr->getX();
61  auto srcY = srcPtr->getY();
62  auto refPos = wcs.skyToPixel(refPtr->getCoord());
63  auto refX = refPos[0];
64  auto refY = refPos[1];
65  val.push_back(::hypot(srcX - refX, srcY - refY));
66  }
67  return afw::math::makeStatistics(val, flags, sctrl);
68 }

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

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

◆ 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]wcsWCS describing pixel to sky transformation
[in]matchListlist 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]flagswhat to calculate; OR constants such as lsst::afw::math::MEAN, MEANCLIP, STDDEV, MEDIAN, defined in lsst/afw/math/Statitics.h's Property enum
[in]sctrlstatistics configuration

Definition at line 71 of file makeMatchStatistics.cc.

73  {
74  if (matchList.empty()) {
75  throw LSST_EXCEPT(pexExcept::RuntimeError, "matchList is empty");
76  }
78  val.reserve(matchList.size());
79 
80  for (auto const& match : matchList) {
81  auto refPtr = match.first;
82  auto srcPtr = match.second;
83  auto refCoord = refPtr->getCoord();
84  auto srcCoord = wcs.pixelToSky(srcPtr->getCentroid());
85  auto angSep = refCoord.separation(srcCoord);
86  val.push_back(angSep.asRadians());
87  }
88  return afw::math::makeStatistics(val, flags, sctrl);
89 }

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

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

◆ 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.

448  {
449  auto coordKey = afwTable::CoordKey(posRefCat.getSchema()["coord"]);
450  ProxyVector r;
451  r.reserve(posRefCat.size());
452  for (afwTable::SimpleCatalog::const_iterator posRefPtr = posRefCat.begin(); posRefPtr != posRefCat.end();
453  ++posRefPtr) {
454  r.push_back(RecordProxy(posRefPtr, tanWcs.skyToPixel(posRefPtr->get(coordKey))));
455  }
456  return r;
457 }

◆ 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.

437  {
438  ProxyVector r;
439  r.reserve(sourceCat.size());
440  for (afwTable::SourceCatalog::const_iterator sourcePtr = sourceCat.begin(); sourcePtr != sourceCat.end();
441  ++sourcePtr) {
442  r.push_back(
443  RecordProxy(sourcePtr, tanWcs.skyToPixel(distortedWcs.pixelToSky(sourcePtr->getCentroid()))));
444  }
445  return r;
446 }

◆ 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]sipForwardMapping from pixel coordinates to intermediate world coordinates.
[in]sipReverseMapping from intermediate world coordinates to pixel coordinates.
[in]skyOriginICRS position of the gnomonic projection that maps sky coordinates to intermediate world coordinates (CRVAL).
Exceptions
pex::exceptions::InvalidParameterErrorif the forward and reverse SIP transforms have different CRPIX values or CD matrices.

Definition at line 148 of file SipTransform.cc.

150  {
151  if (!sipForward.getPixelOrigin().asEigen().isApprox(sipReverse.getPixelOrigin().asEigen())) {
152  std::ostringstream oss;
153  oss << "SIP forward and reverse transforms have inconsistent CRPIX: " << sipForward.getPixelOrigin()
154  << " != " << sipReverse.getPixelOrigin();
155  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, oss.str());
156  }
157  if (!sipForward.getCdMatrix().getMatrix().isApprox(sipReverse.getCdMatrix().getMatrix())) {
158  std::ostringstream oss;
159  oss << "SIP forward and reverse transforms have inconsistent CD matrix: " << sipForward.getCdMatrix()
160  << "\n!=\n"
161  << sipReverse.getCdMatrix();
162  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, oss.str());
163  }
164  Eigen::MatrixXd sipA(ndarray::asEigenMatrix(sipForward.getPoly().getXCoeffs()));
165  Eigen::MatrixXd sipB(ndarray::asEigenMatrix(sipForward.getPoly().getYCoeffs()));
166  Eigen::MatrixXd sipAP(ndarray::asEigenMatrix(sipReverse.getPoly().getXCoeffs()));
167  Eigen::MatrixXd sipBP(ndarray::asEigenMatrix(sipReverse.getPoly().getYCoeffs()));
168 
169  return afw::geom::makeTanSipWcs(sipForward.getPixelOrigin(), skyOrigin,
170  sipForward.getCdMatrix().getMatrix(), sipA, sipB, sipAP, sipBP);
171 }

◆ 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]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 459 of file matchOptimisticB.cc.

463  {
464  control.validate();
465  if (posRefCat.empty()) {
466  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "no entries in posRefCat");
467  }
468  if (sourceCat.empty()) {
469  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "no entries in sourceCat");
470  }
471  if (posRefBegInd < 0) {
472  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd < 0");
473  }
474  if (static_cast<size_t>(posRefBegInd) >= posRefCat.size()) {
475  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd too big");
476  }
477  double const maxRotationRad = geom::degToRad(control.maxRotationDeg);
478 
479  // Create an undistorted Wcs to project everything with
480  // We'll anchor it at the center of the area.
481  geom::Extent2D srcCenter(0, 0);
482  for (auto iter = sourceCat.begin(); iter != sourceCat.end(); ++iter) {
483  srcCenter += geom::Extent2D(iter->getCentroid());
484  }
485  srcCenter /= sourceCat.size();
486 
487  sphgeom::Vector3d refCenter(0, 0, 0);
488  for (auto iter = posRefCat.begin(); iter != posRefCat.end(); ++iter) {
489  refCenter += iter->getCoord().getVector();
490  }
491  refCenter /= posRefCat.size();
492 
493  auto tanWcs =
494  afw::geom::makeSkyWcs(geom::Point2D(srcCenter), geom::SpherePoint(refCenter), wcs.getCdMatrix());
495 
496  ProxyVector posRefProxyCat = makeProxies(posRefCat, *tanWcs);
497  ProxyVector sourceProxyCat = makeProxies(sourceCat, wcs, *tanWcs);
498 
499  // sourceSubCat contains at most the numBrightStars brightest sources, sorted by decreasing flux
500  ProxyVector sourceSubCat =
501  selectPoint(sourceProxyCat, sourceCat.getSchema().find<double>(control.sourceFluxField).key,
502  control.numBrightStars);
503 
504  // posRefSubCat skips the initial posRefBegInd brightest reference objects and contains
505  // at most the next len(sourceSubCat) + 25 brightest reference objects, sorted by decreasing flux
506  ProxyVector posRefSubCat =
507  selectPoint(posRefProxyCat, posRefCat.getSchema().find<double>(control.refFluxField).key,
508  sourceSubCat.size() + 25, posRefBegInd);
509  if (verbose) {
510  std::cout << "Catalog sizes: " << sourceSubCat.size() << " " << posRefSubCat.size() << std::endl;
511  }
512 
513  // Construct a list of pairs of position reference stars sorted by increasing separation
514  std::vector<ProxyPair> posRefPairList;
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++) {
518  posRefPairList.push_back(ProxyPair(posRefSubCat[i], posRefSubCat[j]));
519  }
520  }
521  std::sort(posRefPairList.begin(), posRefPairList.end(), cmpPair);
522 
523  // Construct a list of pairs of sources sorted by increasing separation
524  std::vector<ProxyPair> sourcePairList;
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]));
529  }
530  }
531  std::sort(sourcePairList.begin(), sourcePairList.end(), cmpPair);
532 
534  afwTable::ReferenceMatchVector matPairSave;
536 
537  size_t const fullShapeSize = control.numPointsForShape - 1; // Max size of shape array
538  for (size_t ii = 0; ii < sourcePairList.size(); ii++) {
539  ProxyPair p = sourcePairList[ii];
540 
542  searchPair(posRefPairList, p, control.matchingAllowancePix, maxRotationRad);
543 
544  // If candidate pairs are found
545  if (q.size() != 0) {
546  std::vector<ProxyPair> srcMatPair;
547  std::vector<ProxyPair> catMatPair;
548 
549  // Go through candidate pairs
550  for (size_t l = 0; l < q.size(); l++) {
551  // sign matters, so don't use deltaAng; no need to wrap because
552  // the result is used with deltaAng later
553  double dpa = p.pa - q[l].pa;
554 
555  srcMatPair.clear();
556  catMatPair.clear();
557 
558  srcMatPair.push_back(p);
559  catMatPair.push_back(q[l]);
560 
561  if (verbose) {
562  std::cout << "p dist: " << p.distance << " pa: " << p.pa << std::endl;
563  std::cout << "q dist: " << q[l].distance << " pa: " << q[l].pa << std::endl;
564  }
565 
566  for (size_t k = 0; k < sourceSubCat.size(); k++) {
567  if (p.first == sourceSubCat[k] || p.second == sourceSubCat[k]) continue;
568 
569  ProxyPair pp(p.first, sourceSubCat[k]);
570 
571  std::vector<ProxyPair>::iterator r = searchPair3(
572  posRefPairList, pp, q[l], control.matchingAllowancePix, dpa, maxRotationRad);
573  if (r != posRefPairList.end()) {
574  srcMatPair.push_back(pp);
575  catMatPair.push_back(*r);
576  if (verbose) {
577  std::cout << " p dist: " << pp.distance << " pa: " << pp.pa << std::endl;
578  std::cout << " r dist: " << (*r).distance << " pa: " << (*r).pa << std::endl;
579  }
580  if (srcMatPair.size() == fullShapeSize) {
581  break;
582  }
583  }
584  }
585 
586  bool goodMatch = false;
587  if (srcMatPair.size() == fullShapeSize) {
588  goodMatch = true;
589  for (size_t k = 1; k < catMatPair.size(); k++) {
590  if (catMatPair[0].first != catMatPair[k].first) {
591  goodMatch = false;
592  break;
593  }
594  }
595  }
596 
597  if (goodMatch && srcMatPair.size() == fullShapeSize) {
598  ProxyVector srcMat;
599  ProxyVector catMat;
600 
601  srcMat.push_back(srcMatPair[0].first);
602  catMat.push_back(catMatPair[0].first);
603  for (size_t k = 0; k < srcMatPair.size(); k++) {
604  srcMat.push_back(srcMatPair[k].second);
605  catMat.push_back(catMatPair[k].second);
606  }
607 
608  boost::shared_array<double> coeff = polyfit(1, srcMat, catMat);
609 
610  if (verbose) {
611  for (size_t k = 0; k < srcMat.size(); k++) {
612  std::cout << "circle(" << srcMat[k].getX() << "," << srcMat[k].getY()
613  << ",10) # color=green" << std::endl;
614  std::cout << "circle(" << catMat[k].getX() << "," << catMat[k].getY()
615  << ",10) # color=red" << std::endl;
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;
622  }
623  }
624 
625  double a = coeff[1];
626  double b = coeff[2];
627  double c = coeff[4];
628  double d = coeff[5];
629  geom::Angle const theta(std::acos((a * b + c * d) /
630  (std::sqrt(a * a + c * c) * std::sqrt(b * b + d * d))),
631  geom::radians);
632  if (verbose) {
633  std::cout << "Linear fit from match:" << std::endl;
634  std::cout << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
635  std::cout << coeff[3] << " " << coeff[4] << " " << coeff[5] << std::endl;
636  std::cout << "Determinant (max " << control.maxDeterminant << "): ";
637  std::cout << coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1. << std::endl;
638  std::cout << "Angle between axes (deg; allowed 90 +/- ";
639  std::cout << control.allowedNonperpDeg << "): ";
640  std::cout << theta.asDegrees() << std::endl;
641  }
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) {
646  if (verbose) {
647  std::cout << "Bad; continuing" << std::endl;
648  }
649  continue;
650  } else {
651  double x0, y0, x1, y1;
652  int num = 0;
653  srcMat.clear();
654  catMat.clear();
655  for (size_t i = 0; i < sourceSubCat.size(); i++) {
656  x0 = sourceSubCat[i].getX();
657  y0 = sourceSubCat[i].getY();
658  transform(1, coeff, x0, y0, &x1, &y1);
659  auto refObjDist =
660  searchNearestPoint(posRefSubCat, x1, y1, control.matchingAllowancePix);
661  if (refObjDist.first != posRefSubCat.end()) {
662  num++;
663  srcMat.push_back(sourceSubCat[i]);
664  catMat.push_back(*refObjDist.first);
665  if (verbose) {
666  std::cout << "Match: " << x0 << "," << y0 << " --> " << x1 << "," << y1
667  << " <==> " << refObjDist.first->getX() << ","
668  << refObjDist.first->getY() << std::endl;
669  }
670  }
671  }
672  if (num <= control.numPointsForShape) {
673  // Can get matrix = 0,0,0,0; everything matches a single catalog object
674  if (verbose) {
675  std::cout << "Insufficient initial matches; continuing" << std::endl;
676  }
677  continue;
678  }
679  coeff = polyfit(1, srcMat, catMat);
680  if (verbose) {
681  std::cout << "Coefficients from initial matching:" << std::endl;
682  for (size_t i = 0; i < 6; ++i) {
683  std::cout << coeff[i] << " ";
684  }
685  std::cout << std::endl;
686  }
687 
688  matPair = FinalVerify(coeff, posRefProxyCat, sourceProxyCat,
689  control.matchingAllowancePix, verbose);
690  if (verbose) {
691  std::cout << "Number of matches: " << matPair.size() << " vs "
692  << control.minMatchedPairs << std::endl;
693  }
694  if (matPair.size() <= static_cast<std::size_t>(control.minMatchedPairs)) {
695  if (verbose) {
696  std::cout << "Insufficient final matches; continuing" << std::endl;
697  }
698  if (matPair.size() > matPairSave.size()) {
699  matPairSave = matPair;
700  }
701  continue;
702  } else {
703  if (verbose) {
704  std::cout << "Finish" << std::endl;
705  }
706  matPairCand.push_back(matPair);
707  if (matPairCand.size() == 3) {
708  goto END;
709  }
710  }
711  }
712  }
713  }
714  }
715  }
716 
717 END:
718  if (matPairCand.size() == 0) {
719  return matPairSave;
720  } else {
721  size_t nmatch = matPairCand[0].size();
722  afwTable::ReferenceMatchVector matPairRet = matPairCand[0];
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];
727  }
728  }
729  return matPairRet;
730  }
731 }

◆ PYBIND11_MODULE() [1/5]

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

Definition at line 47 of file makeMatchStatistics.cc.

47  {
48  py::module::import("lsst.afw.math");
49 
50 
51  declareMakeMatchStatistics<afw::table::ReferenceMatch>(mod);
52  declareMakeMatchStatistics<afw::table::SourceMatch>(mod);
53 }

◆ PYBIND11_MODULE() [2/5]

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

Definition at line 91 of file matchOptimisticB.cc.

91  {
92  declareRecordProxy(mod);
93  declareProxyPair(mod);
94  declareMatchOptimisticBControl(mod);
95 
96  mod.def("makeProxies",
97  (ProxyVector(*)(afw::table::SourceCatalog const &, afw::geom::SkyWcs const &,
98  afw::geom::SkyWcs const &)) &
100  "sourceCat"_a, "distortedWcs"_a, "tanWcs"_a);
101  mod.def("makeProxies",
102  (ProxyVector(*)(afw::table::SimpleCatalog const &, afw::geom::SkyWcs const &)) & makeProxies,
103  "posRefCat"_a, "tanWcs"_a);
104 
105  mod.def("matchOptimisticB", &matchOptimisticB, "posRefCat"_a, "sourceCat"_a, "control"_a, "wcs"_a,
106  "posRefBegInd"_a = 0, "verbose"_a = false);
107 }

◆ PYBIND11_MODULE() [3/5]

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

Definition at line 102 of file polynomialTransform.cc.

102  {
103  declarePolynomialTransform(mod);
104  declareScaledPolynomialTransform(mod);
105 
106  mod.def("compose",
107  (PolynomialTransform(*)(geom::AffineTransform const &, PolynomialTransform const &)) & compose,
108  "t1"_a, "t2"_a);
109  mod.def("compose",
110  (PolynomialTransform(*)(PolynomialTransform const &, geom::AffineTransform const &)) & compose,
111  "t1"_a, "t2"_a);
112 }

◆ PYBIND11_MODULE() [4/5]

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

Definition at line 71 of file scaledPolynomialTransformFitter.cc.

71  {
72  declareOutlierRejectionControl(mod);
73  declareScaledPolynomialTransformFitter(mod);
74 }

◆ PYBIND11_MODULE() [5/5]

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

Definition at line 104 of file sipTransform.cc.

104  {
105  declareSipTransformBase(mod);
106  declareSipForwardTransform(mod);
107  declareSipReverseTransform(mod);
108 
109  mod.def("makeWcs", makeWcs, "sipForward"_a, "sipReverse"_a, "skyOrigin"_a);
110  mod.def("transformWcsPixels", transformWcsPixels, "wcs"_a, "s"_a);
111  mod.def("rotateWcsPixelsBy90", rotateWcsPixelsBy90, "wcs"_a, "nQuarter"_a, "dimensions"_a);
112 }

◆ 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]wcsOriginal SkyWcs to be rotated.
[in]nQuarterNumber of 90 degree rotations (positive is counterclockwise).
[in]dimensionsWidth and height of the image.

Definition at line 179 of file SipTransform.cc.

180  {
181  geom::Extent2D offset;
182  switch (nQuarter % 4) {
183  case 0:
184  offset = geom::Extent2D(0, 0);
185  break;
186  case 1:
187  offset = geom::Extent2D(dimensions.getY() - 1, 0);
188  break;
189  case 2:
190  offset = geom::Extent2D(dimensions - geom::Extent2I(1, 1));
191  break;
192  case 3:
193  offset = geom::Extent2D(0, dimensions.getX() - 1);
194  break;
195  }
196  auto rot = geom::LinearTransform::makeRotation(nQuarter * 90.0 * geom::degrees);
197  return transformWcsPixels(wcs, geom::AffineTransform(rot, offset));
198 }

◆ 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]wcsOriginal SkyWcs object.
[in]sAffineTransform 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.

174  {
175  auto affineTransform22 = afw::geom::makeTransform(s);
176  return afw::geom::makeModifiedWcs(*affineTransform22->inverted(), wcs, true);
177 }
lsst::geom::degrees
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
lsst::meas::astrom::PolynomialTransform::compose
compose
lsst::meas::astrom::ProxyPair::second
RecordProxy second
Definition: matchOptimisticB.h:61
lsst::geom::LinearTransform::makeRotation
static LinearTransform makeRotation(Angle t) noexcept
Definition: LinearTransform.h:102
lsst::afw::table._match.second
second
Definition: _match.py:78
std::fabs
T fabs(T... args)
std::vector::reserve
T reserve(T... args)
wcs
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
lsst::afw::table::CoordKey
A FunctorKey used to get or set celestial coordinates from a pair of lsst::geom::Angle keys.
Definition: aggregates.h:210
std::vector< double >
lsst::meas::astrom::makeProxies
ProxyVector makeProxies(afw::table::SourceCatalog const &sourceCat, afw::geom::SkyWcs const &distortedWcs, afw::geom::SkyWcs const &tanWcs)
Definition: matchOptimisticB.cc:436
std::vector::size
T size(T... args)
lsst::afw::table._match.first
first
Definition: _match.py:76
lsst::meas::astrom::transformWcsPixels
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.
Definition: SipTransform.cc:173
lsst::afw::table::SortedCatalogT::const_iterator
Base::const_iterator const_iterator
Definition: SortedCatalog.h:50
lsst::afw::geom::makeModifiedWcs
std::shared_ptr< SkyWcs > makeModifiedWcs(TransformPoint2ToPoint2 const &pixelTransform, SkyWcs const &wcs, bool modifyActualPixels)
Create a new SkyWcs whose pixels are transformed by pixelTransform, as described below.
Definition: SkyWcs.cc:491
val
ImageT val
Definition: CR.cc:146
lsst::afw::image::X
@ X
Definition: ImageUtils.h:36
std::sort
T sort(T... args)
lsst::meas::astrom::ProxyPair::distance
double distance
Definition: matchOptimisticB.h:62
std::sqrt
T sqrt(T... args)
lsst::afw::table::SimpleCatalog
SortedCatalogT< SimpleRecord > SimpleCatalog
Definition: fwd.h:79
std::vector::clear
T clear(T... args)
lsst::afw::table::SourceCatalog
SortedCatalogT< SourceRecord > SourceCatalog
Definition: fwd.h:85
std::vector::push_back
T push_back(T... args)
lsst::meas::astrom::ProxyPair
Definition: matchOptimisticB.h:59
lsst::afw::math::makeStatistics
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:354
lsst::geom::AffineTransform
An affine coordinate transformation consisting of a linear transformation and an offset.
Definition: AffineTransform.h:75
std::hypot
T hypot(T... args)
coeff
table::Key< table::Array< double > > coeff
Definition: PsfexPsf.cc:362
lsst::afw::geom::makeTanSipWcs
std::shared_ptr< SkyWcs > makeTanSipWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Construct a TAN-SIP SkyWcs with forward SIP distortion terms and an iterative inverse.
Definition: SkyWcs.cc:543
std::cout
z
double z
Definition: Match.cc:44
lsst::afw::image::Y
@ Y
Definition: ImageUtils.h:36
lsst::meas::astrom::matchOptimisticB
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.
Definition: matchOptimisticB.cc:459
dimensions
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
lsst::meas::astrom::ProxyPair::pa
double pa
Definition: matchOptimisticB.h:63
pixel
table::PointKey< int > pixel
Definition: DeltaFunctionKernel.cc:101
result
py::object result
Definition: _schema.cc:429
b
table::Key< int > b
Definition: TransmissionCurve.cc:467
std::acos
T acos(T... args)
lsst::geom::degToRad
constexpr double degToRad(double x) noexcept
Definition: Angle.h:51
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst::meas::astrom::ProxyPair::first
RecordProxy first
Definition: matchOptimisticB.h:60
std::ostringstream
STL class.
lsst::meas::astrom::RecordProxy
A wrapper around a SimpleRecord or SourceRecord that allows us to record a pixel position in a way th...
Definition: matchOptimisticB.h:22
lsst::meas::astrom::makeWcs
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.
Definition: SipTransform.cc:148
std::endl
T endl(T... args)
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
a
table::Key< int > a
Definition: TransmissionCurve.cc:466
std::vector::begin
T begin(T... args)
lsst::geom::Point< double, 2 >
lsst::geom::Angle
A class representing an angle.
Definition: Angle.h:127
lsst::geom::radians
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
lsst::afw::geom::makeSkyWcs
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:526
lsst::geom::Extent2D
Extent< double, 2 > Extent2D
Definition: Extent.h:400
std::vector::empty
T empty(T... args)
transform
table::Key< int > transform
Definition: TransformMap.cc:299
std::ostringstream::str
T str(T... args)
std::size_t
std::vector::end
T end(T... args)
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::meas::astrom::rotateWcsPixelsBy90
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 cente...
Definition: SipTransform.cc:179
lsst::afw::geom::makeTransform
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
Definition: transformFactory.cc:154
lsst::geom::Extent< double, 2 >
astshim.fitsChanContinued.iter
def iter(self)
Definition: fitsChanContinued.py:88
lsst::pex::exceptions::RuntimeError
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
m
int m
Definition: SpanSet.cc:49
lsst::meas::astrom::ProxyVector
std::vector< RecordProxy > ProxyVector
Definition: matchOptimisticB.h:52