LSSTApplications  11.0-24-g0a022a1,14.0+77,15.0,15.0+1
LSSTDataManagementBasePackage
astrometry_net.cc
Go to the documentation of this file.
1 // -*- lsst-C++ -*-
2 
3 #include <sstream>
4 #include <utility>
5 #include <vector>
6 
7 #include "lsst/pex/exceptions.h"
9 
10 namespace lsst {
11 namespace meas {
12 namespace extensions {
13 namespace astrometryNet {
14 
15 namespace {
16 
17 struct timer_baton {
18  solver_t* s;
19  double timelimit;
20 };
21 
22 static time_t timer_callback(void* baton) {
23  struct timer_baton* tt = static_cast<struct timer_baton*>(baton);
24  solver_t* solver = tt->s;
25  if (solver->timeused > tt->timelimit)
26  solver->quit_now = 1;
27  return 1;
28 }
29 
30 } // namespace <anonymous>
31 
32 MultiIndex::MultiIndex(std::string const & filepath) : _multiindex(multiindex_new(filepath.c_str())) {
33  if (!_multiindex) {
35  os << "Could not read multi-index star file " << filepath;
37  }
38 }
39 
40 void MultiIndex::addIndex(std::string const & filepath, bool metadataOnly) {
41  int const flags = metadataOnly ? INDEX_ONLY_LOAD_METADATA : 0;
42  if (multiindex_add_index(_multiindex.get(), filepath.c_str(), flags)) {
44  os << "Failed to read multiindex from \"" << filepath << "\""
45  << " with meatadaOnly=" << metadataOnly;
47  }
48 }
49 
50 int MultiIndex::isWithinRange(double ra, double dec, double radius_deg) {
51  return index_is_within_range(multiindex_get(_multiindex.get(), 0), ra, dec, radius_deg);
52 }
53 
55  if (multiindex_reload_starkd(_multiindex.get())) {
57  os << "Failed to reload multi-index star file " << getName();
59  }
60 }
61 
62 
63 Solver::Solver() : _solver(solver_new()) {}
64 
66  // Working around a bug in Astrometry.net: doesn't take ownership of the field.
67  // unseemly familiarity with the innards... but valgrind-clean.
68  starxy_free(_solver->fieldxy);
69  _solver->fieldxy = NULL;
70 }
71 
74  lsst::afw::coord::Coord const &ctrCoord,
75  lsst::afw::geom::Angle const &radius,
76  const char* idCol,
77  std::vector<std::string> const& filterNameList,
78  std::vector<std::string> const& magColList,
79  std::vector<std::string> const& magErrColList,
80  const char* starGalCol,
81  const char* varCol,
82  bool uniqueIds)
83 {
84  if ((filterNameList.size() != magColList.size()) || (filterNameList.size() != magErrColList.size())) {
86  "Filter name, mag column, and mag error column vectors must be the same length.");
87  }
88  std::vector<detail::MagColInfo> magColInfoList;
89  for (size_t i=0; i<filterNameList.size(); ++i) {
91  mc.filterName = filterNameList[i];
92  mc.magCol = magColList[i];
93  mc.magErrCol = magErrColList[i];
94  magColInfoList.push_back(mc);
95  }
96  return detail::getCatalogImpl(inds, ctrCoord, radius,
97  idCol, magColInfoList, starGalCol, varCol, uniqueIds);
98 }
99 
101  // Gather solve stats...
102  auto qa = std::make_shared<daf::base::PropertyList>();
103  // FIXME -- Ticket #1875 prevents dotted-names from working with toString().
104  qa->set("meas_astrom*an*n_tried", _solver->numtries);
105  qa->set("meas_astrom*an*n_matched", _solver->nummatches);
106  qa->set("meas_astrom*an*n_scaleok", _solver->numscaleok);
107  qa->set("meas_astrom*an*n_cxdxcut", _solver->num_cxdx_skipped);
108  qa->set("meas_astrom*an*n_meanxcut", _solver->num_meanx_skipped);
109  qa->set("meas_astrom*an*n_radeccut", _solver->num_radec_skipped);
110  qa->set("meas_astrom*an*n_scalecut", _solver->num_abscale_skipped);
111  qa->set("meas_astrom*an*n_verified", _solver->num_verified);
112  qa->set("meas_astrom*an*time_used", _solver->timeused);
113  qa->set("meas_astrom*an*best_logodds", _solver->best_logodds);
114  if (_solver->best_index) {
115  index_t* ind = _solver->best_index;
116  qa->set("meas_astrom*an*best_index*id", ind->indexid);
117  qa->set("meas_astrom*an*best_index*hp", ind->healpix);
118  qa->set("meas_astrom*an*best_index*nside", ind->hpnside);
119  qa->set("meas_astrom*an*best_index*name", std::string(ind->indexname));
120  }
121  if (_solver->have_best_match) {
122  MatchObj* mo = &(_solver->best_match);
123  std::string s = boost::str(boost::format("%i") % mo->star[0]);
124  for (int i=1; i<mo->dimquads; i++)
125  s = s + boost::str(boost::format(", %i") % mo->star[i]);
126  qa->set("meas_astrom*an*best_match*starinds", s);
127  qa->set("meas_astrom*an*best_match*coderr", std::sqrt(mo->code_err));
128  qa->set("meas_astrom*an*best_match*nmatch", mo->nmatch);
129  qa->set("meas_astrom*an*best_match*ndistract", mo->ndistractor);
130  qa->set("meas_astrom*an*best_match*nconflict", mo->nconflict);
131  qa->set("meas_astrom*an*best_match*nfield", mo->nfield);
132  qa->set("meas_astrom*an*best_match*nindex", mo->nindex);
133  qa->set("meas_astrom*an*best_match*nbest", mo->nbest);
134  qa->set("meas_astrom*an*best_match*logodds", mo->logodds);
135  qa->set("meas_astrom*an*best_match*parity", mo->parity ? 0 : 1);
136  qa->set("meas_astrom*an*best_match*nobjs", mo->objs_tried);
137  }
138  return qa;
139 }
140 
142  MatchObj* match = solver_get_best_match(_solver.get());
143  if (!match)
145  tan_t* wcs = &(match->wcstan);
146 
147  afw::geom::Point2D crpix(wcs->crpix[0], wcs->crpix[1]);
148  auto const crval = afw::coord::IcrsCoord(wcs->crval[0] * afw::geom::degrees,
149  wcs->crval[1] * afw::geom::degrees);
150  Eigen::Matrix2d cdMatrix;
151  for (int i = 0; i < 2; ++i) {
152  for (int j = 0; j < 2; ++j) {
153  cdMatrix(i, j) = wcs->cd[i][j];
154  }
155  }
156  return afw::geom::makeSkyWcs(crpix, crval, cdMatrix);
157 }
158 
159 void Solver::run(double cpulimit) {
160  solver_log_params(_solver.get());
161  struct timer_baton tt;
162  if (cpulimit > 0.) {
163  tt.s = _solver.get();
164  tt.timelimit = cpulimit;
165  _solver->userdata = &tt;
166  _solver->timer_callback = timer_callback;
167  }
168  solver_run(_solver.get());
169  if (cpulimit > 0.) {
170  _solver->timer_callback = NULL;
171  _solver->userdata = NULL;
172  }
173 }
174 
182  for (std::vector<index_t*>::iterator pind = inds.begin();
183  pind != inds.end(); ++pind) {
184  detail::IndexManager man(*pind);
185 // printf("Checking index \"%s\"\n", man.index->indexname);
186  if (_solver->use_radec) {
187  double ra,dec,radius;
188  xyzarr2radecdeg(_solver->centerxyz, &ra, &dec);
189  radius = distsq2deg(_solver->r2);
190  if (!index_is_within_range(man.index, ra, dec, radius)) {
191  //printf("Not within RA,Dec range\n");
192  continue;
193  }
194  }
195  // qlo,qhi in arcsec
196  double qlo, qhi;
197  solver_get_quad_size_range_arcsec(_solver.get(), &qlo, &qhi);
198  if (!index_overlaps_scale_range(man.index, qlo, qhi)) {
199 // printf("Not within quad scale range\n");
200  continue;
201  }
202 // printf("Adding index.\n");
203  if (index_reload(man.index)) {
205  "Failed to index_reload() an astrometry_net_data index file -- out of file descriptors?");
206  }
207 
208  solver_add_index(_solver.get(), man.index);
209  }
210 }
211 
212 void Solver::setImageSize(int width, int height) {
213  solver_set_field_bounds(_solver.get(), 0, width, 0, height);
214  double hi = hypot(width, height);
215  double lo = 0.1 * std::min(width, height);
216  solver_set_quad_size_range(_solver.get(), lo, hi);
217 }
218 
219 void Solver::setStars(lsst::afw::table::SourceCatalog const & srcs, int x0, int y0) {
220  // convert to Astrometry.net "starxy_t"
221  starxy_free(_solver->fieldxy);
222  const size_t N = srcs.size();
223  starxy_t *starxy = starxy_new(N, true, false);
224  for (size_t i=0; i<N; ++i) {
225  double const x = srcs[i].getX();
226  double const y = srcs[i].getY();
227  double const flux = srcs[i].getPsfFlux();
228  starxy_set(starxy, i, x - x0, y - y0);
229  starxy_set_flux(starxy, i, flux);
230  }
231  // Sort the array
232  starxy_sort_by_flux(starxy);
233 
234  starxy_free(solver_get_field(_solver.get()));
235  solver_free_field(_solver.get());
236  solver_set_field(_solver.get(), starxy);
237  solver_reset_field_size(_solver.get());
238  // Find field boundaries and precompute kdtree
239  solver_preprocess_field(_solver.get());
240 }
241 
242 
244  lsst::afw::coord::IcrsCoord icrs = coord.toIcrs();
245  return lsst::afw::geom::Angle(healpix_distance_to_radec(hp, nside, icrs.getLongitude().asDegrees(),
246  icrs.getLatitude().asDegrees(), NULL),
248 }
249 
250 }}}} // namespace lsst::meas::extensions::astrometryNet
std::string magCol
name of magnitude column
Definition: utils.h:19
lsst::afw::table::SimpleCatalog getCatalogImpl(std::vector< index_t *> inds, lsst::afw::coord::Coord const &ctrCoord, lsst::afw::geom::Angle const &radius, const char *idCol, std::vector< MagColInfo > const &magColInfoList, const char *starGalCol, const char *varCol, bool uniqueIds=true)
Implementation for index_t::getCatalog method.
Definition: utils.cc:59
solver_t * s
lsst::afw::geom::Angle getLongitude() const
The main access method for the longitudinal coordinate.
Definition: Coord.h:166
list extensions
Definition: conf.py:22
lsst::afw::geom::Angle Angle
Definition: misc.h:36
std::string magErrCol
name of magnitude sigma column
Definition: utils.h:20
std::unique_ptr< multiindex_t, _Deleter > _multiindex
AngleUnit constexpr degrees
constant with units of degrees
Definition: Angle.h:88
MultiIndex(std::string const &filepath)
Construct a MultiIndex from an astrometry.net multi-index file.
virtual IcrsCoord toIcrs() const
Convert ourself to ICRS: RA, Dec (basically J2000)
Definition: Coord.cc:572
tbl::Key< int > wcs
T end(T... args)
std::shared_ptr< lsst::afw::geom::SkyWcs > getWcs()
lsst::afw::geom::Angle getLatitude() const
The main access method for the latitudinal coordinate.
Definition: Coord.h:173
std::unique_ptr< solver_t, _Deleter > _solver
STL class.
T min(T... args)
Include files required for standard LSST Exception handling.
double timelimit
T push_back(T... args)
A class representing an angle.
Definition: Angle.h:102
A base class for image defects.
Definition: cameraGeom.dox:3
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:63
RAII manager for astrometry.net indices.
Definition: utils.h:32
lsst::afw::coord::IcrsCoord IcrsCoord
Definition: misc.h:35
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:134
T str(T... args)
double x
std::shared_ptr< lsst::daf::base::PropertyList > getSolveStats() const
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type and message and optionally other arguments (dependent on the ty...
Definition: Exception.h:46
STL class.
T begin(T... args)
void addIndex(std::string const &filepath, bool metadataOnly)
Add an index read from a file.
T c_str(T... args)
T sqrt(T... args)
lsst::afw::geom::Angle healpixDistance(int hp, int nside, lsst::afw::coord::Coord const &coord)
Calculate the distance from coordinates to a healpix.
constexpr double asDegrees() const noexcept
Return an Angle&#39;s value in degrees.
Definition: Angle.h:144
This is the base class for spherical coordinates.
Definition: Coord.h:63
A class to handle Icrs coordinates (inherits from Coord)
Definition: Coord.h:291
void setStars(lsst::afw::table::SourceCatalog const &srcs, int x0, int y0)
lsst::afw::table::SimpleCatalog getCatalog(std::vector< index_t *> inds, lsst::afw::coord::Coord const &ctrCoord, lsst::afw::geom::Angle const &radius, const char *idCol, std::vector< std::string > const &filterNameList, std::vector< std::string > const &magColList, std::vector< std::string > const &magErrColList, const char *starGalCol, const char *varCol, bool uniqueIds=true)
Load reference objects in a region of the sky described by a center coordinate and a radius...
int isWithinRange(double ra, double dec, double radius_deg)
Is this multi-index in range of the specified cone?
void addIndices(std::vector< index_t *> inds)
Add indices to the solver.