Review for DM-3257 (Port flux.scaled from HSC)

include/lsst/meas/base.h

@@ -1,7 +1,7 @@
 // -*- lsst-c++ -*-
 /*
  * LSST Data Management System
- * Copyright 2008-2015 LSST Corporation.
+ * Copyright 2008-2016 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -40,6 +40,7 @@
 #include "lsst/meas/base/NaiveCentroid.h"
 #include "lsst/meas/base/PeakLikelihoodFlux.h"
 #include "lsst/meas/base/ApertureFlux.h"
+#include "lsst/meas/base/ScaledApertureFlux.h"
 #include "lsst/meas/base/CircularApertureFlux.h"
 
 #endif // !LSST_MEAS_base_h_INCLUDED

include/lsst/meas/base/ApertureFlux.h

@@ -1,7 +1,7 @@
 // -*- lsst-c++ -*-
 /*
  * LSST Data Management System
- * Copyright 2008-2015 AURA/LSST.
+ * Copyright 2008-2016 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -24,6 +24,7 @@
 #ifndef LSST_MEAS_BASE_ApertureFlux_h_INCLUDED
 #define LSST_MEAS_BASE_ApertureFlux_h_INCLUDED
 
+#include "boost/array.hpp"
 #include "lsst/pex/config.h"
 #include "lsst/afw/image/Exposure.h"
 #include "lsst/afw/table/arrays.h"
@@ -218,6 +219,11 @@ class ApertureFluxAlgorithm : public SimpleAlgorithm {
      */
     static std::string makeFieldPrefix(std::string const & name, double radius);
 
+    /**
+     *  Return the flag definitions which apply to aperture flux measurements.
+     */
+    static boost::array<FlagDefinition,ApertureFluxAlgorithm::N_FLAGS> const & getFlagDefinitions();
+
 protected:
 
     void copyResultToRecord(Result const & result, afw::table::SourceRecord & record, int index) const;

include/lsst/meas/base/ScaledApertureFlux.h

@@ -0,0 +1,101 @@
+// -*- lsst-c++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2008-2016 AURA/LSST.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program.  If not,
+ * see <http://www.lsstcorp.org/LegalNotices/>.
+ */
+
+#ifndef LSST_MEAS_BASE_ScaledApertureFlux_h_INCLUDED
+#define LSST_MEAS_BASE_ScaledApertureFlux_h_INCLUDED
+
+#include "lsst/afw/table.h"
+#include "lsst/afw/image/Exposure.h"
+#include "lsst/meas/base/Algorithm.h"
+#include "lsst/meas/base/FluxUtilities.h"
+#include "lsst/meas/base/FlagHandler.h"
+#include "lsst/meas/base/InputUtilities.h"
+#include "lsst/pex/config.h"
+
+namespace lsst { namespace meas { namespace base {
+
+class ScaledApertureFluxControl {
+public:
+    LSST_CONTROL_FIELD(
+        shiftKernel, std::string,
+        "Warping kernel used to shift Sinc photometry coefficients to different center positions"
+    );
+    LSST_CONTROL_FIELD(scale, double, "Scaling factor of PSF FWHM for aperture radius.");
+
+    // The default scaling factor is chosen such that scaled aperture
+    // magnitudes are expected to be equal to Kron magnitudes, based on
+    // measurements performed by Stephen Gwyn on WIRCam. See:
+    // http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/wirwolf/docs/proc.html#photcal
+    // http://www.cfht.hawaii.edu/fr/news/UM2013/presentations/Session10-SGwyn.pdf
+    ScaledApertureFluxControl() : shiftKernel("lanczos5"), scale(3.14) {}
+};
+
+
+/**
+ *  @brief Measure the flux in an aperture scaled to the PSF.
+ *
+ *  This algorithm performs a sinc aperture flux measurement where they size
+ *  of the aperture is determined by multiplying the FWHM of the PSF by the
+ *  scaling factor specified in the algorithm configuration.
+ */
+class ScaledApertureFluxAlgorithm : public SimpleAlgorithm {
+public:
+
+    typedef ScaledApertureFluxControl Control;
+    typedef ApertureFluxResult Result;
+
+    ScaledApertureFluxAlgorithm(Control const & control, std::string const & name,
+                                afw::table::Schema & schema);
+
+    /**
+     *  Measure the scaled aperture flux on the given image.
+     *
+     *  Python plugins will delegate to this method.
+     *
+     *  @param[in,out] record      Record used to save outputs and retrieve positions.
+     *  @param[in]     exposure    Image to be measured.
+     */
+    virtual void measure(
+        afw::table::SourceRecord & measRecord,
+        afw::image::Exposure<float> const & exposure
+    ) const override;
+
+    virtual void fail(afw::table::SourceRecord & measRecord, MeasurementError * error=NULL) const override;
+
+private:
+    Control _ctrl;
+    FluxResultKey _fluxResultKey;
+    FlagHandler _flagHandler;
+    SafeCentroidExtractor _centroidExtractor;
+};
+
+class ScaledApertureFluxTransform : public FluxTransform {
+public:
+    typedef ScaledApertureFluxControl Control;
+    ScaledApertureFluxTransform(Control const & ctrl, std::string const & name,
+                                afw::table::SchemaMapper & mapper);
+};
+
+}}} // namespace lsst::meas::base
+
+#endif // !LSST_MEAS_BASE_ScaledApertureFlux_h_INCLUDED

