LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
FastFinder.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 <algorithm>
26 
27 #include "lsst/log/Log.h"
28 #include "lsst/jointcal/BaseStar.h"
30 
31 namespace {
32 LOG_LOGGER _log = LOG_GET("jointcal.FastFinder");
33 }
34 
35 namespace lsst {
36 namespace jointcal {
37 
38 FastFinder::FastFinder(const BaseStarList &list, const unsigned nXSlice)
39  : baselist(list), count(list.size()), stars(count), nslice(nXSlice), index(nslice + 1) {
40  if (count == 0) return;
41 
42  // fill "stars"
43  unsigned j = 0;
44  for (auto const &ci : list) {
45  stars[j] = ci;
46  ++j;
47  }
48 
49  sort(stars.begin(), stars.end(),
50  [](const stars_element &E1, const stars_element &E2) { return (E1->x < E2->x); });
51 
52  xmin = stars[0]->x;
53  xmax = stars[count - 1]->x;
55  if (xmin == xmax) nslice = 1;
56 
57  // the x size of each slice:
58  xstep = (xmax - xmin) / nslice;
59 
60  // fill the index array with the first star beyond the slice limit.
61  index[0] = 0; // first
62  unsigned istar = 0;
63  for (unsigned islice = 1; islice < nslice; ++islice) {
64  double xend = xmin + (islice)*xstep;
65  while (istar < count && stars[istar]->x < xend) ++istar;
66  index[islice] = istar;
67  }
68  index[nslice] = count; // last
69  for (unsigned islice = 0; islice < nslice; ++islice) {
70  sort(stars.begin() + index[islice], stars.begin() + index[islice + 1],
71  [](const stars_element &E1, const stars_element &E2) {
72  return (E1->y < E2->y);
73  }); // sort each slice in y.
74  }
75  // dump();
76 }
77 
78 void FastFinder::dump() const {
79  for (unsigned i = 0; i < count; ++i) {
80  stars[i]->dump();
81  }
82 }
83 
85  bool (*SkipIt)(const BaseStar &)) const {
86  if (count == 0) return nullptr;
87  FastFinder::Iterator it = beginScan(where, maxDist);
88  if (*it == nullptr) return nullptr;
90  double minDist2 = maxDist * maxDist;
91  for (; *it != nullptr; ++it) {
92  if (SkipIt && SkipIt(**it)) continue;
93  double dist2 = where.computeDist2(**it);
94  if (dist2 < minDist2) {
95  pbest = *it;
96  minDist2 = dist2;
97  }
98  }
99  return pbest;
100 }
101 
104  bool (*SkipIt)(const BaseStar &)) const {
105  closest = nullptr;
106  if (count == 0) return nullptr;
107  FastFinder::Iterator it = beginScan(where, maxDist);
108  if (*it == nullptr) return nullptr;
109  std::shared_ptr<const BaseStar> pbest1; // closest
110  std::shared_ptr<const BaseStar> pbest2; // second closest
111  double minDist1_2 = maxDist * maxDist;
112  double minDist2_2 = maxDist * maxDist;
113  for (; *it != nullptr; ++it) {
114  if (SkipIt && SkipIt(**it)) continue;
115  double dist2 = where.computeDist2(**it);
116  if (dist2 < minDist1_2) {
117  pbest2 = pbest1;
118  minDist2_2 = minDist1_2;
119  pbest1 = *it;
120  minDist1_2 = dist2;
121  } else if (dist2 < minDist2_2) {
122  pbest2 = *it;
123  minDist2_2 = dist2;
124  }
125  }
126  closest = pbest1;
127  return pbest2;
128 }
129 
130 /* It is by no means clear the the 2 following routines are actually needed.
131  It is nor clear to me (P.A) why they are different... but they really are.
132 */
133 /* Locate the last position (in the sorted array) between begin and
134  end that lies before yVal.*/
136  if (begin == stars.end() || begin == end) return stars.end();
137  int span = end - begin - 1;
138  while (span > 1) {
139  int half_span = span / 2;
140  pstar middle = begin + half_span;
141  if ((*middle)->y < yVal) {
142  begin += half_span;
143  span -= half_span;
144  } else {
145  span -= (span - half_span);
146  }
147  }
148  return begin;
149 }
150 
151 /* Locate the first position (in the sorted array) between begin and
152  end that lies beyond yVal.*/
154  if (begin == stars.end()) return stars.end();
155  int span = end - begin - 1;
156  while (span > 1) {
157  int half_span = span / 2;
158  pstar middle = end - half_span;
159  if ((*middle)->y > yVal) {
160  end -= half_span;
161  span -= half_span;
162  } else {
163  span -= (span - half_span);
164  }
165  }
166  return end - 1;
167 }
168 
169 void FastFinder::findRangeInSlice(const int iSlice, const double yStart, const double yEnd, pstar &start,
170  pstar &end) const {
171  start = locateYStart(stars.begin() + index[iSlice], stars.begin() + index[iSlice + 1], yStart);
172  end = locateYEnd(start, stars.begin() + index[iSlice + 1], yEnd);
173 }
174 
175 FastFinder::Iterator FastFinder::beginScan(const Point &where, double maxDist) const {
176  return FastFinder::Iterator(*this, where, maxDist);
177 }
178 
180 
181 Iterator::Iterator(const FastFinder &F, const Point &where, double maxDist)
182  : finder(F), null_value(F.stars.end()) {
183  current = pend = null_value; // does not iterate
184  int startSlice = 0;
185  if (finder.xstep != 0) // means we have several slices
186  {
187  startSlice = std::max(0, int((where.x - maxDist - finder.xmin) / finder.xstep));
188  /* obviously, endSlice (and starSlice) can be negative.
189  This is why slice indices are "int" rather than "unsigned". */
190  endSlice = std::min(int(finder.nslice), int((where.x + maxDist - finder.xmin) / finder.xstep) + 1);
191  } else {
192  startSlice = 0;
193  endSlice = 1;
194  }
195  // beyond limits:
196  if (startSlice >= int(finder.nslice) || endSlice < 0) return;
197  // we are inside in x, so, we setup the y range:
198  yStart = where.y - maxDist;
199  yEnd = where.y + maxDist;
200  /* rather than initializing here, we step back one
201  slice and let "++" do its job */
202  currentSlice = startSlice - 1; // again, this requires "int" slices
203  ++(*this);
204 }
205 
207  if (current != null_value) return *current;
208  return nullptr;
209 }
210 
212  if (current != pend) {
213  current++;
214  } else
215  do {
216  currentSlice++;
217  if (currentSlice >= endSlice) {
219  return;
220  }
222  } while (current == null_value);
223  check();
224 }
225 
227  if (current != null_value &&
228  (current < finder.stars.begin() || current >= finder.stars.begin() + finder.count)) {
229  LOGLS_ERROR(_log, "Error in FastFinder " << *current << " " << *(finder.stars.begin()) << ' '
230  << *(finder.stars.begin() + finder.count));
231  }
232 }
233 } // namespace jointcal
234 } // namespace lsst
#define LOG_LOGGER
Definition: Log.h:703
A point in a plane.
Definition: Point.h:36
decltype(stars) typedef ::const_iterator pstar
Definition: FastFinder.h:70
std::vector< std::shared_ptr< const BaseStar > > stars
Definition: FastFinder.h:64
Iterator(const FastFinder &f, const Point &where, double maxDist)
Definition: FastFinder.cc:181
void dump() const
mostly for debugging
Definition: FastFinder.cc:78
pstar locateYEnd(pstar begin, pstar end, double yVal) const
Definition: FastFinder.cc:153
LSST DM logging module built on log4cxx.
double x
coordinate
Definition: Point.h:41
T min(T... args)
The base class for handling stars. Used by all matching routines.
Definition: BaseStar.h:50
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:84
FastFinder::Iterator Iterator
Definition: FastFinder.cc:179
std::lists of Stars.
Definition: StarList.h:58
A base class for image defects.
Iterator beginScan(const Point &where, double maxDist) const
Definition: FastFinder.cc:175
Iterator meant to traverse objects within some limiting distance.
Definition: FastFinder.h:90
FastFinder(const BaseStarList &list, const unsigned nXSlice=100)
Constructor.
Definition: FastFinder.cc:38
pstar locateYStart(pstar begin, pstar end, double yVal) const
Definition: FastFinder.cc:135
STL class.
T max(T... args)
double x
T count(T... args)
void findRangeInSlice(const int iSlice, const double yStart, const double yEnd, pstar &start, pstar &end) const
Definition: FastFinder.cc:169
This is an auxillary class for matching objects from starlists.
Definition: FastFinder.h:54
std::vector< unsigned > index
Definition: FastFinder.h:66
decltype(stars) typedef ::value_type stars_element
Definition: FastFinder.h:69
std::shared_ptr< const BaseStar > secondClosest(const Point &where, const double maxDist, std::shared_ptr< const BaseStar > &closest, bool(*SkipIt)(const BaseStar &)=nullptr) const
Definition: FastFinder.cc:102
double computeDist2(const Point &other) const
distance squared to other
Definition: Point.h:55
T sort(T... args)
stars_element operator*() const
Definition: FastFinder.cc:206
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
Fast locator in starlists.
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
Definition: Log.h:668
int end