LSSTApplications  20.0.0
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.");
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 
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 
491  SourceCatalog const &);
493 } // namespace table
494 } // namespace afw
495 } // namespace lsst
lsst::geom::degrees
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
schema
table::Schema schema
Definition: Amplifier.cc:115
lsst::afw::table::CatalogT::end
iterator end()
Definition: Catalog.h:397
LOG_LOGGER
#define LOG_LOGGER
Definition: Log.h:703
lsst::afw::table::MatchControl::symmetricMatch
bool symmetricMatch
"Produce symmetric matches (default: true):\n" "i.e. if (s1, s2, d) is reported, then so is (s2,...
Definition: Match.h:53
lsst::afw::table::matchRaDec
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
std::shared_ptr< RecordT >
lsst::afw::table._match.second
second
Definition: _match.py:78
Match.h
std::cos
T cos(T... args)
std::asin
T asin(T... args)
lsst::afw::table::unpackMatches
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
std::vector
STL class.
lsst::pex::exceptions::RangeError
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
dec
double dec
Definition: Match.cc:41
lsst::afw::table._match.first
first
Definition: _match.py:76
lsst::afw
Definition: imageAlgorithm.dox:1
LOG_GET
#define LOG_GET(logger)
Definition: Log.h:75
std::unique_ptr::get
T get(T... args)
lsst::afw::table::Schema
Defines the fields and offsets for a table.
Definition: Schema.h:50
lsst::afw::table::packMatches
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
std::sort
T sort(T... args)
std::sqrt
T sqrt(T... args)
lsst::afw::table::SimpleCatalog
SortedCatalogT< SimpleRecord > SimpleCatalog
Definition: fwd.h:79
src
std::shared_ptr< RecordT > src
Definition: Match.cc:48
lsst::afw::table::SourceCatalog
SortedCatalogT< SourceRecord > SourceCatalog
Definition: fwd.h:85
std::vector::push_back
T push_back(T... args)
lsst::afw::table::matchXy
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
lsst::afw::table::CatalogT::size
size_type size() const
Return the number of elements in the catalog.
Definition: Catalog.h:408
LOGLS_WARN
#define LOGLS_WARN(logger, message)
Definition: Log.h:648
lsst::geom::Angle::asRadians
constexpr double asRadians() const noexcept
Return an Angle's value in radians.
Definition: Angle.h:166
std::isnan
T isnan(T... args)
Angle.h
z
double z
Definition: Match.cc:44
lsst::afw::table._source.SourceCatalog
Definition: _source.py:33
lsst::afw::table::CatalogT::getSchema
Schema getSchema() const
Return the schema associated with the catalog's table.
Definition: Catalog.h:117
lsst::afw::table::Key
A class used as a handle to a particular field in a table.
Definition: fwd.h:45
lsst::afw::table::CatalogIterator
Iterator class for CatalogT.
Definition: Catalog.h:39
lsst::afw::table::MatchControl::findOnlyClosest
bool findOnlyClosest
"Return only the closest match if more than one is found " "(default: true)" ;
Definition: Match.h:50
std::int64_t
lsst::afw::table::Match
Lightweight representation of a geometric match between two records.
Definition: fwd.h:101
result
py::object result
Definition: _schema.cc:429
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::sin
T sin(T... args)
y
double y
Definition: Match.cc:43
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
lsst.pipe.tasks.insertFakes.radius
radius
Definition: insertFakes.py:288
x
double x
Definition: Match.cc:42
lsst::afw::table::SourceMatch
Match< SourceRecord, SourceRecord > SourceMatch
Definition: fwd.h:105
lsst::afw::table::MatchControl::includeMismatches
bool includeMismatches
"Include failed matches (i.e. one 'match' is NULL) " "(default: false)" ;
Definition: Match.h:56
lsst::geom::Angle
A class representing an angle.
Definition: Angle.h:127
lsst::geom::radians
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
lsst::afw::table::MatchControl
Pass parameters to algorithms that match list of sources.
Definition: Match.h:45
lsst::afw::table::BaseRecord::set
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
lsst::afw::table::CatalogT::begin
iterator begin()
Iterator access.
Definition: Catalog.h:396
LSST_MATCH_RADEC
#define LSST_MATCH_RADEC(RTYPE, C1, C2)
Definition: Match.cc:289
std::unique_ptr
STL class.
lsst::afw::table::CatalogT< BaseRecord >
Log.h
LSST DM logging module built on log4cxx.
lsst::afw::table::SortedCatalogT
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Definition: fwd.h:63
exceptions.h