LSST Applications g0aad566f14+98bdd3b178,g0dd44d6229+e0c2448475,g2079a07aa2+86d27d4dc4,g2305ad1205+28192cebb5,g2bbee38e9b+047b288a59,g337abbeb29+047b288a59,g33d1c0ed96+047b288a59,g3a166c0a6a+047b288a59,g3d1719c13e+98bdd3b178,g487adcacf7+08349b74a1,g4be5004598+98bdd3b178,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+742623b69a,g63cd9335cc+585e252eca,g858d7b2824+98bdd3b178,g88963caddf+0cb8e002cc,g88c60e1ba4+10d1129650,g99cad8db69+26979ce5a5,g9ddcbc5298+9a081db1e4,ga1e77700b3+a912195c07,gae0086650b+585e252eca,gb0e22166c9+60f28cb32d,gb1fda7adca+90d62c6f66,gb3a676b8dc+b4feba26a1,gb4b16eec92+950562e714,gba4ed39666+c2a2e4ac27,gbb8dafda3b+2fe9a6a885,gc120e1dc64+0553034b2f,gc28159a63d+047b288a59,gc3e9b769f7+f8a3f20ab6,gcf0d15dbbd+e0c2448475,gdaeeff99f8+f9a426f77a,gdb0af172c8+f61b8374e1,ge79ae78c31+047b288a59,w.2024.20
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
33
34namespace lsst {
35namespace afw {
36namespace table {
37namespace {
38
39template <typename RecordT>
40struct 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
51template <typename Record1, typename Record2>
52bool operator<(RecordPos<Record1> const &s1, RecordPos<Record2> const &s2) {
53 return (s1.dec < s2.dec);
54}
55
56struct 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
72template <typename Cat>
73size_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("lsst.afw.table.matchRaDec", "At least one source had ra or dec equal to NaN");
94 }
95 return n;
96}
97
98template size_t makeRecordPositions(SimpleCatalog const &, RecordPos<SimpleRecord> *);
99template size_t makeRecordPositions(SourceCatalog const &, RecordPos<SourceRecord> *);
100
101template <typename Cat1, typename Cat2>
102bool 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
108template <typename Cat>
109bool 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
126double 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
139lsst::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
146template <typename Cat1, typename Cat2>
148 Cat2 const &cat2,
149 lsst::geom::Angle radius,
150 bool closest) {
152 mc.findOnlyClosest = closest;
153
154 return matchRaDec(cat1, cat2, radius, mc);
155}
156
157template <typename Cat1, typename Cat2>
159 Cat2 const &cat2,
160 lsst::geom::Angle radius,
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
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 {
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(
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
238template <typename Cat>
240 lsst::geom::Angle radius,
241 bool symmetric) {
243 mc.symmetricMatch = symmetric;
244
245 return matchRaDec(cat, radius, mc);
246}
247
248template <typename Cat>
250 lsst::geom::Angle radius,
251 MatchControl const &mc) {
253 std::vector<MatchT> matches;
254
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();
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) {
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
300 mc.findOnlyClosest = closest;
301
302 return matchXy(cat1, cat2, radius, mc);
303}
304
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) {
378 }
379 }
380 return matches;
381}
382
383SourceMatchVector matchXy(SourceCatalog const &cat, double radius, bool symmetric) {
385 mc.symmetricMatch = symmetric;
386
387 return matchXy(cat, radius, mc);
388}
389
390SourceMatchVector matchXy(SourceCatalog const &cat, double radius, MatchControl const &mc) {
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
431template <typename Record1, typename Record2>
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 using Iter = typename std::vector<Match<Record1, Record2>>::const_iterator;
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
454template <typename Cat1, typename Cat2>
456 Cat1 const &first,
457 Cat2 const &second) {
458 LOG_LOGGER tableLog = LOG_GET("lsst.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 {
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 {
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
py::object result
Definition _schema.cc:429
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition Log.h:659
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition Log.h:75
#define LOG_LOGGER
Definition Log.h:714
std::shared_ptr< RecordT > src
Definition Match.cc:48
double z
Definition Match.cc:44
#define LSST_MATCH_RADEC(RTYPE, C1, C2)
Definition Match.cc:228
double dec
Definition Match.cc:41
int y
Definition SpanSet.cc:48
table::Schema schema
Definition python.h:134
T asin(T... args)
Tag types used to declare specialized field types.
Definition misc.h:31
size_type size() const
Return the number of elements in the catalog.
Definition Catalog.h:412
iterator begin()
Iterator access.
Definition Catalog.h:400
Schema getSchema() const
Return the schema associated with the catalog's table.
Definition Catalog.h:117
Pass parameters to algorithms that match list of sources.
Definition Match.h:45
Defines the fields and offsets for a table.
Definition Schema.h:51
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
typename Base::const_iterator const_iterator
A class representing an angle.
Definition Angle.h:128
constexpr double asRadians() const noexcept
Return an Angle's value in radians.
Definition Angle.h:173
Reports invalid arguments.
Definition Runtime.h:66
Reports when the result of an operation cannot be represented by the destination type.
Definition Runtime.h:115
T cos(T... args)
T isnan(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
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
Match< SourceRecord, SourceRecord > SourceMatch
Definition fwd.h:105
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
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
AngleUnit constexpr degrees
constant with units of degrees
Definition Angle.h:110
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
T push_back(T... args)
T resize(T... args)
T sin(T... args)
T sort(T... args)
T sqrt(T... args)