Skip to content

Commit

Permalink
Use geom rather than afwGeom where appropriate
Browse files Browse the repository at this point in the history
  • Loading branch information
@timj
timj committed Jul 16, 2019
1 parent d0f904b commit 1ba1e73
Show file tree
Hide file tree
Showing 3 changed files with 44 additions and 42 deletions.
23 changes: 12 additions & 11 deletions include/lsst/meas/extensions/photometryKron.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include <cmath>

#include "lsst/pex/config.h"
#include "lsst/geom.h"
#include "lsst/afw/image/Exposure.h"
#include "lsst/afw/geom/SkyWcs.h"
#include "lsst/meas/base/Algorithm.h"
Expand Down Expand Up @@ -138,9 +139,9 @@ class KronFluxAlgorithm : public base::SimpleAlgorithm {
void _applyForced(
afw::table::SourceRecord & source,
afw::image::Exposure<float> const & exposure,
afw::geom::Point2D const & center,
geom::Point2D const & center,
afw::table::SourceRecord const & reference,
afw::geom::AffineTransform const & refToMeas
geom::AffineTransform const & refToMeas
) const;

PTR(KronAperture) _fallbackRadius(afw::table::SourceRecord& source, double const R_K_psf,
Expand All @@ -158,20 +159,20 @@ class KronFluxAlgorithm : public base::SimpleAlgorithm {

class KronAperture {
public:
KronAperture(afw::geom::Point2D const& center, afw::geom::ellipses::BaseCore const& core,
KronAperture(geom::Point2D const& center, afw::geom::ellipses::BaseCore const& core,
float radiusForRadius=std::nanf("")) :
_center(center),
_axes(core),
_radiusForRadius(radiusForRadius)
{}

explicit KronAperture(afw::table::SourceRecord const& source, float radiusForRadius=std::nanf("")) :
_center(afw::geom::Point2D(source.getX(), source.getY())),
_center(geom::Point2D(source.getX(), source.getY())),
_axes(source.getShape()),
_radiusForRadius(radiusForRadius)
{}

KronAperture(afw::table::SourceRecord const& reference, afw::geom::AffineTransform const& refToMeas,
KronAperture(afw::table::SourceRecord const& reference, geom::AffineTransform const& refToMeas,
double radius, float radiusForRadius=std::nanf("")) :
_center(refToMeas(reference.getCentroid())),
_axes(getKronAxes(reference.getShape(), refToMeas.getLinear(), radius)),
Expand All @@ -183,7 +184,7 @@ class KronAperture {
double getY() const { return _center.getY(); }
float getRadiusForRadius() const { return _radiusForRadius; }

afw::geom::Point2D const& getCenter() const { return _center; }
geom::Point2D const& getCenter() const { return _center; }

afw::geom::ellipses::Axes & getAxes() { return _axes; }

Expand All @@ -198,7 +199,7 @@ class KronAperture {
static PTR(KronAperture) determineRadius(
ImageT const& image, ///< Image to measure
afw::geom::ellipses::Axes axes, ///< Shape of aperture
afw::geom::Point2D const& center, ///< Centre of source
geom::Point2D const& center, ///< Centre of source
KronFluxControl const& ctrl ///< control the algorithm
);

Expand All @@ -211,8 +212,8 @@ class KronAperture {
) const;

/// Transform a Kron Aperture to a different frame
PTR(KronAperture) transform(afw::geom::AffineTransform const& trans) const {
afw::geom::Point2D const center = trans(getCenter());
PTR(KronAperture) transform(geom::AffineTransform const& trans) const {
geom::Point2D const center = trans(getCenter());
afw::geom::ellipses::Axes const axes(*getAxes().transform(trans.getLinear()).copy());
return std::make_shared<KronAperture>(center, axes);
}
Expand All @@ -221,12 +222,12 @@ class KronAperture {
static
afw::geom::ellipses::Axes getKronAxes(
afw::geom::ellipses::Axes const& shape,
afw::geom::LinearTransform const& transformation,
geom::LinearTransform const& transformation,
double const radius
);

private:
afw::geom::Point2D const _center; // Center of aperture
geom::Point2D const _center; // Center of aperture
afw::geom::ellipses::Axes _axes; // Ellipse defining aperture shape
float _radiusForRadius; // Radius used to estimate the Kron radius
};
Expand Down
38 changes: 19 additions & 19 deletions src/KronPhotometry.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,16 +26,16 @@
#include <functional>
#include "boost/math/constants/constants.hpp"
#include "lsst/pex/exceptions.h"
#include "lsst/afw/geom/Point.h"
#include "lsst/afw/geom/Box.h"
#include "lsst/geom/Point.h"
#include "lsst/geom/Box.h"
#include "lsst/afw/geom/SpanSet.h"
#include "lsst/afw/image/Exposure.h"
#include "lsst/afw/table/Source.h"
#include "lsst/afw/math/Integrate.h"
#include "lsst/afw/math/FunctionLibrary.h"
#include "lsst/afw/math/KernelFunctions.h"
#include "lsst/afw/detection/Psf.h"
#include "lsst/afw/geom/AffineTransform.h"
#include "lsst/geom/AffineTransform.h"
#include "lsst/afw/geom/ellipses.h"
#include "lsst/meas/base.h"
#include "lsst/meas/base/ApertureFlux.h"
Expand Down Expand Up @@ -83,7 +83,7 @@ template <typename MaskedImageT>
void reset(afw::detection::Footprint const&) {}

/// @brief method called for each pixel by applyFunctor
void operator()(afw::geom::Point2I const & pos,
void operator()(geom::Point2I const & pos,
typename MaskedImageT::Image::Pixel const & ival,
typename MaskedImageT::Variance::Pixel const & vval) {
_sum += ival;
Expand Down Expand Up @@ -116,7 +116,7 @@ template <typename MaskedImageT, typename WeightImageT>
class FootprintFindMoment {
public:
FootprintFindMoment(MaskedImageT const& mimage, ///< The image the source lives in
afw::geom::Point2D const& center, // center of the object
geom::Point2D const& center, // center of the object
double const ab, // axis ratio
double const theta // rotation of ellipse +ve from x axis
) : _xcen(center.getX()), _ycen(center.getY()),
Expand All @@ -139,7 +139,7 @@ class FootprintFindMoment {
#endif

MaskedImageT const& mimage = this->getImage();
afw::geom::Box2I const& bbox(foot.getBBox());
geom::Box2I const& bbox(foot.getBBox());
int const x0 = bbox.getMinX(), y0 = bbox.getMinY(), x1 = bbox.getMaxX(), y1 = bbox.getMaxY();

if (x0 < _imageX0 || y0 < _imageY0 ||
Expand All @@ -154,7 +154,7 @@ class FootprintFindMoment {
}

/// @brief method called for each pixel by applyFunctor
void operator()(afw::geom::Point2I const & pos, typename MaskedImageT::Image::Pixel const & ival) {
void operator()(geom::Point2I const & pos, typename MaskedImageT::Image::Pixel const & ival) {
double x = static_cast<double>(pos.getX());
double y = static_cast<double>(pos.getY());
double const dx = x - _xcen;
Expand All @@ -177,7 +177,7 @@ class FootprintFindMoment {
*/

double const eR = 0.38259771140356325; // <r> for a single square pixel, about the centre
r = ::hypot(r, eR*(1 + ::hypot(dx, dy)/afw::geom::ROOT2));
r = ::hypot(r, eR*(1 + ::hypot(dx, dy)/geom::ROOT2));
}
#endif

Expand Down Expand Up @@ -220,7 +220,7 @@ class FootprintFindMoment {

afw::geom::ellipses::Axes KronAperture::getKronAxes(
afw::geom::ellipses::Axes const& shape,
afw::geom::LinearTransform const& transformation,
geom::LinearTransform const& transformation,
double const radius
)
{
Expand All @@ -233,7 +233,7 @@ template<typename ImageT>
PTR(KronAperture) KronAperture::determineRadius(
ImageT const& image,
afw::geom::ellipses::Axes axes,
afw::geom::Point2D const& center,
geom::Point2D const& center,
KronFluxControl const& ctrl
)
{
Expand All @@ -258,7 +258,7 @@ PTR(KronAperture) KronAperture::determineRadius(
//
afw::detection::Footprint foot(afw::geom::SpanSet::fromShape(
afw::geom::ellipses::Ellipse(axes, center)));
afw::geom::Box2I bbox = !smoothImage ?
geom::Box2I bbox = !smoothImage ?
foot.getBBox() :
kernel.growBBox(foot.getBBox()); // the smallest bbox needed to convolve with Kernel
bbox.clip(image.getBBox());
Expand Down Expand Up @@ -323,22 +323,22 @@ std::pair<double, double> photometer(
LSST_EXCEPT_ADD(e, (boost::format("Measuring Kron flux for object at (%.3f, %.3f);"
" aperture radius %g,%g theta %g")
% aperture.getCenter().getX() % aperture.getCenter().getY()
% axes.getA() % axes.getB() % afw::geom::radToDeg(axes.getTheta())).str());
% axes.getA() % axes.getB() % geom::radToDeg(axes.getTheta())).str());
throw e;
}
}


double calculatePsfKronRadius(
CONST_PTR(afw::detection::Psf) const& psf, // PSF to measure
afw::geom::Point2D const& center, // Centroid of source on parent image
geom::Point2D const& center, // Centroid of source on parent image
double smoothingSigma=0.0 // Gaussian sigma of smoothing applied
)
{
assert(psf);
double const radius = psf->computeShape(center).getDeterminantRadius();
// For a Gaussian N(0, sigma^2), the Kron radius is sqrt(pi/2)*sigma
return ::sqrt(afw::geom::PI/2)*::hypot(radius, std::max(0.0, smoothingSigma));
return ::sqrt(geom::PI/2)*::hypot(radius, std::max(0.0, smoothingSigma));
}

template<typename ImageT>
Expand Down Expand Up @@ -436,9 +436,9 @@ void KronFluxAlgorithm::_applyAperture(
void KronFluxAlgorithm::_applyForced(
afw::table::SourceRecord & source,
afw::image::Exposure<float> const & exposure,
afw::geom::Point2D const & center,
geom::Point2D const & center,
afw::table::SourceRecord const & reference,
afw::geom::AffineTransform const & refToMeas
geom::AffineTransform const & refToMeas
) const
{
float const radius = reference.get(reference.getSchema().find<float>(_ctrl.refRadiusName).key);
Expand All @@ -453,7 +453,7 @@ void KronFluxAlgorithm::measure(
afw::table::SourceRecord & source,
afw::image::Exposure<float> const& exposure
) const {
afw::geom::Point2D center = _centroidExtractor(source, _flagHandler);
geom::Point2D center = _centroidExtractor(source, _flagHandler);

// Did we hit a condition that fundamentally prevented measuring the Kron flux?
// Such conditions include hitting the edge of the image and bad input shape, but not low signal-to-noise.
Expand Down Expand Up @@ -563,7 +563,7 @@ void KronFluxAlgorithm::measureForced(
afw::table::SourceRecord const & refRecord,
afw::geom::SkyWcs const & refWcs
) const {
afw::geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
auto xytransform = afw::geom::makeWcsPairTransform(refWcs, *exposure.getWcs());
_applyForced(measRecord, exposure, center, refRecord,
linearizeTransform(*xytransform, refRecord.getCentroid())
Expand Down Expand Up @@ -600,7 +600,7 @@ PTR(KronAperture) KronFluxAlgorithm::_fallbackRadius(afw::table::SourceRecord& s
template PTR(KronAperture) KronAperture::determineRadius<afw::image::MaskedImage<TYPE> >( \
afw::image::MaskedImage<TYPE> const&, \
afw::geom::ellipses::Axes, \
afw::geom::Point2D const&, \
geom::Point2D const&, \
KronFluxControl const& \
); \
template std::pair<double, double> KronAperture::measureFlux<afw::image::MaskedImage<TYPE> >( \
Expand Down
25 changes: 13 additions & 12 deletions tests/test_kron.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
import lsst.utils.tests
import lsst.afw.detection as afwDetection
import lsst.afw.display as afwDisplay
import lsst.geom as geom
import lsst.afw.geom as afwGeom
import lsst.afw.geom.ellipses as afwEllipses
import lsst.afw.image as afwImage
Expand Down Expand Up @@ -82,7 +83,7 @@ def makeGalaxy(width, height, flux, a, b, theta, dx=0.0, dy=0.0, xy0=None, xcen=

if val < 0:
val = 0
gal[afwGeom.Point2I(x, y), afwImage.LOCAL] = val/nsample**2
gal[geom.Point2I(x, y), afwImage.LOCAL] = val/nsample**2

ii += val
iuu += val*u**2
Expand All @@ -93,10 +94,10 @@ def makeGalaxy(width, height, flux, a, b, theta, dx=0.0, dy=0.0, xy0=None, xcen=

exp = afwImage.makeExposure(afwImage.makeMaskedImage(gal))
exp.getMaskedImage().getVariance().set(1.0)
scale = 1.0e-4*afwGeom.degrees
scale = 1.0e-4*geom.degrees
cdMatrix = afwGeom.makeCdMatrix(scale=scale, flipX=True)
exp.setWcs(afwGeom.makeSkyWcs(crpix=afwGeom.Point2D(0.0, 0.0),
crval=afwGeom.SpherePoint(0.0, 0.0, afwGeom.degrees),
exp.setWcs(afwGeom.makeSkyWcs(crpix=geom.Point2D(0.0, 0.0),
crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
cdMatrix=cdMatrix))
# add a dummy Psf. The new SdssCentroid needs one
exp.setPsf(afwDetection.GaussianPsf(11, 11, 0.01))
Expand Down Expand Up @@ -196,7 +197,7 @@ def makeAndMeasure(self, measureKron, a, b, theta, dx=0.0, dy=0.0, nsigma=6, kfa
makeImage = True
if makeImage:
self.objImg = makeGalaxy(self.width, self.height, self.flux, a, b, theta, dx=dx, dy=dy,
xy0=afwGeom.Point2I(10, 10), xcen=xcen, ycen=ycen)
xy0=geom.Point2I(10, 10), xcen=xcen, ycen=ycen)
if display:
disp = afwDisplay.Display(frame=0)
disp.mtv(self.objImg, title=self._testMethodName + ": %g %g" % (a, b))
Expand Down Expand Up @@ -231,7 +232,7 @@ def measureKron(self, objImg, xcen, ycen, nsigma, kfac, nIterForRadius, disp=Non
#
# Now measure things
#
center = afwGeom.Point2D(xcen, ycen)
center = geom.Point2D(xcen, ycen)
msConfig = makeMeasurementConfig(False, nsigma, nIterForRadius, kfac)
source = measureFree(objImg, center, msConfig)
algMeta = source.getTable().getMetadata()
Expand Down Expand Up @@ -292,7 +293,7 @@ def measureKronInPython(self, objImg, xcen, ycen, nsigma, kfac, nIterForRadius,
#
# Note: this code was converted to the new meas_framework, but is not exercised.
msConfig = makeMeasurementConfig(False, nsigma, nIterForRadius, kfac)
center = afwGeom.Point2D(xcen, ycen)
center = geom.Point2D(xcen, ycen)
source = self.measureFree(objImg, center, msConfig)
algMeta = source.getTable().getMetadata()
self.assertTrue(algMeta.exists('ext_photometryKron_KronFlux_nRadiusForFlux'))
Expand All @@ -315,7 +316,7 @@ def measureKronInPython(self, objImg, xcen, ycen, nsigma, kfac, nIterForRadius,
# Get footprint
#
ellipse = afwEllipses.Ellipse(afwEllipses.Axes(nsigma*a, nsigma*b, theta),
afwGeom.Point2D(xcen, ycen))
geom.Point2D(xcen, ycen))
fpEllipse = afwDetection.Footprint(ellipse)

sumI = 0.0
Expand Down Expand Up @@ -527,7 +528,7 @@ def testForced(self):
b = a*axisRatio
for theta in (0, 30, 45):
width, height = 256, 256
center = afwGeom.Point2D(0.5*width, 0.5*height)
center = geom.Point2D(0.5*width, 0.5*height)
original = makeGalaxy(width, height, 1000.0, a, b, theta)
msConfig = makeMeasurementConfig(forced=False, kfac=kfac)
source = measureFree(original, center, msConfig)
Expand All @@ -536,16 +537,16 @@ def testForced(self):
if source.get("ext_photometryKron_KronFlux_flag"):
continue

angleList = [val*afwGeom.degrees for val in (45, 90)]
angleList = [val*geom.degrees for val in (45, 90)]
scaleList = [1.0, 0.5]
offsetList = [(1.23, 4.56), (12.3, 45.6)]

for angle, scale, offset in itertools.product(angleList, scaleList, offsetList):
dx, dy = offset
pixelScale = original.getWcs().getPixelScale()*scale
cdMatrix = afwGeom.makeCdMatrix(scale=pixelScale, orientation=angle, flipX=True)
wcs = afwGeom.makeSkyWcs(crpix=afwGeom.Point2D(dx, dy),
crval=afwGeom.SpherePoint(0.0, 0.0, afwGeom.degrees),
wcs = afwGeom.makeSkyWcs(crpix=geom.Point2D(dx, dy),
crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
cdMatrix=cdMatrix)

warped = warper.warpExposure(wcs, original)
Expand Down

0 comments on commit 1ba1e73

Please sign in to comment.