LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
ListMatch.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * This file is part of jointcal.
4 *
5 * Developed for the LSST Data Management System.
6 * This product includes software developed by the LSST Project
7 * (https://www.lsst.org).
8 * See the COPYRIGHT file at the top-level directory of this distribution
9 * for details of code ownership.
10 *
11 * This program is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see <https://www.gnu.org/licenses/>.
23 */
24
25#include <iostream>
26#include <cmath>
27#include <list>
28#include <memory>
29#include <algorithm>
30#ifndef M_PI
31#define M_PI 3.14159265358979323846 /* pi */
32#endif
33
34#include "lsst/log/Log.h"
42
43namespace {
44LOG_LOGGER _log = LOG_GET("lsst.jointcal.ListMatch");
45}
46
47namespace lsst {
48namespace jointcal {
49
50// cuts.. limits, etc for combinatorial match
51
52/* a Segment is a pair of stars form the same image. it is used for matching starlists */
53
54struct Segment {
55 /* data */
56 double r, dx, dy;
58 int s1rank;
59
60 /* constructor (could set last argument to identity by default) */
63 s1rank = star1Rank;
64 s1 = std::move(star1);
65 s2 = std::move(star2);
66 Point P1 = transform.apply(*s1);
67 Point P2 = transform.apply(*s2);
68 dx = P2.x - P1.x;
69 dy = P2.y - P1.y;
70 r = sqrt(dx * dx + dy * dy);
71 }
72
73 /* arg(other/(*this)) if considered as complex(dx,dy) */
74 double relativeAngle(Segment *other) const {
75 return atan2(other->dx * dy - dx * other->dy, dx * other->dx + dy * other->dy);
76 }
77
78 friend std::ostream &operator<<(std::ostream &stream, const Segment &segment) {
79 stream << " dx " << segment.dx << " dy " << segment.dy << " r " << segment.r << std::endl;
80 return stream;
81 }
82};
83
84class SegmentList : public std::list<Segment> {
85public:
86 // SegmentList(const BaseStarList &list, const int nStar);
87 SegmentList(const BaseStarList &list, int nStar,
89};
90
91using SegmentIterator = std::list<Segment>::iterator;
92using SegmentCIterator = std::list<Segment>::const_iterator;
93
94static bool DecreasingLength(const Segment &first, const Segment &second) { return (first.r > second.r); }
95
97 BaseStarCIterator siStop;
98
99 /* find the fence */
100 siStop = list.begin();
101 int limit = std::min(nStars, int(list.size())) - 1; // -1 because test happens after incrementation
102 for (int count = 0; count < limit; count++) ++siStop;
103
104 // iterate on star pairs
105 int rank = 0;
106 for (auto si1 = list.begin(); si1 != siStop; ++si1, rank++)
107 for (auto si2 = siStop; si2 != si1; --si2) {
108 push_back(Segment(*si1, *si2, rank, transform));
109 }
110 sort(DecreasingLength); /* allows a break in loops */
111}
112
113//#include <pair>
114
115struct SegmentPair : public std::pair<Segment *, Segment *> {
117};
118
120using SegmentPairListIterator = SegmentPairList::iterator;
121using SegmentPairListCIterator = SegmentPairList::const_iterator;
122
123static std::unique_ptr<StarMatchList> MatchListExtract(const SegmentPairList &pairList, int rank1, int rank2,
125 /* first Select in the segment pairs list the ones which make use of star rank1 in segment1
126 and star s2 in segment2 */
127
129
130 for (const auto & a_pair : pairList) {
131 if (a_pair.first->s1rank != rank1 || a_pair.second->s1rank != rank2) continue;
132 /* now we store as star matches both ends of segment pairs ,
133 but only once the beginning of segments because they all have the same,
134 given the selection 3 lines above */
135 if (matchList->size() == 0)
136 matchList->push_back(StarMatch(transform.apply(*(a_pair.first->s1)), *(a_pair.second->s1),
137 a_pair.first->s1, a_pair.second->s1));
138 /* always store the match at end */
139 matchList->push_back(StarMatch(transform.apply(*(a_pair.first->s2)), *(a_pair.second->s2),
140 a_pair.first->s2, a_pair.second->s2));
141 }
142 return matchList;
143}
144
145static bool DecreasingQuality(const std::unique_ptr<StarMatchList> &first,
146 const std::unique_ptr<StarMatchList> &second) {
147 int idiff = first->size() - second->size();
148 if (idiff != 0)
149 return (idiff > 0);
150 else
151 return (first->getDist2() < second->getDist2());
152}
153
154/* many matching solutions (StarMatchList) will be compared. Store them in a SolList : */
155
157
158/* This one searches a general transformation by histogramming the relative size and orientation
159of star pairs ( Segment's) built from the 2 lists */
160
161static std::unique_ptr<StarMatchList> ListMatchupRotShift_Old(BaseStarList &list1, BaseStarList &list2,
163 const MatchConditions &conditions) {
164 SegmentList sList1(list1, conditions.nStarsList1, transform);
165 SegmentList sList2(list2, conditions.nStarsList2, AstrometryTransformIdentity());
166
167 /* choose the binning of the histogram so that
168 1: ratio = 1 and rotation angle = n * (pi/2) are bin centers. since
169 the angle is computed using atan2, its range is [-pi,pi],
170 and the histogram range is [-pi-eps, pi-eps], so
171 if (angle>pi- angleOffset) angle -= 2*pi before filling. */
172 int nBinsR = 21;
173 int nBinsAngle = 180; /* can be divided by 4 */
174 double angleOffset = M_PI / nBinsAngle;
175 double minRatio = conditions.minSizeRatio();
176 double maxRatio = conditions.maxSizeRatio();
177 Histo2d histo(nBinsR, minRatio, maxRatio, nBinsAngle, -M_PI - angleOffset, M_PI - angleOffset);
178
179 SegmentIterator segi1, segi2;
180 Segment *seg1, *seg2;
181 double ratio, angle;
182 for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
183 seg1 = &(*segi1);
184 if (seg1->r == 0) continue;
185 for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
186 seg2 = &(*segi2);
187 /* if one considers the 2 segments as complex numbers z1 and z2, ratio=mod(z1/z2) and angle =
188 * arg(z1/z2) */
189 /* I did not put a member function in Segment to compute both because we apply a cut on ratio
190 before actually
191 computing the angle (which involves a call to atan2 (expensive)) */
192 ratio = seg2->r / seg1->r;
193 if (ratio > maxRatio) continue;
194 if (ratio < minRatio) break; /* use th fact that segment lists are sorted by decresing length */
195 angle = seg1->relativeAngle(seg2);
196 if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
197 histo.fill(ratio, angle);
198 }
199 }
200 double binr, bina;
201 histo.binWidth(binr, bina);
202
203 SolList Solutions;
204 /* now we want to find in the (r,theta) bins that have the highest counts, the star pair
205 (one in l1, one in list2) that contribute to the largest number of segment pairs in this bin :
206 so, we histogram a couple of integer that uniquely defines the stars, for the segment pairs
207 that contribute to the maximum bin. We choose to histogram the rank of s1 of segment 1
208 versus the rank of s1 for segment 2 */
209
210 for (int i = 0; i < conditions.maxTrialCount; ++i) {
211 double ratioMax, angleMax;
212 double maxContent = histo.maxBin(ratioMax, angleMax);
213 histo.fill(ratioMax, angleMax, -maxContent);
214
215 if (conditions.printLevel >= 1)
216 LOGLS_DEBUG(_log, " valMax " << maxContent << " ratio " << ratioMax << " angle " << angleMax);
217
218 minRatio = ratioMax - binr / 2;
219 maxRatio = ratioMax + binr / 2;
220 double minAngle = angleMax - bina / 2;
221 double maxAngle = angleMax + bina / 2;
222 SegmentPairList pairList;
223 Histo2d historank(conditions.nStarsList1, 0., conditions.nStarsList1, conditions.nStarsList2, 0.,
224 conditions.nStarsList2);
225 /* reloop on segment pairs to select the ones in this specific bin */
226
227 for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
228 seg1 = &(*segi1);
229 if (seg1->r == 0) continue;
230 for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
231 seg2 = &(*segi2);
232 ratio = seg2->r / seg1->r;
233 if (ratio > maxRatio) continue;
234 if (ratio < minRatio)
235 break; /* use the fact that segment lists are sorted by decresing length */
236 angle = seg1->relativeAngle(seg2);
237 if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
238 if (angle < minAngle || angle > maxAngle) continue;
239 pairList.push_back(SegmentPair(seg1, seg2)); /* store the match */
240 historank.fill(seg1->s1rank + 0.5, seg2->s1rank + 0.5);
241 }
242 }
243 for (int iteration = 0; iteration < conditions.maxTrialCount; iteration++) {
244 double dr1, dr2;
245 double maxval = historank.maxBin(dr1, dr2);
246 /* set this bin to zero so that next iteration will find next maximum */
247 historank.fill(dr1, dr2, -maxval);
248 auto a_list = MatchListExtract(pairList, int(dr1), int(dr2), AstrometryTransformIdentity());
249 a_list->refineTransform(conditions.nSigmas); // mandatory for the sorting fields to be filled
250 Solutions.push_back(std::move(a_list));
251 }
252 } /* end of loop on (r,theta) bins */
253 Solutions.sort(DecreasingQuality);
255 best.swap(*Solutions.begin());
256 /* remove the first one from the list */
257 Solutions.pop_front();
258 if (conditions.printLevel >= 1) {
259 LOGLS_DEBUG(_log, "Best solution " << best->computeResidual() << " npairs " << best->size());
260 LOGLS_DEBUG(_log, *(best->getTransform()));
261 LOGLS_DEBUG(_log, "Chi2 " << best->getChi2() << ',' << " Number of solutions " << Solutions.size());
262 }
263 return best;
264}
265
266/* this matching routine searches brutally a match between lists in
267 the 4 parameter space: size ratio, rotation angle, x and y
268 shifts. This is done by histogramming where combinations of four
269 objets (2 on each list) fall in this 4 parameter space.
270
271 One trick is that rather than using actual offsets, we histogram
272 object indices of the combination:
273*/
274
275static std::unique_ptr<StarMatchList> ListMatchupRotShift_New(BaseStarList &list1, BaseStarList &list2,
276 const AstrometryTransform &transform,
277 const MatchConditions &conditions) {
278 if (list1.size() <= 4 || list2.size() <= 4) {
279 LOGL_FATAL(_log, "ListMatchupRotShift_New : (at least) one of the lists is too short.");
280 return nullptr;
281 }
282
283 SegmentList sList1(list1, conditions.nStarsList1, transform);
284 SegmentList sList2(list2, conditions.nStarsList2, AstrometryTransformIdentity());
285
286 /* choose the binning of the histogram so that
287 1: ratio = 1 and rotation angle = n * (pi/2) are bin centers. since
288 the angle is computed using atan2, its range is [-pi,pi],
289 and the histogram range is [-pi-eps, pi-eps], so
290 if (angle>pi- angleOffset) angle -= 2*pi before filling. */
291 int nBinsR = 21;
292 int nBinsAngle = 180; /* can be divided by 4 */
293 double angleOffset = M_PI / nBinsAngle;
294 double minRatio = conditions.minSizeRatio();
295 double maxRatio = conditions.maxSizeRatio();
296 SparseHisto4d histo(nBinsR, minRatio, maxRatio, nBinsAngle, -M_PI - angleOffset, M_PI - angleOffset,
297 conditions.nStarsList1, 0., conditions.nStarsList1, conditions.nStarsList2, 0.,
298 conditions.nStarsList2, sList1.size() * sList2.size());
299
300 SegmentIterator segi1, segi2;
301 Segment *seg1, *seg2;
302 double ratio, angle;
303
304 for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
305 seg1 = &(*segi1);
306 if (seg1->r == 0) continue;
307 for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
308 seg2 = &(*segi2);
309 /* if one considers the 2 segments as complex numbers z1 and z2, ratio=mod(z1/z2) and angle =
310 * arg(z1/z2) */
311 /* I did not put a member function in Segment to compute both because we apply a cut on ratio
312 before actually
313 computing the angle (which involves a call to atan2 (expensive)) */
314 ratio = seg2->r / seg1->r;
315 if (ratio > maxRatio) continue;
316 if (ratio < minRatio) break; /* use th fact that segment lists are sorted by decresing length */
317 angle = seg1->relativeAngle(seg2);
318 if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
319 histo.fill(ratio, angle, seg1->s1rank + 0.5, seg2->s1rank + 0.5);
320 }
321 }
322
323 SolList Solutions;
324 /* now we find the highest bins of the histogram, and recover the original objects.
325 This involves actually re-looping on the combinations, but it is much
326 faster that the original histogram filling loop, since we only compute
327 angle and ratio for Segments that have the right first object
328 */
329
330 int oldMaxContent = 0;
331
332 for (int i = 0; i < 4 * conditions.maxTrialCount;
333 ++i) // leave a limit to make avoid (almost) infinite loops
334 {
335 double pars[4];
336 int maxContent = histo.maxBin(pars);
337 if (maxContent == 0) break;
338 if (conditions.printLevel >= 1) {
339 LOGLS_DEBUG(_log, "ValMax " << maxContent << " ratio " << pars[0] << " angle " << pars[1]);
340 }
341 histo.zeroBin(pars);
342 if (i > 0) { /* the match possibilities come out in a random order when they have the same content.
343 so, we stop investigating guesses when the content goes down AND the requested search
344 depth
345 (maxTrialCount) is reached */
346 if (maxContent < oldMaxContent && i >= conditions.maxTrialCount) break;
347 }
348 oldMaxContent = maxContent;
349 /* reloop on segment pairs to select the ones in this specific bin */
350 int rank1L1 = int(pars[2]);
351 int rank1L2 = int(pars[3]);
352 double minAngle, maxAngle;
353 histo.binLimits(pars, 0, minRatio, maxRatio);
354 histo.binLimits(pars, 1, minAngle, maxAngle);
355
356 std::unique_ptr<StarMatchList> a_list(new StarMatchList);
357
358 for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
359 seg1 = &(*segi1);
360 if (seg1->s1rank != rank1L1) continue;
361 if (seg1->r == 0) continue;
362 for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
363 seg2 = &(*segi2);
364 if (seg2->s1rank != rank1L2) continue;
365 // push in the list the match corresponding to end number 1 of segments
366 if (a_list->empty())
367 a_list->push_back(StarMatch(*(seg1->s1), *(seg2->s1), seg1->s1, seg2->s1));
368 ratio = seg2->r / seg1->r;
369 if (ratio > maxRatio) continue;
370 if (ratio < minRatio)
371 break; /* use the fact that segment lists are sorted by decresing length */
372 angle = seg1->relativeAngle(seg2);
373 if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
374 if (angle < minAngle || angle > maxAngle) continue;
375 /* here we have 2 segments which have the right
376 - length ratio
377 - relative angle
378 - first objects (objects on the end number 1).
379 The objects on the end number 2 are the actual matches : */
380 a_list->push_back(StarMatch(*(seg1->s2), *(seg2->s2), seg1->s2, seg2->s2));
381 }
382 }
383
384 // a basic check for sanity of the algorithm :
385
386 if (int(a_list->size()) != maxContent + 1) {
387 LOGLS_ERROR(_log, "There is an internal inconsistency in ListMatchupRotShift.");
388 LOGLS_ERROR(_log, "maxContent = " << maxContent);
389 LOGLS_ERROR(_log, "matches->size() = " << a_list->size());
390 }
391 a_list->refineTransform(conditions.nSigmas);
392 Solutions.push_back(std::move(a_list));
393 }
394
395 if (Solutions.empty()) {
396 LOGLS_ERROR(_log, "Error In ListMatchup : not a single pair match.");
397 LOGLS_ERROR(_log, "Probably, the relative scale of lists is not within bounds.");
398 LOGLS_ERROR(_log, "min/max ratios: " << minRatio << ' ' << maxRatio);
399 return nullptr;
400 }
401
402 Solutions.sort(DecreasingQuality);
404 best.swap(*Solutions.begin());
405 /* remove the first one from the list */
406 Solutions.pop_front();
407 if (conditions.printLevel >= 1) {
408 LOGLS_INFO(_log, "Best solution " << best->computeResidual() << " npairs " << best->size());
409 LOGLS_INFO(_log, *(best->getTransform()));
410 LOGLS_INFO(_log, "Chi2 " << best->getChi2() << ", Number of solutions " << Solutions.size());
411 }
412 return best;
413}
414
415static std::unique_ptr<StarMatchList> ListMatchupRotShift(BaseStarList &list1, BaseStarList &list2,
416 const AstrometryTransform &transform,
417 const MatchConditions &conditions) {
418 if (conditions.algorithm == 1)
419 return ListMatchupRotShift_Old(list1, list2, transform, conditions);
420 else
421 return ListMatchupRotShift_New(list1, list2, transform, conditions);
422}
423
425 const MatchConditions &conditions) {
426 list1.fluxSort();
427 list2.fluxSort();
428
429 return ListMatchupRotShift(list1, list2, AstrometryTransformIdentity(), conditions);
430}
431
433 const MatchConditions &conditions) {
434 list1.fluxSort();
435 list2.fluxSort();
436
437 AstrometryTransformLinear flip(0, 0, 1, 0, 0, -1);
438 std::unique_ptr<StarMatchList> flipped(ListMatchupRotShift(list1, list2, flip, conditions));
440 ListMatchupRotShift(list1, list2, AstrometryTransformIdentity(), conditions));
441 if (!flipped || !unflipped) return std::unique_ptr<StarMatchList>(nullptr);
442 if (conditions.printLevel >= 1) {
443 LOGLS_DEBUG(_log,
444 "unflipped Residual " << unflipped->computeResidual() << " nused " << unflipped->size());
445 LOGLS_DEBUG(_log, "flipped Residual " << flipped->computeResidual() << " nused " << flipped->size());
446 }
447 if (DecreasingQuality(flipped, unflipped)) {
448 if (conditions.printLevel >= 1) LOGL_DEBUG(_log, "Keeping flipped solution.");
449 // One should NOT apply the flip to the result because the matchlist
450 // (even the flipped one) contains the actual coordinates of stars.
451 // MatchListExtract is always called with AstrometryTransformIdentity() as last parameter
452 return flipped;
453 } else {
454 if (conditions.printLevel >= 1) LOGL_DEBUG(_log, "Keeping unflipped solution.");
455 return unflipped;
456 }
457}
458
459#ifdef STORAGE
460// timing : 2.5 s for l1 of 1862 objects and l2 of 2617 objects
462 const BaseStarList &list2,
463 const AstrometryTransform &transform,
464 double maxShift) {
465 int ncomb = list1.size() * list2.size();
466 if (!ncomb) return nullptr;
467 int nx;
468 if (ncomb > 10000)
469 nx = 100;
470 else
471 nx = (int)sqrt(ncomb);
472
473 Histo2d histo(nx, -maxShift, maxShift, nx, -maxShift, maxShift);
474
475 BaseStarCIterator s1, s2;
476 double x1, y1;
477 for (s1 = list1.begin(); s1 != list1.end(); ++s1) {
478 transform.apply((*s1)->x, (*s1)->y, x1, y1);
479 for (s2 = list2.begin(); s2 != list2.end(); ++s2) {
480 histo.fill((*s2)->x - x1, (*s2)->y - y1);
481 }
482 }
483 double dx = 0, dy = 0;
484 histo.maxBin(dx, dy);
485 return std::unique_ptr<AstrometryTransformLinear>(new AstrometryTransformLinearShift(dx, dy));
486}
487#endif /*STORAGE*/
488
489// timing : 140 ms for l1 of 1862 objects and l2 of 2617 objects (450 MHz, "-O4") maxShift = 200.
491 const BaseStarList &list2,
493 double maxShift, double binSize) {
494 int nx;
495 if (binSize == 0) {
496 int ncomb = list1.size() * list2.size();
497 if (ncomb > 10000)
498 nx = 100;
499 else
500 nx = (int)sqrt(double(ncomb));
501 if (!ncomb) return std::unique_ptr<AstrometryTransformLinear>(nullptr);
502 } else
503 nx = int(2 * maxShift / binSize + 0.5);
504
505 Histo2d histo(nx, -maxShift, maxShift, nx, -maxShift, maxShift);
506 double binSizeNew = 2 * maxShift / nx;
507
509 FastFinder finder(list2);
510 double x1, y1;
511 for (s1 = list1.begin(); s1 != list1.end(); ++s1) {
512 transform.apply((*s1)->x, (*s1)->y, x1, y1);
513 FastFinder::Iterator it = finder.beginScan(Point(x1, y1), maxShift);
514 while (*it) {
515 auto s2 = *it;
516 histo.fill(s2->x - x1, s2->y - y1);
517 ++it;
518 }
519 }
520 SolList Solutions;
521 for (int i = 0; i < 4; ++i) {
522 double dx = 0, dy = 0;
523 double count = histo.maxBin(dx, dy);
524 histo.fill(dx, dy, -count); // zero the maxbin
525 AstrometryTransformLinearShift shift(dx, dy);
526 auto newGuess = compose(shift, transform);
527 auto raw_matches = listMatchCollect(list1, list2, newGuess.get(), binSizeNew);
529 raw_matches->applyTransform(*matches, &transform);
530 matches->setTransformOrder(1);
531 matches->refineTransform(3.);
532 Solutions.push_back(std::move(matches));
533 }
534 Solutions.sort(DecreasingQuality);
536 new AstrometryTransformLinear(*std::const_pointer_cast<AstrometryTransformLinear>(
537 std::dynamic_pointer_cast<const AstrometryTransformLinear>(
538 Solutions.front()->getTransform()))));
539 return best;
540}
541
542#ifdef STORAGE
543
544// this is the old fashioned way...
545
546std::unique_ptr<StarMatchList> listMatchCollect_Slow(const BaseStarList &list1, const BaseStarList &list2,
547 const AstrometryTransform *guess, const double maxDist) {
548 std::unique_ptr<StarMatchList> matches(new StarMatchList);
549 /****** Collect ***********/
550 for (BaseStarCIterator si = list1.begin(); si != list1.end(); ++si) {
551 const Point *p1 = (*si);
552 const Point p2 = guess->apply(*p1);
553 const BaseStar *neighbour = list2.findClosest(p2);
554 if (!neighbour) continue;
555 double distance = p2.Distance(*neighbour);
556 if (distance < maxDist) {
557 matches->push_back(StarMatch(*p1, *neighbour, *si, neighbour));
558 // assign the distance, since we have it in hand:
559 matches->back().distance = distance;
560 }
561 }
562 return matches;
563}
564#endif
565
566// here is the real active routine:
567
569 const AstrometryTransform *guess, const double maxDist) {
571 /****** Collect ***********/
572 FastFinder finder(list2);
573 for (const auto & si : list1) {
574 auto p1 = si;
575 Point p2 = guess->apply(*p1);
576 auto neighbour = finder.findClosest(p2, maxDist);
577 if (!neighbour) continue;
578 double distance = p2.Distance(*neighbour);
579 if (distance < maxDist) {
580 matches->push_back(StarMatch(*p1, *neighbour, p1, neighbour));
581 // assign the distance, since we have it in hand:
582 matches->back().distance = distance;
583 }
584 }
585 matches->setTransform(guess);
586
587 return matches;
588}
589
590#ifdef STORAGE
591// unused
593std::unique_ptr<StarMatchList> CollectAndFit(const BaseStarList &list1, const BaseStarList &list2,
594 const AstrometryTransform *guess, const double maxDist) {
595 const AstrometryTransform *bestTransform = guess;
597 while (true) {
598 auto m = listMatchCollect(list1, list2, bestTransform, maxDist);
599 m->setTransform(bestTransform);
600 m->refineTransform(3.);
601 LOGLS_INFO(_log, "Iterating: resid " << m->computeResidual() << " size " << m->size());
602 if (!prevMatch ||
603 (prevMatch && m->computeResidual() < prevMatch->computeResidual() * 0.999 && m->Chi2() > 0)) {
604 prevMatch.swap(m);
605 bestTransform = prevMatch->Transform();
606 } else {
607 break;
608 }
609 }
610 return prevMatch;
611}
612#endif
613
615 const double maxDist) {
617 FastFinder finder(list2);
618 for (const auto & si : list1) {
619 auto p1 = si;
620 auto neighbour = finder.findClosest(*p1, maxDist);
621 if (!neighbour) continue;
622 double distance = p1->Distance(*neighbour);
623 if (distance < maxDist) {
624 matches->push_back(StarMatch(*p1, *neighbour, p1, neighbour));
625 // assign the distance, since we have it in hand:
626 matches->back().distance = distance;
627 }
628 }
629
630 matches->setTransform(std::make_shared<AstrometryTransformIdentity>());
631
632 return matches;
633}
634
635static bool is_transform_ok(const StarMatchList *match, double pixSizeRatio2, const size_t nmin) {
636 if ((fabs(fabs(std::dynamic_pointer_cast<const AstrometryTransformLinear>(match->getTransform())
637 ->determinant()) -
638 pixSizeRatio2) /
639 pixSizeRatio2 <
640 0.2) &&
641 (match->size() > nmin))
642 return true;
643 LOGL_ERROR(_log, "transform is not ok!");
644 match->printTransform();
645 return false;
646}
647
648// utility to check current transform difference
649static double transform_diff(const BaseStarList &List, const AstrometryTransform *T1,
650 const AstrometryTransform *T2) {
651 double diff2 = 0;
652 FatPoint tf1;
653 Point tf2;
654 int count = 0;
655 for (const auto & it : List) {
656 const BaseStar &s = *it;
657 T1->transformPosAndErrors(s, tf1);
658 T2->apply(s, tf2);
659 double dx = tf1.x - tf2.x;
660 double dy = tf1.y - tf2.y;
661 diff2 += (tf1.vy * dx * dx + tf1.vx * dy * dy - 2 * tf1.vxy * dx * dy) /
662 (tf1.vx * tf1.vy - tf1.vxy * tf1.vxy);
663 count++;
664 }
665 if (count) return diff2 / double(count);
666 return 0;
667}
668
669static double median_distance(const StarMatchList *match, const AstrometryTransform *transform) {
670 size_t nstars = match->size();
671 std::vector<double> resid(nstars);
672 auto ir = resid.begin();
673 for (auto it = match->begin(); it != match->end(); ++it, ++ir)
674 *ir = sqrt(transform->apply(it->point1).computeDist2(it->point2));
675 sort(resid.begin(), resid.end());
676 return (nstars & 1) ? resid[nstars / 2] : (resid[nstars / 2 - 1] + resid[nstars / 2]) * 0.5;
677}
678
680 const BaseStarList &List2,
681 const MatchConditions &conditions) {
682 BaseStarList list1, list2;
683 List1.copyTo(list1);
684 list1.fluxSort();
685 List2.copyTo(list2);
686 list2.fluxSort();
687
688 LOGLS_INFO(_log, "listMatchCombinatorial: find match between " << list1.size() << " and " << list2.size()
689 << " stars...");
690 auto match = matchSearchRotShiftFlip(list1, list2, conditions);
691 double pixSizeRatio2 = std::pow(conditions.sizeRatio, 2);
692 size_t nmin =
693 std::min(size_t(10), size_t(std::min(List1.size(), List2.size()) * conditions.minMatchRatio));
694
696 if (is_transform_ok(match.get(), pixSizeRatio2, nmin))
697 transform = match->getTransform()->clone();
698 else {
699 LOGL_ERROR(_log, "listMatchCombinatorial: direct transform failed, trying reverse");
700 match = matchSearchRotShiftFlip(list2, list1, conditions);
701 if (is_transform_ok(match.get(), pixSizeRatio2, nmin))
702 transform = match->inverseTransform();
703 else {
704 LOGL_FATAL(_log, "FAILED");
705 }
706 }
707
708 if (transform) {
709 LOGL_INFO(_log, "FOUND");
710 if (conditions.printLevel >= 1) {
711 LOGL_DEBUG(_log, " listMatchCombinatorial: found the following transform.");
712 LOGLS_DEBUG(_log, *transform);
713 }
714 } else
715 LOGL_ERROR(_log, "listMatchCombinatorial: failed to find a transform");
716 return transform;
717}
718
721 const int maxOrder) {
722 if (!transform) {
724 }
725
726 // some hard-coded constants that could go in a param file
727 const double brightDist = 2.; // distance in pixels in a match
728 const double fullDist = 4.; // distance in pixels in a match between entire lists
729 const double nSigmas = 3.; // k-sigma clipping on residuals
730 const size_t nStars = 500; // max number of bright stars to fit
731
732 int order = 1;
733 size_t nstarmin = 3;
734
735 BaseStarList list1, list2;
736 List1.copyTo(list1);
737 list1.fluxSort();
738 list1.cutTail(nStars);
739 List2.copyTo(list2);
740 list2.fluxSort();
741 list2.cutTail(nStars);
742
743 auto fullMatch = listMatchCollect(List1, List2, transform.get(), fullDist);
744 auto brightMatch = listMatchCollect(list1, list2, transform.get(), brightDist);
745 double curChi2 = computeChi2(*brightMatch, *transform) / brightMatch->size();
746
747 LOGLS_INFO(_log, "listMatchRefine: start: med.resid " << median_distance(fullMatch.get(), transform.get())
748 << " #match " << fullMatch->size());
749
750 do { // loop on transform order on full list of stars
751 auto curTransform = brightMatch->getTransform()->clone();
752 unsigned iter = 0;
753 double transDiff;
754 do { // loop on transform diff only on bright stars
755 brightMatch->setTransformOrder(order);
756 brightMatch->refineTransform(nSigmas);
757 transDiff = transform_diff(list1, brightMatch->getTransform().get(), curTransform.get());
758 curTransform = brightMatch->getTransform()->clone();
759 brightMatch = listMatchCollect(list1, list2, curTransform.get(), brightDist);
760 } while (brightMatch->size() > nstarmin && transDiff > 0.05 && ++iter < 5);
761
762 double prevChi2 = curChi2;
763 curChi2 = computeChi2(*brightMatch, *curTransform) / brightMatch->size();
764
765 fullMatch = listMatchCollect(List1, List2, curTransform.get(), fullDist);
766 LOGLS_INFO(_log, "listMatchRefine: order " << order << " med.resid "
767 << median_distance(fullMatch.get(), curTransform.get())
768 << " #match " << fullMatch->size());
769 if (((prevChi2 - curChi2) > 0.01 * curChi2) && curChi2 > 0) {
770 LOGLS_INFO(_log, " listMatchRefine: order " << order << " was a better guess.");
771 transform = brightMatch->getTransform()->clone();
772 }
773 nstarmin = brightMatch->getTransform()->getNpar();
774 } while (++order <= maxOrder);
775
776 return transform;
777}
778} // namespace jointcal
779} // namespace lsst
Fast locator in starlists.
table::Key< double > angle
#define M_PI
Definition ListMatch.cc:31
Combinatorial searches for linear transformations to go from list1 to list2.
LSST DM logging module built on log4cxx.
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition Log.h:75
#define LOGL_INFO(logger, message...)
Log a info-level message using a varargs/printf style interface.
Definition Log.h:531
#define LOGLS_INFO(logger, message)
Log a info-level message using an iostream-based interface.
Definition Log.h:639
#define LOG_LOGGER
Definition Log.h:714
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
Definition Log.h:679
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
Definition Log.h:563
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
Definition Log.h:619
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition Log.h:515
#define LOGL_FATAL(logger, message...)
Log a fatal-level message using a varargs/printf style interface.
Definition Log.h:579
int m
Definition SpanSet.cc:48
pairs of points
table::Key< int > transform
T atan2(T... args)
T begin(T... args)
a virtual (interface) class for geometric transformations.
virtual void apply(double xIn, double yIn, double &xOut, double &yOut) const =0
A do-nothing transformation. It anyway has dummy routines to mimick a AstrometryTransform.
implements the linear transformations (6 real coefficients).
just here to provide a specialized constructor, and fit.
Iterator meant to traverse objects within some limiting distance.
Definition FastFinder.h:90
This is an auxillary class for matching objects from starlists.
Definition FastFinder.h:54
std::shared_ptr< const BaseStar > findClosest(const Point &where, double maxDist, bool(*SkipIt)(const BaseStar &)=nullptr) const
Find the closest with some rejection capability.
Definition FastFinder.cc:83
Iterator beginScan(const Point &where, double maxDist) const
double maxBin(double &x, double &y) const
Definition Histo2d.cc:80
void binWidth(double &Hdx, double &Hdy) const
Definition Histo2d.h:44
void fill(float x, float y, float weight=1.)
Definition Histo2d.cc:75
A point in a plane.
Definition Point.h:37
double x
coordinate
Definition Point.h:42
double Distance(const Point &other) const
Definition Point.h:51
SegmentList(const BaseStarList &list, int nStar, const AstrometryTransform &transform=AstrometryTransformIdentity())
Definition ListMatch.cc:96
std::lists of Stars.
Definition StarList.h:58
void cutTail(int nKeep)
cuts the end of the std::list
Definition StarList.cc:48
void copyTo(StarList< Star > &copy) const
clears copy and makes a copy of the std::list to copy
Definition StarList.cc:68
void fluxSort()
a model routine to sort the std::list
Definition StarList.cc:42
A hanger for star associations.
Definition StarMatch.h:55
T count(T... args)
T distance(T... args)
T end(T... args)
T endl(T... args)
T fabs(T... args)
T front(T... args)
T min(T... args)
T move(T... args)
std::unique_ptr< AstrometryTransform > listMatchRefine(const BaseStarList &list1, const BaseStarList &list2, std::unique_ptr< AstrometryTransform > transform, int maxOrder=3)
Definition ListMatch.cc:719
std::list< Segment >::const_iterator SegmentCIterator
Definition ListMatch.cc:92
BaseStarList::const_iterator BaseStarCIterator
Definition BaseStar.h:123
std::list< std::unique_ptr< StarMatchList > > SolList
Definition ListMatch.cc:156
std::unique_ptr< AstrometryTransform > listMatchCombinatorial(const BaseStarList &list1, const BaseStarList &list2, const MatchConditions &conditions=MatchConditions())
Definition ListMatch.cc:679
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
std::unique_ptr< StarMatchList > listMatchCollect(const BaseStarList &list1, const BaseStarList &list2, const AstrometryTransform *guess, double maxDist)
assembles star matches.
Definition ListMatch.cc:568
StarList< BaseStar > BaseStarList
Definition BaseStar.h:121
double computeChi2(const StarMatchList &L, const AstrometryTransform &transform)
the actual chi2
Definition StarMatch.cc:244
SegmentPairList::const_iterator SegmentPairListCIterator
Definition ListMatch.cc:121
std::unique_ptr< StarMatchList > matchSearchRotShift(BaseStarList &list1, BaseStarList &list2, const MatchConditions &conditions)
searches a geometrical transformation that goes from list1 to list2.
Definition ListMatch.cc:424
std::unique_ptr< StarMatchList > matchSearchRotShiftFlip(BaseStarList &list1, BaseStarList &list2, const MatchConditions &conditions)
same as above but searches also a flipped solution.
Definition ListMatch.cc:432
std::list< Segment >::iterator SegmentIterator
Definition ListMatch.cc:91
std::unique_ptr< AstrometryTransformLinear > listMatchupShift(const BaseStarList &list1, const BaseStarList &list2, const AstrometryTransform &transform, double maxShift, double binSize=0)
searches for a 2 dimensional shift using a very crude histogram method.
Definition ListMatch.cc:490
SegmentPairList::iterator SegmentPairListIterator
Definition ListMatch.cc:120
STL namespace.
T pop_front(T... args)
T pow(T... args)
Segment push_back(Segment ... args)
T size(T... args)
Segment sort(Segment ... args)
T sqrt(T... args)
Parameters to be provided to combinatorial searches.
Definition ListMatch.h:40
std::shared_ptr< const BaseStar > s1
Definition ListMatch.cc:57
double relativeAngle(Segment *other) const
Definition ListMatch.cc:74
Segment(std::shared_ptr< const BaseStar > star1, std::shared_ptr< const BaseStar > star2, const int star1Rank, const AstrometryTransform &transform)
Definition ListMatch.cc:61
friend std::ostream & operator<<(std::ostream &stream, const Segment &segment)
Definition ListMatch.cc:78
std::shared_ptr< const BaseStar > s2
Definition ListMatch.cc:57
SegmentPair(Segment *f, Segment *s)
Definition ListMatch.cc:116
T swap(T... args)
table::Key< int > order