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