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
fieldImpl.cc
Go to the documentation of this file.
1// -*- lsst-C++ -*-
2#include <cstring>
3#include "lsst/meas/extensions/psfex/Field.hh"
4#undef PI
5#include "lsst/geom/Point.h"
7
8extern "C" {
9 #include "globals.h"
10}
11
12namespace lsst { namespace meas { namespace extensions { namespace psfex {
13
14Field::Field(std::string const& ident) :
15 impl(NULL), _isInitialized(false)
16{
17 QCALLOC(impl, fieldstruct, 1);
18 impl->next = 0;
19
20 strcpy(impl->catname, ident.c_str());
21 impl->rcatname = impl->catname;
22#if 0
23 strncpy(impl->rtcatname, impl->rcatname, sizeof(impl->rtcatname) - 1);
24 strncpy(impl->ident, "??", sizeof(impl->ident) - 1);
25#elif 1
26 if (!(impl->rcatname = strrchr(impl->catname, '/'))) {
27 impl->rcatname = impl->catname;
28 } else {
29 ++impl->rcatname;
30 }
31
32 strncpy(impl->rtcatname, impl->rcatname, sizeof(impl->rtcatname) - 1);
33 {
34 char *pstr=strrchr(impl->rtcatname, '.');
35 if (pstr) {
36 *pstr = '\0';
37 }
38 }
39
40 strncpy(impl->ident, "??", sizeof(impl->ident) - 1);
41#endif
42
43 impl->ndet = 0;
44 impl->psf = NULL;
45 impl->wcs = NULL;
46
47 _finalize();
48}
49
50Field::~Field()
51{
52 for (int i = 0; i != impl->next; ++i) {
53 free(impl->wcs[i]); // psfex's wcs isn't quite the same as ours ...
54 impl->wcs[i] = NULL; // ... so don't let psfex free it
55 }
56 field_end(impl);
57 impl = NULL;
58}
59
60/************************************************************************************************************/
61
62void
63Field::_finalize(bool force)
64{
65 if (force || !_isInitialized) {
67 _isInitialized = true;
68 }
69}
70
71/************************************************************************************************************/
72
74Field::getPsfs() const
75{
76 if (_psfs.empty()) {
77 _psfs.reserve(impl->next);
78 for (int i = 0; i != impl->next; ++i) {
79 _psfs.push_back(Psf(impl->psf[i]));
80 }
81 }
82
83 return _psfs;
84}
85
86void
87Field::addExt(lsst::afw::geom::SkyWcs const& wcs_,
88 int const naxis1, int const naxis2,
89 int const nobj)
90{
91 QREALLOC(impl->psf, psfstruct *, impl->next + 1);
92 impl->psf[impl->next] = 0;
93 QREALLOC(impl->wcs, wcsstruct *, impl->next + 1);
94 impl->wcs[impl->next] = 0;
95 /*
96 * We're going to fake psfex's wcsstruct object. We only need enough of it for field_locate
97 */
98 QMALLOC(impl->wcs[impl->next], wcsstruct, 1);
99 wcsstruct *wcs = impl->wcs[impl->next];
100
101 wcs->naxis = 2;
102 wcs->naxisn[0] = naxis1;
103 wcs->naxisn[1] = naxis2;
104
105 auto const crval = wcs_.getSkyOrigin();
106 // crpix using the FITS standard = pixel origin using the LSST standard + 1
107 auto const crpix = wcs_.getPixelOrigin() + geom::Extent2D(1, 1);
108 auto const cdMatrix = wcs_.getCdMatrix();
109 std::string const cunit("DEG");
110 auto metadata = wcs_.getFitsMetadata();
111 for (int i = 0; i != wcs->naxis; ++i) {
112 auto ifits = i + 1;
113 auto ctype = metadata->getAsString("CTYPE" + std::to_string(ifits));
114 strncpy(wcs->ctype[i], ctype.c_str(), ctype.size() + 1);
115 strncpy(wcs->cunit[i], cunit.c_str(), cunit.size() + 1);
116 wcs->crpix[i] = crpix[i];
117 wcs->crval[i] = crval[i].asDegrees();
118 wcs->cdelt[i] = 1.0; // scale is in the CD matrix (is this even needed?)
119 wcs->crder[i] = 0;
120 wcs->csyer[i] = 0;
121 }
122 for (int i = 0, k = 0; i < 2; ++i) {
123 for (int j = 0; j < 2; ++j, ++k) {
124 wcs->cd[k] = cdMatrix(i, j);
125 }
126 }
127 wcs->lng = 0;
128 wcs->lat = 1;
129 wcs->equinox = 2000;
130
131 auto center = wcs_.pixelToSky(geom::Point2D(0.5*naxis1, 0.5*naxis2));
132 wcs->wcsscalepos[0] = center.getLongitude().asDegrees();
133 wcs->wcsscalepos[1] = center.getLatitude().asDegrees();
134
135 double maxradius = 0.0; // Maximum distance to wcsscalepos
136 for (int x = 0; x <= 1; ++x) {
137 for (int y = 0; y <= 1; ++y) {
138 geom::Point2D point(x*naxis1, y*naxis2); // Corner
139 double const radius = center.separation(wcs_.pixelToSky(point)).asDegrees();
140 if (radius > maxradius) {
141 maxradius = radius;
142 }
143 }
144 }
145 wcs->wcsmaxradius = maxradius;
146
147 impl->ndet += nobj;
148
149 ++impl->next;
150}
151
152}}}}
table::PointKey< double > crpix
Definition OldWcs.cc:129
table::PointKey< double > crval
Definition OldWcs.cc:128
int y
Definition SpanSet.cc:48
T c_str(T... args)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition SkyWcs.h:117
Eigen::Matrix2d getCdMatrix(lsst::geom::Point2D const &pixel) const
Get the 2x2 CD matrix at the specified pixel position.
Definition SkyWcs.cc:193
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Definition SkyWcs.h:334
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
Definition SkyWcs.cc:182
lsst::geom::Point2D getPixelOrigin() const
Get the pixel origin, in pixels, using the LSST convention.
Definition SkyWcs.h:215
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
Definition SkyWcs.cc:213
Tag types used to declare specialized field types.
Definition misc.h:31
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
Definition Angle.h:176
Angle getLongitude() const noexcept
The longitude of this point.
Extent< double, 2 > Extent2D
Definition Extent.h:400
T reserve(T... args)
T strcpy(T... args)
T strncpy(T... args)
T strrchr(T... args)
T to_string(T... args)