52template std::shared_ptr<meas::extensions::psfex::PsfexPsf>
62namespace afw = lsst::afw;
65 lsst::meas::extensions::psfex::Psf
const& psf,
67 ) :
ImagePsf(), _averagePosition(averagePosition),
69 _comp(psf.impl->npix),
70 _context(psf.impl->poly->ndim)
73 _poly = poly_copy(psf.impl->poly);
75 _pixstep = psf.impl->pixstep;
77 std::copy(psf.impl->size, psf.impl->size + psf.impl->dim, _size.begin());
79 std::copy(psf.impl->comp, psf.impl->comp + psf.impl->npix, _comp.begin());
81 for (
int i = 0; i != psf.impl->poly->ndim; ++i) {
82 _context[i].first = psf.impl->contextoffset[i];
83 _context[i].second = psf.impl->contextscale[i];
88 _averagePosition(
geom::Point2I(0, 0)),
114 return _context.size();
120 double pos[MAXCONTEXT];
121 int const ndim = _context.size();
124 str(boost::format(
"Only spatial variation (ndim == 2) is supported; saw %d")
130 position = _averagePosition;
133 for (
int i = 0; i < ndim; ++i) {
134 pos[i] = (position[i] - _context[i].first)/_context[i].second;
137 poly_func(_poly, pos);
139 int const w = _size[0], h = _size[1];
144 const int nbasis = _size.size() > 2 ? _size[2] : 1;
148 float const vigstep = 1/_pixstep;
149 float const dx = 0.0, dy = 0.0;
159 for (
int i = 0; i != nbasis; ++i) {
163 vignet_resample(
const_cast<float *
>(&_comp[i*w*h]), w, h,
164 &sampledBasis[0], sampleW, sampleH,
165 -dx*vigstep, -dy*vigstep, vigstep, 1.0);
170 float *pl = &sampledBasis[0];
171 for (
int y = 0; y != sampleH; ++y) {
172 for (
int x = 0; x != sampleW; ++x) {
207 int const w = _size[0], h = _size[1];
208 int sampleW =
static_cast<int>(w*_pixstep);
209 int sampleH =
static_cast<int>(h*_pixstep);
212 if (sampleW % 2 == 0) sampleW += 1;
213 if (sampleH % 2 == 0) sampleH += 1;
215 float dx = center[0] -
static_cast<int>(center[0]);
216 float dy = center[1] -
static_cast<int>(center[1]);
218 if (dx > 0.5) dx -= 1.0;
219 if (dy > 0.5) dy -= 1.0;
223 static_cast<int>(center[1] - dy + 0.5) - sampleH/2),
234 double pos[MAXCONTEXT];
235 int const ndim = _context.size();
238 str(boost::format(
"Only spatial variation (ndim == 2) is supported; saw %d")
243 for (
int i = 0; i < ndim; ++i) {
244 pos[i] = (position[i] - _context[i].first)/_context[i].second;
247 poly_func(_poly, pos);
249 int const w = _size[0], h = _size[1];
251 const int nbasis = _size.size() > 2 ? _size[2] : 1;
254 int const npix = w*h;
255 for (
int i = 0; i != nbasis; ++i) {
256 float *pl = &fullresIm[0];
257 float const fac = _poly->basis[i];
258 float const *ppc = &_comp[i*w*h];
260 for (
int j = 0; j != npix; ++j) {
268 float const vigstep = 1/_pixstep;
269 float dx = center[0] -
static_cast<int>(center[0]);
270 float dy = center[1] -
static_cast<int>(center[1]);
271 if (dx > 0.5) dx -= 1.0;
272 if (dy > 0.5) dy -= 1.0;
276 geom::Box2I bbox = _doComputeBBox(position, center);
284 vignet_resample(&fullresIm[0], w, h,
285 &sampledIm[0], sampleW, sampleH,
286 -dx*vigstep, -dy*vigstep, vigstep, 1.0);
288 float *pl = &sampledIm[0];
289 float const sum =
std::accumulate(pl, pl + sampleW*sampleH,
static_cast<float>(0));
290 for (
int y = 0; y != sampleH; ++y) {
291 for (
int x = 0; x != sampleW; ++x) {
292 (*im)(x, y) = *pl++/sum;
308class PsfexPsfSchema1 {
312 ndim(schema.addField<int>(
"ndim",
"Number of elements in group")),
313 ngroup(schema.addField<int>(
"ngroup",
"Number of elements in degree")),
314 ncoeff(schema.addField<int>(
"ncoeff",
"Number of coefficients")),
316 _size_size(schema.addField<int>(
"_size_size",
"Size of _size array")),
317 _comp_size(schema.addField<int>(
"_comp_size",
"Size of _comp array")),
318 _context_size(schema.addField<int>(
"_context_size",
"Size of _context array")),
320 averagePosition(
afw::table::PointKey<double>::addFields(schema,
"averagePosition",
"average position of stars used to make the PSF",
"pixel")),
321 _pixstep(schema.addField<float>(
"_pixstep",
"oversampling",
"pixel"))
326 table::Schema schema;
329 table::Key<int> ndim;
330 table::Key<int> ngroup;
331 table::Key<int> ncoeff;
333 table::Key<int> _size_size;
334 table::Key<int> _comp_size;
335 table::Key<int> _context_size;
337 table::PointKey<double> averagePosition;
338 table::Key<float> _pixstep;
341class PsfexPsfSchema2 {
343 PsfexPsfSchema2(
int const ndim,
int const ngroup,
int const ncoeff,
344 int size_size,
int comp_size,
int context_size) :
347 group(schema.addField<table::Array<int> >(
"group",
"Groups (of coefficients?)", ndim)),
348 degree(schema.addField<table::Array<int> >(
"degree",
"Degree in each group", ngroup)),
349 basis(schema.addField<table::Array<double> >(
"basis",
"Values of the basis functions", ncoeff)),
350 coeff(schema.addField<table::Array<double> >(
"coeff",
"Polynomial coefficients", ncoeff)),
351 _size(schema.addField<table::Array<int> >(
"_size",
"PSF dimensions", size_size)),
352 _comp(schema.addField<table::Array<float> >(
"_comp",
"Complete pixel data", comp_size)),
353 _context_first( schema.addField<table::Array<double> >(
"_context_first",
354 "Offset to apply to context data", context_size)),
355 _context_second(schema.addField<table::Array<double> >(
"_context_second",
356 "Scale to apply to context data", context_size))
361 table::Schema schema;
363 table::Key<table::Array<int> > group;
364 table::Key<table::Array<int> > degree;
365 table::Key<table::Array<double> > basis;
366 table::Key<table::Array<double> > coeff;
368 table::Key<table::Array<int> > _size;
369 table::Key<table::Array<float> > _comp;
370 table::Key<table::Array<double> > _context_first;
371 table::Key<table::Array<double> > _context_second;
374std::string getPsfexPsfPersistenceName() {
return "PsfexPsf"; }
390 int ndim, ngroup, ncoeff;
391 int size_size, comp_size, context_size;
393 PsfexPsfSchema1
const & keys = PsfexPsfSchema1();
396 table::BaseRecord
const & record = catalogs[0].front();
399 ndim = record.get(keys.ndim);
400 ngroup = record.get(keys.ngroup);
401 ncoeff = record.get(keys.ncoeff);
403 result->_averagePosition = record.get(keys.averagePosition);
404 result->_pixstep = record.get(keys._pixstep);
406 size_size = record.get(keys._size_size);
407 comp_size = record.get(keys._comp_size);
408 context_size = record.get(keys._context_size);
412 PsfexPsfSchema2
const keys(ndim, ngroup, ncoeff,
413 size_size, comp_size, context_size);
417 table::BaseRecord
const & record = catalogs[1].front();
422 int const *begin = record.getElement(keys.group);
423 group.
assign(begin, begin + ndim);
425 for (
int i = 0; i != ndim; ++i) {
431 int const *begin = record.getElement(keys.degree);
432 degree.
assign(begin, begin + ngroup);
434 result->_poly = poly_init(&group[0], group.
size(), °ree[0], degree.
size());
438 double const *begin = record.getElement(keys.basis);
439 std::copy(begin, begin + ncoeff, result->_poly->basis);
442 double const *begin = record.getElement(keys.coeff);
443 std::copy(begin, begin + ncoeff, result->_poly->coeff);
447 int const *begin = record.getElement(keys._size);
448 result->_size.assign(begin, begin + size_size);
451 float const *begin = record.getElement(keys._comp);
452 result->_comp.assign(begin, begin + comp_size);
455 double const *begin1 = record.getElement(keys._context_first);
456 double const *begin2 = record.getElement(keys._context_second);
457 result->_context.resize(context_size);
458 for (
int i = 0; i != context_size; ++i) {
459 result->_context[i].first = begin1[i];
460 result->_context[i].second = begin2[i];
474 detail::PsfexPsfFactory registration(getPsfexPsfPersistenceName());
487 PsfexPsfSchema1
const keys;
492 record->set(keys.ndim, _poly->ndim);
493 record->set(keys.ngroup, _poly->ngroup);
494 record->set(keys.ncoeff, _poly->ncoeff);
496 record->set(keys.averagePosition, _averagePosition);
497 record->set(keys._pixstep, _pixstep);
498 record->set(keys._size_size, _size.size());
499 record->set(keys._comp_size, _comp.size());
500 record->set(keys._context_size, _context.size());
506 PsfexPsfSchema2
const keys(_poly->ndim, _poly->ngroup, _poly->ncoeff,
507 _size.size(), _comp.size(), _context.size());
512 int *begin = record->getElement(keys.group);
513 std::copy(_poly->group, _poly->group + _poly->ndim, begin);
516 int *begin = record->getElement(keys.degree);
517 std::copy(_poly->degree, _poly->degree + _poly->ngroup, begin);
520 double *begin = record->getElement(keys.basis);
521 std::copy(_poly->basis, _poly->basis + _poly->ncoeff, begin);
524 double *begin = record->getElement(keys.coeff);
525 std::copy(_poly->coeff, _poly->coeff + _poly->ncoeff, begin);
529 int *begin = record->getElement(keys._size);
530 std::copy(_size.begin(), _size.end(), begin);
533 float *begin = record->getElement(keys._comp);
534 std::copy(_comp.begin(), _comp.end(), begin);
537 double *begin1 = record->getElement(keys._context_first);
538 double *begin2 = record->getElement(keys._context_second);
539 for (
unsigned int i = 0; i != _context.size(); ++i) {
540 begin1[i] = _context[i].first;
541 begin2[i] = _context[i].second;
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Utilities for persisting KernelPsf and subclasses thereof.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
image::Image< Pixel > Image
Image type returned by computeImage.
Describe the colour of a source.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
io::CatalogVector CatalogVector
io::InputArchive InputArchive
PersistableFactory(std::string const &name)
Constructor for the factory.
An integer coordinate rectangle.
int getHeight() const noexcept
int getWidth() const noexcept
ImagePsf(bool isFixed=false)
Represent a PSF as a linear combination of PSFEX (== Karhunen-Loeve) basis functions.
std::shared_ptr< lsst::afw::math::LinearCombinationKernel const > getKernel(lsst::geom::Point2D=lsst::geom::Point2D(std::numeric_limits< double >::quiet_NaN())) const
Return the PSF's basis functions as a spatially-invariant LinearCombinationKernel with unit weights.
void write(lsst::afw::table::io::OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
lsst::geom::Box2I doComputeBBox(lsst::geom::Point2D const &position, lsst::afw::image::Color const &color) const override
Compute the bbox of the kernel image at the specified position/color.
PsfexPsf(lsst::meas::extensions::psfex::Psf const &psf, lsst::geom::Point2D const &averagePosition=lsst::geom::Point2D())
Constructor for a PsfexPsf.
std::shared_ptr< lsst::afw::detection::Psf::Image > doComputeImage(lsst::geom::Point2D const &position, lsst::afw::image::Color const &color) const override
Compute an image of the Psf at the specified position/colour, at pixel position in the output image.
int getNdim() const
Return the number of dependency parameters in the psfex polynomial fit.
std::shared_ptr< lsst::afw::detection::Psf > clone() const override
Polymorphic deep copy; should usually be unnecessary as Psfs are immutable.x.
virtual std::shared_ptr< lsst::afw::detection::Psf::Image > _doComputeImage(lsst::geom::Point2D const &position, lsst::afw::image::Color const &color, lsst::geom::Point2D const ¢er) const
Compute an image of the Psf at the specified position/colour, at pixel position in the output image.
std::shared_ptr< lsst::afw::detection::Psf::Image > doComputeKernelImage(lsst::geom::Point2D const &position, lsst::afw::image::Color const &color) const override
Compute an image of the Psf at the specified position/colour, at pixel (0.0, 0.0) in the output image...
std::string getPersistenceName() const override
Name used in table persistence.
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
std::string getPythonModule() const override
The python module name (for use in table persistence)
PsfexPsfFactory(std::string const &name)
virtual std::shared_ptr< table::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const
Construct a new object from the given InputArchive and vector of catalogs.
Reports invalid arguments.
Reports errors in the logical structure of the program.
std::vector< std::shared_ptr< Kernel > > KernelList
CatalogT< BaseRecord > BaseCatalog
Point< double, 2 > Point2D
Extent< int, 2 > Extent2I