LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
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
18#include "lsst/pex/exceptions.h"
19#include "lsst/geom/Angle.h"
20#include "lsst/geom/Extent.h"
21#include "lsst/geom/Point.h"
25
27namespace afwTable = lsst::afw::table;
28
29namespace {
30using namespace lsst::meas::astrom;
31
32afwTable::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
35struct refIdTag {}; // tag for retrieval by reference object ID
36struct insertionOrderTag {}; // tag for retrieval by insertion order
37typedef boost::multi_index_container<
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
55inline double absDeltaAngle(double ang1, double ang2) { return std::fmod(std::fabs(ang1 - ang2), M_PI * 2); }
56
57bool cmpPair(ProxyPair const &a, ProxyPair const &b) { return a.distance > b.distance; }
58
64struct 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
90ProxyVector 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};
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
102std::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
141void 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
157boost::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
272std::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
306void 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
351afwTable::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
399namespace lsst {
400namespace meas {
401namespace 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
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
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))),
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 }
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
717END:
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:128
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
Definition: Angle.h:176
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 empty(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:52
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:109
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