LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
30
31namespace {
32LOG_LOGGER _log = LOG_GET("lsst.jointcal.FastFinder");
33}
34
35namespace lsst {
36namespace jointcal {
37
38FastFinder::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
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.*/
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 auto 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.*/
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 auto 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
168void 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
174FastFinder::Iterator FastFinder::beginScan(const Point &where, double maxDist) const {
175 return FastFinder::Iterator(*this, where, maxDist);
176}
177
179
180Iterator::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
int end
Fast locator in starlists.
LSST DM logging module built on log4cxx.
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition Log.h:75
#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
T begin(T... args)
The base class for handling stars. Used by all matching routines.
Definition BaseStar.h:51
Iterator meant to traverse objects within some limiting distance.
Definition FastFinder.h:90
Iterator(const FastFinder &f, const Point &where, double maxDist)
This is an auxillary class for matching objects from starlists.
Definition FastFinder.h:54
decltype(stars)::const_iterator pstar
Definition FastFinder.h:70
decltype(stars)::value_type stars_element
Definition FastFinder.h:69
pstar locateYStart(pstar begin, pstar end, double yVal) const
void findRangeInSlice(int iSlice, double yStart, double yEnd, pstar &start, pstar &end) const
void print(std::ostream &out) const
mostly for debugging
Definition FastFinder.cc:77
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
FastFinder(const BaseStarList &list, unsigned nXSlice=100)
Constructor.
Definition FastFinder.cc:38
std::shared_ptr< const BaseStar > secondClosest(const Point &where, double maxDist, std::shared_ptr< const BaseStar > &closest, bool(*SkipIt)(const BaseStar &)=nullptr) const
Iterator beginScan(const Point &where, double maxDist) const
std::vector< unsigned > index
Definition FastFinder.h:66
std::vector< std::shared_ptr< const BaseStar > > stars
Definition FastFinder.h:64
pstar locateYEnd(pstar begin, pstar end, double yVal) const
A point in a plane.
Definition Point.h:37
std::lists of Stars.
Definition StarList.h:58
T count(T... args)
daf::base::PropertyList * list
Definition fits.cc:932
T max(T... args)
T min(T... args)
T sort(T... args)