Skip to content
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

Jump to bottom

DM-38544: Allow getCutouts to extend off the edge of chips #698

Merged
merged 5 commits into from
Jul 18, 2023
Merged

DM-38544: Allow getCutouts to extend off the edge of chips #698

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c59dd30
Add pixel-space version of getCutout
parejkoj Jul 7, 2023
3b12cb9
clang-format changes
parejkoj Jul 7, 2023
cf58252
Fix typo
parejkoj Jul 7, 2023
b7b57ed
Box-based cutouts
parejkoj Jul 7, 2023
b650b91
Normalize getCutout docstrings
parejkoj Jul 17, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
50 changes: 42 additions & 8 deletions include/lsst/afw/image/Exposure.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ class Exposure {
* Return a subimage corresponding to the given box.
*
* @param bbox Bounding box of the subimage returned.
* @param origin Origin bbox is rleative to; PARENT accounts for xy0, LOCAL does not.
* @param origin Origin bbox is relative to; PARENT accounts for xy0, LOCAL does not.
* @return A subimage view into this.
*
* This method is wrapped as __getitem__ in Python.
Expand Down Expand Up @@ -294,9 +294,7 @@ class Exposure {
}

/// Set the Exposure's filter information
void setFilter(std::shared_ptr<FilterLabel const> filterLabel) {
_info->setFilter(filterLabel);
}
void setFilter(std::shared_ptr<FilterLabel const> filterLabel) { _info->setFilter(filterLabel); }

/// Set the Exposure's PhotoCalib object
void setPhotoCalib(std::shared_ptr<PhotoCalib const> photoCalib) { _info->setPhotoCalib(photoCalib); }
Expand Down Expand Up @@ -415,8 +413,8 @@ class Exposure {
/**
* Return an Exposure that is a small cutout of the original.
*
* @param center desired center of cutout (in RA and Dec)
* @param size width and height (in that order) of cutout in pixels
* @param center Desired center of cutout in RA and Dec.
* @param size Width and height (in that order) of cutout in pixels
*
* @return An Exposure of the requested size centered on `center` to within
* half a pixel in either dimension. Pixels past the edge of the original
Expand All @@ -425,11 +423,47 @@ class Exposure {
*
* @throws lsst::pex::exceptions::LogicError Thrown if this Exposure does
* not have a WCS.
* @throws lsst::pex::exceptions::InvalidParameterError Thrown if ``center``
* falls outside this Exposure or if ``size`` is not a valid size.
* @throws lsst::pex::exceptions::InvalidParameterError Thrown if `center`
* falls outside this Exposure or if `size` is not a valid size.
*/
Exposure getCutout(lsst::geom::SpherePoint const& center, lsst::geom::Extent2I const& size) const;

/**
* Return an Exposure that is a small cutout of the original.
*
* This is distinguished from subset() by allowing the cutout to extend off the image, and
* padding with NaN and mask NO_DATA in that case.
*
* @param center Desired center of cutout in image pixels.
* @param size Width and height (in that order) of cutout in pixels.
*
* @return An Exposure of the requested size centered on `center` to within
* half a pixel in either dimension. Pixels past the edge of the original
* exposure will equal @ref math::edgePixel(image::detail::MaskedImage_tag)
* "math::edgePixel<MaskedImageT>".
*
* @throws lsst::pex::exceptions::InvalidParameterError Thrown if `center`
* falls outside this Exposure or if `size` is not a valid size.
*/
Exposure getCutout(lsst::geom::Point2D const& center, lsst::geom::Extent2I const& size) const;

/**
* Return an Exposure that is a small cutout of the original.
*
* This is distinguished from subset() by allowing the box to extend off the image, and padding with
* NaN and mask NO_DATA in that case.
parejkoj marked this conversation as resolved.
Show resolved Hide resolved
*
* @param box Pixel box to cut from this exposure.
*
* @return An Exposure of the requested box to within half a pixel in either dimension. Pixels past the
* edge of the original exposure will equal @ref math::edgePixel(image::detail::MaskedImage_tag)
* "math::edgePixel<MaskedImageT>", with their Mask set to NO_DATA.
*
* @throws lsst::pex::exceptions::InvalidParameterError Thrown if the center of the box
* falls outside this Exposure.
*/
Exposure getCutout(lsst::geom::Box2I const& box) const;
Copy link
Member

Choose a reason for hiding this comment

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

The Box2I overload is not documented. Did you mean to use @{ and @}? (Possibly irrelevant, given my comment below.)

Copy link
Member

Choose a reason for hiding this comment

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

The Box2D overload worries me, given the unintuitive way conversions between Box2I and Box2D work. Is there a need for it? At the very least, I think you'd need to document how the conversion from Box2D to Box2I is done (I think your implementation is equivalent to EXPAND?).

Copy link
Member

Choose a reason for hiding this comment

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

For me the question is whether there is something extra getCutout can do with a Box2D that it couldn't do with a Box2I (e.g. set some metadata field on return). If it's just a convenience for some way of getting a Box2I from a Box2D and then calling getCutout, then I don't think we want the overload.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thanks, yes, I meant to wrap it in @{ @}.

I added the Box2D version because that seemed like the first thing a user would want (cutout of a generic box), while the Box2I version is what you'd do with e.g. footprint.getBBox(). I figured it would be helpful for users to not have to remember to do geom.ceil(box2d), as I think it's challenging for users to remember how to convert between our various primitives. I could add the box2d version as just a python (or pybind11) interface, if that would help?

Copy link
Member

Choose a reason for hiding this comment

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

I could add the box2d version as just a python (or pybind11) interface, if that would help?

I don't think it would. The issue is that the correct conversion from Box2D to Box2I is context-dependent no matter which language you call it in. I assume Box2I(Box2D, EdgeHandlingEnum) is accessible from Python (possibly with a string replacing the enum)? If so, requiring that users call it themselves might be the least of the evils.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I've removed the Box2D overload. We'll see if any users want that functionality.


private:
void _readFits(fits::Fits& fitsfile, lsst::geom::Box2I const& bbox, ImageOrigin origin,
bool conformMasks);
Expand Down
16 changes: 13 additions & 3 deletions python/lsst/afw/image/_exposure/_exposure.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,9 +137,9 @@ PyExposure<PixelT> declareExposure(lsst::utils::python::WrapperCollection &wrapp

cls.def("subset", &ExposureT::subset, "bbox"_a, "origin"_a = PARENT);

cls.def("writeFits", (void (ExposureT::*)(std::string const &) const) & ExposureT::writeFits);
cls.def("writeFits", (void(ExposureT::*)(std::string const &) const) & ExposureT::writeFits);
cls.def("writeFits",
(void (ExposureT::*)(fits::MemFileManager &) const) & ExposureT::writeFits);
(void(ExposureT::*)(fits::MemFileManager &) const) & ExposureT::writeFits);
cls.def("writeFits", [](ExposureT &self, fits::Fits &fits) { self.writeFits(fits); });

cls.def(
Expand Down Expand Up @@ -172,7 +172,17 @@ PyExposure<PixelT> declareExposure(lsst::utils::python::WrapperCollection &wrapp
cls.def_static("readFits", (ExposureT(*)(std::string const &))ExposureT::readFits);
cls.def_static("readFits", (ExposureT(*)(fits::MemFileManager &))ExposureT::readFits);

cls.def("getCutout", &ExposureT::getCutout, "center"_a, "size"_a);
cls.def("getCutout",
py::overload_cast<lsst::geom::SpherePoint const &, lsst::geom::Extent2I const &>(
&ExposureT::getCutout, py::const_),
"center"_a, "size"_a);
cls.def("getCutout",
py::overload_cast<lsst::geom::Point2D const &, lsst::geom::Extent2I const &>(
&ExposureT::getCutout, py::const_),
"center"_a, "size"_a);
cls.def("getCutout",
py::overload_cast<lsst::geom::Box2I const &>(&ExposureT::getCutout, py::const_),
"box"_a);
});
}
} // namespace
Expand Down
18 changes: 15 additions & 3 deletions src/image/Exposure.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,9 +208,15 @@ Exposure<ImageT, MaskT, VarianceT> Exposure<ImageT, MaskT, VarianceT>::getCutout
}
lsst::geom::Point2D pixelCenter = getWcs()->skyToPixel(center);

