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
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"
6 #include "lsst/afw/geom/SkyWcs.h"
7 
8 extern "C" {
9  #include "globals.h"
10 }
11 
12 namespace lsst { namespace meas { namespace extensions { namespace psfex {
13 
14 Field::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 
50 Field::~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 
62 void
63 Field::_finalize(bool force)
64 {
65  if (force || !_isInitialized) {
66  field_init_finalize(impl);
67  _isInitialized = true;
68  }
69 }
70 
71 /************************************************************************************************************/
72 
74 Field::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 
86 void
87 Field::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 }}}}
double x
table::PointKey< double > crpix
Definition: OldWcs.cc:129
table::PointKey< double > crval
Definition: OldWcs.cc:128
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
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
Extent< double, 2 > Extent2D
Definition: Extent.h:400
A base class for image defects.
T reserve(T... args)
T strcpy(T... args)
T strncpy(T... args)
T strrchr(T... args)
T to_string(T... args)