Skip to content

DM-55220: add accessors for template provenance to CoaddPsf #479

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
TallJimbo merged 2 commits into from
Jun 30, 2026
Merged

DM-55220: add accessors for template provenance to CoaddPsf #479

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
54aa070
Fix bad const/shared_ptr non-usage in CoaddPsf.
TallJimbo Jun 16, 2026
d5a9214
Add accessors for input-image tract and patch on CoaddPsf.
TallJimbo Jun 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 28 additions & 11 deletions include/lsst/meas/algorithms/CoaddPsf.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @param[in] warpingKernelName Name of warping kernel
* @param[in] cacheSize Warping kernel cache size
*/
explicit CoaddPsf(afw::table::ExposureCatalog const& catalog, afw::geom::SkyWcs const& coaddWcs,
explicit CoaddPsf(afw::table::ExposureCatalog const& catalog,
std::shared_ptr<afw::geom::SkyWcs const> coaddWcs,
std::string const& weightFieldName = "weight",
std::string const& warpingKernelName = "lanczos3", int cacheSize = 10000);

Expand All @@ -90,12 +91,12 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @param[in] weightFieldName Field name that contains the weight of the exposure in the coadd;
* defaults to "weight".
*/
CoaddPsf(afw::table::ExposureCatalog const& catalog, afw::geom::SkyWcs const& coaddWcs,
CoaddPsf(afw::table::ExposureCatalog const& catalog, std::shared_ptr<afw::geom::SkyWcs const> coaddWcs,
CoaddPsfControl const& ctrl, std::string const& weightFieldName = "weight")
: CoaddPsf(catalog, coaddWcs, weightFieldName, ctrl.warpingKernelName, ctrl.cacheSize) {}

explicit CoaddPsf(afw::table::ExposureCatalog const& catalog, ///< Unpersisted catalog
afw::geom::SkyWcs const& coaddWcs, ///< WCS for the coadd
std::shared_ptr<afw::geom::SkyWcs const> coaddWcs, ///< WCS for the coadd
geom::Point2D const& averagePosition, ///< Default position for accessors
std::string const& warpingKernelName = "lanczos3", ///< Warping kernel name
int cacheSize = 10000 ///< Kernel cache size
Expand All @@ -115,7 +116,7 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
geom::Point2D getAveragePosition() const override { return _averagePosition; }

/// Return the Wcs of the coadd (defines the coordinate system of the Psf).
afw::geom::SkyWcs getCoaddWcs() { return _coaddWcs; }
std::shared_ptr<afw::geom::SkyWcs const> getCoaddWcs() const { return _coaddWcs; }

/// Return the number of component Psfs in this CoaddPsf
int getComponentCount() const;
Expand All @@ -127,7 +128,7 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @returns Corresponding Psf.
* @throws RangeError Index of component is out of range.
*/
std::shared_ptr<afw::detection::Psf const> getPsf(int index);
std::shared_ptr<afw::detection::Psf const> getPsf(int index) const;

/**
* Get the Wcs of the component image at index.
Expand All @@ -136,7 +137,7 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @returns Corresponding Wcs.
* @throws RangeError Index of component is out of range.
*/
afw::geom::SkyWcs getWcs(int index);
std::shared_ptr<afw::geom::SkyWcs const> getWcs(int index) const;

/**
* Get the weight of the component image at index.
Expand All @@ -145,7 +146,7 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @returns Corresponding weight.
* @throws RangeError Index of component is out of range.
*/
double getWeight(int index);
double getWeight(int index) const;

/**
* Get the exposure ID of the component image at index.
Expand All @@ -154,7 +155,7 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @returns Corresponding exposure ID.
* @throws RangeError Index of component is out of range.
*/
afw::table::RecordId getId(int index);
afw::table::RecordId getId(int index) const;

/**
* Get the bounding box (in component image Pixel coordinates) of the component image at index.
Expand All @@ -163,7 +164,7 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @returns Corresponding bounding box.
* @throws RangeError Index of component is out of range.
*/
geom::Box2I getBBox(int index);
geom::Box2I getBBox(int index) const;

/**
* Get the validPolygon (in component image Pixel coordinates) of the component image at index.
Expand All @@ -172,7 +173,21 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag
* @returns Corresponding validPolygon.
* @throws RangeError Index of component is out of range.
*/
std::shared_ptr<afw::geom::polygon::Polygon const> getValidPolygon(int index);
std::shared_ptr<afw::geom::polygon::Polygon const> getValidPolygon(int index) const;

/**
* Get the tract ID of the component image at the given index.
*
* @throws lsst::pex::exceptions::NotFoundError There is no tract column.
*/
int getTract(int index) const;

/**
* Get the patch ID of the component image at the given index.
*
* @throws lsst::pex::exceptions::NotFoundError There is no patch column.
*/
int getPatch(int index) const;
Comment on lines +183 to +190

@isullivan isullivan Jun 22, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see how getTract and getPatch actually get set. src/CoadPsf.cc initializes with a minimal schema, which does not include patch or tract columns, so it looks like the exposure catalog won't actually have tract or patch columns.

@TallJimbo TallJimbo Jun 23, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are a few different constructors for CoaddPsf, and you're right about most of them. But there's one - the one we actually use in GetTemplateTask - that provides more control over the schema.

From the perspective of CoaddPsf, this usage was not expected or intentional. But it happens to be where we've stuffed a lot of useful information in all of the DP2 files. I am looking forward to being able to replace it with something more intentional in the future when we make GetTemplateTask write an lsst.images output directly.


/**
* @brief Return true if the CoaddPsf persistable (always true).
Expand Down Expand Up @@ -206,8 +221,10 @@ class CoaddPsf : public afw::table::io::PersistableFacade<CoaddPsf>, public Imag

private:
afw::table::ExposureCatalog _catalog;
afw::geom::SkyWcs _coaddWcs;
std::shared_ptr<afw::geom::SkyWcs const> _coaddWcs;
afw::table::Key<double> _weightKey;
afw::table::Key<int> _tractKey;
afw::table::Key<int> _patchKey;
geom::Point2D _averagePosition;
std::string _warpingKernelName; // could be removed if we could get this from _warpingControl (#2949)
std::shared_ptr<afw::math::WarpingControl const> _warpingControl;
Expand Down
8 changes: 5 additions & 3 deletions python/lsst/meas/algorithms/coaddPsf/coaddPsf.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,15 +49,15 @@ void declareCoaddPsf(lsst::cpputils::python::WrapperCollection &wrappers) {

auto clsCoaddPsf = wrappers.wrapType(PyCoaddPsf(wrappers.module, "CoaddPsf"), [](auto &mod, auto &cls) {
/* Constructors */
cls.def(py::init<afw::table::ExposureCatalog const &, afw::geom::SkyWcs const &,
cls.def(py::init<afw::table::ExposureCatalog const &, std::shared_ptr<afw::geom::SkyWcs>,
std::string const &, std::string const &, int>(),
"catalog"_a, "coaddWcs"_a, "weightFieldName"_a = "weight",
"warpingKernelName"_a = "lanczos3", "cacheSize"_a = 10000);
cls.def(py::init<afw::table::ExposureCatalog const &, afw::geom::SkyWcs const &,
cls.def(py::init<afw::table::ExposureCatalog const &, std::shared_ptr<afw::geom::SkyWcs>,
geom::Point2D const &, std::string const &, int>(),
"catalog"_a, "coaddWcs"_a, "averagePosition"_a,
"warpingKernelName"_a = "lanczos3", "cacheSize"_a = 10000);
cls.def(py::init<afw::table::ExposureCatalog const &, afw::geom::SkyWcs const &,
cls.def(py::init<afw::table::ExposureCatalog const &, std::shared_ptr<afw::geom::SkyWcs>,
CoaddPsfControl const &, std::string const &>(),
"catalog"_a, "coaddWcs"_a, "ctrl"_a, "weightFieldName"_a = "weight");

Expand All @@ -72,6 +72,8 @@ void declareCoaddPsf(lsst::cpputils::python::WrapperCollection &wrappers) {
cls.def("getId", &CoaddPsf::getId);
cls.def("getBBox", &CoaddPsf::getBBox);
cls.def("getValidPolygon", &CoaddPsf::getValidPolygon);
cls.def("getTract", &CoaddPsf::getTract);
cls.def("getPatch", &CoaddPsf::getPatch);
cls.def("isPersistable", &CoaddPsf::isPersistable);
});

Expand Down
70 changes: 51 additions & 19 deletions src/CoaddPsf.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,8 @@ geom::Point2D computeAveragePosition(afw::table::ExposureCatalog const &catalog,

} // namespace

CoaddPsf::CoaddPsf(afw::table::ExposureCatalog const &catalog, afw::geom::SkyWcs const &coaddWcs,
CoaddPsf::CoaddPsf(afw::table::ExposureCatalog const &catalog,
std::shared_ptr<afw::geom::SkyWcs const> coaddWcs,
std::string const &weightFieldName, std::string const &warpingKernelName, int cacheSize)
: _coaddWcs(coaddWcs),
_warpingKernelName(warpingKernelName),
Expand All @@ -165,17 +166,29 @@ CoaddPsf::CoaddPsf(afw::table::ExposureCatalog const &catalog, afw::geom::SkyWcs
record->assign(*i, mapper);
_catalog.push_back(record);
}
_averagePosition = computeAveragePosition(_catalog, _coaddWcs, _weightKey);
_averagePosition = computeAveragePosition(_catalog, *_coaddWcs, _weightKey);
}

CoaddPsf::CoaddPsf(afw::table::ExposureCatalog const &catalog, afw::geom::SkyWcs const &coaddWcs,
CoaddPsf::CoaddPsf(afw::table::ExposureCatalog const &catalog,
std::shared_ptr<afw::geom::SkyWcs const> coaddWcs,
geom::Point2D const &averagePosition, std::string const &warpingKernelName, int cacheSize)
: _catalog(catalog),
_coaddWcs(coaddWcs),
_weightKey(_catalog.getSchema()["weight"]),
_tractKey(),
_patchKey(),
_averagePosition(averagePosition),
_warpingKernelName(warpingKernelName),
_warpingControl(new afw::math::WarpingControl(warpingKernelName, "", cacheSize)) {}
_warpingControl(new afw::math::WarpingControl(warpingKernelName, "", cacheSize)) {
try {
_tractKey = _catalog.getSchema()["tract"];
} catch (pex::exceptions::NotFoundError &) {
}
try {
_patchKey = _catalog.getSchema()["patch"];
} catch (pex::exceptions::NotFoundError &) {
}
}

std::shared_ptr<afw::detection::Psf> CoaddPsf::clone() const { return std::make_shared<CoaddPsf>(*this); }

Expand Down Expand Up @@ -227,7 +240,7 @@ void addToImage(std::shared_ptr<afw::image::Image<double>> image,
} // anonymous

geom::Box2I CoaddPsf::doComputeBBox(geom::Point2D const &ccdXY, afw::image::Color const &color) const {
afw::table::ExposureCatalog subcat = _catalog.subsetContaining(ccdXY, _coaddWcs, true);
afw::table::ExposureCatalog subcat = _catalog.subsetContaining(ccdXY, *_coaddWcs, true);
if (subcat.empty()) {
throw LSST_EXCEPT(
lsst::afw::detection::InvalidPsfError,
Expand All @@ -238,7 +251,7 @@ geom::Box2I CoaddPsf::doComputeBBox(geom::Point2D const &ccdXY, afw::image::Colo
geom::Box2I ret;
for (auto const &exposureRecord : subcat) {
// compute transform from exposure pixels to coadd pixels
auto exposureToCoadd = afw::geom::makeWcsPairTransform(*exposureRecord.getWcs(), _coaddWcs);
auto exposureToCoadd = afw::geom::makeWcsPairTransform(*exposureRecord.getWcs(), *_coaddWcs);
WarpedPsf warpedPsf = WarpedPsf(exposureRecord.getPsf(), exposureToCoadd, _warpingControl);
geom::Box2I componentBBox = warpedPsf.computeBBox(ccdXY, color);
ret.include(componentBBox);
Expand All @@ -250,7 +263,7 @@ geom::Box2I CoaddPsf::doComputeBBox(geom::Point2D const &ccdXY, afw::image::Colo
std::shared_ptr<afw::detection::Psf::Image>
CoaddPsf::doComputeKernelImage(geom::Point2D const &ccdXY, afw::image::Color const &color) const {
// Get the subset of exposures which contain our coordinate within their validPolygons.
afw::table::ExposureCatalog subcat = _catalog.subsetContaining(ccdXY, _coaddWcs, true);
afw::table::ExposureCatalog subcat = _catalog.subsetContaining(ccdXY, *_coaddWcs, true);
if (subcat.empty()) {
throw LSST_EXCEPT(
lsst::afw::detection::InvalidPsfError,
Expand All @@ -267,7 +280,7 @@ CoaddPsf::doComputeKernelImage(geom::Point2D const &ccdXY, afw::image::Color con

for (auto const &exposureRecord : subcat) {
// compute transform from exposure pixels to coadd pixels
auto exposureToCoadd = afw::geom::makeWcsPairTransform(*exposureRecord.getWcs(), _coaddWcs);
auto exposureToCoadd = afw::geom::makeWcsPairTransform(*exposureRecord.getWcs(), *_coaddWcs);
std::shared_ptr<afw::image::Image<double>> componentImg;
try {
WarpedPsf warpedPsf = WarpedPsf(exposureRecord.getPsf(), exposureToCoadd, _warpingControl);
Expand Down Expand Up @@ -295,48 +308,68 @@ CoaddPsf::doComputeKernelImage(geom::Point2D const &ccdXY, afw::image::Color con

int CoaddPsf::getComponentCount() const { return _catalog.size(); }

std::shared_ptr<afw::detection::Psf const> CoaddPsf::getPsf(int index) {
std::shared_ptr<afw::detection::Psf const> CoaddPsf::getPsf(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
return _catalog[index].getPsf();
}

afw::geom::SkyWcs CoaddPsf::getWcs(int index) {
std::shared_ptr<afw::geom::SkyWcs const> CoaddPsf::getWcs(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
return *_catalog[index].getWcs();
return _catalog[index].getWcs();
}

std::shared_ptr<afw::geom::polygon::Polygon const> CoaddPsf::getValidPolygon(int index) {
std::shared_ptr<afw::geom::polygon::Polygon const> CoaddPsf::getValidPolygon(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
return _catalog[index].getValidPolygon();
}

double CoaddPsf::getWeight(int index) {
double CoaddPsf::getWeight(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
return _catalog[index].get(_weightKey);
}

afw::table::RecordId CoaddPsf::getId(int index) {
afw::table::RecordId CoaddPsf::getId(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
return _catalog[index].getId();
}

geom::Box2I CoaddPsf::getBBox(int index) {
geom::Box2I CoaddPsf::getBBox(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
return _catalog[index].getBBox();
}

int CoaddPsf::getTract(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
if (!_tractKey.isValid()) {
throw LSST_EXCEPT(pex::exceptions::NotFoundError, "CoaddPsf has no tract column");
}
return _catalog[index][_tractKey];
}

int CoaddPsf::getPatch(int index) const {
if (index < 0 || index >= getComponentCount()) {
throw LSST_EXCEPT(pex::exceptions::RangeError, "index of CoaddPsf component out of range");
}
if (!_patchKey.isValid()) {
throw LSST_EXCEPT(pex::exceptions::NotFoundError, "CoaddPsf has no patch column");
}
return _catalog[index][_patchKey];
}

// ---------- Persistence -----------------------------------------------------------------------------------

// For persistence of CoaddPsf, we have two catalogs: the first has just one record, and contains
Expand Down Expand Up @@ -388,7 +421,7 @@ class CoaddPsf::Factory : public afw::table::io::PersistableFactory {
afw::table::BaseRecord const &record1 = catalogs.front().front();
return std::shared_ptr<CoaddPsf>(
new CoaddPsf(afw::table::ExposureCatalog::readFromArchive(archive, catalogs.back()),
*archive.get<afw::geom::SkyWcs>(record1.get(keys1.coaddWcs)),
archive.get<afw::geom::SkyWcs>(record1.get(keys1.coaddWcs)),
record1.get(keys1.averagePosition), record1.get(keys1.warpingKernelName),
record1.get(keys1.cacheSize)));
}
Expand All @@ -413,7 +446,7 @@ class CoaddPsf::Factory : public afw::table::io::PersistableFactory {
} catch (pex::exceptions::NotFoundError &) {
}
auto averagePos = computeAveragePosition(internalCat, *coaddWcs, weightKey);
return std::shared_ptr<CoaddPsf>(new CoaddPsf(internalCat, *coaddWcs, averagePos));
return std::shared_ptr<CoaddPsf>(new CoaddPsf(internalCat, coaddWcs, averagePos));
}

Factory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
Expand All @@ -435,8 +468,7 @@ void CoaddPsf::write(OutputArchiveHandle &handle) const {
CoaddPsfPersistenceHelper const &keys1 = CoaddPsfPersistenceHelper::get();
afw::table::BaseCatalog cat1 = handle.makeCatalog(keys1.schema);
std::shared_ptr<afw::table::BaseRecord> record1 = cat1.addNew();
auto coaddWcsPtr = std::make_shared<afw::geom::SkyWcs>(_coaddWcs);
record1->set(keys1.coaddWcs, handle.put(coaddWcsPtr));
record1->set(keys1.coaddWcs, handle.put(_coaddWcs));
record1->set(keys1.cacheSize, _warpingControl->getCacheSize());
record1->set(keys1.averagePosition, _averagePosition);
record1->set(keys1.warpingKernelName, _warpingKernelName);
Expand Down
Loading