if (!lsst::geom::Box2D(getBBox()).contains(pixelCenter)) {
return getCutout(pixelCenter, size);
}

template <typename ImageT, typename MaskT, typename VarianceT>
Exposure<ImageT, MaskT, VarianceT> Exposure<ImageT, MaskT, VarianceT>::getCutout(
lsst::geom::Point2D const &center, lsst::geom::Extent2I const &size) const {
if (!lsst::geom::Box2D(getBBox()).contains(center)) {
std::stringstream buffer;
buffer << "Point " << center << " lies at pixel " << pixelCenter << ", which lies outside Exposure "
buffer << "Point " << center << " lies at pixel " << center << ", which lies outside Exposure "
<< getBBox();
throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, buffer.str());
}
Expand All @@ -219,7 +225,7 @@ Exposure<ImageT, MaskT, VarianceT> Exposure<ImageT, MaskT, VarianceT>::getCutout
buffer << "Cannot create bounding box with dimensions " << size;
throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, buffer.str());
}
lsst::geom::Box2I bbox = lsst::geom::Box2I::makeCenteredBox(pixelCenter, size);
lsst::geom::Box2I bbox = lsst::geom::Box2I::makeCenteredBox(center, size);

// cutout must have independent ExposureInfo
auto copyInfo = std::make_shared<ExposureInfo>(*getInfo());
Expand All @@ -232,6 +238,12 @@ Exposure<ImageT, MaskT, VarianceT> Exposure<ImageT, MaskT, VarianceT>::getCutout
return cutout;
}

