LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
matchOptimisticB.cc
Go to the documentation of this file.
1 #include <algorithm>
2 #include <cmath>
3 #include <fstream>
4 #include <iostream>
5 #include <memory>
6 #include <utility>
7 #include <vector>
8 
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/sphgeom/Vector3d.h"
18 #include "lsst/pex/exceptions.h"
19 #include "lsst/geom/Angle.h"
20 #include "lsst/geom/Extent.h"
21 #include "lsst/geom/Point.h"
22 #include "lsst/geom/SpherePoint.h"
23 #include "lsst/afw/geom/SkyWcs.h"
25 
27 namespace afwTable = lsst::afw::table;
28 
29 namespace {
30 using namespace lsst::meas::astrom;
31 
32 afwTable::RecordId getRefId(ProxyPair const &proxyPair) { return proxyPair.first.record->getId(); }
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<boost::multi_index::tag<insertionOrderTag> >,
41  boost::multi_index::ordered_unique< // index by reference object ID
42  boost::multi_index::tag<refIdTag>,
43  boost::multi_index::global_fun<ProxyPair const &, afwTable::RecordId, &getRefId> > > >
44  MultiIndexedProxyPairList;
45 
46 // Algorithm is based on V.Tabur 2007, PASA, 24, 189-198
47 // "Fast Algorithms for Matching CCD Images to a Stellar Catalogue"
48 
55 inline double absDeltaAngle(double ang1, double ang2) { return std::fmod(std::fabs(ang1 - ang2), M_PI * 2); }
56 
57 bool cmpPair(ProxyPair const &a, ProxyPair const &b) { return a.distance > b.distance; }
58 
64 struct CompareProxyFlux {
65  bool operator()(RecordProxy const &a, RecordProxy const &b) const {
66  double aFlux = a.record->get(key);
67  double bFlux = b.record->get(key);
68  if (std::isnan(aFlux)) {
69  aFlux = 0.0;
70  }
71  if (std::isnan(bFlux)) {
72  bFlux = 0.0;
73  }
74  return aFlux > bFlux;
75  }
76 
78 };
79 
90 ProxyVector selectPoint(ProxyVector const &a, afwTable::Key<double> const &key, std::size_t num,
91  std::size_t startInd = 0) {
92  if (startInd >= a.size()) {
93  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "startInd too big");
94  }
95  CompareProxyFlux cmp = {key};
96  ProxyVector b(a);
97  std::sort(b.begin(), b.end(), cmp);
98  std::size_t const endInd = std::min(startInd + num, b.size());
99  return ProxyVector(b.begin() + startInd, b.begin() + endInd);
100 }
101 
102 std::vector<ProxyPair> searchPair(std::vector<ProxyPair> const &a, ProxyPair const &p, double e,
103  double e_dpa) {
105 
106  for (size_t i = 0; i < a.size(); i++) {
107  double dd = std::fabs(a[i].distance - p.distance);
108  double dpa = absDeltaAngle(a[i].pa, p.pa);
109  if (dd < e && dpa < e_dpa) {
110  v.push_back(a[i]);
111  }
112  }
113 
114  return v;
115 }
116 
118  ProxyPair const &q, double e, double dpa, double e_dpa = 0.02) {
120  double dd_min = 1.E+10;
121  // double dpa_min = e_dpa;
122 
123  for (auto i = a.begin(); i != a.end(); ++i) {
124  double dd = std::fabs(i->distance - p.distance);
125 #if 1
126  if (dd < e && absDeltaAngle(p.pa, i->pa - dpa) < e_dpa && dd < dd_min && (i->first == q.first)) {
127  dd_min = dd;
128  idx = i;
129  }
130 #else
131  if (dd < e && absDeltaAngle(p.pa, i->pa - dpa) < dpa_min) {
132  dpa_min = std::fabs(p.pa - i->pa - dpa);
133  idx = i;
134  }
135 #endif
136  }
137 
138  return idx;
139 }
140 
141 void transform(int order, boost::shared_array<double> const &coeff, double x, double y, double *xn,
142  double *yn) {
143  int ncoeff = (order + 1) * (order + 2) / 2;
144  *xn = 0.0;
145  *yn = 0.0;
146  int n = 0;
147  for (int i = 0; i <= order; i++) {
148  for (int k = 0; k <= i; k++) {
149  int j = i - k;
150  *xn += coeff[n] * pow(x, j) * pow(y, k);
151  *yn += coeff[n + ncoeff] * pow(x, j) * pow(y, k);
152  n++;
153  }
154  }
155 }
156 
157 boost::shared_array<double> polyfit(int order, ProxyVector const &img, ProxyVector const &posRefCat) {
158  int ncoeff = (order + 1) * (order + 2) / 2;
159  std::unique_ptr<int[]> xorder(new int[ncoeff]);
160  std::unique_ptr<int[]> yorder(new int[ncoeff]);
161 
162  int n = 0;
163  for (int i = 0; i <= order; i++) {
164  for (int k = 0; k <= i; k++) {
165  int j = i - k;
166  xorder[n] = j;
167  yorder[n] = k;
168  n++;
169  }
170  }
171 
172  std::unique_ptr<int[]> flag(new int[img.size()]);
173  for (size_t k = 0; k < img.size(); k++) {
174  flag[k] = 1;
175  }
176 
177  std::unique_ptr<double[]> a_data(new double[ncoeff * ncoeff]);
178  std::unique_ptr<double[]> b_data(new double[ncoeff]);
179  std::unique_ptr<double[]> c_data(new double[ncoeff]);
180 
181  boost::shared_array<double> coeff(new double[ncoeff * 2]);
182 
183  for (int loop = 0; loop < 1; loop++) {
184  for (int i = 0; i < ncoeff; i++) {
185  for (int j = 0; j < ncoeff; j++) {
186  a_data[i * ncoeff + j] = 0.0;
187  for (size_t k = 0; k < img.size(); k++) {
188  if (flag[k] == 1) {
189  a_data[i * ncoeff + j] +=
190  pow(img[k].getX(), xorder[i]) * pow(img[k].getY(), yorder[i]) *
191  pow(img[k].getX(), xorder[j]) * pow(img[k].getY(), yorder[j]);
192  }
193  }
194  }
195  b_data[i] = c_data[i] = 0.0;
196  for (unsigned int k = 0; k < img.size(); k++) {
197  if (flag[k] == 1) {
198  b_data[i] += pow(img[k].getX(), xorder[i]) * pow(img[k].getY(), yorder[i]) *
199  posRefCat[k].getX();
200  c_data[i] += pow(img[k].getX(), xorder[i]) * pow(img[k].getY(), yorder[i]) *
201  posRefCat[k].getY();
202  }
203  }
204  }
205 
206  gsl_matrix_view a = gsl_matrix_view_array(a_data.get(), ncoeff, ncoeff);
207  gsl_vector_view b = gsl_vector_view_array(b_data.get(), ncoeff);
208  gsl_vector_view c = gsl_vector_view_array(c_data.get(), ncoeff);
209 
210  std::shared_ptr<gsl_vector> x(gsl_vector_alloc(ncoeff), gsl_vector_free);
211  std::shared_ptr<gsl_vector> y(gsl_vector_alloc(ncoeff), gsl_vector_free);
212 
213  int s;
214 
215  std::shared_ptr<gsl_permutation> p(gsl_permutation_alloc(ncoeff), gsl_permutation_free);
216 
217  gsl_linalg_LU_decomp(&a.matrix, p.get(), &s);
218  gsl_linalg_LU_solve(&a.matrix, p.get(), &b.vector, x.get());
219  gsl_linalg_LU_solve(&a.matrix, p.get(), &c.vector, y.get());
220 
221  for (int i = 0; i < ncoeff; i++) {
222  coeff[i] = x->data[i];
223  coeff[i + ncoeff] = y->data[i];
224  }
225 
226  double S, Sx, Sy, Sxx, Syy;
227  S = Sx = Sy = Sxx = Syy = 0.0;
228  for (size_t k = 0; k < img.size(); k++) {
229  if (flag[k] == 1) {
230  double x0 = img[k].getX();
231  double y0 = img[k].getY();
232  double x1, y1;
233  transform(order, coeff, x0, y0, &x1, &y1);
234  S += 1.;
235  Sx += (x1 - posRefCat[k].getX());
236  Sxx += (x1 - posRefCat[k].getX()) * (x1 - posRefCat[k].getX());
237  Sy += (y1 - posRefCat[k].getY());
238  Syy += (y1 - posRefCat[k].getY()) * (y1 - posRefCat[k].getY());
239  }
240  }
241  double x_sig = std::sqrt((Sxx - Sx * Sx / S) / S);
242  double y_sig = std::sqrt((Syy - Sy * Sy / S) / S);
243  // std::cout << x_sig << " " << y_sig << std::endl;
244 
245  for (size_t k = 0; k < img.size(); k++) {
246  double x0 = img[k].getX();
247  double y0 = img[k].getY();
248  double x1, y1;
249  transform(order, coeff, x0, y0, &x1, &y1);
250  if (std::fabs(x1 - posRefCat[k].getX()) > 2. * x_sig ||
251  std::fabs(y1 - posRefCat[k].getY()) > 2. * y_sig) {
252  flag[k] = 0;
253  }
254  }
255  }
256 
257  return coeff;
258 }
259 
272 std::pair<ProxyVector::const_iterator, double> searchNearestPoint(ProxyVector const &posRefCat, double x,
273  double y, double matchingAllowancePix) {
274  auto minDistSq = matchingAllowancePix * matchingAllowancePix;
275  auto foundPtr = posRefCat.end();
276  for (auto posRefPtr = posRefCat.begin(); posRefPtr != posRefCat.end(); ++posRefPtr) {
277  auto const dx = posRefPtr->getX() - x;
278  auto const dy = posRefPtr->getY() - y;
279  auto const distSq = dx * dx + dy * dy;
280  if (distSq < minDistSq) {
281  foundPtr = posRefPtr;
282  minDistSq = distSq;
283  }
284  }
285  return std::make_pair(foundPtr, std::sqrt(minDistSq));
286 }
287 
306 void addNearestMatch(MultiIndexedProxyPairList &proxyPairList, boost::shared_array<double> coeff,
307  ProxyVector const &posRefCat, RecordProxy const &source, double matchingAllowancePix) {
308  double x1, y1;
309  auto x0 = source.getX();
310  auto y0 = source.getY();
311  transform(1, coeff, x0, y0, &x1, &y1);
312  auto refObjDist = searchNearestPoint(posRefCat, x1, y1, matchingAllowancePix);
313  if (refObjDist.first == posRefCat.end()) {
314  // no reference object sufficiently close; do nothing
315  return;
316  }
317  auto existingMatch = proxyPairList.get<refIdTag>().find(refObjDist.first->record->getId());
318  if (existingMatch == proxyPairList.get<refIdTag>().end()) {
319  // reference object not used before; add new entry
320  auto proxyPair = ProxyPair(*refObjDist.first, source);
321  proxyPairList.get<refIdTag>().insert(proxyPair);
322  return;
323 #if 0
324  } else {
325  // This code keeps the closest match. It is disabled because it appears to be unnecessary
326  // and may be undesirable, since the image may have very faint sources which may give
327  // false matches. Note that the calling algorithm checks for matches starting with the
328  // brightest source and works down in source brightness.
329  if (existingMatch->distance <= refObjDist.second) {
330  // reference object used before and old source is closer; do nothing
331  return;
332  } else {
333  // reference object used before, but new source is closer; delete old entry and add new
334  proxyPairList.get<refIdTag>().erase(existingMatch);
335  auto proxyPair = ProxyPair(*refObjDist.first, source);
336  proxyPairList.get<refIdTag>().insert(proxyPair);
337  return;
338  }
339 #endif
340  }
341 }
342 
351 afwTable::ReferenceMatchVector FinalVerify(boost::shared_array<double> coeff, ProxyVector const &posRefCat,
352  ProxyVector const &sourceCat, double matchingAllowancePix,
353  bool verbose) {
354  MultiIndexedProxyPairList proxyPairList;
355 
356  for (auto sourcePtr = sourceCat.begin(); sourcePtr != sourceCat.end(); ++sourcePtr) {
357  addNearestMatch(proxyPairList, coeff, posRefCat, *sourcePtr, matchingAllowancePix);
358  }
359  int order = 1;
360  if (proxyPairList.size() > 5) {
361  for (auto j = 0; j < 100; j++) {
362  auto prevNumMatches = proxyPairList.size();
363  auto srcMat = ProxyVector();
364  auto catMat = ProxyVector();
365  srcMat.reserve(proxyPairList.size());
366  catMat.reserve(proxyPairList.size());
367  for (auto matchPtr = proxyPairList.get<refIdTag>().begin();
368  matchPtr != proxyPairList.get<refIdTag>().end(); ++matchPtr) {
369  catMat.push_back(matchPtr->first);
370  srcMat.push_back(matchPtr->second);
371  }
372  coeff = polyfit(order, srcMat, catMat);
373  proxyPairList.clear();
374 
375  for (ProxyVector::const_iterator sourcePtr = sourceCat.begin(); sourcePtr != sourceCat.end();
376  ++sourcePtr) {
377  addNearestMatch(proxyPairList, coeff, posRefCat, *sourcePtr, matchingAllowancePix);
378  }
379  if (proxyPairList.size() == prevNumMatches) {
380  break;
381  }
382  }
383  }
384  auto matPair = afwTable::ReferenceMatchVector();
385  matPair.reserve(proxyPairList.size());
386  for (auto proxyPairIter = proxyPairList.get<insertionOrderTag>().begin();
387  proxyPairIter != proxyPairList.get<insertionOrderTag>().end(); ++proxyPairIter) {
388  matPair.push_back(afwTable::ReferenceMatch(
389  proxyPairIter->first.record,
390  std::static_pointer_cast<afwTable::SourceRecord>(proxyPairIter->second.record),
391  proxyPairIter->distance));
392  }
393 
394  return matPair;
395 }
396 
397 } // anonymous namespace
398 
399 namespace lsst {
400 namespace meas {
401 namespace astrom {
402 
404  if (refFluxField.empty()) {
405  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "refFluxField must be specified");
406  }
407  if (sourceFluxField.empty()) {
408  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "sourceFluxField must be specified");
409  }
410  if (numBrightStars <= 0) {
411  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "numBrightStars must be positive");
412  }
413  if (minMatchedPairs < 0) {
414  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "minMatchedPairs must not be negative");
415  }
416  if (matchingAllowancePix <= 0) {
417  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "matchingAllowancePix must be positive");
418  }
419  if (maxOffsetPix <= 0) {
420  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "maxOffsetPix must be positive");
421  }
422  if (maxRotationDeg <= 0) {
423  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "maxRotationRad must be positive");
424  }
425  if (allowedNonperpDeg <= 0) {
426  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "allowedNonperpDeg must be positive");
427  }
428  if (numPointsForShape <= 0) {
429  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "numPointsForShape must be positive");
430  }
431  if (maxDeterminant <= 0) {
432  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "maxDeterminant must be positive");
433  }
434 }
435 
436 ProxyVector makeProxies(afwTable::SourceCatalog const &sourceCat, afw::geom::SkyWcs const &distortedWcs,
437  afw::geom::SkyWcs const &tanWcs) {
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 }
447 
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 }
458 
460  afwTable::SourceCatalog const &sourceCat,
461  MatchOptimisticBControl const &control,
462  afw::geom::SkyWcs const &wcs, int posRefBegInd,
463  bool verbose) {
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 }
732 
733 } // namespace astrom
734 } // namespace meas
735 } // namespace lsst
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
#define M_PI
Definition: ListMatch.cc:31
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
int y
Definition: SpanSet.cc:48
table::Key< int > transform
table::Key< int > b
table::Key< int > a
This file declares a class for representing vectors in ℝ³.
T acos(T... args)
T begin(T... args)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Definition: SkyWcs.h:334
lsst::geom::Point2D skyToPixel(lsst::geom::SpherePoint const &sky) const
Compute pixel position(s) from sky position(s)
Definition: SkyWcs.h:349
size_type size() const
Return the number of elements in the catalog.
Definition: Catalog.h:413
bool empty() const
Return true if the catalog has no records.
Definition: Catalog.h:410
iterator begin()
Iterator access.
Definition: Catalog.h:401
Schema getSchema() const
Return the schema associated with the catalog's table.
Definition: Catalog.h:118
A FunctorKey used to get or set celestial coordinates from a pair of lsst::geom::Angle keys.
Definition: aggregates.h:210
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition: Schema.cc:467
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Definition: SortedCatalog.h:42
typename Base::const_iterator const_iterator
Definition: SortedCatalog.h:50
A class representing an angle.
Definition: Angle.h:127
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
Definition: Angle.h:169
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
Reports invalid arguments.
Definition: Runtime.h:66
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
T clear(T... args)
T end(T... args)
T endl(T... args)
T fabs(T... args)
T find(T... args)
T fmod(T... args)
T isnan(T... args)
T make_pair(T... args)
T min(T... args)
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224
def erase(frame=None)
Definition: ds9.py:96
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
std::vector< ReferenceMatch > ReferenceMatchVector
Definition: fwd.h:108
Extent< double, 2 > Extent2D
Definition: Extent.h:400
constexpr double degToRad(double x) noexcept
Definition: Angle.h:51
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
std::vector< RecordProxy > ProxyVector
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.
ProxyVector makeProxies(afw::table::SourceCatalog const &sourceCat, afw::geom::SkyWcs const &distortedWcs, afw::geom::SkyWcs const &tanWcs)
A base class for image defects.
T pow(T... args)
T push_back(T... args)
T reserve(T... args)
T size(T... args)
T sort(T... args)
T sqrt(T... args)
table::Key< int > ncoeff
Definition: PsfexPsf.cc:327
table::Key< table::Array< double > > coeff
Definition: PsfexPsf.cc:362
Lightweight representation of a geometric match between two records.
Definition: Match.h:67
std::string refFluxField
"name of flux field in reference catalog" ;
double maxRotationDeg
"maximum allowed frame rotation (deg)" ;
double matchingAllowancePix
"maximum allowed distance between reference objects and sources (pixels)" ;
int numPointsForShape
"number of points in a matching shape" ;
int minMatchedPairs
"minimum number of matches" ;
int numBrightStars
"maximum number of bright reference stars to use" ;
std::string sourceFluxField
"name of flux field in source catalog" ;
double maxOffsetPix
"maximum allowed frame translation (pixels)" ;
double allowedNonperpDeg
"allowed non-perpendicularity of x and y axes (deg)" ;
A wrapper around a SimpleRecord or SourceRecord that allows us to record a pixel position in a way th...
std::shared_ptr< afw::table::SimpleRecord > record
table::Key< int > order