python/lsst/meas/base/plugins.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 #
 # LSST Data Management System
-# Copyright 2008-2015 AURA/LSST.
+# Copyright 2008-2016 AURA/LSST.
 #
 # This product includes software developed by the
 # LSST Project (http://www.lsst.org/).
@@ -67,6 +67,8 @@
 wrapSimpleAlgorithm(bl.PixelFlagsAlgorithm, Control=bl.PixelFlagsControl, executionOrder=BasePlugin.FLUX_ORDER)
 wrapSimpleAlgorithm(bl.SdssShapeAlgorithm, Control=bl.SdssShapeControl,
                 TransformClass=bl.SdssShapeTransform, executionOrder=BasePlugin.SHAPE_ORDER)
+wrapSimpleAlgorithm(bl.ScaledApertureFluxAlgorithm, Control=bl.ScaledApertureFluxControl,
+                TransformClass=bl.ScaledApertureFluxTransform, executionOrder=BasePlugin.FLUX_ORDER)
 
 wrapSimpleAlgorithm(bl.CircularApertureFluxAlgorithm, needsMetadata=True, Control=bl.ApertureFluxControl,
                     TransformClass=bl.ApertureFluxTransform, executionOrder=BasePlugin.FLUX_ORDER)

python/lsst/meas/base/pluginsLib.i