template <typename ImageT, typename MaskT, typename VarianceT>
Exposure<ImageT, MaskT, VarianceT> Exposure<ImageT, MaskT, VarianceT>::getCutout(
lsst::geom::Box2I const &box) const {
return getCutout(box.getCenter(), box.getDimensions());
}

// Explicit instantiations
/// @cond
template class Exposure<std::uint16_t>;
Expand Down
41 changes: 40 additions & 1 deletion tests/test_exposure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -714,7 +714,10 @@ def testTicket2861(self):
exposure3 = afwImage.ExposureF(tmpFile)
self.assertIsNotNone(exposure3.getInfo().getCoaddInputs())

def testGetCutout(self):
def testGetCutoutSky(self):
"""Test we can get cutouts in sky coordinates, so long as there is a
valid WCS.
"""
wcs = self.smallExposure.getWcs()

dimensions = [lsst.geom.Extent2I(100, 50), lsst.geom.Extent2I(15, 15), lsst.geom.Extent2I(0, 10),
Expand Down Expand Up @@ -746,6 +749,42 @@ def testGetCutout(self):
with self.assertRaises(pexExcept.InvalidParameterError, msg=msg):
self.smallExposure.getCutout(cutoutCenter, cutoutSize)

def testGetCutoutPixel(self):
"""Test that we can get cutouts in pixel coordinates, even if the
extent is off the edge of the image, even if there is no WCS.
"""
dimensions = [lsst.geom.Extent2I(100, 50), lsst.geom.Extent2I(15, 15), lsst.geom.Extent2I(0, 10),
lsst.geom.Extent2I(25, 30), lsst.geom.Extent2I(15, -5),
2*self.exposureMiOnly.getDimensions()]
locations = [("center", lsst.geom.Box2D(self.exposureMiOnly.getBBox()).getCenter()),
("edge", lsst.geom.Point2D(0, 0)),
("rounding test", lsst.geom.Point2D(0.2, 0.7)),
("just inside", lsst.geom.Point2D(-0.5 + 1e-4, -0.5 + 1e-4)),
# These two should raise; center must be within image box.
("just outside", lsst.geom.Point2D(-0.5 - 1e-4, -0.5 - 1e-4)),
("outside", lsst.geom.Point2D(-1000, -1000))]
for cutoutSize in dimensions:
for label, cutoutCenter in locations:
msg = 'Cutout size = %s, location = %s' % (cutoutSize, label)
if "outside" not in label and all(cutoutSize.gt(0)):
cutout = self.exposureMiOnly.getCutout(cutoutCenter, cutoutSize)
self._checkCutoutPixels(
cutout,
self._getValidCorners(self.exposureMiOnly.getBBox(), cutout.getBBox()),
msg)

# Same result even if there is a wcs.
cutoutWithWcs = self.smallExposure.getCutout(cutoutCenter, cutoutSize)
self.assertMaskedImagesEqual(cutout.maskedImage, cutoutWithWcs.maskedImage)

# Getting a cutout with a bbox should produce the same result.
box = lsst.geom.Box2I.makeCenteredBox(cutoutCenter, lsst.geom.Extent2I(cutoutSize))
cutoutBox2I = self.exposureMiOnly.getCutout(box)
self.assertMaskedImagesEqual(cutout.maskedImage, cutoutBox2I.maskedImage)
else:
with self.assertRaises(pexExcept.InvalidParameterError, msg=msg):
self.exposureMiOnly.getCutout(cutoutCenter, cutoutSize)
parejkoj marked this conversation as resolved.
Show resolved Hide resolved

def testGetConvexPolygon(self):
"""Test the convex polygon."""
# Check that we do not have a convex polygon for the plain exposure.
Expand Down