33#include "boost/iterator/iterator_adaptor.hpp"
34#include "boost/iterator/transform_iterator.hpp"
35#include "ndarray/eigen.h"
72 explicit AvgPosItem(
double wx_ = 0.0,
double wy_ = 0.0,
double w_ = 0.0) :
wx(wx_),
wy(wy_),
w(w_) {}
78 bool operator<(AvgPosItem
const &other)
const {
return w < other.w; }
80 AvgPosItem &
operator+=(AvgPosItem
const &other) {
87 AvgPosItem &
operator-=(AvgPosItem
const &other) {
94 friend AvgPosItem
operator+(AvgPosItem a, AvgPosItem
const &
b) {
return a +=
b; }
96 friend AvgPosItem
operator-(AvgPosItem a, AvgPosItem
const &
b) {
return a -=
b; }
100 afw::geom::SkyWcs
const &
coaddWcs, afw::table::Key<double> weightKey) {
101 afw::table::Key<int> goodPixKey;
103 goodPixKey = catalog.getSchema()[
"goodpix"];
104 }
catch (pex::exceptions::NotFoundError &) {
110 AvgPosItem item(p.getX(), p.getY(), i->get(weightKey));
111 if (goodPixKey.isValid()) {
112 item.w *= i->get(goodPixKey);
116 items.push_back(item);
126 for (std::vector<AvgPosItem>::iterator iter =
items.begin();
128 if (iter ==
items.end()) {
133 "Could not find a valid average position for CoaddPsf");
153 mapper.addMapping(goodPixKey,
true);
160 _weightKey =
mapper.addMapping(weightKey, weightField);
165 record->assign(*i,
mapper);
168 _averagePosition = computeAveragePosition(_catalog, _coaddWcs, _weightKey);
175 _weightKey(_catalog.getSchema()[
"weight"]),
192 for (
unsigned int i = 0; i < imgVector.size(); i++) {
207 assert(imgVector.size() == weightVector.
size());
208 for (
unsigned int i = 0; i < imgVector.size(); i++) {
210 double weight = weightVector[i];
211 double sum = ndarray::asEigenMatrix(componentImg->getArray()).sum();
218 overlap.clip(
image->getBBox());
223 targetSubImage.scaledPlus(
weight / sum, cSubImage);
231 if (subcat.
empty()) {
234 (boost::format(
"Cannot compute BBox at point %s; no input images at that point.") % ccdXY)
239 for (
auto const &exposureRecord : subcat) {
242 WarpedPsf warpedPsf =
WarpedPsf(exposureRecord.getPsf(), exposureToCoadd, _warpingControl);
254 if (subcat.
empty()) {
257 (boost::format(
"Cannot compute CoaddPsf at point %s; no input images at that point.") % ccdXY)
260 double weightSum = 0.0;
268 for (
auto const &exposureRecord : subcat) {
273 WarpedPsf warpedPsf =
WarpedPsf(exposureRecord.getPsf(), exposureToCoadd, _warpingControl);
276 LSST_EXCEPT_ADD(exc, (boost::format(
"Computing WarpedPsf kernel image for id=%d") %
277 exposureRecord.getId())
282 weightSum += exposureRecord.get(_weightKey);
283 weightVector.
push_back(exposureRecord.get(_weightKey));
291 addToImage(
image, imgVector, weightVector);
302 return _catalog[index].getPsf();
309 return *_catalog[index].getWcs();
316 return _catalog[index].getValidPolygon();
323 return _catalog[index].
get(_weightKey);
330 return _catalog[index].getId();
337 return _catalog[index].getBBox();
349class CoaddPsfPersistenceHelper {
351 afw::table::Schema schema;
357 static CoaddPsfPersistenceHelper
const &get() {
358 static CoaddPsfPersistenceHelper
const instance;
363 CoaddPsfPersistenceHelper()
365 coaddWcs(schema.addField<
int>(
"coaddwcs",
"archive ID of the coadd's WCS")),
366 cacheSize(schema.addField<
int>(
"cachesize",
"size of the warping cache")),
368 schema,
"avgpos",
"PSF accessors default position",
"pixel")),
370 schema.addField<
std::string>(
"warpingkernelname",
"warping kernel name", 32)) {}
379 if (catalogs.size() == 1u) {
383 return readV0(archive, catalogs);
386 CoaddPsfPersistenceHelper
const &keys1 = CoaddPsfPersistenceHelper::get();
392 record1.
get(keys1.averagePosition), record1.
get(keys1.warpingKernelName),
393 record1.
get(keys1.cacheSize)));
405 auto coaddWcs = internalCat.back().getWcs();
406 internalCat.pop_back();
412 weightKey = internalCat.getSchema()[
"weight"];
415 auto averagePos = computeAveragePosition(internalCat, *
coaddWcs, weightKey);
424std::string getCoaddPsfPersistenceName() {
return "CoaddPsf"; }
426CoaddPsf::Factory registration(getCoaddPsfPersistenceName());
435 CoaddPsfPersistenceHelper
const &keys1 = CoaddPsfPersistenceHelper::get();
438 auto coaddWcsPtr = std::make_shared<afw::geom::SkyWcs>(_coaddWcs);
439 record1->set(keys1.coaddWcs, handle.
put(coaddWcsPtr));
440 record1->set(keys1.cacheSize, _warpingControl->getCacheSize());
441 record1->set(keys1.averagePosition, _averagePosition);
442 record1->set(keys1.warpingKernelName, _warpingKernelName);
std::vector< SchemaItem< Flag > > * items
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
An exception thrown when we have an invalid PSF.
lsst::geom::Box2I computeBBox(lsst::geom::Point2D position, image::Color color=image::Color()) const
Return the bounding box of the image returned by computeKernelImage()
std::shared_ptr< Image > computeKernelImage(lsst::geom::Point2D position, image::Color color=image::Color(), ImageOwnerEnum owner=COPY) const
Return an Image of the PSF, in a form suitable for convolution.
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Describe the colour of a source.
A class to represent a 2-dimensional array of pixels.
Tag types used to declare specialized field types.
Base class for all records.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
std::shared_ptr< RecordT > const get(size_type i) const
Return a pointer to the record at index i.
size_type size() const
Return the number of elements in the catalog.
bool empty() const
Return true if the catalog has no records.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
std::shared_ptr< Table > getTable() const
Return the table associated with the catalog.
void push_back(Record const &r)
Add a copy of the given record to the end of the catalog.
Custom catalog class for ExposureRecord/Table.
typename Base::const_iterator const_iterator
static ExposureCatalogT readFromArchive(io::InputArchive const &archive, BaseCatalog const &catalog)
Convenience input function for Persistables that contain an ExposureCatalog.
void writeToArchive(io::OutputArchiveHandle &handle, bool ignoreUnpersistable=true) const
Convenience output function for Persistables that contain an ExposureCatalog.
ExposureCatalogT subsetContaining(lsst::geom::SpherePoint const &coord, bool includeValidPolygon=false) const
Return a shallow subset of the catalog with only those records that contain the given point.
static Schema makeMinimalSchema()
Return a minimal schema for Exposure tables and records.
A mapping between the keys of two Schemas, used to copy data between them.
A vector of catalogs used by Persistable.
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.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
A base class for factory classes used to reconstruct objects from records.
An integer coordinate rectangle.
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
std::shared_ptr< afw::table::io::Persistable > readV0(InputArchive const &archive, CatalogVector const &catalogs) const
virtual std::shared_ptr< afw::table::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const
Construct a new object from the given InputArchive and vector of catalogs.
Factory(std::string const &name)
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy. Usually unnecessary, as Psfs are immutable.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
double getWeight(int index)
Get the weight of the component image at index.
std::shared_ptr< afw::detection::Psf const > getPsf(int index)
Get the Psf of the component image at index.
afw::table::RecordId getId(int index)
Get the exposure ID of the component image at index.
geom::Box2I getBBox(int index)
Get the bounding box (in component image Pixel coordinates) of the component image at index.
std::shared_ptr< afw::detection::Psf::Image > doComputeKernelImage(geom::Point2D const &ccdXY, afw::image::Color const &color) const override
These virtual member functions are private, not protected, because we only want derived classes to im...
std::shared_ptr< afw::geom::polygon::Polygon const > getValidPolygon(int index)
Get the validPolygon (in component image Pixel coordinates) of the component image at index.
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
CoaddPsf(afw::table::ExposureCatalog const &catalog, afw::geom::SkyWcs const &coaddWcs, std::string const &weightFieldName="weight", std::string const &warpingKernelName="lanczos3", int cacheSize=10000)
Main constructors for CoaddPsf.
geom::Box2I doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const override
int getComponentCount() const
Return the number of component Psfs in this CoaddPsf.
afw::geom::SkyWcs getWcs(int index)
Get the Wcs of the component image at index.
A Psf class that maps an arbitrary Psf through a coordinate transformation.
Reports errors in the logical structure of the program.
Reports attempts to access elements using an invalid key.
Reports when the result of an operation cannot be represented by the destination type.
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
std::shared_ptr< Image< PixelT > > operator+(Image< PixelT > const &img, ImageSlice< PixelT > const &slc)
Overload operator+()
Image< LhsPixelT > & operator+=(Image< LhsPixelT > &lhs, Image< RhsPixelT > const &rhs)
Add lhs to Image rhs (i.e. pixel-by-pixel addition) where types are different.
Image< LhsPixelT > & operator-=(Image< LhsPixelT > &lhs, Image< RhsPixelT > const &rhs)
Subtract lhs from Image rhs (i.e. pixel-by-pixel subtraction) where types are different.
std::shared_ptr< Image< PixelT > > operator-(Image< PixelT > const &img, ImageSlice< PixelT > const &slc)
Overload operator-()
ExposureCatalogT< ExposureRecord > ExposureCatalog
Point< double, 2 > Point2D
geom::Box2I getOverallBBox(std::vector< std::shared_ptr< afw::image::Image< double > > > const &imgVector)
afw::table::Key< double > weight
afw::table::Key< int > coaddWcs
afw::table::PointKey< double > averagePosition
afw::table::Key< std::string > warpingKernelName
table::Key< int > cacheSize