LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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"
31 #include "lsst/jointcal/CcdImage.h"
36 
37 #include "lsst/pex/exceptions.h"
39 
40 #include <memory>
41 #include <string>
42 #include <iostream>
43 
44 namespace {
45 // Append the keys of this map into a comma-separated string.
46 template <typename KeyType, typename ValueType>
47 void outputMapKeys(std::map<KeyType, ValueType> const &map, std::ostream &os) {
48  bool first = true;
49  os << "[";
50  for (auto const &i : map) {
51  if (first)
52  first = false;
53  else
54  os << ", ";
55  os << i.first;
56  }
57  os << "]";
58 }
59 } // namespace
60 
61 namespace lsst {
62 namespace jointcal {
63 
65  CcdImageList const &ccdImageList, std::shared_ptr<ProjectionHandler const> projectionHandler,
66  int chipOrder, int visitOrder)
67  : AstrometryModel(LOG_GET("jointcal.ConstrainedAstrometryModel")),
68  _skyToTangentPlane(projectionHandler) {
69  // keep track of which chip we want to hold fixed (the one closest to the middle of the focal plane)
70  double minRadius2 = std::numeric_limits<double>::infinity();
71  CcdIdType constrainedChip = -1;
72 
73  // first loop to initialize all visit and chip transforms.
74  for (auto &ccdImage : ccdImageList) {
75  const CcdImage &im = *ccdImage;
76  auto visit = im.getVisit();
77  auto chip = im.getCcdId();
78  auto visitp = _visitMap.find(visit);
79  if (visitp == _visitMap.end()) {
80  _visitMap[visit] = std::make_shared<SimplePolyMapping>(AstrometryTransformLinear(),
81  AstrometryTransformPolynomial(visitOrder));
82  }
83  auto chipp = _chipMap.find(chip);
84  if (chipp == _chipMap.end()) {
85  auto center = ccdImage->getDetector()->getCenter(afw::cameraGeom::FOCAL_PLANE);
86  double radius2 = std::pow(center.getX(), 2) + std::pow(center.getY(), 2);
87  if (radius2 < minRadius2) {
88  minRadius2 = radius2;
89  constrainedChip = chip;
90  }
91 
92  auto pixelsToFocal =
94  Frame const &frame = im.getImageFrame();
95  // construct the chip transform by approximating the pixel->Focal afw::geom::Transform.
97  AstrometryTransformPolynomial(pixelsToFocal, frame, chipOrder);
99  _chipMap[chip] = std::make_shared<SimplePolyMapping>(shiftAndNormalize,
100  pol * shiftAndNormalize.inverted());
101  }
102  }
103 
104  // Hold the "central" chip map fixed and don't fit it, to remove a degeneracy.
105  _chipMap.at(constrainedChip)->setToBeFit(false);
106 
107  // now, second loop to set the mappings of the CCdImages
108  for (auto &ccdImage : ccdImageList) {
109  const CcdImage &im = *ccdImage;
110  auto visit = im.getVisit();
111  auto chip = im.getCcdId();
112 
113  // check that the chip_indexed part was indeed assigned
114  // (i.e. the reference visit was complete)
115  if (_chipMap.find(chip) == _chipMap.end()) {
116  LOGLS_WARN(_log, "Chip " << chip << " is missing in the reference exposure, expect troubles.");
118  _chipMap[chip] =
119  std::make_shared<SimplePolyMapping>(norm, AstrometryTransformPolynomial(chipOrder));
120  }
121  _mappings[ccdImage->getHashKey()] =
122  std::make_unique<ChipVisitAstrometryMapping>(_chipMap[chip], _visitMap[visit]);
123  }
124  LOGLS_INFO(_log, "Got " << _chipMap.size() << " chip mappings and " << _visitMap.size()
125  << " visit mappings; holding chip " << constrainedChip << " fixed ("
126  << getTotalParameters() << " total parameters).");
127  LOGLS_DEBUG(_log, "CcdImage map has " << _mappings.size() << " mappings, with "
128  << _mappings.bucket_count() << " buckets and a load factor of "
129  << _mappings.load_factor());
130 }
131 
133  return findMapping(ccdImage);
134 }
135 
141  Eigen::Index firstIndex) {
142  Eigen::Index index = firstIndex;
143  if (whatToFit.find("Distortions") == std::string::npos) {
144  LOGLS_ERROR(_log, "assignIndices was called and Distortions is *not* in whatToFit");
145  return 0;
146  }
147  // if we get here "Distortions" is in whatToFit
148  _fittingChips = (whatToFit.find("DistortionsChip") != std::string::npos);
149  _fittingVisits = (whatToFit.find("DistortionsVisit") != std::string::npos);
150  // If nothing more than "Distortions" is specified, it means all:
151  if ((!_fittingChips) && (!_fittingVisits)) {
152  _fittingChips = _fittingVisits = true;
153  }
154  if (_fittingChips)
155  for (auto &i : _chipMap) {
156  i.second->setIndex(index);
157  index += i.second->getNpar();
158  }
159  if (_fittingVisits)
160  for (auto &i : _visitMap) {
161  i.second->setIndex(index);
162  index += i.second->getNpar();
163  }
164  // Tell the mappings which derivatives they will have to fill:
165  for (auto &i : _mappings) {
166  i.second->setWhatToFit(_fittingChips, _fittingVisits);
167  }
168  return index;
169 }
170 
171 void ConstrainedAstrometryModel::offsetParams(Eigen::VectorXd const &delta) {
172  if (_fittingChips)
173  for (auto &i : _chipMap) {
174  auto mapping = i.second.get();
175  mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
176  }
177  if (_fittingVisits)
178  for (auto &i : _visitMap) {
179  auto mapping = i.second.get();
180  mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
181  }
182 }
183 
185  for (auto i = _visitMap.begin(); i != _visitMap.end(); ++i) i->second->freezeErrorTransform();
186  for (auto i = _chipMap.begin(); i != _chipMap.end(); ++i) i->second->freezeErrorTransform();
187 }
188 
190  auto chipp = _chipMap.find(chip);
191  if (chipp == _chipMap.end()) {
192  std::stringstream errMsg;
193  errMsg << "No such chipId: " << chip << " among ";
194  outputMapKeys(_chipMap, errMsg);
195  std::cout << std::endl;
196  throw pexExcept::InvalidParameterError(errMsg.str());
197  }
198  return chipp->second->getTransform();
199 }
200 
201 // Array of visits involved in the solution.
204  res.reserve(_visitMap.size());
205  for (auto i = _visitMap.begin(); i != _visitMap.end(); ++i) res.push_back(i->first);
206  return res;
207 }
208 
210  auto visitp = _visitMap.find(visit);
211  if (visitp == _visitMap.end()) {
212  std::stringstream errMsg;
213  errMsg << "No such visitId: " << visit << " among ";
214  outputMapKeys(_visitMap, errMsg);
215  std::cout << std::endl;
216  throw pexExcept::InvalidParameterError(errMsg.str());
217  }
218  return visitp->second->getTransform();
219 }
220 
222  std::size_t total = 0;
223  for (auto &i : _chipMap) {
224  total += i.second->getNpar();
225  }
226  for (auto &i : _visitMap) {
227  total += i.second->getNpar();
228  }
229  return total;
230 }
231 
233  auto proj = std::dynamic_pointer_cast<const TanRaDecToPixel>(getSkyToTangentPlane(ccdImage));
234  jointcal::Point tangentPoint(proj->getTangentPoint());
235 
236  auto imageFrame = ccdImage.getImageFrame();
237  auto pixelsToFocal = getChipTransform(ccdImage.getCcdId()).toAstMap(imageFrame);
238  jointcal::Frame focalBox = getChipTransform(ccdImage.getCcdId()).apply(imageFrame, false);
239  auto focalToIwc = getVisitTransform(ccdImage.getVisit()).toAstMap(focalBox);
240 
241  ast::Frame pixelFrame(2, "Domain=PIXELS");
242  ast::Frame focalFrame(2, "Domain=FOCAL");
243  ast::Frame iwcFrame(2, "Domain=IWC");
244 
245  // make a basic SkyWcs and extract the IWC portion
246  auto iwcToSkyWcs = afw::geom::makeSkyWcs(
247  geom::Point2D(0, 0), geom::SpherePoint(tangentPoint.x, tangentPoint.y, geom::degrees),
249  auto iwcToSkyMap = iwcToSkyWcs->getFrameDict()->getMapping("PIXELS", "SKY");
250  auto skyFrame = iwcToSkyWcs->getFrameDict()->getFrame("SKY");
251 
252  ast::FrameDict frameDict(pixelFrame);
253  frameDict.addFrame("PIXELS", *pixelsToFocal, focalFrame);
254  frameDict.addFrame("FOCAL", *focalToIwc, iwcFrame);
255  frameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
256  return std::make_shared<afw::geom::SkyWcs>(frameDict);
257 }
258 
260  out << "Constrained Astrometry Model (" << _mappings.size() << " composite mappings; " << _chipMap.size()
261  << " sensor mappings, " << _visitMap.size() << " visit mappings):" << std::endl;
262  out << *_skyToTangentPlane << std::endl;
263  out << "Sensor to sky transforms:" << std::endl;
264  for (auto &i : _mappings) {
265  out << i.first << std::endl;
266  out << *(i.second) << std::endl;
267  }
268 }
269 
270 AstrometryMapping *ConstrainedAstrometryModel::findMapping(CcdImage const &ccdImage) const {
271  auto i = _mappings.find(ccdImage.getHashKey());
272  if (i == _mappings.end())
274  "ConstrainedAstrometryModel cannot find CcdImage " + ccdImage.getName());
275  return i->second.get();
276 }
277 
278 } // namespace jointcal
279 } // namespace lsst
lsst::geom::degrees
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
lsst::jointcal::CcdImage::getHashKey
CcdImageKey getHashKey() const
Definition: CcdImage.h:152
std::string
STL class.
std::shared_ptr
STL class.
lsst::jointcal::ConstrainedAstrometryModel::ConstrainedAstrometryModel
ConstrainedAstrometryModel(CcdImageList const &ccdImageList, std::shared_ptr< ProjectionHandler const > projectionHandler, int chipOrder, int visitOrder)
Definition: ConstrainedAstrometryModel.cc:64
lsst::jointcal::ConstrainedAstrometryModel::getVisits
std::vector< VisitIdType > getVisits() const
Access to array of visits involved in the solution.
Definition: ConstrainedAstrometryModel.cc:202
std::list< std::shared_ptr< CcdImage > >
CcdImage.h
lsst::jointcal::ConstrainedAstrometryModel::getTotalParameters
std::size_t getTotalParameters() const override
Return the total number of parameters in this model.
Definition: ConstrainedAstrometryModel.cc:221
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:194
lsst::jointcal::CcdIdType
int CcdIdType
Definition: CcdImage.h:49
AstrometryTransform.h
std::vector::reserve
T reserve(T... args)
Eigenstuff.h
lsst::jointcal::AstrometryTransformLinear
implements the linear transformations (6 real coefficients).
Definition: AstrometryTransform.h:426
std::vector
STL class.
std::string::find
T find(T... args)
lsst::jointcal::ConstrainedAstrometryModel::freezeErrorTransform
void freezeErrorTransform() override
From there on, measurement errors are propagated using the current transforms (and no longer evolve).
Definition: ConstrainedAstrometryModel.cc:184
lsst::afw::table._match.first
first
Definition: _match.py:76
std::stringstream
STL class.
LOG_GET
#define LOG_GET(logger)
Definition: Log.h:75
ast::FrameDict
A FrameSet whose frames can be referenced by domain name.
Definition: FrameDict.h:67
LOGLS_INFO
#define LOGLS_INFO(logger, message)
Definition: Log.h:628
ast::FrameDict::addFrame
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
LOGLS_ERROR
#define LOGLS_ERROR(logger, message)
Definition: Log.h:668
geom.h
geom.h
astshim.h
lsst::jointcal::CcdImage::getName
std::string getName() const
Return the _name that identifies this ccdImage.
Definition: CcdImage.h:79
std::vector::push_back
T push_back(T... args)
lsst::jointcal::Point::y
double y
Definition: Point.h:41
LOGLS_WARN
#define LOGLS_WARN(logger, message)
Definition: Log.h:648
lsst::jointcal::ConstrainedAstrometryModel::assignIndices
Eigen::Index assignIndices(std::string const &whatToFit, Eigen::Index firstIndex) override
Positions the various parameter sets into the parameter vector, starting at firstIndex.
Definition: ConstrainedAstrometryModel.cc:140
ast::Frame
Frame is used to represent a coordinate system.
Definition: Frame.h:157
lsst::jointcal::CcdImage::getImageFrame
Frame const & getImageFrame() const
Frame in pixels.
Definition: CcdImage.h:191
lsst::jointcal::AstrometryTransform::apply
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const =0
std::numeric_limits::infinity
T infinity(T... args)
lsst::jointcal::AstrometryModel
Interface between AstrometryFit and the combinations of Mappings from pixels to some tangent plane (a...
Definition: AstrometryModel.h:50
ProjectionHandler.h
std::ostream
STL class.
lsst::jointcal::AstrometryTransformLinear::inverted
AstrometryTransformLinear inverted() const
returns the inverse: T1 = T2.inverted();
Definition: AstrometryTransform.cc:1263
lsst::jointcal::CcdImage
Handler of an actual image from a single CCD.
Definition: CcdImage.h:64
lsst::jointcal::ConstrainedAstrometryModel::getMapping
AstrometryMapping const * getMapping(CcdImage const &) const override
Mapping associated to a given CcdImage.
Definition: ConstrainedAstrometryModel.cc:132
lsst::jointcal::ConstrainedAstrometryModel::getVisitTransform
AstrometryTransform const & getVisitTransform(VisitIdType const &visit) const
Access to mappings.
Definition: ConstrainedAstrometryModel.cc:209
lsst::jointcal::normalizeCoordinatesTransform
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
Definition: AstrometryTransform.cc:833
lsst::jointcal::AstrometryModel::_log
LOG_LOGGER _log
lsst.logging instance, to be created by a subclass so that messages have consistent name.
Definition: AstrometryModel.h:116
std::map
STL class.
lsst::jointcal
Definition: Associations.h:49
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LOGLS_DEBUG
#define LOGLS_DEBUG(logger, message)
Definition: Log.h:608
os
std::ostream * os
Definition: Schema.cc:746
lsst::jointcal::CcdImage::getVisit
VisitIdType getVisit() const
returns visit ID
Definition: CcdImage.h:148
lsst::afw::cameraGeom::FOCAL_PLANE
CameraSys const FOCAL_PLANE
Focal plane coordinates: Position on a 2-d planar approximation to the focal plane (x,...
Definition: CameraSys.cc:30
AstrometryModel.h
std::endl
T endl(T... args)
lsst::jointcal::CcdImage::getCcdId
CcdIdType getCcdId() const
returns ccd ID
Definition: CcdImage.h:145
lsst::jointcal::ConstrainedAstrometryModel::makeSkyWcs
std::shared_ptr< afw::geom::SkyWcs > makeSkyWcs(CcdImage const &ccdImage) const override
Make a SkyWcs that contains this model.
Definition: ConstrainedAstrometryModel.cc:232
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
lsst::afw::geom::makeCdMatrix
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:138
lsst::geom::Point< double, 2 >
lsst::jointcal::ConstrainedAstrometryModel::getSkyToTangentPlane
const std::shared_ptr< AstrometryTransform const > getSkyToTangentPlane(CcdImage const &ccdImage) const override
The mapping of sky coordinates (i.e.
Definition: ConstrainedAstrometryModel.h:113
lsst::jointcal::AstrometryMapping
virtual class needed in the abstraction of the distortion model
Definition: AstrometryMapping.h:39
lsst::jointcal::Point
A point in a plane.
Definition: Point.h:36
lsst::afw::geom::makeSkyWcs
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:526
lsst::jointcal::CcdImage::getDetector
std::shared_ptr< afw::cameraGeom::Detector > getDetector() const
Definition: CcdImage.h:150
ConstrainedAstrometryModel.h
std::stringstream::str
T str(T... args)
lsst::pex::exceptions
Definition: Exception.h:37
lsst::jointcal::AstrometryTransformPolynomial
Polynomial transformation class.
Definition: AstrometryTransform.h:280
std::size_t
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::jointcal::AstrometryTransform
a virtual (interface) class for geometric transformations.
Definition: AstrometryTransform.h:65
lsst::jointcal::Point::x
double x
coordinate
Definition: Point.h:41
lsst::jointcal::VisitIdType
int VisitIdType
Definition: CcdImage.h:48
lsst.pipe.tasks.assembleCoadd.visit
visit
Definition: assembleCoadd.py:336
lsst::jointcal::AstrometryTransform::toAstMap
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.
Definition: AstrometryTransform.h:189
lsst::jointcal::ConstrainedAstrometryModel::offsetParams
void offsetParams(Eigen::VectorXd const &Delta) override
Dispaches the offsets after a fit step into the actual locations of parameters.
Definition: ConstrainedAstrometryModel.cc:171
lsst::afw::cameraGeom::PIXELS
CameraSysPrefix const PIXELS
Pixel coordinates: Nominal position on the entry surface of a given detector (x, y unbinned pixels).
Definition: CameraSys.cc:34
StarMatch.h
pairs of points
Log.h
LSST DM logging module built on log4cxx.
lsst::jointcal::ConstrainedAstrometryModel::print
void print(std::ostream &out) const override
Print a string representation of the contents of this mapping, for debugging.
Definition: ConstrainedAstrometryModel.cc:259
exceptions.h
lsst::jointcal::Frame
rectangle with sides parallel to axes.
Definition: Frame.h:38
lsst::jointcal::ConstrainedAstrometryModel::getChipTransform
AstrometryTransform const & getChipTransform(CcdIdType const chip) const
Access to mappings.
Definition: ConstrainedAstrometryModel.cc:189
std::pow
T pow(T... args)