LSSTApplications  19.0.0-14-gb0260a2+72efe9b372,20.0.0+7927753e06,20.0.0+8829bf0056,20.0.0+995114c5d2,20.0.0+b6f4b2abd1,20.0.0+bddc4f4cbe,20.0.0-1-g253301a+8829bf0056,20.0.0-1-g2b7511a+0d71a2d77f,20.0.0-1-g5b95a8c+7461dd0434,20.0.0-12-g321c96ea+23efe4bbff,20.0.0-16-gfab17e72e+fdf35455f6,20.0.0-2-g0070d88+ba3ffc8f0b,20.0.0-2-g4dae9ad+ee58a624b3,20.0.0-2-g61b8584+5d3db074ba,20.0.0-2-gb780d76+d529cf1a41,20.0.0-2-ged6426c+226a441f5f,20.0.0-2-gf072044+8829bf0056,20.0.0-2-gf1f7952+ee58a624b3,20.0.0-20-geae50cf+e37fec0aee,20.0.0-25-g3dcad98+544a109665,20.0.0-25-g5eafb0f+ee58a624b3,20.0.0-27-g64178ef+f1f297b00a,20.0.0-3-g4cc78c6+e0676b0dc8,20.0.0-3-g8f21e14+4fd2c12c9a,20.0.0-3-gbd60e8c+187b78b4b8,20.0.0-3-gbecbe05+48431fa087,20.0.0-38-ge4adf513+a12e1f8e37,20.0.0-4-g97dc21a+544a109665,20.0.0-4-gb4befbc+087873070b,20.0.0-4-gf910f65+5d3db074ba,20.0.0-5-gdfe0fee+199202a608,20.0.0-5-gfbfe500+d529cf1a41,20.0.0-6-g64f541c+d529cf1a41,20.0.0-6-g9a5b7a1+a1cd37312e,20.0.0-68-ga3f3dda+5fca18c6a4,20.0.0-9-g4aef684+e18322736b,w.2020.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 }
76 
77 void FastFinder::print(std::ostream &out) const {
78  for (unsigned i = 0; i < count; ++i) {
79  stars[i]->print(out);
80  }
81 }
82 
84  bool (*SkipIt)(const BaseStar &)) const {
85  if (count == 0) return nullptr;
86  FastFinder::Iterator it = beginScan(where, maxDist);
87  if (*it == nullptr) return nullptr;
89  double minDist2 = maxDist * maxDist;
90  for (; *it != nullptr; ++it) {
91  if (SkipIt && SkipIt(**it)) continue;
92  double dist2 = where.computeDist2(**it);
93  if (dist2 < minDist2) {
94  pbest = *it;
95  minDist2 = dist2;
96  }
97  }
98  return pbest;
99 }
100 
103  bool (*SkipIt)(const BaseStar &)) const {
104  closest = nullptr;
105  if (count == 0) return nullptr;
106  FastFinder::Iterator it = beginScan(where, maxDist);
107  if (*it == nullptr) return nullptr;
108  std::shared_ptr<const BaseStar> pbest1; // closest
109  std::shared_ptr<const BaseStar> pbest2; // second closest
110  double minDist1_2 = maxDist * maxDist;
111  double minDist2_2 = maxDist * maxDist;
112  for (; *it != nullptr; ++it) {
113  if (SkipIt && SkipIt(**it)) continue;
114  double dist2 = where.computeDist2(**it);
115  if (dist2 < minDist1_2) {
116  pbest2 = pbest1;
117  minDist2_2 = minDist1_2;
118  pbest1 = *it;
119  minDist1_2 = dist2;
120  } else if (dist2 < minDist2_2) {
121  pbest2 = *it;
122  minDist2_2 = dist2;
123  }
124  }
125  closest = pbest1;
126  return pbest2;
127 }
128 
129 /* It is by no means clear the the 2 following routines are actually needed.
130  It is nor clear to me (P.A) why they are different... but they really are.
131 */
132 /* Locate the last position (in the sorted array) between begin and
133  end that lies before yVal.*/
134 FastFinder::pstar FastFinder::locateYStart(pstar begin, pstar end, double yVal) const {
135  if (begin == stars.end() || begin == end) return stars.end();
136  int span = end - begin - 1;
137  while (span > 1) {
138  int half_span = span / 2;
139  pstar middle = begin + half_span;
140  if ((*middle)->y < yVal) {
141  begin += half_span;
142  span -= half_span;
143  } else {
144  span -= (span - half_span);
145  }
146  }
147  return begin;
148 }
149 
150 /* Locate the first position (in the sorted array) between begin and
151  end that lies beyond yVal.*/
152 FastFinder::pstar FastFinder::locateYEnd(pstar begin, pstar end, double yVal) const {
153  if (begin == stars.end()) return stars.end();
154  int span = end - begin - 1;
155  while (span > 1) {
156  int half_span = span / 2;
157  pstar middle = end - half_span;
158  if ((*middle)->y > yVal) {
159  end -= half_span;
160  span -= half_span;
161  } else {
162  span -= (span - half_span);
163  }
164  }
165  return end - 1;
166 }
167 
168 void FastFinder::findRangeInSlice(const int iSlice, const double yStart, const double yEnd, pstar &start,
169  pstar &end) const {
170  start = locateYStart(stars.begin() + index[iSlice], stars.begin() + index[iSlice + 1], yStart);
171  end = locateYEnd(start, stars.begin() + index[iSlice + 1], yEnd);
172 }
173 
174 FastFinder::Iterator FastFinder::beginScan(const Point &where, double maxDist) const {
175  return FastFinder::Iterator(*this, where, maxDist);
176 }
177 
179 
180 Iterator::Iterator(const FastFinder &F, const Point &where, double maxDist)
181  : finder(F), null_value(F.stars.end()) {
182  current = pend = null_value; // does not iterate
183  int startSlice = 0;
184  if (finder.xstep != 0) // means we have several slices
185  {
186  startSlice = std::max(0, int((where.x - maxDist - finder.xmin) / finder.xstep));
187  /* obviously, endSlice (and starSlice) can be negative.
188  This is why slice indices are "int" rather than "unsigned". */
189  endSlice = std::min(int(finder.nslice), int((where.x + maxDist - finder.xmin) / finder.xstep) + 1);
190  } else {
191  startSlice = 0;
192  endSlice = 1;
193  }
194  // beyond limits:
195  if (startSlice >= int(finder.nslice) || endSlice < 0) return;
196  // we are inside in x, so, we setup the y range:
197  yStart = where.y - maxDist;
198  yEnd = where.y + maxDist;
199  /* rather than initializing here, we step back one
200  slice and let "++" do its job */
201  currentSlice = startSlice - 1; // again, this requires "int" slices
202  ++(*this);
203 }
204 
206  if (current != null_value) return *current;
207  return nullptr;
208 }
209 
211  if (current != pend) {
212  current++;
213  } else
214  do {
215  currentSlice++;
216  if (currentSlice >= endSlice) {
218  return;
219  }
221  } while (current == null_value);
222  check();
223 }
224 
226  if (current != null_value &&
227  (current < finder.stars.begin() || current >= finder.stars.begin() + finder.count)) {
228  LOGLS_ERROR(_log, "Error in FastFinder " << *current << " " << *(finder.stars.begin()) << ' '
229  << *(finder.stars.begin() + finder.count));
230  }
231 }
232 } // namespace jointcal
233 } // namespace lsst
lsst::jointcal::FastFinder::locateYEnd
pstar locateYEnd(pstar begin, pstar end, double yVal) const
Definition: FastFinder.cc:152
lsst::jointcal::FastFinder::Iterator::check
void check() const
Definition: FastFinder.cc:225
LOG_LOGGER
#define LOG_LOGGER
Definition: Log.h:703
lsst::jointcal::FastFinder::count
unsigned count
Definition: FastFinder.h:58
lsst::jointcal::FastFinder::findRangeInSlice
void findRangeInSlice(const int iSlice, const double yStart, const double yEnd, pstar &start, pstar &end) const
Definition: FastFinder.cc:168
lsst::jointcal::FastFinder::beginScan
Iterator beginScan(const Point &where, double maxDist) const
Definition: FastFinder.cc:174
lsst::jointcal::FastFinder::print
void print(std::ostream &out) const
mostly for debugging
Definition: FastFinder.cc:77
std::shared_ptr
STL class.
std::list
STL class.
lsst::jointcal::FastFinder::xstep
double xstep
Definition: FastFinder.h:67
lsst::jointcal::FastFinder::secondClosest
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:101
lsst::jointcal::FastFinder::index
std::vector< unsigned > index
Definition: FastFinder.h:66
lsst::jointcal::FastFinder::Iterator
Iterator meant to traverse objects within some limiting distance.
Definition: FastFinder.h:90
LOG_GET
#define LOG_GET(logger)
Definition: Log.h:75
lsst::jointcal::FastFinder::Iterator::current
pstar current
Definition: FastFinder.h:97
LOGLS_ERROR
#define LOGLS_ERROR(logger, message)
Definition: Log.h:668
std::sort
T sort(T... args)
lsst::jointcal::FastFinder::Iterator::yStart
double yStart
Definition: FastFinder.h:94
lsst::jointcal::FastFinder::stars
std::vector< std::shared_ptr< const BaseStar > > stars
Definition: FastFinder.h:64
end
int end
Definition: BoundedField.cc:105
lsst::jointcal::Point::y
double y
Definition: Point.h:41
lsst::jointcal::FastFinder::stars_element
decltype(stars) typedef ::value_type stars_element
Definition: FastFinder.h:69
lsst::jointcal::FastFinder::Iterator::operator++
void operator++()
Definition: FastFinder.cc:210
lsst::jointcal::Point::computeDist2
double computeDist2(const Point &other) const
distance squared to other
Definition: Point.h:55
std::ostream
STL class.
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::jointcal::FastFinder::Iterator::null_value
pstar null_value
Definition: FastFinder.h:98
lsst::jointcal::StarList
std::lists of Stars.
Definition: StarList.h:58
lsst::jointcal::FastFinder::Iterator::endSlice
int endSlice
Definition: FastFinder.h:93
lsst::jointcal::Iterator
FastFinder::Iterator Iterator
Definition: FastFinder.cc:178
lsst::jointcal::FastFinder::Iterator::yEnd
double yEnd
Definition: FastFinder.h:94
lsst::jointcal::FastFinder::Iterator::Iterator
Iterator(const FastFinder &f, const Point &where, double maxDist)
Definition: FastFinder.cc:180
lsst::jointcal::FastFinder::Iterator::operator*
stars_element operator*() const
Definition: FastFinder.cc:205
lsst::jointcal::FastFinder::findClosest
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
lsst::jointcal::FastFinder::pstar
decltype(stars) typedef ::const_iterator pstar
Definition: FastFinder.h:70
lsst::jointcal
Definition: Associations.h:49
lsst::jointcal::FastFinder::xmax
double xmax
Definition: FastFinder.h:67
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::jointcal::FastFinder::locateYStart
pstar locateYStart(pstar begin, pstar end, double yVal) const
Definition: FastFinder.cc:134
std::min
T min(T... args)
lsst::jointcal::FastFinder
This is an auxillary class for matching objects from starlists.
Definition: FastFinder.h:54
BaseStar.h
std::begin
T begin(T... args)
std::count
T count(T... args)
lsst::jointcal::Point
A point in a plane.
Definition: Point.h:36
FastFinder.h
Fast locator in starlists.
lsst::jointcal::FastFinder::Iterator::finder
const FastFinder & finder
Definition: FastFinder.h:92
lsst::jointcal::FastFinder::xmin
double xmin
Definition: FastFinder.h:67
lsst::jointcal::BaseStar
The base class for handling stars. Used by all matching routines.
Definition: BaseStar.h:50
lsst::jointcal::FastFinder::Iterator::currentSlice
int currentSlice
Definition: FastFinder.h:93
lsst::jointcal::Point::x
double x
coordinate
Definition: Point.h:41
std::max
T max(T... args)
lsst::jointcal::FastFinder::nslice
unsigned nslice
Definition: FastFinder.h:65
Log.h
LSST DM logging module built on log4cxx.
lsst::jointcal::FastFinder::Iterator::pend
pstar pend
Definition: FastFinder.h:97
lsst::jointcal::FastFinder::FastFinder
FastFinder(const BaseStarList &list, const unsigned nXSlice=100)
Constructor.
Definition: FastFinder.cc:38