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-32411: Make WarpedPsf delegate to source PSF's computeImage. #280

Merged
merged 4 commits into from
May 26, 2022
Merged

DM-32411: Make WarpedPsf delegate to source PSF's computeImage. #280

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6d3f65b
Minor cleanups in CoaddPsf.
TallJimbo May 2, 2022
d4a0a94
Make WarpedPsf delegate to source PSF's computeImage.
TallJimbo May 2, 2022
9474d70
Rewrite WarpedPsf's image bbox calculation to ensure odd dimensions.
TallJimbo May 23, 2022
4620c95
Round down to ensure odd size
arunkannawadi May 23, 2022
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
6 changes: 5 additions & 1 deletion src/CoaddPsf.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,6 +198,8 @@ geom::Box2I getOverallBBox(std::vector<std::shared_ptr<afw::image::Image<double>

// Read all the images from the Image Vector and add them to image

namespace {

void addToImage(std::shared_ptr<afw::image::Image<double>> image,
std::vector<std::shared_ptr<afw::image::Image<double>>> const &imgVector,
std::vector<double> const &weightVector) {
Expand All @@ -221,6 +223,8 @@ 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);
if (subcat.empty()) {
Expand All @@ -244,7 +248,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 expoures which contain our coordinate within their validPolygons.
// Get the subset of exposures which contain our coordinate within their validPolygons.
afw::table::ExposureCatalog subcat = _catalog.subsetContaining(ccdXY, _coaddWcs, true);
if (subcat.empty()) {
throw LSST_EXCEPT(
Expand Down
120 changes: 74 additions & 46 deletions src/WarpedPsf.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,14 +50,6 @@ namespace algorithms {

namespace {

inline double min4(double a, double b, double c, double d) {
return std::min(std::min(a, b), std::min(c, d));
}

inline double max4(double a, double b, double c, double d) {
return std::max(std::max(a, b), std::max(c, d));
}

// TODO: make this routine externally callable and more generic using templates
// (also useful in e.g. math/offsetImage.cc)
std::shared_ptr<afw::detection::Psf::Image> zeroPadImage(afw::detection::Psf::Image const &im, int xPad,
Expand All @@ -75,8 +67,6 @@ std::shared_ptr<afw::detection::Psf::Image> zeroPadImage(afw::detection::Psf::Im
}

geom::Box2I computeBBoxFromTransform(geom::Box2I const bbox, geom::AffineTransform const &t) {
static const int dst_padding = 0;

// This is the maximum scaling I can imagine we'd want - it's approximately what you'd
// get from trying to coadd 2"-pixel data (e.g. 2MASS) onto a 0.01"-pixel grid (e.g.
// from JWST). Anything beyond that is probably a bug elsewhere in the code, and
Expand All @@ -87,29 +77,25 @@ geom::Box2I computeBBoxFromTransform(geom::Box2I const bbox, geom::AffineTransfo
throw LSST_EXCEPT(pex::exceptions::RangeError, "Unexpectedly large transform passed to WarpedPsf");
}

// min/max coordinate values in input image
int in_xlo = bbox.getMinX();
int in_xhi = bbox.getMinX() + bbox.getWidth() - 1;
int in_ylo = bbox.getMinY();
int in_yhi = bbox.getMinY() + bbox.getHeight() - 1;

// corners of output image
geom::Point2D c00 = t(geom::Point2D(in_xlo, in_ylo));
geom::Point2D c01 = t(geom::Point2D(in_xlo, in_yhi));
geom::Point2D c10 = t(geom::Point2D(in_xhi, in_ylo));
geom::Point2D c11 = t(geom::Point2D(in_xhi, in_yhi));

//
// bounding box for output image
//
int out_xlo = floor(min4(c00.getX(), c01.getX(), c10.getX(), c11.getX())) - dst_padding;
int out_ylo = floor(min4(c00.getY(), c01.getY(), c10.getY(), c11.getY())) - dst_padding;
int out_xhi = ceil(max4(c00.getX(), c01.getX(), c10.getX(), c11.getX())) + dst_padding;
int out_yhi = ceil(max4(c00.getY(), c01.getY(), c10.getY(), c11.getY())) + dst_padding;
// Floating point version of input bbox (expanded to include entire pixels).
geom::Box2D in_box_fp(bbox);
// Floating point version of output bbox (to be populated as bbox of
// transformed points.
geom::Box2D out_box_fp;
auto const in_corners = in_box_fp.getCorners();
for (auto const & in_corner : in_corners) {
auto out_corner = t(in_corner);
out_box_fp.include(out_corner);
}

geom::Box2I ret = geom::Box2I(geom::Point2I(out_xlo, out_ylo),
geom::Extent2I(out_xhi - out_xlo + 1, out_yhi - out_ylo + 1));
return ret;
// We want to guarantee that the output bbox has odd dimensions, so instead
// of using the Box2I converting constructor directly, we round the center
// point of the floating point box and dilate by its half-dimensions.
geom::Extent2I out_half_dims = geom::floor(0.5*out_box_fp.getDimensions());
geom::Box2I out_box;
geom::Point2I out_center(out_box_fp.getCenter());
out_box.include(out_center);
return out_box.dilatedBy(out_half_dims);
}

/**
Expand Down Expand Up @@ -148,6 +134,52 @@ std::shared_ptr<afw::detection::Psf::Image> warpAffine(afw::detection::Psf::Imag
return ret;
}

// A helper struct for logic and intermediate results used by both
// doComputeKernelImage and doComputeBBox.
//
// This linearizes the transform from destination to source coordinates, then
// computes both a source position that corresponds to a rounded version of the
// destination position and a transform from source to destination that also
// undoes the rounding.
struct PreparedTransforms {
arunkannawadi marked this conversation as resolved.
Show resolved Hide resolved

static PreparedTransforms compute(
afw::detection::Psf const & src_psf,
afw::geom::TransformPoint2ToPoint2 const & distortion,
geom::PointD const & dest_position
) {
// Linearize the transform from destination coordinates to source
// coordinates to avoid expensive WCS calls.
geom::AffineTransform dest_to_src = afw::geom::linearizeTransform(
*distortion.inverted(),
dest_position
);
// Instead of calling computeKernelImage on the source PSFs, we want to
// call computeImage with the source coordinate version of an
// integer-valued destination coordinate; that gives the lower-level PSF
// model the responsibility of doing the shifting from integer to
// fractional source coordinates, and it may be able to do a better job of
// that than we could do here (by e.g. using an internal analytic or
// oversampled model). So here's that integer-valued destination position
// and the corresponding source position.
geom::Point2D rounded_dest_position(
std::round(dest_position.getX()),
std::round(dest_position.getY())
);
auto src_position = dest_to_src(rounded_dest_position);
// Now we want to warp by the inverse of src_to_dest, but we want the
// image to be centered at (0, 0), not rounded_dest_position, so we
// compose that trivial shift with the dest_to_src transform.
auto src_to_dest = geom::AffineTransform(dest_to_src.inverted());
auto unround = geom::AffineTransform(-geom::Extent2D(rounded_dest_position));
auto src_to_dest_unrounded = unround * src_to_dest;
return PreparedTransforms{src_position, src_to_dest_unrounded};
}

geom::Point2D src_position;
geom::AffineTransform src_to_dest_unrounded;
};

} // namespace

WarpedPsf::WarpedPsf(std::shared_ptr<afw::detection::Psf const> undistortedPsf,
Expand Down Expand Up @@ -201,14 +233,13 @@ std::shared_ptr<afw::detection::Psf> WarpedPsf::resized(int width, int height) c

std::shared_ptr<afw::detection::Psf::Image> WarpedPsf::doComputeKernelImage(
geom::Point2D const &position, afw::image::Color const &color) const {
geom::AffineTransform t = afw::geom::linearizeTransform(*_distortion->inverted(), position);
geom::Point2D tp = t(position);

std::shared_ptr<Image> im = _undistortedPsf->computeKernelImage(tp, color);

// Go to the warped coordinate system with 'p' at the origin
auto srcToDest = geom::AffineTransform(t.inverted().getLinear());
std::shared_ptr<afw::detection::Psf::Psf::Image> ret = warpAffine(*im, srcToDest, *_warpingControl);
auto prepared_transforms = PreparedTransforms::compute(*_undistortedPsf, *_distortion, position);
std::shared_ptr<Image> src_im = _undistortedPsf->computeImage(prepared_transforms.src_position, color);
std::shared_ptr<afw::detection::Psf::Psf::Image> ret = warpAffine(
*src_im,
prepared_transforms.src_to_dest_unrounded,
*_warpingControl
);

double normFactor = 1.0;
//
Expand All @@ -230,12 +261,9 @@ std::shared_ptr<afw::detection::Psf::Image> WarpedPsf::doComputeKernelImage(
}

geom::Box2I WarpedPsf::doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const {
geom::AffineTransform t = afw::geom::linearizeTransform(*_distortion->inverted(), position);
geom::Point2D tp = t(position);
geom::Box2I bboxUndistorted = _undistortedPsf->computeBBox(tp, color);
geom::Box2I ret =
computeBBoxFromTransform(bboxUndistorted, geom::AffineTransform(t.inverted().getLinear()));
return ret;
auto prepared_transforms = PreparedTransforms::compute(*_undistortedPsf, *_distortion, position);
auto src_bbox = _undistortedPsf->computeImageBBox(prepared_transforms.src_position, color);
return computeBBoxFromTransform(src_bbox, prepared_transforms.src_to_dest_unrounded);
}

namespace {
Expand Down