LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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"
35 #include "lsst/jointcal/BaseStar.h"
38 #include "lsst/jointcal/Histo2d.h"
39 #include "lsst/jointcal/Histo4d.h"
42 
43 namespace {
44 LOG_LOGGER _log = LOG_GET("lsst.jointcal.ListMatch");
45 }
46 
47 namespace lsst {
48 namespace 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 
54 struct 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(*star1);
67  Point P2 = transform.apply(*star2);
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) {
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 
84 class SegmentList : public std::list<Segment> {
85 public:
86  // SegmentList(const BaseStarList &list, const int nStar);
87  SegmentList(const BaseStarList &list, const int nStar,
89 };
90 
93 
94 static 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 
115 struct SegmentPair : public std::pair<Segment *, Segment *> {
116  SegmentPair(Segment *f, Segment *s) : std::pair<Segment *, Segment *>(f, s){};
117 };
118 
120 typedef SegmentPairList::iterator SegmentPairListIterator;
121 typedef SegmentPairList::const_iterator SegmentPairListCIterator;
122 
123 static 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 (SegmentPairListCIterator spi = pairList.begin(); spi != pairList.end(); spi++) {
131  const SegmentPair &a_pair = *spi;
132  if (a_pair.first->s1rank != rank1 || a_pair.second->s1rank != rank2) continue;
133  /* now we store as star matches both ends of segment pairs ,
134  but only once the beginning of segments because they all have the same,
135  given the selection 3 lines above */
136  if (matchList->size() == 0)
137  matchList->push_back(StarMatch(transform.apply(*(a_pair.first->s1)), *(a_pair.second->s1),
138  a_pair.first->s1, a_pair.second->s1));
139  /* always store the match at end */
140  matchList->push_back(StarMatch(transform.apply(*(a_pair.first->s2)), *(a_pair.second->s2),
141  a_pair.first->s2, a_pair.second->s2));
142  }
143  return matchList;
144 }
145 
146 static bool DecreasingQuality(const std::unique_ptr<StarMatchList> &first,
148  int idiff = first->size() - second->size();
149  if (idiff != 0)
150  return (idiff > 0);
151  else
152  return (first->getDist2() < second->getDist2());
153 }
154 
155 /* many matching solutions (StarMatchList) will be compared. Store them in a SolList : */
156 
158 
159 /* This one searches a general transformation by histogramming the relative size and orientation
160 of star pairs ( Segment's) built from the 2 lists */
161 
162 static std::unique_ptr<StarMatchList> ListMatchupRotShift_Old(BaseStarList &list1, BaseStarList &list2,
164  const MatchConditions &conditions) {
165  SegmentList sList1(list1, conditions.nStarsList1, transform);
166  SegmentList sList2(list2, conditions.nStarsList2, AstrometryTransformIdentity());
167 
168  /* choose the binning of the histogram so that
169  1: ratio = 1 and rotation angle = n * (pi/2) are bin centers. since
170  the angle is computed using atan2, its range is [-pi,pi],
171  and the histogram range is [-pi-eps, pi-eps], so
172  if (angle>pi- angleOffset) angle -= 2*pi before filling. */
173  int nBinsR = 21;
174  int nBinsAngle = 180; /* can be divided by 4 */
175  double angleOffset = M_PI / nBinsAngle;
176  double minRatio = conditions.minSizeRatio();
177  double maxRatio = conditions.maxSizeRatio();
178  Histo2d histo(nBinsR, minRatio, maxRatio, nBinsAngle, -M_PI - angleOffset, M_PI - angleOffset);
179 
180  SegmentIterator segi1, segi2;
181  Segment *seg1, *seg2;
182  double ratio, angle;
183  for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
184  seg1 = &(*segi1);
185  if (seg1->r == 0) continue;
186  for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
187  seg2 = &(*segi2);
188  /* if one considers the 2 segments as complex numbers z1 and z2, ratio=mod(z1/z2) and angle =
189  * arg(z1/z2) */
190  /* I did not put a member function in Segment to compute both because we apply a cut on ratio
191  before actually
192  computing the angle (which involves a call to atan2 (expensive)) */
193  ratio = seg2->r / seg1->r;
194  if (ratio > maxRatio) continue;
195  if (ratio < minRatio) break; /* use th fact that segment lists are sorted by decresing length */
196  angle = seg1->relativeAngle(seg2);
197  if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
198  histo.fill(ratio, angle);
199  }
200  }
201  double binr, bina;
202  histo.binWidth(binr, bina);
203 
204  SolList Solutions;
205  /* now we want to find in the (r,theta) bins that have the highest counts, the star pair
206  (one in l1, one in list2) that contribute to the largest number of segment pairs in this bin :
207  so, we histogram a couple of integer that uniquely defines the stars, for the segment pairs
208  that contribute to the maximum bin. We choose to histogram the rank of s1 of segment 1
209  versus the rank of s1 for segment 2 */
210 
211  for (int i = 0; i < conditions.maxTrialCount; ++i) {
212  double ratioMax, angleMax;
213  double maxContent = histo.maxBin(ratioMax, angleMax);
214  histo.fill(ratioMax, angleMax, -maxContent);
215 
216  if (conditions.printLevel >= 1)
217  LOGLS_DEBUG(_log, " valMax " << maxContent << " ratio " << ratioMax << " angle " << angleMax);
218 
219  minRatio = ratioMax - binr / 2;
220  maxRatio = ratioMax + binr / 2;
221  double minAngle = angleMax - bina / 2;
222  double maxAngle = angleMax + bina / 2;
223  SegmentPairList pairList;
224  Histo2d historank(conditions.nStarsList1, 0., conditions.nStarsList1, conditions.nStarsList2, 0.,
225  conditions.nStarsList2);
226  /* reloop on segment pairs to select the ones in this specific bin */
227 
228  for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
229  seg1 = &(*segi1);
230  if (seg1->r == 0) continue;
231  for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
232  seg2 = &(*segi2);
233  ratio = seg2->r / seg1->r;
234  if (ratio > maxRatio) continue;
235  if (ratio < minRatio)
236  break; /* use the fact that segment lists are sorted by decresing length */
237  angle = seg1->relativeAngle(seg2);
238  if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
239  if (angle < minAngle || angle > maxAngle) continue;
240  pairList.push_back(SegmentPair(seg1, seg2)); /* store the match */
241  historank.fill(seg1->s1rank + 0.5, seg2->s1rank + 0.5);
242  }
243  }
244  for (int iteration = 0; iteration < conditions.maxTrialCount; iteration++) {
245  double dr1, dr2;
246  double maxval = historank.maxBin(dr1, dr2);
247  /* set this bin to zero so that next iteration will find next maximum */
248  historank.fill(dr1, dr2, -maxval);
249  auto a_list = MatchListExtract(pairList, int(dr1), int(dr2), AstrometryTransformIdentity());
250  a_list->refineTransform(conditions.nSigmas); // mandatory for the sorting fields to be filled
251  Solutions.push_back(std::move(a_list));
252  }
253  } /* end of loop on (r,theta) bins */
254  Solutions.sort(DecreasingQuality);
256  best.swap(*Solutions.begin());
257  /* remove the first one from the list */
258  Solutions.pop_front();
259  if (conditions.printLevel >= 1) {
260  LOGLS_DEBUG(_log, "Best solution " << best->computeResidual() << " npairs " << best->size());
261  LOGLS_DEBUG(_log, *(best->getTransform()));
262  LOGLS_DEBUG(_log, "Chi2 " << best->getChi2() << ',' << " Number of solutions " << Solutions.size());
263  }
264  return best;
265 }
266 
267 /* this matching routine searches brutally a match between lists in
268  the 4 parameter space: size ratio, rotation angle, x and y
269  shifts. This is done by histogramming where combinations of four
270  objets (2 on each list) fall in this 4 parameter space.
271 
272  One trick is that rather than using actual offsets, we histogram
273  object indices of the combination:
274 */
275 
276 static std::unique_ptr<StarMatchList> ListMatchupRotShift_New(BaseStarList &list1, BaseStarList &list2,
277  const AstrometryTransform &transform,
278  const MatchConditions &conditions) {
279  if (list1.size() <= 4 || list2.size() <= 4) {
280  LOGL_FATAL(_log, "ListMatchupRotShift_New : (at least) one of the lists is too short.");
281  return nullptr;
282  }
283 
284  SegmentList sList1(list1, conditions.nStarsList1, transform);
285  SegmentList sList2(list2, conditions.nStarsList2, AstrometryTransformIdentity());
286 
287  /* choose the binning of the histogram so that
288  1: ratio = 1 and rotation angle = n * (pi/2) are bin centers. since
289  the angle is computed using atan2, its range is [-pi,pi],
290  and the histogram range is [-pi-eps, pi-eps], so
291  if (angle>pi- angleOffset) angle -= 2*pi before filling. */
292  int nBinsR = 21;
293  int nBinsAngle = 180; /* can be divided by 4 */
294  double angleOffset = M_PI / nBinsAngle;
295  double minRatio = conditions.minSizeRatio();
296  double maxRatio = conditions.maxSizeRatio();
297  SparseHisto4d histo(nBinsR, minRatio, maxRatio, nBinsAngle, -M_PI - angleOffset, M_PI - angleOffset,
298  conditions.nStarsList1, 0., conditions.nStarsList1, conditions.nStarsList2, 0.,
299  conditions.nStarsList2, sList1.size() * sList2.size());
300 
301  SegmentIterator segi1, segi2;
302  Segment *seg1, *seg2;
303  double ratio, angle;
304 
305  for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
306  seg1 = &(*segi1);
307  if (seg1->r == 0) continue;
308  for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
309  seg2 = &(*segi2);
310  /* if one considers the 2 segments as complex numbers z1 and z2, ratio=mod(z1/z2) and angle =
311  * arg(z1/z2) */
312  /* I did not put a member function in Segment to compute both because we apply a cut on ratio
313  before actually
314  computing the angle (which involves a call to atan2 (expensive)) */
315  ratio = seg2->r / seg1->r;
316  if (ratio > maxRatio) continue;
317  if (ratio < minRatio) break; /* use th fact that segment lists are sorted by decresing length */
318  angle = seg1->relativeAngle(seg2);
319  if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
320  histo.fill(ratio, angle, seg1->s1rank + 0.5, seg2->s1rank + 0.5);
321  }
322  }
323 
324  SolList Solutions;
325  /* now we find the highest bins of the histogram, and recover the original objects.
326  This involves actually re-looping on the combinations, but it is much
327  faster that the original histogram filling loop, since we only compute
328  angle and ratio for Segments that have the right first object
329  */
330 
331  int oldMaxContent = 0;
332 
333  for (int i = 0; i < 4 * conditions.maxTrialCount;
334  ++i) // leave a limit to make avoid (almost) infinite loops
335  {
336  double pars[4];
337  int maxContent = histo.maxBin(pars);
338  if (maxContent == 0) break;
339  if (conditions.printLevel >= 1) {
340  LOGLS_DEBUG(_log, "ValMax " << maxContent << " ratio " << pars[0] << " angle " << pars[1]);
341  }
342  histo.zeroBin(pars);
343  if (i > 0) { /* the match possibilities come out in a random order when they have the same content.
344  so, we stop investigating guesses when the content goes down AND the requested search
345  depth
346  (maxTrialCount) is reached */
347  if (maxContent < oldMaxContent && i >= conditions.maxTrialCount) break;
348  }
349  oldMaxContent = maxContent;
350  /* reloop on segment pairs to select the ones in this specific bin */
351  int rank1L1 = int(pars[2]);
352  int rank1L2 = int(pars[3]);
353  double minAngle, maxAngle;
354  histo.binLimits(pars, 0, minRatio, maxRatio);
355  histo.binLimits(pars, 1, minAngle, maxAngle);
356 
357  std::unique_ptr<StarMatchList> a_list(new StarMatchList);
358 
359  for (segi1 = sList1.begin(); segi1 != sList1.end(); ++segi1) {
360  seg1 = &(*segi1);
361  if (seg1->s1rank != rank1L1) continue;
362  if (seg1->r == 0) continue;
363  for (segi2 = sList2.begin(); segi2 != sList2.end(); ++segi2) {
364  seg2 = &(*segi2);
365  if (seg2->s1rank != rank1L2) continue;
366  // push in the list the match corresponding to end number 1 of segments
367  if (a_list->size() == 0)
368  a_list->push_back(StarMatch(*(seg1->s1), *(seg2->s1), seg1->s1, seg2->s1));
369  ratio = seg2->r / seg1->r;
370  if (ratio > maxRatio) continue;
371  if (ratio < minRatio)
372  break; /* use the fact that segment lists are sorted by decresing length */
373  angle = seg1->relativeAngle(seg2);
374  if (angle > M_PI - angleOffset) angle -= 2. * M_PI;
375  if (angle < minAngle || angle > maxAngle) continue;
376  /* here we have 2 segments which have the right
377  - length ratio
378  - relative angle
379  - first objects (objects on the end number 1).
380  The objects on the end number 2 are the actual matches : */
381  a_list->push_back(StarMatch(*(seg1->s2), *(seg2->s2), seg1->s2, seg2->s2));
382  }
383  }
384 
385  // a basic check for sanity of the algorithm :
386 
387  if (int(a_list->size()) != maxContent + 1) {
388  LOGLS_ERROR(_log, "There is an internal inconsistency in ListMatchupRotShift.");
389  LOGLS_ERROR(_log, "maxContent = " << maxContent);
390  LOGLS_ERROR(_log, "matches->size() = " << a_list->size());
391  }
392  a_list->refineTransform(conditions.nSigmas);
393  Solutions.push_back(std::move(a_list));
394  }
395 
396  if (Solutions.size() == 0) {
397  LOGLS_ERROR(_log, "Error In ListMatchup : not a single pair match.");
398  LOGLS_ERROR(_log, "Probably, the relative scale of lists is not within bounds.");
399  LOGLS_ERROR(_log, "min/max ratios: " << minRatio << ' ' << maxRatio);
400  return nullptr;
401  }
402 
403  Solutions.sort(DecreasingQuality);
405  best.swap(*Solutions.begin());
406  /* remove the first one from the list */
407  Solutions.pop_front();
408  if (conditions.printLevel >= 1) {
409  LOGLS_INFO(_log, "Best solution " << best->computeResidual() << " npairs " << best->size());
410  LOGLS_INFO(_log, *(best->getTransform()));
411  LOGLS_INFO(_log, "Chi2 " << best->getChi2() << ", Number of solutions " << Solutions.size());
412  }
413  return best;
414 }
415 
416 static std::unique_ptr<StarMatchList> ListMatchupRotShift(BaseStarList &list1, BaseStarList &list2,
417  const AstrometryTransform &transform,
418  const MatchConditions &conditions) {
419  if (conditions.algorithm == 1)
420  return ListMatchupRotShift_Old(list1, list2, transform, conditions);
421  else
422  return ListMatchupRotShift_New(list1, list2, transform, conditions);
423 }
424 
426  const MatchConditions &conditions) {
427  list1.fluxSort();
428  list2.fluxSort();
429 
430  return ListMatchupRotShift(list1, list2, AstrometryTransformIdentity(), conditions);
431 }
432 
434  const MatchConditions &conditions) {
435  list1.fluxSort();
436  list2.fluxSort();
437 
438  AstrometryTransformLinear flip(0, 0, 1, 0, 0, -1);
439  std::unique_ptr<StarMatchList> flipped(ListMatchupRotShift(list1, list2, flip, conditions));
441  ListMatchupRotShift(list1, list2, AstrometryTransformIdentity(), conditions));
442  if (!flipped || !unflipped) return std::unique_ptr<StarMatchList>(nullptr);
443  if (conditions.printLevel >= 1) {
444  LOGLS_DEBUG(_log,
445  "unflipped Residual " << unflipped->computeResidual() << " nused " << unflipped->size());
446  LOGLS_DEBUG(_log, "flipped Residual " << flipped->computeResidual() << " nused " << flipped->size());
447  }
448  if (DecreasingQuality(flipped, unflipped)) {
449  if (conditions.printLevel >= 1) LOGL_DEBUG(_log, "Keeping flipped solution.");
450  // One should NOT apply the flip to the result because the matchlist
451  // (even the flipped one) contains the actual coordinates of stars.
452  // MatchListExtract is always called with AstrometryTransformIdentity() as last parameter
453  return flipped;
454  } else {
455  if (conditions.printLevel >= 1) LOGL_DEBUG(_log, "Keeping unflipped solution.");
456  return unflipped;
457  }
458 }
459 
460 #ifdef STORAGE
461 // timing : 2.5 s for l1 of 1862 objects and l2 of 2617 objects
463  const BaseStarList &list2,
464  const AstrometryTransform &transform,
465  double maxShift) {
466  int ncomb = list1.size() * list2.size();
467  if (!ncomb) return nullptr;
468  int nx;
469  if (ncomb > 10000)
470  nx = 100;
471  else
472  nx = (int)sqrt(ncomb);
473 
474  Histo2d histo(nx, -maxShift, maxShift, nx, -maxShift, maxShift);
475 
476  BaseStarCIterator s1, s2;
477  double x1, y1;
478  for (s1 = list1.begin(); s1 != list1.end(); ++s1) {
479  transform.apply((*s1)->x, (*s1)->y, x1, y1);
480  for (s2 = list2.begin(); s2 != list2.end(); ++s2) {
481  histo.fill((*s2)->x - x1, (*s2)->y - y1);
482  }
483  }
484  double dx = 0, dy = 0;
485  histo.maxBin(dx, dy);
486  return std::unique_ptr<AstrometryTransformLinear>(new AstrometryTransformLinearShift(dx, dy));
487 }
488 #endif /*STORAGE*/
489 
490 // timing : 140 ms for l1 of 1862 objects and l2 of 2617 objects (450 MHz, "-O4") maxShift = 200.
492  const BaseStarList &list2,
494  double maxShift, double binSize) {
495  int nx;
496  if (binSize == 0) {
497  int ncomb = list1.size() * list2.size();
498  if (ncomb > 10000)
499  nx = 100;
500  else
501  nx = (int)sqrt(double(ncomb));
502  if (!ncomb) return std::unique_ptr<AstrometryTransformLinear>(nullptr);
503  } else
504  nx = int(2 * maxShift / binSize + 0.5);
505 
506  Histo2d histo(nx, -maxShift, maxShift, nx, -maxShift, maxShift);
507  double binSizeNew = 2 * maxShift / nx;
508 
510  FastFinder finder(list2);
511  double x1, y1;
512  for (s1 = list1.begin(); s1 != list1.end(); ++s1) {
513  transform.apply((*s1)->x, (*s1)->y, x1, y1);
514  FastFinder::Iterator it = finder.beginScan(Point(x1, y1), maxShift);
515  while (*it) {
516  auto s2 = *it;
517  histo.fill(s2->x - x1, s2->y - y1);
518  ++it;
519  }
520  }
521  SolList Solutions;
522  for (int i = 0; i < 4; ++i) {
523  double dx = 0, dy = 0;
524  double count = histo.maxBin(dx, dy);
525  histo.fill(dx, dy, -count); // zero the maxbin
526  AstrometryTransformLinearShift shift(dx, dy);
527  auto newGuess = compose(shift, transform);
528  auto raw_matches = listMatchCollect(list1, list2, newGuess.get(), binSizeNew);
530  raw_matches->applyTransform(*matches, &transform);
531  matches->setTransformOrder(1);
532  matches->refineTransform(3.);
533  Solutions.push_back(std::move(matches));
534  }
535  Solutions.sort(DecreasingQuality);
537  new AstrometryTransformLinear(*std::const_pointer_cast<AstrometryTransformLinear>(
538  std::dynamic_pointer_cast<const AstrometryTransformLinear>(
539  Solutions.front()->getTransform()))));
540  return best;
541 }
542 
543 #ifdef STORAGE
544 
545 // this is the old fashioned way...
546 
547 std::unique_ptr<StarMatchList> listMatchCollect_Slow(const BaseStarList &list1, const BaseStarList &list2,
548  const AstrometryTransform *guess, const double maxDist) {
549  std::unique_ptr<StarMatchList> matches(new StarMatchList);
550  /****** Collect ***********/
551  for (BaseStarCIterator si = list1.begin(); si != list1.end(); ++si) {
552  const Point *p1 = (*si);
553  const Point p2 = guess->apply(*p1);
554  const BaseStar *neighbour = list2.findClosest(p2);
555  if (!neighbour) continue;
556  double distance = p2.Distance(*neighbour);
557  if (distance < maxDist) {
558  matches->push_back(StarMatch(*p1, *neighbour, *si, neighbour));
559  // assign the distance, since we have it in hand:
560  matches->back().distance = distance;
561  }
562  }
563  return matches;
564 }
565 #endif
566 
567 // here is the real active routine:
568 
570  const AstrometryTransform *guess, const double maxDist) {
572  /****** Collect ***********/
573  FastFinder finder(list2);
574  for (BaseStarCIterator si = list1.begin(); si != list1.end(); ++si) {
575  auto p1 = (*si);
576  Point p2 = guess->apply(*p1);
577  auto neighbour = finder.findClosest(p2, maxDist);
578  if (!neighbour) continue;
579  double distance = p2.Distance(*neighbour);
580  if (distance < maxDist) {
581  matches->push_back(StarMatch(*p1, *neighbour, p1, neighbour));
582  // assign the distance, since we have it in hand:
583  matches->back().distance = distance;
584  }
585  }
586  matches->setTransform(guess);
587 
588  return matches;
589 }
590 
591 #ifdef STORAGE
592 // unused
594 std::unique_ptr<StarMatchList> CollectAndFit(const BaseStarList &list1, const BaseStarList &list2,
595  const AstrometryTransform *guess, const double maxDist) {
596  const AstrometryTransform *bestTransform = guess;
598  while (true) {
599  auto m = listMatchCollect(list1, list2, bestTransform, maxDist);
600  m->setTransform(bestTransform);
601  m->refineTransform(3.);
602  LOGLS_INFO(_log, "Iterating: resid " << m->computeResidual() << " size " << m->size());
603  if (!prevMatch ||
604  (prevMatch && m->computeResidual() < prevMatch->computeResidual() * 0.999 && m->Chi2() > 0)) {
605  prevMatch.swap(m);
606  bestTransform = prevMatch->Transform();
607  } else {
608  break;
609  }
610  }
611  return prevMatch;
612 }
613 #endif
614 
616  const double maxDist) {
618  FastFinder finder(list2);
619  for (BaseStarCIterator si = list1.begin(); si != list1.end(); ++si) {
620  auto p1 = (*si);
621  auto neighbour = finder.findClosest(*p1, maxDist);
622  if (!neighbour) continue;
623  double distance = p1->Distance(*neighbour);
624  if (distance < maxDist) {
625  matches->push_back(StarMatch(*p1, *neighbour, p1, neighbour));
626  // assign the distance, since we have it in hand:
627  matches->back().distance = distance;
628  }
629  }
630 
631  matches->setTransform(std::make_shared<AstrometryTransformIdentity>());
632 
633  return matches;
634 }
635 
636 static bool is_transform_ok(const StarMatchList *match, double pixSizeRatio2, const size_t nmin) {
637  if ((fabs(fabs(std::dynamic_pointer_cast<const AstrometryTransformLinear>(match->getTransform())
638  ->determinant()) -
639  pixSizeRatio2) /
640  pixSizeRatio2 <
641  0.2) &&
642  (match->size() > nmin))
643  return true;
644  LOGL_ERROR(_log, "transform is not ok!");
645  match->printTransform();
646  return false;
647 }
648 
649 // utility to check current transform difference
650 static double transform_diff(const BaseStarList &List, const AstrometryTransform *T1,
651  const AstrometryTransform *T2) {
652  double diff2 = 0;
653  FatPoint tf1;
654  Point tf2;
655  int count = 0;
656  for (BaseStarCIterator it = List.begin(); it != List.end(); ++it) {
657  const BaseStar &s = **it;
658  T1->transformPosAndErrors(s, tf1);
659  T2->apply(s, tf2);
660  double dx = tf1.x - tf2.x;
661  double dy = tf1.y - tf2.y;
662  diff2 += (tf1.vy * dx * dx + tf1.vx * dy * dy - 2 * tf1.vxy * dx * dy) /
663  (tf1.vx * tf1.vy - tf1.vxy * tf1.vxy);
664  count++;
665  }
666  if (count) return diff2 / double(count);
667  return 0;
668 }
669 
670 static double median_distance(const StarMatchList *match, const AstrometryTransform *transform) {
671  size_t nstars = match->size();
672  std::vector<double> resid(nstars);
674  for (auto it = match->begin(); it != match->end(); ++it, ++ir)
675  *ir = sqrt(transform->apply(it->point1).computeDist2(it->point2));
676  sort(resid.begin(), resid.end());
677  return (nstars & 1) ? resid[nstars / 2] : (resid[nstars / 2 - 1] + resid[nstars / 2]) * 0.5;
678 }
679 
681  const BaseStarList &List2,
682  const MatchConditions &conditions) {
683  BaseStarList list1, list2;
684  List1.copyTo(list1);
685  list1.fluxSort();
686  List2.copyTo(list2);
687  list2.fluxSort();
688 
689  LOGLS_INFO(_log, "listMatchCombinatorial: find match between " << list1.size() << " and " << list2.size()
690  << " stars...");
691  auto match = matchSearchRotShiftFlip(list1, list2, conditions);
692  double pixSizeRatio2 = std::pow(conditions.sizeRatio, 2);
693  size_t nmin =
694  std::min(size_t(10), size_t(std::min(List1.size(), List2.size()) * conditions.minMatchRatio));
695 
697  if (is_transform_ok(match.get(), pixSizeRatio2, nmin))
698  transform = match->getTransform()->clone();
699  else {
700  LOGL_ERROR(_log, "listMatchCombinatorial: direct transform failed, trying reverse");
701  match = matchSearchRotShiftFlip(list2, list1, conditions);
702  if (is_transform_ok(match.get(), pixSizeRatio2, nmin))
703  transform = match->inverseTransform();
704  else {
705  LOGL_FATAL(_log, "FAILED");
706  }
707  }
708 
709  if (transform) {
710  LOGL_INFO(_log, "FOUND");
711  if (conditions.printLevel >= 1) {
712  LOGL_DEBUG(_log, " listMatchCombinatorial: found the following transform.");
713  LOGLS_DEBUG(_log, *transform);
714  }
715  } else
716  LOGL_ERROR(_log, "listMatchCombinatorial: failed to find a transform");
717  return transform;
718 }
719 
722  const int maxOrder) {
723  if (!transform) {
724  return std::unique_ptr<AstrometryTransform>(nullptr);
725  }
726 
727  // some hard-coded constants that could go in a param file
728  const double brightDist = 2.; // distance in pixels in a match
729  const double fullDist = 4.; // distance in pixels in a match between entire lists
730  const double nSigmas = 3.; // k-sigma clipping on residuals
731  const size_t nStars = 500; // max number of bright stars to fit
732 
733  int order = 1;
734  size_t nstarmin = 3;
735 
736  BaseStarList list1, list2;
737  List1.copyTo(list1);
738  list1.fluxSort();
739  list1.cutTail(nStars);
740  List2.copyTo(list2);
741  list2.fluxSort();
742  list2.cutTail(nStars);
743 
744  auto fullMatch = listMatchCollect(List1, List2, transform.get(), fullDist);
745  auto brightMatch = listMatchCollect(list1, list2, transform.get(), brightDist);
746  double curChi2 = computeChi2(*brightMatch, *transform) / brightMatch->size();
747 
748  LOGLS_INFO(_log, "listMatchRefine: start: med.resid " << median_distance(fullMatch.get(), transform.get())
749  << " #match " << fullMatch->size());
750 
751  do { // loop on transform order on full list of stars
752  auto curTransform = brightMatch->getTransform()->clone();
753  unsigned iter = 0;
754  double transDiff;
755  do { // loop on transform diff only on bright stars
756  brightMatch->setTransformOrder(order);
757  brightMatch->refineTransform(nSigmas);
758  transDiff = transform_diff(list1, brightMatch->getTransform().get(), curTransform.get());
759  curTransform = brightMatch->getTransform()->clone();
760  brightMatch = listMatchCollect(list1, list2, curTransform.get(), brightDist);
761  } while (brightMatch->size() > nstarmin && transDiff > 0.05 && ++iter < 5);
762 
763  double prevChi2 = curChi2;
764  curChi2 = computeChi2(*brightMatch, *curTransform) / brightMatch->size();
765 
766  fullMatch = listMatchCollect(List1, List2, curTransform.get(), fullDist);
767  LOGLS_INFO(_log, "listMatchRefine: order " << order << " med.resid "
768  << median_distance(fullMatch.get(), curTransform.get())
769  << " #match " << fullMatch->size());
770  if (((prevChi2 - curChi2) > 0.01 * curChi2) && curChi2 > 0) {
771  LOGLS_INFO(_log, " listMatchRefine: order " << order << " was a better guess.");
772  transform = brightMatch->getTransform()->clone();
773  }
774  nstarmin = brightMatch->getTransform()->getNpar();
775  } while (++order <= maxOrder);
776 
777  return transform;
778 }
779 } // namespace jointcal
780 } // 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(const double xIn, const 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, const 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
Definition: FastFinder.cc:174
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, const int nStar, const AstrometryTransform &transform=AstrometryTransformIdentity())
Definition: ListMatch.cc:96
std::lists of Stars.
Definition: StarList.h:58
void cutTail(const 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:54
T count(T... args)
T end(T... args)
T endl(T... args)
T front(T... args)
T min(T... args)
T move(T... args)
StarList< BaseStar > BaseStarList
Definition: BaseStar.h:121
std::list< std::unique_ptr< StarMatchList > > SolList
Definition: ListMatch.cc:157
std::unique_ptr< StarMatchList > listMatchCollect(const BaseStarList &list1, const BaseStarList &list2, const AstrometryTransform *guess, const double maxDist)
assembles star matches.
Definition: ListMatch.cc:569
std::unique_ptr< AstrometryTransform > listMatchCombinatorial(const BaseStarList &list1, const BaseStarList &list2, const MatchConditions &conditions=MatchConditions())
Definition: ListMatch.cc:680
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
std::list< Segment >::iterator SegmentIterator
Definition: ListMatch.cc:91
std::list< Segment >::const_iterator SegmentCIterator
Definition: ListMatch.cc:92
SegmentPairList::const_iterator SegmentPairListCIterator
Definition: ListMatch.cc:121
SegmentPairList::iterator SegmentPairListIterator
Definition: ListMatch.cc:120
BaseStarList::const_iterator BaseStarCIterator
Definition: BaseStar.h:123
std::unique_ptr< AstrometryTransform > listMatchRefine(const BaseStarList &list1, const BaseStarList &list2, std::unique_ptr< AstrometryTransform > transform, const int maxOrder=3)
Definition: ListMatch.cc:720
std::list< SegmentPair > SegmentPairList
Definition: ListMatch.cc:119
double computeChi2(const StarMatchList &L, const AstrometryTransform &transform)
the actual chi2
Definition: StarMatch.cc:246
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:425
std::unique_ptr< StarMatchList > matchSearchRotShiftFlip(BaseStarList &list1, BaseStarList &list2, const MatchConditions &conditions)
same as above but searches also a flipped solution.
Definition: ListMatch.cc:433
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:491
A base class for image defects.
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
friend std::ostream & operator<<(std::ostream &stream, const Segment &segment)
Definition: ListMatch.cc:78
Segment(std::shared_ptr< const BaseStar > star1, std::shared_ptr< const BaseStar > star2, const int star1Rank, const AstrometryTransform &transform)
Definition: ListMatch.cc:61
std::shared_ptr< const BaseStar > s2
Definition: ListMatch.cc:57
double relativeAngle(Segment *other)
Definition: ListMatch.cc:74
SegmentPair(Segment *f, Segment *s)
Definition: ListMatch.cc:116
T swap(T... args)
table::Key< int > order