LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
Match.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010, 2016 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include <algorithm>
26 #include <cmath>
27 #include <memory>
28 
29 #include "lsst/pex/exceptions.h"
30 #include "lsst/log/Log.h"
31 #include "lsst/geom/Angle.h"
32 #include "lsst/afw/table/Match.h"
33 
34 namespace lsst {
35 namespace afw {
36 namespace table {
37 namespace {
38 
39 template <typename RecordT>
40 struct RecordPos {
41  double dec;
42  double x;
43  double y;
44  double z;
45  // JFB removed extra pointer here; this may have performance implications, but hopefully not
46  // significant ones. BaseCatalog iterators yield temporary BaseRecord shared_ptrs, so storing
47  // their address was no longer an option.
49 };
50 
51 template <typename Record1, typename Record2>
52 bool operator<(RecordPos<Record1> const &s1, RecordPos<Record2> const &s2) {
53  return (s1.dec < s2.dec);
54 }
55 
56 struct CmpRecordPtr {
57  bool operator()(std::shared_ptr<SourceRecord> const s1, std::shared_ptr<SourceRecord> const s2) {
58  return s1->getY() < s2->getY();
59  }
60 };
61 
72 template <typename Cat>
73 size_t makeRecordPositions(Cat const &cat, RecordPos<typename Cat::Record> *positions) {
74  size_t n = 0;
75  Key<lsst::geom::Angle> raKey = Cat::Table::getCoordKey().getRa();
76  Key<lsst::geom::Angle> decKey = Cat::Table::getCoordKey().getDec();
77  for (typename Cat::const_iterator i(cat.begin()), e(cat.end()); i != e; ++i) {
78  lsst::geom::Angle ra = i->get(raKey);
79  lsst::geom::Angle dec = i->get(decKey);
80  if (std::isnan(ra.asRadians()) || std::isnan(dec.asRadians())) {
81  continue;
82  }
83  double cosDec = std::cos(dec);
84  positions[n].dec = dec.asRadians();
85  positions[n].x = std::cos(ra) * cosDec;
86  positions[n].y = std::sin(ra) * cosDec;
87  positions[n].z = std::sin(dec);
88  positions[n].src = i;
89  ++n;
90  }
91  std::sort(positions, positions + n);
92  if (n < cat.size()) {
93  LOGLS_WARN("afw.table.matchRaDec", "At least one source had ra or dec equal to NaN");
94  }
95  return n;
96 }
97 
98 template size_t makeRecordPositions(SimpleCatalog const &, RecordPos<SimpleRecord> *);
99 template size_t makeRecordPositions(SourceCatalog const &, RecordPos<SourceRecord> *);
100 
101 template <typename Cat1, typename Cat2>
102 bool doSelfMatchIfSame(std::vector<Match<typename Cat1::Record, typename Cat2::Record> > &result,
103  Cat1 const &cat1, Cat2 const &cat2, lsst::geom::Angle radius) {
104  // types are different, so the catalogs are never the same.
105  return false;
106 }
107 
108 template <typename Cat>
109 bool doSelfMatchIfSame(std::vector<Match<typename Cat::Record, typename Cat::Record> > &result,
110  Cat const &cat1, Cat const &cat2, lsst::geom::Angle radius) {
111  if (&cat1 == &cat2) {
112  result = matchRaDec(cat1, radius);
113  return true;
114  }
115  return false;
116 }
117 
126 double toUnitSphereDistanceSquared(lsst::geom::Angle theta) noexcept {
127  return 2. * (1. - std::cos(theta.asRadians()));
128  // == 4.0 * pow(std::sin(0.5 * theta.asRadians()), 2.0)
129 }
130 
139 lsst::geom::Angle fromUnitSphereDistanceSquared(double d2) noexcept {
140  // acos(1 - 0.5*d2) doesn't require sqrt but isn't as precise for small d2
141  return 2.0 * std::asin(0.5 * std::sqrt(d2)) * lsst::geom::radians;
142 }
143 
144 } // namespace
145 
146 template <typename Cat1, typename Cat2>
148  Cat2 const &cat2,
150  bool closest) {
151  MatchControl mc;
152  mc.findOnlyClosest = closest;
153 
154  return matchRaDec(cat1, cat2, radius, mc);
155 }
156 
157 template <typename Cat1, typename Cat2>
159  Cat2 const &cat2,
161  MatchControl const &mc) {
163  std::vector<MatchT> matches;
164 
165  if (doSelfMatchIfSame(matches, cat1, cat2, radius)) return matches;
166 
167  if (radius < 0.0 || (radius > (45. * lsst::geom::degrees))) {
168  throw LSST_EXCEPT(pex::exceptions::RangeError, "match radius out of range (0 to 45 degrees)");
169  }
170  if (cat1.size() == 0 || cat2.size() == 0) {
171  return matches;
172  }
173  // setup match parameters
174  double const d2Limit = toUnitSphereDistanceSquared(radius);
175 
176  // Build position lists
177  size_t len1 = cat1.size();
178  size_t len2 = cat2.size();
179 
180  typedef RecordPos<typename Cat1::Record> Pos1;
181  typedef RecordPos<typename Cat2::Record> Pos2;
182  std::unique_ptr<Pos1[]> pos1(new Pos1[len1]);
183  std::unique_ptr<Pos2[]> pos2(new Pos2[len2]);
184  len1 = makeRecordPositions(cat1, pos1.get());
185  len2 = makeRecordPositions(cat2, pos2.get());
187 
188  for (size_t i = 0, start = 0; i < len1; ++i) {
189  double minDec = pos1[i].dec - radius.asRadians();
190  while (start < len2 && pos2[start].dec < minDec) {
191  ++start;
192  }
193  if (start == len2) {
194  break;
195  }
196  double maxDec = pos1[i].dec + radius.asRadians();
197  size_t closestIndex = -1; // Index of closest match (if any)
198  double d2Include = d2Limit; // Squared distance for inclusion of match
199  bool found = false; // Found anything?
200  size_t nMatches = 0; // Number of matches
201  for (size_t j = start; j < len2 && pos2[j].dec <= maxDec; ++j) {
202  double dx = pos1[i].x - pos2[j].x;
203  double dy = pos1[i].y - pos2[j].y;
204  double dz = pos1[i].z - pos2[j].z;
205  double d2 = dx * dx + dy * dy + dz * dz;
206  if (d2 < d2Include) {
207  if (mc.findOnlyClosest) {
208  d2Include = d2;
209  closestIndex = j;
210  found = true;
211  } else {
212  matches.push_back(MatchT(pos1[i].src, pos2[j].src, fromUnitSphereDistanceSquared(d2)));
213  }
214  ++nMatches;
215  }
216  }
217  if (mc.includeMismatches && nMatches == 0) {
218  matches.push_back(MatchT(pos1[i].src, nullRecord, NAN));
219  }
220  if (mc.findOnlyClosest && found) {
221  matches.push_back(
222  MatchT(pos1[i].src, pos2[closestIndex].src, fromUnitSphereDistanceSquared(d2Include)));
223  }
224  }
225  return matches;
226 }
227 
228 #define LSST_MATCH_RADEC(RTYPE, C1, C2) \
229  template RTYPE matchRaDec(C1 const &, C2 const &, lsst::geom::Angle, bool); \
230  template RTYPE matchRaDec(C1 const &, C2 const &, lsst::geom::Angle, MatchControl const &)
231 
235 
236 #undef LSST_MATCH_RADEC
237 
238 template <typename Cat>
241  bool symmetric) {
242  MatchControl mc;
243  mc.symmetricMatch = symmetric;
244 
245  return matchRaDec(cat, radius, mc);
246 }
247 
248 template <typename Cat>
251  MatchControl const &mc) {
253  std::vector<MatchT> matches;
254 
255  if (radius < 0.0 || radius > (45.0 * lsst::geom::degrees)) {
256  throw LSST_EXCEPT(pex::exceptions::RangeError, "match radius out of range (0 to 45 degrees)");
257  }
258  if (cat.size() == 0) {
259  return matches;
260  }
261  // setup match parameters
262  double const d2Limit = toUnitSphereDistanceSquared(radius);
263 
264  // Build position list
265  size_t len = cat.size();
266  typedef RecordPos<typename Cat::Record> Pos;
267  std::unique_ptr<Pos[]> pos(new Pos[len]);
268  len = makeRecordPositions(cat, pos.get());
269 
270  for (size_t i = 0; i < len; ++i) {
271  double maxDec = pos[i].dec + radius.asRadians();
272  for (size_t j = i + 1; j < len && pos[j].dec <= maxDec; ++j) {
273  double dx = pos[i].x - pos[j].x;
274  double dy = pos[i].y - pos[j].y;
275  double dz = pos[i].z - pos[j].z;
276  double d2 = dx * dx + dy * dy + dz * dz;
277  if (d2 < d2Limit) {
278  lsst::geom::Angle d = fromUnitSphereDistanceSquared(d2);
279  matches.push_back(MatchT(pos[i].src, pos[j].src, d));
280  if (mc.symmetricMatch) {
281  matches.push_back(MatchT(pos[j].src, pos[i].src, d));
282  }
283  }
284  }
285  }
286  return matches;
287 }
288 
289 #define LSST_MATCH_RADEC(RTYPE, C) \
290  template RTYPE matchRaDec(C const &, lsst::geom::Angle, bool); \
291  template RTYPE matchRaDec(C const &, lsst::geom::Angle, MatchControl const &)
292 
295 
296 #undef LSST_MATCH_RADEC
297 
298 SourceMatchVector matchXy(SourceCatalog const &cat1, SourceCatalog const &cat2, double radius, bool closest) {
299  MatchControl mc;
300  mc.findOnlyClosest = closest;
301 
302  return matchXy(cat1, cat2, radius, mc);
303 }
304 
305 SourceMatchVector matchXy(SourceCatalog const &cat1, SourceCatalog const &cat2, double radius,
306  MatchControl const &mc) {
307  if (&cat1 == &cat2) {
308  return matchXy(cat1, radius);
309  }
310  // setup match parameters
311  double const r2 = radius * radius;
312 
313  // copy and sort array of pointers on y
314  size_t len1 = cat1.size();
315  size_t len2 = cat2.size();
319  size_t n = 0;
320  for (SourceCatalog::const_iterator i(cat1.begin()), e(cat1.end()); i != e; ++i) {
321  if (std::isnan(i->getX()) || std::isnan(i->getY())) {
322  continue;
323  }
324  pos1[n] = i;
325  ++n;
326  }
327  len1 = n;
328  n = 0;
329  for (SourceCatalog::const_iterator i(cat2.begin()), e(cat2.end()); i != e; ++i) {
330  if (std::isnan(i->getX()) || std::isnan(i->getY())) {
331  continue;
332  }
333  pos2[n] = i;
334  ++n;
335  }
336  len2 = n;
337 
338  std::sort(pos1.get(), pos1.get() + len1, CmpRecordPtr());
339  std::sort(pos2.get(), pos2.get() + len2, CmpRecordPtr());
340 
341  SourceMatchVector matches;
342  for (size_t i = 0, start = 0; i < len1; ++i) {
343  double y = pos1[i]->getY();
344  double minY = y - radius;
345  while (start < len2 && pos2[start]->getY() < minY) {
346  ++start;
347  }
348  if (start == len2) {
349  break;
350  }
351  double x = pos1[i]->getX();
352  double maxY = y + radius;
353  double y2;
354  size_t closestIndex = -1; // Index of closest match (if any)
355  double r2Include = r2; // Squared radius for inclusion of match
356  bool found = false; // Found anything?
357  size_t nMatches = 0; // Number of matches
358  for (size_t j = start; j < len2 && (y2 = pos2[j]->getY()) <= maxY; ++j) {
359  double dx = x - pos2[j]->getX();
360  double dy = y - y2;
361  double d2 = dx * dx + dy * dy;
362  if (d2 < r2Include) {
363  if (mc.findOnlyClosest) {
364  r2Include = d2;
365  closestIndex = j;
366  found = true;
367  } else {
368  matches.push_back(SourceMatch(pos1[i], pos2[j], std::sqrt(d2)));
369  }
370  ++nMatches;
371  }
372  }
373  if (mc.includeMismatches && nMatches == 0) {
374  matches.push_back(SourceMatch(pos1[i], nullRecord, NAN));
375  }
376  if (mc.findOnlyClosest && found) {
377  matches.push_back(SourceMatch(pos1[i], pos2[closestIndex], std::sqrt(r2Include)));
378  }
379  }
380  return matches;
381 }
382 
383 SourceMatchVector matchXy(SourceCatalog const &cat, double radius, bool symmetric) {
384  MatchControl mc;
385  mc.symmetricMatch = symmetric;
386 
387  return matchXy(cat, radius, mc);
388 }
389 
391  // setup match parameters
392  double const r2 = radius * radius;
393 
394  // copy and sort array of pointers on y
395  size_t len = cat.size();
397  size_t n = 0;
398  for (SourceCatalog::const_iterator i(cat.begin()), e(cat.end()); i != e; ++i) {
399  if (std::isnan(i->getX()) || std::isnan(i->getY())) {
400  continue;
401  }
402  pos[n] = i;
403  ++n;
404  }
405  len = n;
406 
407  std::sort(pos.get(), pos.get() + len, CmpRecordPtr());
408 
409  SourceMatchVector matches;
410  for (size_t i = 0; i < len; ++i) {
411  double x = pos[i]->getX();
412  double y = pos[i]->getY();
413  double maxY = y + radius;
414  double y2;
415  for (size_t j = i + 1; j < len && (y2 = pos[j]->getY()) <= maxY; ++j) {
416  double dx = x - pos[j]->getX();
417  double dy = y - y2;
418  double d2 = dx * dx + dy * dy;
419  if (d2 < r2) {
420  double d = std::sqrt(d2);
421  matches.push_back(SourceMatch(pos[i], pos[j], d));
422  if (mc.symmetricMatch) {
423  matches.push_back(SourceMatch(pos[j], pos[i], d));
424  }
425  }
426  }
427  }
428  return matches;
429 }
430 
431 template <typename Record1, typename Record2>
433  Schema schema;
434  Key<RecordId> outKey1 = schema.addField<RecordId>("first", "ID for first source record in match.");
435  Key<RecordId> outKey2 = schema.addField<RecordId>("second", "ID for second source record in match.");
436  Key<double> keyD = schema.addField<double>("distance", "Distance between matches sources.");
437  BaseCatalog result(schema);
438  result.getTable()->preallocate(matches.size());
439  result.reserve(matches.size());
440  typedef typename std::vector<Match<Record1, Record2> >::const_iterator Iter;
441  for (Iter i = matches.begin(); i != matches.end(); ++i) {
442  std::shared_ptr<BaseRecord> record = result.addNew();
443  record->set(outKey1, i->first->getId());
444  record->set(outKey2, i->second->getId());
445  record->set(keyD, i->distance);
446  }
447  return result;
448 }
449 
450 template BaseCatalog packMatches(SimpleMatchVector const &);
452 template BaseCatalog packMatches(SourceMatchVector const &);
453 
454 template <typename Cat1, typename Cat2>
456  Cat1 const &first,
457  Cat2 const &second) {
458  LOG_LOGGER tableLog = LOG_GET("afw.table");
459  Key<RecordId> inKey1 = matches.getSchema()["first"];
460  Key<RecordId> inKey2 = matches.getSchema()["second"];
461  Key<double> keyD = matches.getSchema()["distance"];
462  if (!first.isSorted() || !second.isSorted())
464  "Catalogs passed to unpackMatches must be sorted.");
467  result.resize(matches.size());
468  typename std::vector<MatchT>::iterator j = result.begin();
469  for (BaseCatalog::const_iterator i = matches.begin(); i != matches.end(); ++i, ++j) {
470  typename Cat1::const_iterator k1 = first.find(i->get(inKey1));
471  typename Cat2::const_iterator k2 = second.find(i->get(inKey2));
472  if (k1 != first.end()) {
473  j->first = k1;
474  } else {
475  LOGLS_WARN(tableLog,
476  "Persisted match record with ID " << i->get(inKey1) << " not found in catalog 1.");
477  }
478  if (k2 != second.end()) {
479  j->second = k2;
480  } else {
481  LOGLS_WARN(tableLog,
482  "Persisted match record with ID " << i->get(inKey2) << " not found in catalog 2.");
483  }
484  j->distance = i->get(keyD);
485  }
486  return result;
487 }
488 
489 template SimpleMatchVector unpackMatches(BaseCatalog const &, SimpleCatalog const &, SimpleCatalog const &);
491  SourceCatalog const &);
492 template SourceMatchVector unpackMatches(BaseCatalog const &, SourceCatalog const &, SourceCatalog const &);
493 } // namespace table
494 } // namespace afw
495 } // namespace lsst
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:648
Defines the fields and offsets for a table.
Definition: Schema.h:50
bool findOnlyClosest
"Return only the closest match if more than one is found " "(default: true)" ;
Definition: Match.h:50
double dec
Definition: Match.cc:41
#define LOG_LOGGER
Definition: Log.h:703
constexpr double asRadians() const noexcept
Return an Angle&#39;s value in radians.
Definition: Angle.h:166
Schema getSchema() const
Return the schema associated with the catalog&#39;s table.
Definition: Catalog.h:117
double z
Definition: Match.cc:44
AngleUnit constexpr radians
constant with units of radians
Definition: Angle.h:108
A class representing an angle.
Definition: Angle.h:127
T resize(T... args)
std::vector< Match< typename Cat1::Record, typename Cat2::Record > > unpackMatches(BaseCatalog const &matches, Cat1 const &cat1, Cat2 const &cat2)
Reconstruct a MatchVector from a BaseCatalog representation of the matches and a pair of catalogs...
Definition: Match.cc:455
LSST DM logging module built on log4cxx.
T sin(T... args)
Pass parameters to algorithms that match list of sources.
Definition: Match.h:45
T push_back(T... args)
AngleUnit constexpr degrees
constant with units of degrees
Definition: Angle.h:109
A base class for image defects.
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:63
std::shared_ptr< RecordT > src
Definition: Match.cc:48
iterator end()
Iterator access.
Definition: Catalog.h:397
Lightweight representation of a geometric match between two records.
Definition: fwd.h:101
std::vector< SimpleMatch > SimpleMatchVector
Definition: fwd.h:107
table::Schema schema
Definition: Amplifier.cc:115
Iterator class for CatalogT.
Definition: Catalog.h:38
T cos(T... args)
double x
Definition: Match.cc:42
BaseCatalog packMatches(std::vector< Match< Record1, Record2 > > const &matches)
Return a table representation of a MatchVector that can be used to persist it.
Definition: Match.cc:432
double y
Definition: Match.cc:43
T asin(T... args)
T get(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
STL class.
STL class.
A class used as a handle to a particular field in a table.
Definition: fwd.h:45
T begin(T... args)
#define LSST_MATCH_RADEC(RTYPE, C1, C2)
Definition: Match.cc:289
SortedCatalogT< SourceRecord > SourceCatalog
Definition: fwd.h:85
Reports invalid arguments.
Definition: Runtime.h:66
T isnan(T... args)
std::vector< Match< typename Cat1::Record, typename Cat2::Record > > matchRaDec(Cat1 const &cat1, Cat2 const &cat2, lsst::geom::Angle radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
Definition: Match.cc:158
SourceMatchVector matchXy(SourceCatalog const &cat1, SourceCatalog const &cat2, double radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
Definition: Match.cc:305
size_type size() const
Return the number of elements in the catalog.
Definition: Catalog.h:408
T sort(T... args)
T sqrt(T... args)
bool includeMismatches
"Include failed matches (i.e. one &#39;match&#39; is NULL) " "(default: false)" ;
Definition: Match.h:56
iterator begin()
Iterator access.
Definition: Catalog.h:396
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
SortedCatalogT< SimpleRecord > SimpleCatalog
Definition: fwd.h:79
Match< SourceRecord, SourceRecord > SourceMatch
Definition: fwd.h:105
py::object result
Definition: _schema.cc:429
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
Key< T > addField(Field< T > const &field, bool doReplace=false)
Add a new field to the Schema, and return the associated Key.
Definition: Schema.cc:668
bool symmetricMatch
"Produce symmetric matches (default: true):\n" "i.e. if (s1, s2, d) is reported, then so is (s2...
Definition: Match.h:53
std::vector< SourceMatch > SourceMatchVector
Definition: fwd.h:109