@@ -1,7 +1,7 @@
 // -*- lsst-c++ -*-
 /*
  * LSST Data Management System
- * Copyright 2008-2015 AURA/LSST.
+ * Copyright 2008-2016 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -39,6 +39,7 @@
 %convertConfig(lsst::meas::base, GaussianFluxTransform)
 %convertConfig(lsst::meas::base, PeakLikelihoodFluxTransform)
 %convertConfig(lsst::meas::base, PsfFluxTransform)
+%convertConfig(lsst::meas::base, ScaledApertureFluxTransform)
 
 %convertConfig(lsst::meas::base, NaiveCentroidTransform)
 %convertConfig(lsst::meas::base, GaussianCentroidTransform)
@@ -61,6 +62,9 @@
 %feature("notabstract") lsst::meas::base::GaussianFluxAlgorithm;
 %include "lsst/meas/base/GaussianFlux.h"
 
+%feature("notabstract") lsst::meas::base::ScaledApertureFluxAlgorithm;
+%include "lsst/meas/base/ScaledApertureFlux.h"
+
 // centroid algorithms
 
 %feature("notabstract") lsst::meas::base::GaussianCentroidAlgorithm;

src/ApertureFlux.cc

@@ -1,7 +1,7 @@
 // -*- lsst-c++ -*-
 /*
  * LSST Data Management System
- * Copyright 2008-2015 AURA/LSST.
+ * Copyright 2008-2016 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -23,7 +23,6 @@
 
 #include <numeric>
 
-#include "boost/array.hpp"
 #include "boost/algorithm/string/replace.hpp"
 
 #include "ndarray/eigen.h"
@@ -44,21 +43,15 @@ ApertureFluxControl::ApertureFluxControl() : radii(10), maxSincRadius(10.0), shi
     std::copy(defaultRadii.begin(), defaultRadii.end(), radii.begin());
 }
 
-namespace {
-
-boost::array<FlagDefinition,ApertureFluxAlgorithm::N_FLAGS> const & getFlagDefinitions() {
+boost::array<FlagDefinition,ApertureFluxAlgorithm::N_FLAGS> const & ApertureFluxAlgorithm::getFlagDefinitions() {
     static boost::array<FlagDefinition,ApertureFluxAlgorithm::N_FLAGS> flagDefs = {{
-        {"flag", "flag set if aperture failed for any reason"},
-        {"flag_apertureTruncated", "flag set if aperture did not fit within the measurement image"},
-        {"flag_sincCoeffsTruncated",
-         "flag set if the full sinc coefficient image for aperture %d did not "
-         "fit within the measurement image"}
+        {"flag", "general failure flag"},
+        {"flag_apertureTruncated", "aperture did not fit within measurement image"},
+        {"flag_sincCoeffsTruncated", "full sinc coefficient image did not fit within measurement image"}
     }};
     return flagDefs;
 }
 
-} // anonymous
-
 std::string ApertureFluxAlgorithm::makeFieldPrefix(std::string const & name, double radius) {
     std::string prefix = (boost::format("%s_%.1f") % name % radius).str();
     return boost::replace_all_copy(prefix, ".", "_");
@@ -71,8 +64,8 @@ ApertureFluxAlgorithm::Keys::Keys(
     flags(
         FlagHandler::addFields(
             schema, prefix,
-            getFlagDefinitions().begin(),
-            getFlagDefinitions().begin() + (isSinc ? 3 : 2)
+            ApertureFluxAlgorithm::getFlagDefinitions().begin(),
+            ApertureFluxAlgorithm::getFlagDefinitions().begin() + (isSinc ? 3 : 2)
         )
     )
 {}
@@ -329,8 +322,9 @@ ApertureFluxTransform::ApertureFluxTransform(
     _ctrl(ctrl)
 {
     for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
-        for (auto flag = getFlagDefinitions().begin();
-            flag < getFlagDefinitions().begin() + (_ctrl.radii[i] <= _ctrl.maxSincRadius ? 3 : 2); flag++) {
+        for (auto flag = ApertureFluxAlgorithm::getFlagDefinitions().begin();
+            flag < ApertureFluxAlgorithm::getFlagDefinitions().begin() +
+                   (_ctrl.radii[i] <= _ctrl.maxSincRadius ? 3 : 2); flag++) {
             mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>((boost::format("%s_%s") %
                               ApertureFluxAlgorithm::makeFieldPrefix(name, _ctrl.radii[i]) %
                               flag->name).str()).key);

src/ScaledApertureFlux.cc

@@ -0,0 +1,95 @@
+// -*- lsst-c++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2008-2016 AURA/LSST.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program.  If not,
+ * see <http://www.lsstcorp.org/LegalNotices/>.
+ */
+
+#include "lsst/afw/table/Source.h"
+#include "lsst/meas/base/ApertureFlux.h"
+#include "lsst/meas/base/ScaledApertureFlux.h"
+#include "lsst/meas/base/SincCoeffs.h"
+#include "lsst/afw/detection/Psf.h"
+
+namespace lsst { namespace meas { namespace base {
+
+ScaledApertureFluxAlgorithm::ScaledApertureFluxAlgorithm(
+    Control const & ctrl,
+    std::string const & name,
+    afw::table::Schema & schema
+) : _ctrl(ctrl),
+    _fluxResultKey(
+        FluxResultKey::addFields(schema, name, "flux derived from PSF-scaled aperture")
+    ),
+    _centroidExtractor(schema, name)
+{
+    _flagHandler = FlagHandler::addFields(schema, name,
+                                          ApertureFluxAlgorithm::getFlagDefinitions().begin(),
+                                          ApertureFluxAlgorithm::getFlagDefinitions().end());
+}
+
+void ScaledApertureFluxAlgorithm::measure(
+    afw::table::SourceRecord & measRecord,
+    afw::image::Exposure<float> const & exposure
+) const {
+    afw::geom::Point2D const center = _centroidExtractor(measRecord, _flagHandler);
+    double const radius = exposure.getPsf()->computeShape(center).getDeterminantRadius();
+    double const fwhm = 2.0*std::sqrt(2.0*std::log(2))*radius;
+    double const size = _ctrl.scale*fwhm;
+    afw::geom::ellipses::Axes const axes(size, size);
+
+    // ApertureFluxAlgorithm::computeSincFlux requires an ApertureFluxControl as an
+    // argument. All that it uses it for is to read the type of warping kernel.
+    ApertureFluxControl apCtrl;
+    apCtrl.shiftKernel = _ctrl.shiftKernel;
+
+    Result result = ApertureFluxAlgorithm::computeSincFlux(exposure.getMaskedImage(),
+                                                           afw::geom::ellipses::Ellipse(axes, center),
+                                                           apCtrl);
+    measRecord.set(_fluxResultKey, result);
+    if (result.getFlag(ApertureFluxAlgorithm::FAILURE)) {
+        _flagHandler.setValue(measRecord, ApertureFluxAlgorithm::FAILURE, true);
+    }
+    if (result.getFlag(ApertureFluxAlgorithm::APERTURE_TRUNCATED)) {
+        _flagHandler.setValue(measRecord, ApertureFluxAlgorithm::APERTURE_TRUNCATED, true);
+    }
+    if (result.getFlag(ApertureFluxAlgorithm::SINC_COEFFS_TRUNCATED)) {
+        _flagHandler.setValue(measRecord, ApertureFluxAlgorithm::SINC_COEFFS_TRUNCATED, true);
+    }
+}
+
+void ScaledApertureFluxAlgorithm::fail(afw::table::SourceRecord & measRecord,
+                                       MeasurementError * error) const {
+    _flagHandler.handleFailure(measRecord, error);
+}
+
+ScaledApertureFluxTransform::ScaledApertureFluxTransform(
+    Control const & ctrl,
+    std::string const & name,
+    afw::table::SchemaMapper & mapper
+) :
+    FluxTransform{name, mapper}
+{
+    for (auto flag = ApertureFluxAlgorithm::getFlagDefinitions().begin() + 1;
+         flag < ApertureFluxAlgorithm::getFlagDefinitions().end(); flag++) {
+        mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag->name).key);
+    }
+}
+
+}}} // namespace lsst::meas::base