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
matchOptimisticB.cc
Go to the documentation of this file.
1 #include <algorithm>
2 #include <cmath>
3 #include <fstream>
4 #include <iostream>
5 #include <utility>
6 #include <vector>
7 
8 #include "boost/scoped_array.hpp"
9 #include "boost/shared_array.hpp"
10 #include "boost/multi_index_container.hpp"
11 #include "boost/multi_index/sequenced_index.hpp"
12 #include "boost/multi_index/ordered_index.hpp"
13 #include "boost/multi_index/global_fun.hpp"
14 
15 #include "gsl/gsl_linalg.h"
16 
17 #include "lsst/utils/ieee.h"
18 #include "lsst/pex/exceptions.h"
19 #include "lsst/afw/image/Wcs.h"
21 #include "lsst/afw/geom/Angle.h"
23 
24 namespace pexExcept = lsst::pex::exceptions;
25 namespace afwTable = lsst::afw::table;
26 
27 namespace {
28  using namespace lsst::meas::astrom;
29 
30  afwTable::RecordId getRefId(ProxyPair const &proxyPair) {
31  return proxyPair.first.record->getId();
32  }
33 
34  // a collection of ReferenceMatch indexed by insertion order and by reference object ID
35  struct refIdTag{}; // tag for retrieval by reference object ID
36  struct insertionOrderTag{}; // tag for retrieval by insertion order
37  typedef boost::multi_index_container<
38  ProxyPair,
39  boost::multi_index::indexed_by< // index by insertion order, so brightest sources first
40  boost::multi_index::sequenced<
41  boost::multi_index::tag<insertionOrderTag>
42  >,
43  boost::multi_index::ordered_unique< // index by reference object ID
44  boost::multi_index::tag<refIdTag>,
45  boost::multi_index::global_fun<
46  ProxyPair const &, afwTable::RecordId , &getRefId
47  > >
48  >
49  > MultiIndexedProxyPairList;
50 
51  // Algorithm is based on V.Tabur 2007, PASA, 24, 189-198
52  // "Fast Algorithms for Matching CCD Images to a Stellar Catalogue"
53 
60  inline double absDeltaAngle(double ang1, double ang2) {
61  return std::fmod(std::fabs(ang1 - ang2), M_PI*2);
62  }
63 
64  bool cmpPair(ProxyPair const &a, ProxyPair const &b) {
65  return a.distance > b.distance;
66  }
67 
73  struct CompareProxyFlux {
74 
75  bool operator()(RecordProxy const & a, RecordProxy const & b) const {
76  double aFlux = a.record->get(key);
77  double bFlux = b.record->get(key);
78  if (lsst::utils::isnan(aFlux)) {
79  aFlux = 0.0;
80  }
81  if (lsst::utils::isnan(bFlux)) {
82  bFlux = 0.0;
83  }
84  return aFlux > bFlux;
85  }
86 
88  };
89 
100  ProxyVector selectPoint(
101  ProxyVector const &a,
102  afwTable::Key<double> const & key,
103  std::size_t num,
104  std::size_t startInd=0
105  ) {
106  if (startInd >= a.size()) {
107  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "startInd too big");
108  }
109  CompareProxyFlux cmp = {key};
110  ProxyVector b(a);
111  std::sort(b.begin(), b.end(), cmp);
112  std::size_t const endInd = std::min(startInd + num, b.size());
113  return ProxyVector(b.begin() + startInd, b.begin() + endInd);
114  }
115 
116  std::vector<ProxyPair> searchPair(
117  std::vector<ProxyPair> const &a,
118  ProxyPair const &p,
119  double e,
120  double e_dpa
121  ) {
122  std::vector<ProxyPair> v;
123 
124  for (size_t i = 0; i < a.size(); i++) {
125  double dd = std::fabs(a[i].distance - p.distance);
126  double dpa = absDeltaAngle(a[i].pa, p.pa);
127  if (dd < e && dpa < e_dpa) {
128  v.push_back(a[i]);
129  }
130  }
131 
132  return v;
133  }
134 
135  std::vector<ProxyPair>::iterator searchPair3(
136  std::vector<ProxyPair> &a,
137  ProxyPair const &p,
138  ProxyPair const &q,
139  double e,
140  double dpa,
141  double e_dpa = 0.02
142  ) {
143  std::vector<ProxyPair>::iterator idx = a.end();
144  double dd_min = 1.E+10;
145  //double dpa_min = e_dpa;
146 
147  for (auto i = a.begin(); i != a.end(); ++i) {
148  double dd = std::fabs(i->distance - p.distance);
149  #if 1
150  if (dd < e &&
151  absDeltaAngle(p.pa, i->pa - dpa) < e_dpa &&
152  dd < dd_min &&
153  (i->first == q.first)) {
154  dd_min = dd;
155  idx = i;
156  }
157  #else
158  if (dd < e &&
159  absDeltaAngle(p.pa, i->pa - dpa) < dpa_min) {
160  dpa_min = std::fabs(p.pa - i->pa - dpa);
161  idx = i;
162  }
163  #endif
164  }
165 
166  return idx;
167  }
168 
169  void transform(
170  int order,
171  boost::shared_array<double> const & coeff,
172  double x,
173  double y,
174  double *xn,
175  double *yn
176  ) {
177  int ncoeff = (order + 1) * (order + 2) / 2;
178  *xn = 0.0;
179  *yn = 0.0;
180  int n = 0;
181  for (int i = 0; i <= order; i++) {
182  for (int k = 0; k <= i; k++) {
183  int j = i - k;
184  *xn += coeff[n] * pow(x, j) * pow(y, k);
185  *yn += coeff[n+ncoeff] * pow(x, j) * pow(y, k);
186  n++;
187  }
188  }
189  }
190 
191  boost::shared_array<double> polyfit(
192  int order,
193  ProxyVector const &img,
194  ProxyVector const &posRefCat
195  ) {
196  int ncoeff = (order + 1) * (order + 2) / 2;
197  boost::scoped_array<int> xorder(new int[ncoeff]);
198  boost::scoped_array<int> yorder(new int[ncoeff]);
199 
200  int n = 0;
201  for (int i = 0; i <= order; i++) {
202  for (int k = 0; k <= i; k++) {
203  int j = i - k;
204  xorder[n] = j;
205  yorder[n] = k;
206  n++;
207  }
208  }
209 
210  boost::scoped_array<int> flag(new int[img.size()]);
211  for (size_t k = 0; k < img.size(); k++) {
212  flag[k] = 1;
213  }
214 
215  boost::scoped_array<double> a_data(new double[ncoeff*ncoeff]);
216  boost::scoped_array<double> b_data(new double[ncoeff]);
217  boost::scoped_array<double> c_data(new double[ncoeff]);
218 
219  boost::shared_array<double> coeff(new double[ncoeff*2]);
220 
221  for (int loop = 0; loop < 1; loop++) {
222  for (int i = 0; i < ncoeff; i++) {
223  for (int j = 0; j < ncoeff; j++) {
224  a_data[i*ncoeff+j] = 0.0;
225  for (size_t k = 0; k < img.size(); k++) {
226  if (flag[k] == 1) {
227  a_data[i*ncoeff+j] += pow(img[k].getX(), xorder[i]) *
228  pow(img[k].getY(), yorder[i]) *
229  pow(img[k].getX(), xorder[j]) *
230  pow(img[k].getY(), yorder[j]);
231  }
232  }
233  }
234  b_data[i] = c_data[i] = 0.0;
235  for (unsigned int k = 0; k < img.size(); k++) {
236  if (flag[k] == 1) {
237  b_data[i] += pow(img[k].getX(), xorder[i]) *
238  pow(img[k].getY(), yorder[i]) *
239  posRefCat[k].getX();
240  c_data[i] += pow(img[k].getX(), xorder[i]) *
241  pow(img[k].getY(), yorder[i]) *
242  posRefCat[k].getY();
243  }
244  }
245  }
246 
247  gsl_matrix_view a = gsl_matrix_view_array(a_data.get(), ncoeff, ncoeff);
248  gsl_vector_view b = gsl_vector_view_array(b_data.get(), ncoeff);
249  gsl_vector_view c = gsl_vector_view_array(c_data.get(), ncoeff);
250 
251  boost::shared_ptr<gsl_vector> x(gsl_vector_alloc(ncoeff), gsl_vector_free);
252  boost::shared_ptr<gsl_vector> y(gsl_vector_alloc(ncoeff), gsl_vector_free);
253 
254  int s;
255 
256  boost::shared_ptr<gsl_permutation> p(gsl_permutation_alloc(ncoeff), gsl_permutation_free);
257 
258  gsl_linalg_LU_decomp(&a.matrix, p.get(), &s);
259  gsl_linalg_LU_solve(&a.matrix, p.get(), &b.vector, x.get());
260  gsl_linalg_LU_solve(&a.matrix, p.get(), &c.vector, y.get());
261 
262  for (int i = 0; i < ncoeff; i++) {
263  coeff[i] = x->data[i];
264  coeff[i+ncoeff] = y->data[i];
265  }
266 
267  double S, Sx, Sy, Sxx, Syy;
268  S = Sx = Sy = Sxx = Syy = 0.0;
269  for (size_t k = 0; k < img.size(); k++) {
270  if (flag[k] == 1) {
271  double x0 = img[k].getX();
272  double y0 = img[k].getY();
273  double x1, y1;
274  transform(order, coeff, x0, y0, &x1, &y1);
275  S += 1.;
276  Sx += (x1 - posRefCat[k].getX());
277  Sxx += (x1 - posRefCat[k].getX()) * (x1 - posRefCat[k].getX());
278  Sy += (y1 - posRefCat[k].getY());
279  Syy += (y1 - posRefCat[k].getY()) * (y1 - posRefCat[k].getY());
280  }
281  }
282  double x_sig = std::sqrt((Sxx - Sx * Sx / S) / S);
283  double y_sig = std::sqrt((Syy - Sy * Sy / S) / S);
284  //std::cout << x_sig << " " << y_sig << std::endl;
285 
286  for (size_t k = 0; k < img.size(); k++) {
287  double x0 = img[k].getX();
288  double y0 = img[k].getY();
289  double x1, y1;
290  transform(order, coeff, x0, y0, &x1, &y1);
291  if (std::fabs(x1-posRefCat[k].getX()) > 2. * x_sig ||
292  std::fabs(y1-posRefCat[k].getY()) > 2. * y_sig) {
293  flag[k] = 0;
294  }
295  }
296 
297  }
298 
299  return coeff;
300  }
301 
314  std::pair<ProxyVector::const_iterator, double> searchNearestPoint(
315  ProxyVector const &posRefCat,
316  double x,
317  double y,
318  double matchingAllowancePix
319  ) {
320  auto minDistSq = matchingAllowancePix*matchingAllowancePix;
321  auto foundPtr = posRefCat.end();
322  for (auto posRefPtr = posRefCat.begin(); posRefPtr != posRefCat.end(); ++posRefPtr) {
323  auto const dx = posRefPtr->getX() - x;
324  auto const dy = posRefPtr->getY() - y;
325  auto const distSq = dx*dx + dy*dy;
326  if (distSq < minDistSq) {
327  foundPtr = posRefPtr;
328  minDistSq = distSq;
329  }
330  }
331  return std::make_pair(foundPtr, std::sqrt(minDistSq));
332  }
333 
352  void addNearestMatch(
353  MultiIndexedProxyPairList &proxyPairList,
354  boost::shared_array<double> coeff,
355  ProxyVector const & posRefCat,
356  RecordProxy const &source,
357  double matchingAllowancePix
358  ) {
359  double x1, y1;
360  auto x0 = source.getX();
361  auto y0 = source.getY();
362  transform(1, coeff, x0, y0, &x1, &y1);
363  auto refObjDist = searchNearestPoint(posRefCat, x1, y1, matchingAllowancePix);
364  if (refObjDist.first == posRefCat.end()) {
365  // no reference object sufficiently close; do nothing
366  return;
367  }
368  auto existingMatch = proxyPairList.get<refIdTag>().find(refObjDist.first->record->getId());
369  if (existingMatch == proxyPairList.get<refIdTag>().end()) {
370  // reference object not used before; add new entry
371  auto proxyPair = ProxyPair(*refObjDist.first, source);
372  proxyPairList.get<refIdTag>().insert(proxyPair);
373  return;
374 #if 0
375  } else {
376  // This code keeps the closest match. It is disabled because it appears to be unnecessary
377  // and may be undesirable, since the image may have very faint sources which may give
378  // false matches. Note that the calling algorithm checks for matches starting with the
379  // brightest source and works down in source brightness.
380  if (existingMatch->distance <= refObjDist.second) {
381  // reference object used before and old source is closer; do nothing
382  return;
383  } else {
384  // reference object used before, but new source is closer; delete old entry and add new
385  proxyPairList.get<refIdTag>().erase(existingMatch);
386  auto proxyPair = ProxyPair(*refObjDist.first, source);
387  proxyPairList.get<refIdTag>().insert(proxyPair);
388  return;
389  }
390 #endif
391  }
392  }
393 
402  afwTable::ReferenceMatchVector FinalVerify(
403  boost::shared_array<double> coeff,
404  ProxyVector const & posRefCat,
405  ProxyVector const & sourceCat,
406  double matchingAllowancePix,
407  bool verbose
408  ) {
409  MultiIndexedProxyPairList proxyPairList;
410 
411  for (auto sourcePtr = sourceCat.begin(); sourcePtr != sourceCat.end(); ++sourcePtr) {
412  addNearestMatch(proxyPairList, coeff, posRefCat, *sourcePtr, matchingAllowancePix);
413  }
414  int order = 1;
415  if (proxyPairList.size() > 5) {
416  for (auto j = 0; j < 100; j++) {
417  auto prevNumMatches = proxyPairList.size();
418  auto srcMat = ProxyVector();
419  auto catMat = ProxyVector();
420  srcMat.reserve(proxyPairList.size());
421  catMat.reserve(proxyPairList.size());
422  for (auto matchPtr = proxyPairList.get<refIdTag>().begin();
423  matchPtr != proxyPairList.get<refIdTag>().end(); ++matchPtr) {
424  catMat.push_back(matchPtr->first);
425  srcMat.push_back(matchPtr->second);
426  }
427  coeff = polyfit(order, srcMat, catMat);
428  proxyPairList.clear();
429 
430  for (ProxyVector::const_iterator sourcePtr = sourceCat.begin();
431  sourcePtr != sourceCat.end(); ++sourcePtr) {
432  addNearestMatch(proxyPairList, coeff, posRefCat, *sourcePtr, matchingAllowancePix);
433  }
434  if (proxyPairList.size() == prevNumMatches) {
435  break;
436  }
437  }
438  }
439  auto matPair = afwTable::ReferenceMatchVector();
440  matPair.reserve(proxyPairList.size());
441  for (auto proxyPairIter = proxyPairList.get<insertionOrderTag>().begin();
442  proxyPairIter != proxyPairList.get<insertionOrderTag>().end(); ++proxyPairIter) {
443  matPair.push_back(afwTable::ReferenceMatch(
444  proxyPairIter->first.record,
445  boost::static_pointer_cast<afwTable::SourceRecord>(proxyPairIter->second.record),
446  proxyPairIter->distance
447  ));
448  }
449 
450  return matPair;
451  }
452 
453 } // anonymous namespace
454 
455 namespace lsst {
456 namespace meas {
457 namespace astrom {
458 
459  void MatchOptimisticBControl::validate() const {
460  if (refFluxField.empty()) {
461  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "refFluxField must be specified");
462  }
463  if (sourceFluxField.empty()) {
464  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "sourceFluxField must be specified");
465  }
466  if (numBrightStars <= 0) {
467  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "numBrightStars must be positive");
468  }
469  if (minMatchedPairs < 0) {
470  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "minMatchedPairs must not be negative");
471  }
472  if (matchingAllowancePix <= 0) {
473  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "matchingAllowancePix must be positive");
474  }
475  if (maxOffsetPix <= 0) {
476  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "maxOffsetPix must be positive");
477  }
478  if (maxRotationDeg <= 0) {
479  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "maxRotationRad must be positive");
480  }
481  if (allowedNonperpDeg <= 0) {
482  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "allowedNonperpDeg must be positive");
483  }
484  if (numPointsForShape <= 0) {
485  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "numPointsForShape must be positive");
486  }
487  if (maxDeterminant <= 0) {
488  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "maxDeterminant must be positive");
489  }
490  }
491 
493  afw::image::Wcs const& distortedWcs,
494  afw::image::Wcs const& tanWcs
495  )
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  }
506 
508  afw::image::Wcs const& tanWcs
509  )
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  }
520 
522  afwTable::SimpleCatalog const &posRefCat,
523  afwTable::SourceCatalog const &sourceCat,
524  MatchOptimisticBControl const &control,
525  afw::image::Wcs const& wcs,
526  int posRefBegInd,
527  bool verbose
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>(
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 
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))),
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  }
804 
805 }}} // namespace lsst::afw::astrom
int y
int iter
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
int numPointsForShape
&quot;number of points in a matching shape&quot; ;
Extent< double, 2 > Extent2D
Definition: Extent.h:358
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g. RA/dec) to pixel positions.
Definition: Wcs.cc:782
#define CONST_PTR(...)
Definition: base.h:47
tbl::Key< int > wcs
double matchingAllowancePix
&quot;maximum allowed distance between reference objects and sources (pixels)&quot; ;
std::vector< ReferenceMatch > ReferenceMatchVector
Definition: fwd.h:97
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
lsst::afw::coord::IcrsCoord IcrsCoord
Definition: misc.h:38
int isnan(T t)
Definition: ieee.h:110
int const x0
Definition: saturated.cc:45
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)
double asDegrees() const
Definition: Angle.h:124
std::vector< RecordProxy > ProxyVector
Extent< double, 3 > Extent3D
Definition: Extent.h:359
A coordinate class intended to represent absolute positions.
Definition: PSF.h:39
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
Schema getSchema() const
Return the schema associated with the catalog&#39;s table.
Definition: Catalog.h:114
Lightweight representation of a geometric match between two records.
Definition: fwd.h:90
Eigen::Matrix2d getCDMatrix() const
Returns the CD matrix.
Definition: Wcs.cc:573
virtual bool hasDistortion() const
Definition: Wcs.h:220
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
double maxOffsetPix
&quot;maximum allowed frame translation (pixels)&quot; ;
int x
double degToRad(long double angleInDegrees)
Definition: RaDecStr.cc:67
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int numBrightStars
&quot;maximum number of bright reference stars to use&quot; ;
std::string refFluxField
&quot;name of flux field in reference catalog&quot; ;
afw::table::Key< double > b
double maxRotationDeg
&quot;maximum allowed frame rotation (deg)&quot; ;
Point< double, 2 > Point2D
Definition: Point.h:286
ProxyVector makeProxies(lsst::afw::table::SourceCatalog const &sourceCat, afw::image::Wcs const &distortedWcs, afw::image::Wcs const &tanWcs)
std::string sourceFluxField
&quot;name of flux field in source catalog&quot; ;
Base::const_iterator const_iterator
Definition: SortedCatalog.h:49
int const y0
Definition: saturated.cc:45
boost::int64_t RecordId
Type used for unique IDs for records.
Definition: misc.h:21
double allowedNonperpDeg
&quot;allowed non-perpendicularity of x and y axes (deg)&quot; ;
int minMatchedPairs
&quot;minimum number of matches&quot; ;
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
boost::shared_ptr< lsst::afw::table::SimpleRecord > record
Include files required for standard LSST Exception handling.
A FunctorKey used to get or set celestial coordiantes from a pair of Angle keys.
Definition: aggregates.h:119
size_type size() const
Return the number of elements in the catalog.
Definition: Catalog.h:401