LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
SimpleAstrometryModel.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 <memory>
26 #include <string>
27 
28 #include "astshim.h"
29 #include "lsst/log/Log.h"
33 #include "lsst/jointcal/CcdImage.h"
35 #include "lsst/pex/exceptions.h"
37 
38 namespace lsst {
39 namespace jointcal {
40 
42  const std::shared_ptr<ProjectionHandler const> projectionHandler,
43  bool initFromWcs, unsigned nNotFit, unsigned order)
44  : AstrometryModel(LOG_GET("lsst.jointcal.SimpleAstrometryModel")),
45  _skyToTangentPlane(projectionHandler)
46 
47 {
48  std::size_t count = 0;
49 
50  for (auto i = ccdImageList.cbegin(); i != ccdImageList.cend(); ++i, ++count) {
51  const CcdImage &im = **i;
52  if (count < nNotFit) {
55  id->setIndex(-1); // non sense, because it has no parameters
56  _myMap[im.getHashKey()] = std::move(id);
57  } else
58  // Given how AssignIndices works, only the SimplePolyMappings
59  // will actually be fitted, as nNotFit requests.
60  {
61  // first check that there are enough measurements for the requested polynomial order.
62  size_t nObj = im.getCatalogForFit().size();
63  if (nObj == 0) {
64  LOGLS_WARN(_log, "Empty catalog from image: " << im.getName());
65  continue;
66  }
68  if (pol.getOrder() > 0) // if not, it cannot be decreased
69  {
70  while (pol.getNpar() > 2 * nObj) {
71  LOGLS_WARN(_log, "Reducing polynomial order from "
72  << pol.getOrder() << ", due to too few sources (" << nObj
73  << " vs. " << pol.getNpar() << " parameters)");
74  pol.setOrder(pol.getOrder() - 1);
75  }
76  }
77  /* We have to center and normalize the coordinates so that
78  the fit matrix is not too ill-conditionned. Basically, x
79  and y in pixels are mapped to [-1,1]. When the
80  transformation of SimplePolyMapping transformation is
81  accessed, the combination of the normalization and the
82  fitted transformation is returned, so that the trick
83  remains hidden
84  */
85  const Frame &frame = im.getImageFrame();
87  if (initFromWcs) {
89  pol = pol * shiftAndNormalize.inverted();
90  }
91  _myMap[im.getHashKey()] =
93  }
94  }
95 }
96 
98  return findMapping(ccdImage);
99 }
100 
101 Eigen::Index SimpleAstrometryModel::assignIndices(std::string const &whatToFit, Eigen::Index firstIndex) {
102  if (whatToFit.find("Distortions") == std::string::npos) {
103  LOGLS_ERROR(_log, "AssignIndices was called and Distortions is *not* in whatToFit.");
104  return 0;
105  }
106  Eigen::Index index = firstIndex;
107  for (auto i = _myMap.begin(); i != _myMap.end(); ++i) {
108  SimplePolyMapping *p = dynamic_cast<SimplePolyMapping *>(&*(i->second));
109  if (!p) continue; // it should be AstrometryTransformIdentity
110  p->setIndex(index);
111  index += p->getNpar();
112  }
113  return index;
114 }
115 
116 void SimpleAstrometryModel::offsetParams(Eigen::VectorXd const &delta) {
117  for (auto &i : _myMap) {
118  auto mapping = i.second.get();
119  mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
120  }
121 }
122 
124  for (auto &i : _myMap) i.second->freezeErrorTransform();
125 }
126 
128  std::size_t total = 0;
129  for (auto &i : _myMap) {
130  total += i.second->getNpar();
131  }
132  return total;
133 }
134 
136  out << "SimpleAstrometryModel: " << _myMap.size() << " mappings" << std::endl;
137  out << *_skyToTangentPlane << std::endl;
138  out << "Sensor to sky transforms:" << std::endl;
139  for (auto &i : _myMap) {
140  out << i.first << std::endl;
141  out << *(i.second) << std::endl;
142  out << std::endl;
143  }
144 }
145 
147  return dynamic_cast<const SimplePolyMapping *>(findMapping(ccdImage))->getTransform();
148 }
149 
151  auto proj = std::dynamic_pointer_cast<const TanRaDecToPixel>(getSkyToTangentPlane(ccdImage));
152  jointcal::Point tangentPoint(proj->getTangentPoint());
153 
154  auto polyMap = getTransform(ccdImage).toAstMap(ccdImage.getImageFrame());
155  ast::Frame pixelFrame(2, "Domain=PIXELS");
156  ast::Frame iwcFrame(2, "Domain=IWC");
157 
158  // make a basic SkyWcs and extract the IWC portion
159  auto iwcToSkyWcs = afw::geom::makeSkyWcs(
160  geom::Point2D(0, 0), geom::SpherePoint(tangentPoint.x, tangentPoint.y, geom::degrees),
162  auto iwcToSkyMap = iwcToSkyWcs->getFrameDict()->getMapping("PIXELS", "SKY");
163  auto skyFrame = iwcToSkyWcs->getFrameDict()->getFrame("SKY");
164 
165  ast::FrameDict frameDict(pixelFrame, *polyMap, iwcFrame);
166  frameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
167  return std::make_shared<afw::geom::SkyWcs>(frameDict);
168 }
169 
170 AstrometryMapping *SimpleAstrometryModel::findMapping(CcdImage const &ccdImage) const {
171  auto i = _myMap.find(ccdImage.getHashKey());
172  if (i == _myMap.end())
174  "SimpleAstrometryModel cannot find CcdImage " + ccdImage.getName());
175  return i->second.get();
176 }
177 
178 } // namespace jointcal
179 } // namespace lsst
table::Key< int > id
Definition: Detector.cc:162
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:659
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
Definition: Log.h:679
T cbegin(T... args)
A FrameSet whose frames can be referenced by domain name.
Definition: FrameDict.h:67
void addFrame(int iframe, Mapping const &map, Frame const &frame) override
Add a new Frame and an associated Mapping to this FrameSet so as to define a new coordinate system,...
Definition: FrameDict.cc:32
Frame is used to represent a coordinate system.
Definition: Frame.h:157
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
virtual class needed in the abstraction of the distortion model
Interface between AstrometryFit and the combinations of Mappings from pixels to some tangent plane (a...
LOG_LOGGER _log
lsst.logging instance, to be created by a subclass so that messages have consistent name.
a virtual (interface) class for geometric transformations.
virtual std::shared_ptr< ast::Mapping > toAstMap(jointcal::Frame const &domain) const
Create an equivalent AST mapping for this transformation, including an analytic inverse if possible.
A do-nothing transformation. It anyway has dummy routines to mimick a AstrometryTransform.
implements the linear transformations (6 real coefficients).
AstrometryTransformLinear inverted() const
returns the inverse: T1 = T2.inverted();
std::size_t getNpar() const override
total number of parameters
void setOrder(std::size_t order)
Sets the polynomial order (the highest sum of exponents of the largest monomial).
std::size_t getOrder() const
Returns the polynomial order.
Handler of an actual image from a single CCD.
Definition: CcdImage.h:64
std::shared_ptr< AstrometryTransform > const getPixelToTangentPlane() const
Definition: CcdImage.h:140
Frame const & getImageFrame() const
Frame in pixels.
Definition: CcdImage.h:190
MeasuredStarList const & getCatalogForFit() const
Gets the catalog to be used for fitting, which may have been cleaned-up.
Definition: CcdImage.h:94
std::string getName() const
Return the _name that identifies this ccdImage.
Definition: CcdImage.h:79
CcdImageKey getHashKey() const
Definition: CcdImage.h:152
rectangle with sides parallel to axes.
Definition: Frame.h:38
A point in a plane.
Definition: Point.h:37
double x
coordinate
Definition: Point.h:42
std::size_t getNpar() const override
Number of parameters in total.
void setIndex(Eigen::Index i)
Set the index of this mapping in the grand fit.
Eigen::Index assignIndices(std::string const &whatToFit, Eigen::Index firstIndex) override
Positions the various parameter sets into the parameter vector, starting at firstIndex.
void offsetParams(Eigen::VectorXd const &delta) override
Offset the parameters by the provided amounts (by -delta).
std::size_t getTotalParameters() const override
Return the total number of parameters in this model.
AstrometryTransform const & getTransform(CcdImage const &ccdImage) const
Access to mappings.
void print(std::ostream &out) const override
Print a string representation of the contents of this mapping, for debugging.
const std::shared_ptr< AstrometryTransform const > getSkyToTangentPlane(CcdImage const &ccdImage) const override
the mapping of sky coordinates (i.e.
const AstrometryMapping * getMapping(CcdImage const &) const override
Mapping associated to a given CcdImage.
std::shared_ptr< afw::geom::SkyWcs > makeSkyWcs(CcdImage const &ccdImage) const override
Make a SkyWcs that contains this model.
SimpleAstrometryModel(CcdImageList const &ccdImageList, const std::shared_ptr< ProjectionHandler const > projectionHandler, bool initFromWCS, unsigned nNotFit=0, unsigned order=3)
Mapping implementation for a polynomial transformation.
Reports invalid arguments.
Definition: Runtime.h:66
T count(T... args)
T cend(T... args)
T endl(T... args)
T find(T... args)
T move(T... args)
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const &scale, lsst::geom::Angle const &orientation=0 *lsst::geom::degrees, bool flipX=false)
Make a WCS CD matrix.
Definition: SkyWcs.cc:133
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
A base class for image defects.
T size(T... args)
table::Key< int > order