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
ConstrainedAstrometryModel.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 "astshim.h"
26#include "lsst/afw/geom.h"
27#include "lsst/geom.h"
28#include "lsst/log/Log.h"
36
37#include "lsst/pex/exceptions.h"
39
40#include <iostream>
41#include <memory>
42#include <string>
43#include <utility>
44
45namespace {
46// Append the keys of this map into a comma-separated string.
47template <typename KeyType, typename ValueType>
48void outputMapKeys(std::map<KeyType, ValueType> const &map, std::ostream &os) {
49 bool first = true;
50 os << "[";
51 for (auto const &i : map) {
52 if (first)
53 first = false;
54 else
55 os << ", ";
56 os << i.first;
57 }
58 os << "]";
59}
60} // namespace
61
62namespace lsst {
63namespace jointcal {
64
66 CcdImageList const &ccdImageList, std::shared_ptr<ProjectionHandler const> projectionHandler,
67 int chipOrder, int visitOrder)
68 : AstrometryModel(LOG_GET("lsst.jointcal.ConstrainedAstrometryModel")),
69 _skyToTangentPlane(std::move(projectionHandler)) {
70 // keep track of which chip we want to hold fixed (the one closest to the middle of the focal plane)
71 double minRadius2 = std::numeric_limits<double>::infinity();
72 CcdIdType constrainedChip = -1;
73
74 // first loop to initialize all visit and chip transforms.
75 for (auto &ccdImage : ccdImageList) {
76 const CcdImage &im = *ccdImage;
77 auto visit = im.getVisit();
78 auto chip = im.getCcdId();
79 auto visitp = _visitMap.find(visit);
80 if (visitp == _visitMap.end()) {
81 _visitMap[visit] = std::make_shared<SimplePolyMapping>(AstrometryTransformLinear(),
83 }
84 auto chipp = _chipMap.find(chip);
85 if (chipp == _chipMap.end()) {
86 auto center = ccdImage->getDetector()->getCenter(afw::cameraGeom::FOCAL_PLANE);
87 double radius2 = std::pow(center.getX(), 2) + std::pow(center.getY(), 2);
88 if (radius2 < minRadius2) {
89 minRadius2 = radius2;
90 constrainedChip = chip;
91 }
92
93 auto pixelsToFocal =
94 im.getDetector()->getTransform(afw::cameraGeom::PIXELS, afw::cameraGeom::FOCAL_PLANE);
95 Frame const &frame = im.getImageFrame();
96 // construct the chip transform by approximating the pixel->Focal afw::geom::Transform.
98 AstrometryTransformPolynomial(pixelsToFocal, frame, chipOrder);
100 _chipMap[chip] = std::make_shared<SimplePolyMapping>(shiftAndNormalize,
101 pol * shiftAndNormalize.inverted());
102 }
103 }
104
105 // Hold the "central" chip map fixed and don't fit it, to remove a degeneracy.
106 _chipMap.at(constrainedChip)->setToBeFit(false);
107
108 // now, second loop to set the mappings of the CCdImages
109 for (auto &ccdImage : ccdImageList) {
110 const CcdImage &im = *ccdImage;
111 auto visit = im.getVisit();
112 auto chip = im.getCcdId();
113
114 // check that the chip_indexed part was indeed assigned
115 // (i.e. the reference visit was complete)
116 if (_chipMap.find(chip) == _chipMap.end()) {
117 LOGLS_WARN(_log, "Chip " << chip << " is missing in the reference exposure, expect troubles.");
119 _chipMap[chip] =
120 std::make_shared<SimplePolyMapping>(norm, AstrometryTransformPolynomial(chipOrder));
121 }
122 _mappings[ccdImage->getHashKey()] =
123 std::make_unique<ChipVisitAstrometryMapping>(_chipMap[chip], _visitMap[visit]);
124 }
125 LOGLS_INFO(_log, "Got " << _chipMap.size() << " chip mappings and " << _visitMap.size()
126 << " visit mappings; holding chip " << constrainedChip << " fixed ("
127 << getTotalParameters() << " total parameters).");
128 LOGLS_DEBUG(_log, "CcdImage map has " << _mappings.size() << " mappings, with "
129 << _mappings.bucket_count() << " buckets and a load factor of "
130 << _mappings.load_factor());
131}
132
134 return findMapping(ccdImage);
135}
136
142 Eigen::Index firstIndex) {
143 Eigen::Index index = firstIndex;
144 if (whatToFit.find("Distortions") == std::string::npos) {
145 LOGLS_ERROR(_log, "assignIndices was called and Distortions is *not* in whatToFit");
146 return 0;
147 }
148 // if we get here "Distortions" is in whatToFit
149 _fittingChips = (whatToFit.find("DistortionsChip") != std::string::npos);
150 _fittingVisits = (whatToFit.find("DistortionsVisit") != std::string::npos);
151 // If nothing more than "Distortions" is specified, it means all:
152 if ((!_fittingChips) && (!_fittingVisits)) {
153 _fittingChips = _fittingVisits = true;
154 }
155 if (_fittingChips)
156 for (auto &i : _chipMap) {
157 i.second->setIndex(index);
158 index += i.second->getNpar();
159 }
160 if (_fittingVisits)
161 for (auto &i : _visitMap) {
162 i.second->setIndex(index);
163 index += i.second->getNpar();
164 }
165 // Tell the mappings which derivatives they will have to fill:
166 for (auto &i : _mappings) {
167 i.second->setWhatToFit(_fittingChips, _fittingVisits);
168 }
169 return index;
170}
171
172void ConstrainedAstrometryModel::offsetParams(Eigen::VectorXd const &delta) {
173 if (_fittingChips)
174 for (auto &i : _chipMap) {
175 auto mapping = i.second.get();
176 mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
177 }
178 if (_fittingVisits)
179 for (auto &i : _visitMap) {
180 auto mapping = i.second.get();
181 mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
182 }
183}
184
186 for (auto & i : _visitMap) i.second->freezeErrorTransform();
187 for (auto & i : _chipMap) i.second->freezeErrorTransform();
188}
189
191 auto chipp = _chipMap.find(chip);
192 if (chipp == _chipMap.end()) {
193 std::stringstream errMsg;
194 errMsg << "No such chipId: " << chip << " among ";
195 outputMapKeys(_chipMap, errMsg);
198 }
199 return chipp->second->getTransform();
200}
201
202// Array of visits involved in the solution.
205 res.reserve(_visitMap.size());
206 for (const auto & i : _visitMap) res.push_back(i.first);
207 return res;
208}
209
211 auto visitp = _visitMap.find(visit);
212 if (visitp == _visitMap.end()) {
213 std::stringstream errMsg;
214 errMsg << "No such visitId: " << visit << " among ";
215 outputMapKeys(_visitMap, errMsg);
218 }
219 return visitp->second->getTransform();
220}
221
223 std::size_t total = 0;
224 for (auto &i : _chipMap) {
225 total += i.second->getNpar();
226 }
227 for (auto &i : _visitMap) {
228 total += i.second->getNpar();
229 }
230 return total;
231}
232
234 auto proj = std::dynamic_pointer_cast<const TanRaDecToPixel>(getSkyToTangentPlane(ccdImage));
235 jointcal::Point tangentPoint(proj->getTangentPoint());
236
237 auto imageFrame = ccdImage.getImageFrame();
238 auto pixelsToFocal = getChipTransform(ccdImage.getCcdId()).toAstMap(imageFrame);
239 jointcal::Frame focalBox = getChipTransform(ccdImage.getCcdId()).apply(imageFrame, false);
240 auto focalToIwc = getVisitTransform(ccdImage.getVisit()).toAstMap(focalBox);
241
242 ast::Frame pixelFrame(2, "Domain=PIXELS");
243 ast::Frame focalFrame(2, "Domain=FOCAL");
244 ast::Frame iwcFrame(2, "Domain=IWC");
245
246 // make a basic SkyWcs and extract the IWC portion
247 auto iwcToSkyWcs = afw::geom::makeSkyWcs(
248 geom::Point2D(0, 0), geom::SpherePoint(tangentPoint.x, tangentPoint.y, geom::degrees),
249 afw::geom::makeCdMatrix(1.0 * geom::degrees, 0 * geom::degrees, true));
250 auto iwcToSkyMap = iwcToSkyWcs->getFrameDict()->getMapping("PIXELS", "SKY");
251 auto skyFrame = iwcToSkyWcs->getFrameDict()->getFrame("SKY");
252
253 ast::FrameDict frameDict(pixelFrame);
254 frameDict.addFrame("PIXELS", *pixelsToFocal, focalFrame);
255 frameDict.addFrame("FOCAL", *focalToIwc, iwcFrame);
256 frameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
257 return std::make_shared<afw::geom::SkyWcs>(frameDict);
258}
259
261 out << "Constrained Astrometry Model (" << _mappings.size() << " composite mappings; " << _chipMap.size()
262 << " sensor mappings, " << _visitMap.size() << " visit mappings):" << std::endl;
263 out << *_skyToTangentPlane << std::endl;
264 out << "Sensor to sky transforms:" << std::endl;
265 for (auto &i : _mappings) {
266 out << i.first << std::endl;
267 out << *(i.second) << std::endl;
268 }
269}
270
272 auto i = _mappings.find(ccdImage.getHashKey());
273 if (i == _mappings.end())
275 "ConstrainedAstrometryModel cannot find CcdImage " + ccdImage.getName());
276 return i->second.get();
277}
278
279} // namespace jointcal
280} // namespace lsst
#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_INFO(logger, message)
Log a info-level message using an iostream-based interface.
Definition Log.h:639
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
Definition Log.h:679
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
Definition Log.h:619
pairs of points
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 void apply(double xIn, double yIn, double &xOut, double &yOut) const =0
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.
implements the linear transformations (6 real coefficients).
AstrometryTransformLinear inverted() const
returns the inverse: T1 = T2.inverted();
Handler of an actual image from a single CCD.
Definition CcdImage.h:64
CcdIdType getCcdId() const
returns ccd ID
Definition CcdImage.h:145
std::shared_ptr< afw::cameraGeom::Detector > getDetector() const
Definition CcdImage.h:150
VisitIdType getVisit() const
returns visit ID
Definition CcdImage.h:148
Frame const & getImageFrame() const
Frame in pixels.
Definition CcdImage.h:190
std::string getName() const
Return the _name that identifies this ccdImage.
Definition CcdImage.h:79
CcdImageKey getHashKey() const
Definition CcdImage.h:152
ConstrainedAstrometryModel(CcdImageList const &ccdImageList, std::shared_ptr< ProjectionHandler const > projectionHandler, int chipOrder, int visitOrder)
std::size_t getTotalParameters() const override
Return the total number of parameters in this model.
void offsetParams(Eigen::VectorXd const &Delta) override
Dispaches the offsets after a fit step into the actual locations of parameters.
void freezeErrorTransform() override
From there on, measurement errors are propagated using the current transforms (and no longer evolve).
void print(std::ostream &out) const override
Print a string representation of the contents of this mapping, for debugging.
AstrometryTransform const & getVisitTransform(VisitIdType const &visit) const
Access to mappings.
AstrometryTransform const & getChipTransform(CcdIdType chip) const
Access to mappings.
const std::shared_ptr< AstrometryTransform const > getSkyToTangentPlane(CcdImage const &ccdImage) const override
The mapping of sky coordinates (i.e.
AstrometryMapping * findMapping(CcdImage const &ccdImage) const override
Return a pointer to the mapping associated with this ccdImage.
std::vector< VisitIdType > getVisits() const
Access to array of visits involved in the solution.
Eigen::Index assignIndices(std::string const &whatToFit, Eigen::Index firstIndex) override
Positions the various parameter sets into the parameter vector, starting at firstIndex.
std::shared_ptr< afw::geom::SkyWcs > makeSkyWcs(CcdImage const &ccdImage) const override
Make a SkyWcs that contains this model.
AstrometryMapping const * getMapping(CcdImage const &) const override
Mapping associated to a given CcdImage.
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
Reports invalid arguments.
Definition Runtime.h:66
T endl(T... args)
T find(T... args)
T infinity(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
AngleUnit constexpr degrees
constant with units of degrees
Definition Angle.h:110
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
long VisitIdType
Definition CcdImage.h:48
STL namespace.
T pow(T... args)
T push_back(T... args)
T reserve(T... args)
T str(T... args)