39 template <
typename RecordT>
51 template <
typename Record1,
typename Record2>
52 bool operator<(RecordPos<Record1>
const &s1, RecordPos<Record2>
const &s2) {
53 return (s1.dec < s2.dec);
58 return s1->getY() < s2->getY();
72 template <
typename Cat>
73 size_t makeRecordPositions(Cat
const &cat, RecordPos<typename Cat::Record> *positions) {
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) {
84 positions[n].dec =
dec.asRadians();
85 positions[n].x =
std::cos(ra) * cosDec;
86 positions[n].y =
std::sin(ra) * cosDec;
93 LOGLS_WARN(
"afw.table.matchRaDec",
"At least one source had ra or dec equal to NaN");
98 template size_t makeRecordPositions(
SimpleCatalog const &, RecordPos<SimpleRecord> *);
99 template size_t makeRecordPositions(
SourceCatalog const &, RecordPos<SourceRecord> *);
101 template <
typename Cat1,
typename Cat2>
102 bool doSelfMatchIfSame(
std::vector<Match<typename Cat1::Record, typename Cat2::Record> > &
result,
108 template <
typename Cat>
109 bool doSelfMatchIfSame(
std::vector<Match<typename Cat::Record, typename Cat::Record> > &
result,
111 if (&cat1 == &cat2) {
127 return 2. * (1. -
std::cos(theta.asRadians()));
146 template <
typename Cat1,
typename Cat2>
157 template <
typename Cat1,
typename Cat2>
165 if (doSelfMatchIfSame(matches, cat1, cat2,
radius))
return matches;
170 if (cat1.size() == 0 || cat2.size() == 0) {
174 double const d2Limit = toUnitSphereDistanceSquared(
radius);
177 size_t len1 = cat1.size();
178 size_t len2 = cat2.size();
180 typedef RecordPos<typename Cat1::Record> Pos1;
181 typedef RecordPos<typename Cat2::Record> Pos2;
184 len1 = makeRecordPositions(cat1, pos1.
get());
185 len2 = makeRecordPositions(cat2, pos2.
get());
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) {
196 double maxDec = pos1[i].dec +
radius.asRadians();
197 size_t closestIndex = -1;
198 double d2Include = d2Limit;
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) {
212 matches.
push_back(MatchT(pos1[i].
src, pos2[j].
src, fromUnitSphereDistanceSquared(d2)));
218 matches.
push_back(MatchT(pos1[i].
src, nullRecord, NAN));
222 MatchT(pos1[i].
src, pos2[closestIndex].
src, fromUnitSphereDistanceSquared(d2Include)));
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 &)
236 #undef LSST_MATCH_RADEC
238 template <
typename Cat>
248 template <
typename Cat>
258 if (cat.size() == 0) {
262 double const d2Limit = toUnitSphereDistanceSquared(
radius);
265 size_t len = cat.size();
266 typedef RecordPos<typename Cat::Record> Pos;
268 len = makeRecordPositions(cat, pos.
get());
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;
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 &)
296 #undef LSST_MATCH_RADEC
307 if (&cat1 == &cat2) {
314 size_t len1 = cat1.size();
315 size_t len2 = cat2.size();
342 for (
size_t i = 0, start = 0; i < len1; ++i) {
343 double y = pos1[i]->getY();
345 while (start < len2 && pos2[start]->getY() < minY) {
351 double x = pos1[i]->getX();
354 size_t closestIndex = -1;
355 double r2Include = r2;
358 for (
size_t j = start; j < len2 && (y2 = pos2[j]->getY()) <= maxY; ++j) {
359 double dx =
x - pos2[j]->getX();
361 double d2 = dx * dx + dy * dy;
362 if (d2 < r2Include) {
395 size_t len = cat.size();
410 for (
size_t i = 0; i < len; ++i) {
411 double x = pos[i]->getX();
412 double y = pos[i]->getY();
415 for (
size_t j = i + 1; j < len && (y2 = pos[j]->getY()) <= maxY; ++j) {
416 double dx =
x - pos[j]->getX();
418 double d2 = dx * dx + dy * dy;
431 template <
typename Record1,
typename Record2>
436 Key<double> keyD =
schema.addField<
double>(
"distance",
"Distance between matches sources.");
438 result.getTable()->preallocate(matches.size());
439 result.reserve(matches.size());
441 for (Iter i = matches.begin(); i != matches.end(); ++i) {
443 record->
set(outKey1, i->first->getId());
444 record->
set(outKey2, i->second->getId());
445 record->
set(keyD, i->distance);
454 template <
typename Cat1,
typename Cat2>
464 "Catalogs passed to unpackMatches must be sorted.");
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()) {
476 "Persisted match record with ID " << i->get(inKey1) <<
" not found in catalog 1.");
482 "Persisted match record with ID " << i->get(inKey2) <<
" not found in catalog 2.");
484 j->distance = i->get(keyD);