LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
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)