LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
StarMatch.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 <iomanip>
27
31#include "algorithm" // for copy
32
33/* TO DO:
34 think about imposing a maximum number of matches that may
35 be discarded in Cleanup.
36*/
37
38namespace lsst {
39namespace jointcal {
40
41static double sq(double x) { return x * x; }
42
44 FatPoint tr;
45 transform.transformPosAndErrors(point1, tr);
46 double vxx = tr.vx + point2.vx;
47 double vyy = tr.vy + point2.vy;
48 double vxy = tr.vxy + point2.vxy;
49 double det = vxx * vyy - vxy * vxy;
50 return (vyy * sq(tr.x - point2.x) + vxx * sq(tr.y - point2.y) -
51 2 * vxy * (tr.x - point2.x) * (tr.y - point2.y)) /
52 det;
53}
54
56 stream << match.point1.x << ' ' << match.point1.y << ' ' << match.point2.x << ' ' << match.point2.y << ' '
57 << match.distance << std::endl;
58 return stream;
59}
60
61std::ostream &operator<<(std::ostream &stream, const StarMatchList &starMatchList) {
62 stream << " number of elements " << starMatchList.size() << std::endl;
63 copy(starMatchList.begin(), starMatchList.end(), std::ostream_iterator<StarMatch>(stream));
64 return stream;
65}
66
67static std::vector<double> chi2_array(const StarMatchList &starMatchList,
68 const AstrometryTransform &transform) {
69 std::vector<double> result(starMatchList.size());
70 unsigned count = 0;
71 for (auto const &it : starMatchList) result[count++] = it.computeChi2(transform);
72 return result;
73}
74
75static unsigned chi2_cleanup(StarMatchList &starMatchList, const double chi2Cut,
76 const AstrometryTransform &transform) {
77 unsigned erased = starMatchList.removeAmbiguities(transform);
78 for (auto smi = starMatchList.begin(); smi != starMatchList.end();) {
79 if (smi->chi2 > chi2Cut) {
80 smi = starMatchList.erase(smi);
81 erased++;
82 } else
83 ++smi;
84 }
85 return erased;
86}
87
92void StarMatchList::refineTransform(double nSigmas) {
93 double cut;
94 unsigned nremoved;
95 if (!_transform) _transform.reset(new AstrometryTransformLinear);
96 do {
97 int nused = size();
98 if (nused <= 2) {
99 _chi2 = -1;
100 break;
101 }
102 _chi2 = _transform->fit(*this);
103 /* convention of the fitted routines :
104 - chi2 = 0 means zero degrees of freedom
105 (this was not enforced in AstrometryTransform{Lin,Quad,Cub} ...)
106 - chi2 = -1 means ndof <0 (and hence no possible fit)
107 --> in either case, refinement is over
108 The fact that chi2 = 0 was not enforced when necessary means
109 that in this (rare) case, we were discarding matches at random....
110 With AstrometryTransformPolynomial::fit, this is no longer the case.
111 */
112 if (_chi2 <= 0) return;
113 unsigned npair = int(size());
114 if (npair == 0) break; // should never happen
115
116 // compute some chi2 statistics
117 std::vector<double> chi2_array(npair);
118 unsigned count = 0;
119 for (auto &starMatch : *this)
120 chi2_array[count++] = starMatch.chi2 = starMatch.computeChi2(*_transform);
121
122 std::sort(chi2_array.begin(), chi2_array.end());
123 double median = (npair & 1) ? chi2_array[npair / 2]
124 : (chi2_array[npair / 2 - 1] + chi2_array[npair / 2]) * 0.5;
125
126 // discard outliers : the cut is understood as a "distance" cut
127 cut = sq(nSigmas) * median;
128 nremoved = chi2_cleanup(*this, cut, *_transform);
129 } while (nremoved);
130 _dist2 = computeDist2(*this, *_transform);
131}
132
133/* not very robust : assumes that we went through Refine just before... */
135 int deno = (2. * size() - _transform->getNpar());
136 return (deno > 0) ? sqrt(_dist2 / deno) : -1; // is -1 a good idea?
137}
138
140 for (auto &smi : *this) smi.setDistance(transform); // c'est compact
141}
142
144 if (!which) return 0;
146 int initial_count = size();
147 if (which & 1) {
150 }
151 if (which & 2) {
154 }
155 return (initial_count - size());
156}
157
159 if (order == 0)
160 setTransform(std::make_shared<AstrometryTransformLinearShift>());
161 else if (order == 1)
162 setTransform(std::make_shared<AstrometryTransformLinear>());
163 else
165 // might consider throwing if order does not make sense (e.g. >10)
166 _order = order;
167}
168
169/* This routine should operate on a copy : refineTransform
170 might shorten the list */
171/* it is not const although it tries not to change anything */
173 if (!_transform) return nullptr;
174
175 auto old_transform = _transform->clone();
176 double old_chi2 = _chi2;
177
178 swap();
179 setTransformOrder(_order);
180 refineTransform(3.); // keep same order
181 auto inverted_transform = _transform->clone();
182 setTransform(old_transform.get());
183 swap();
184 _chi2 = old_chi2;
185
186 return inverted_transform;
187}
188
189void StarMatchList::cutTail(int nKeep) {
190 iterator si;
191 int count = 0;
192 for (si = begin(); si != end() && count < nKeep; ++count, ++si)
193 ;
194 erase(si, end());
195}
196
198 for (auto &starMatch : *this) {
199 starMatch.swap();
200 }
201}
202
203int StarMatchList::recoveredNumber(double mindist) const {
204 int n = 0;
206 for (auto const &starMatch : *this) {
207 if (starMatch.computeDistance(identity) < mindist) n++;
208 }
209 return (n);
210}
211
212void StarMatchList::applyTransform(StarMatchList &transformed, const AstrometryTransform *priorTransform,
213 const AstrometryTransform *posteriorTransform) const {
214 transformed.clear();
216 const AstrometryTransform &T1 = (priorTransform) ? *priorTransform : id;
217 const AstrometryTransform &T2 = (posteriorTransform) ? *posteriorTransform : id;
218
219 for (auto const &starMatch : *this) {
220 FatPoint p1;
221 T1.transformPosAndErrors(starMatch.point1, p1);
222 FatPoint p2;
223 T2.transformPosAndErrors(starMatch.point2, p2);
224 transformed.push_back(StarMatch(p1, p2, starMatch.s1, starMatch.s2));
225 }
226}
227
229 stream << " ================================================================" << std::endl
230 << " Transformation between lists of order " << getTransformOrder() << std::endl
231 << *_transform //<< endl
232 << " Chi2 = " << getChi2() << " Residual = " << computeResidual() << std::endl
233 << " Number in the list = " << size() << std::endl
234 << " ================================================================" << std::endl;
235}
236
237double computeDist2(const StarMatchList &starMatchList, const AstrometryTransform &transform) {
238 double dist2 = 0;
239 for (auto const &starMatch : starMatchList)
240 dist2 += transform.apply(starMatch.point1).computeDist2(starMatch.point2);
241 return dist2;
242}
243
244double computeChi2(const StarMatchList &starMatchList, const AstrometryTransform &transform) {
245 auto chi2s(chi2_array(starMatchList, transform));
246 double chi2 = 0;
247 for (auto const &c : chi2s) chi2 += c;
248 return chi2;
249}
250} // namespace jointcal
251} // namespace lsst
py::object result
Definition: _schema.cc:429
double x
table::Key< int > id
Definition: Detector.cc:162
pairs of points
table::Key< int > transform
T begin(T... args)
a virtual (interface) class for geometric transformations.
virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const
A do-nothing transformation. It anyway has dummy routines to mimick a AstrometryTransform.
implements the linear transformations (6 real coefficients).
A Point with uncertainties.
Definition: FatPoint.h:34
double x
coordinate
Definition: Point.h:42
A hanger for star associations.
Definition: StarMatch.h:55
double computeChi2(const AstrometryTransform &transform) const
returns the chi2 (using errors in the FatPoint's)
Definition: StarMatch.cc:43
FatPoint point2
2 points
Definition: StarMatch.h:61
std::unique_ptr< AstrometryTransform > inverseTransform()
returns the inverse transform (swap, fit(refineTransform) , and swap).
Definition: StarMatch.cc:172
void setDistance(const AstrometryTransform &transform)
Sets the distance (residual) field of all std::list elements. Mandatory before sorting on distances.
Definition: StarMatch.cc:139
void printTransform(std::ostream &stream=std::cout) const
print the matching transformation quality (transform, chi2, residual)
Definition: StarMatch.cc:228
void swap()
swaps elements 1 and 2 of each starmatch in std::list.
Definition: StarMatch.cc:197
double getChi2() const
access to the chi2 of the last call to refineTransform.
Definition: StarMatch.h:171
int getTransformOrder() const
returns the order of the used transform
Definition: StarMatch.h:174
void setTransform(const AstrometryTransform *transform)
sets a transform between the 2 std::lists and deletes the previous or default one....
Definition: StarMatch.h:188
int recoveredNumber(double mindist) const
count the number of elements for which distance is < mindist
Definition: StarMatch.cc:203
void setTransformOrder(int order)
set transform according to the given order.
Definition: StarMatch.cc:158
double computeResidual() const
returns the average 1d Residual (last call to refineTransform)
Definition: StarMatch.cc:134
unsigned removeAmbiguities(const AstrometryTransform &transform, int which=3)
cleans up the std::list of pairs for pairs that share one of their stars, keeping the closest one.
Definition: StarMatch.cc:143
void cutTail(int nKeep)
deletes the tail of the match std::list
Definition: StarMatch.cc:189
void applyTransform(StarMatchList &transformed, const AstrometryTransform *priorTransform, const AstrometryTransform *posteriorTransform=nullptr) const
enables to get a transformed StarMatchList.
Definition: StarMatch.cc:212
void refineTransform(double nSigmas)
removes pairs beyond nSigmas in distance (where the sigma scale is set by the fit) and iterates until...
Definition: StarMatch.cc:92
T clear(T... args)
T copy(T... args)
T count(T... args)
T end(T... args)
T endl(T... args)
StarMatch erase(StarMatch ... args)
bool compareStar1(const StarMatch &one, const StarMatch &two)
Definition: StarMatch.h:116
bool sameStar2(const StarMatch &one, const StarMatch &two)
Definition: StarMatch.h:126
bool sameStar1(const StarMatch &one, const StarMatch &two)
Definition: StarMatch.h:120
bool compareStar2(const StarMatch &one, const StarMatch &two)
Definition: StarMatch.h:122
double computeDist2(const StarMatchList &S, const AstrometryTransform &transform)
sum of distance squared
Definition: StarMatch.cc:237
double computeChi2(const StarMatchList &L, const AstrometryTransform &transform)
the actual chi2
Definition: StarMatch.cc:244
std::ostream & operator<<(std::ostream &stream, AstrometryMapping const &mapping)
A base class for image defects.
T push_back(T... args)
T size(T... args)
T sort(T... args)
T sqrt(T... args)
StarMatch unique(StarMatch ... args)
table::Key< int > order