LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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 }}}}
std::strcpy
T strcpy(T... args)
y
int y
Definition: SpanSet.cc:49
std::strrchr
T strrchr(T... args)
std::string
STL class.
std::vector::reserve
T reserve(T... args)
wcs
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
std::vector
STL class.
SkyWcs.h
lsst::afw::geom::SkyWcs
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
crval
table::PointKey< double > crval
Definition: OldWcs.cc:130
x
double x
Definition: ChebyshevBoundedField.cc:277
crpix
table::PointKey< double > crpix
Definition: OldWcs.cc:131
std::string::c_str
T c_str(T... args)
std::to_string
T to_string(T... args)
lsst::afw::geom::SkyWcs::pixelToSky
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Definition: SkyWcs.h:334
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::afw::geom::SkyWcs::getCdMatrix
Eigen::Matrix2d getCdMatrix(lsst::geom::Point2D const &pixel) const
Get the 2x2 CD matrix at the specified pixel position.
Definition: SkyWcs.cc:198
std::strncpy
T strncpy(T... args)
lsst.pipe.tasks.insertFakes.radius
radius
Definition: insertFakes.py:288
lsst::geom::Point< double, 2 >
lsst::afw::geom::SkyWcs::getFitsMetadata
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
Definition: SkyWcs.cc:218
lsst::geom::Extent2D
Extent< double, 2 > Extent2D
Definition: Extent.h:400
lsst::afw::geom::SkyWcs::getPixelOrigin
lsst::geom::Point2D getPixelOrigin() const
Get the pixel origin, in pixels, using the LSST convention.
Definition: SkyWcs.h:215
Point.h
lsst::afw::geom::SkyWcs::getSkyOrigin
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
Definition: SkyWcs.cc:187