LSST Applications g063fba187b+cac8b7c890,g0f08755f38+6aee506743,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+b4475c5878,g1dcb35cd9c+8f9bc1652e,g20f6ffc8e0+6aee506743,g217e2c1bcf+73dee94bd0,g28da252d5a+1f19c529b9,g2bbee38e9b+3f2625acfc,g2bc492864f+3f2625acfc,g3156d2b45e+6e55a43351,g32e5bea42b+1bb94961c2,g347aa1857d+3f2625acfc,g35bb328faa+a8ce1bb630,g3a166c0a6a+3f2625acfc,g3e281a1b8c+c5dd892a6c,g3e8969e208+a8ce1bb630,g414038480c+5927e1bc1e,g41af890bb2+8a9e676b2a,g7af13505b9+809c143d88,g80478fca09+6ef8b1810f,g82479be7b0+f568feb641,g858d7b2824+6aee506743,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+2903d499ea,gb58c049af0+d64f4d3760,gc28159a63d+3f2625acfc,gcab2d0539d+b12535109e,gcf0d15dbbd+46a3f46ba9,gda6a2b7d83+46a3f46ba9,gdaeeff99f8+1711a396fd,ge79ae78c31+3f2625acfc,gef2f8181fd+0a71e47438,gf0baf85859+c1f95f4921,gfa517265be+6aee506743,gfa999e8aa5+17cd334064,w.2024.51
LSST Data Management Base Package
Loading...
Searching...
No Matches
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#include <utility>
28
29#include "astshim.h"
30#include "lsst/log/Log.h"
36#include "lsst/pex/exceptions.h"
38
39namespace lsst {
40namespace jointcal {
41
44 bool initFromWcs, unsigned nNotFit, unsigned order)
45 : AstrometryModel(LOG_GET("lsst.jointcal.SimpleAstrometryModel")),
46 _skyToTangentPlane(std::move(projectionHandler))
47
48{
50
51 for (auto i = ccdImageList.cbegin(); i != ccdImageList.cend(); ++i, ++count) {
52 const CcdImage &im = **i;
53 if (count < nNotFit) {
56 id->setIndex(-1); // non sense, because it has no parameters
57 _myMap[im.getHashKey()] = std::move(id);
58 } else
59 // Given how AssignIndices works, only the SimplePolyMappings
60 // will actually be fitted, as nNotFit requests.
61 {
62 // first check that there are enough measurements for the requested polynomial order.
63 size_t nObj = im.getCatalogForFit().size();
64 if (nObj == 0) {
65 LOGLS_WARN(_log, "Empty catalog from image: " << im.getName());
66 continue;
67 }
69 if (pol.getOrder() > 0) // if not, it cannot be decreased
70 {
71 while (pol.getNpar() > 2 * nObj) {
72 LOGLS_WARN(_log, "Reducing polynomial order from "
73 << pol.getOrder() << ", due to too few sources (" << nObj
74 << " vs. " << pol.getNpar() << " parameters)");
75 pol.setOrder(pol.getOrder() - 1);
76 }
77 }
78 /* We have to center and normalize the coordinates so that
79 the fit matrix is not too ill-conditionned. Basically, x
80 and y in pixels are mapped to [-1,1]. When the
81 transformation of SimplePolyMapping transformation is
82 accessed, the combination of the normalization and the
83 fitted transformation is returned, so that the trick
84 remains hidden
85 */
86 const Frame &frame = im.getImageFrame();
88 if (initFromWcs) {
90 pol = pol * shiftAndNormalize.inverted();
91 }
92 _myMap[im.getHashKey()] =
94 }
95 }
96}
97
99 return findMapping(ccdImage);
100}
101
102Eigen::Index SimpleAstrometryModel::assignIndices(std::string const &whatToFit, Eigen::Index firstIndex) {
103 if (whatToFit.find("Distortions") == std::string::npos) {
104 LOGLS_ERROR(_log, "AssignIndices was called and Distortions is *not* in whatToFit.");
105 return 0;
106 }
107 Eigen::Index index = firstIndex;
108 for (auto & i : _myMap) {
109 auto *p = dynamic_cast<SimplePolyMapping *>(&*(i.second));
110 if (!p) continue; // it should be AstrometryTransformIdentity
111 p->setIndex(index);
112 index += p->getNpar();
113 }
114 return index;
115}
116
117void SimpleAstrometryModel::offsetParams(Eigen::VectorXd const &delta) {
118 for (auto &i : _myMap) {
119 auto mapping = i.second.get();
120 mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
121 }
122}
123
125 for (auto &i : _myMap) i.second->freezeErrorTransform();
126}
127
129 std::size_t total = 0;
130 for (auto &i : _myMap) {
131 total += i.second->getNpar();
132 }
133 return total;
134}
135
137 out << "SimpleAstrometryModel: " << _myMap.size() << " mappings" << std::endl;
138 out << *_skyToTangentPlane << std::endl;
139 out << "Sensor to sky transforms:" << std::endl;
140 for (auto &i : _myMap) {
141 out << i.first << std::endl;
142 out << *(i.second) << std::endl;
143 out << std::endl;
144 }
145}
146
148 return dynamic_cast<const SimplePolyMapping *>(findMapping(ccdImage))->getTransform();
149}
150
152 auto proj = std::dynamic_pointer_cast<const TanRaDecToPixel>(getSkyToTangentPlane(ccdImage));
153 jointcal::Point tangentPoint(proj->getTangentPoint());
154
155 auto polyMap = getTransform(ccdImage).toAstMap(ccdImage.getImageFrame());
156 ast::Frame pixelFrame(2, "Domain=PIXELS");
157 ast::Frame iwcFrame(2, "Domain=IWC");
158
159 // make a basic SkyWcs and extract the IWC portion
160 auto iwcToSkyWcs = afw::geom::makeSkyWcs(
161 geom::Point2D(0, 0), geom::SpherePoint(tangentPoint.x, tangentPoint.y, geom::degrees),
162 afw::geom::makeCdMatrix(1.0 * geom::degrees, 0 * geom::degrees, true));
163 auto iwcToSkyMap = iwcToSkyWcs->getFrameDict()->getMapping("PIXELS", "SKY");
164 auto skyFrame = iwcToSkyWcs->getFrameDict()->getFrame("SKY");
165
166 ast::FrameDict frameDict(pixelFrame, *polyMap, iwcFrame);
167 frameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
168 return std::make_shared<afw::geom::SkyWcs>(frameDict);
169}
170
172 auto i = _myMap.find(ccdImage.getHashKey());
173 if (i == _myMap.end())
175 "SimpleAstrometryModel cannot find CcdImage " + ccdImage.getName());
176 return i->second.get();
177}
178
179} // namespace jointcal
180} // 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
Class for a simple mapping implementing a generic AstrometryTransform.
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.
SimpleAstrometryModel(CcdImageList const &ccdImageList, std::shared_ptr< ProjectionHandler const > projectionHandler, bool initFromWCS, unsigned nNotFit=0, unsigned order=3)
const std::shared_ptr< AstrometryTransform const > getSkyToTangentPlane(CcdImage const &ccdImage) const override
the mapping of sky coordinates (i.e.
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.
AstrometryMapping * findMapping(CcdImage const &ccdImage) const override
Return a pointer to the mapping associated with this ccdImage.
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.
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
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].
STL namespace.
T size(T... args)
table::Key< int > order