Skip to content

Commit

Permalink
Merge pull request #253 from lsst/tickets/DM-10728
Browse files Browse the repository at this point in the history 
DM-10728: miscellaneous changes in support of PhotoCalib outputs in meas_mosaic
  • Loading branch information
@TallJimbo
TallJimbo committed Jul 28, 2017
2 parents bade6ec + 5a2997e commit 036366c
Show file tree
Hide file tree
Showing 10 changed files with 272 additions and 81 deletions.
34 changes: 17 additions & 17 deletions include/lsst/afw/image/PhotoCalib.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,10 +173,10 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
* (sourceRecord.get('x'), sourceRecord.get('y')) (pixels) to maggies and maggie error.
*
* @param[in] sourceRecord The source record to get instFlux and position from.
* @param[in] instFluxField The instFlux field: Keys of the form "*_instFlux" and "*_instFluxErr" must
* @param[in] instFluxField The instFlux field: Keys of the form "*_flux" and "*_fluxSigma" must
* exist.
* For example: instFluxField = "PsfFlux" -> "PsfFlux_instFlux",
* "PsfFlux_instFluxErr"
* For example: instFluxField = "PsfFlux" -> "PsfFlux_flux",
* "PsfFlux_fluxSigma"
*
* @returns The flux in maggies and error (err) for this source.
*/
Expand All @@ -188,10 +188,10 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
* (sourceCatalog.get('x'), sourceCatalog.get('y')) (pixels) to maggies.
*
* @param[in] sourceCatalog The source catalog to get instFlux and position from.
* @param[in] instFluxField The instFlux field: Keys of the form "*_instFlux" and "*_instFluxErr" must
* @param[in] instFluxField The instFlux field: Keys of the form "*_flux" and "*_fluxSigma" must
* exist.
* For example: instFluxField = "PsfFlux" -> "PsfFlux_instFlux",
* "PsfFlux_instFluxErr"
* For example: instFluxField = "PsfFlux" -> "PsfFlux_flux",
* "PsfFlux_fluxSigma"
*
* @returns The flux in maggies and error (err) for this source.
*/
Expand All @@ -204,12 +204,12 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
* and write the results back to sourceCatalog[outField_mag].
*
* @param[in] sourceCatalog The source catalog to get instFlux and position from.
* @param[in] instFluxField The instFlux field: Keys of the form "*_instFlux" and "*_instFluxErr" must
* @param[in] instFluxField The instFlux field: Keys of the form "*_flux" and "*_fluxSigma" must
* exist.
* For example: instFluxField = "PsfFlux" -> "PsfFlux_instFlux",
* "PsfFlux_instFluxErr"
* For example: instFluxField = "PsfFlux" -> "PsfFlux_flux",
* "PsfFlux_fluxSigma"
* @param[in] outField The field to write the maggies and maggie errors to.
* Keys of the form "*_flux" and "*_fluxErr" must exist in the schema.
* Keys of the form "*_flux" and "*_fluxSigma" must exist in the schema.
*
* @warning Not implemented yet: See DM-10155.
*/
Expand Down Expand Up @@ -253,10 +253,10 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
* (sourceRecord.get('x'), sourceRecord.get('y')) (pixels) to AB magnitude.
*
* @param[in] sourceRecord The source record to get instFlux and position from.
* @param[in] instFluxField The instFlux field: Keys of the form "*_instFlux" and "*_instFluxErr" must
* @param[in] instFluxField The instFlux field: Keys of the form "*_flux" and "*_fluxSigma" must
* exist.
* For example: instFluxField = "PsfFlux" -> "PsfFlux_instFlux",
* "PsfFlux_instFluxErr"
* "PsfFlux_fluxSigma"
*
* @returns The magnitude and magnitude error for this source.
*/
Expand All @@ -268,10 +268,10 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
* (sourceCatalog.get('x'), sourceCatalog.get('y')) (pixels) to AB magnitudes.
*
* @param[in] sourceCatalog The source catalog to get instFlux and position from.
* @param[in] instFluxField The instFlux field: Keys of the form "*_instFlux" and "*_instFluxErr" must
* @param[in] instFluxField The instFlux field: Keys of the form "*_flux" and "*_fluxSigma" must
* exist.
* For example: instFluxField = "PsfFlux" -> "PsfFlux_instFlux",
* "PsfFlux_instFluxErr"
* "PsfFlux_fluxSigma"
*
* @returns The magnitudes and magnitude errors for the sources.
*/
Expand All @@ -284,12 +284,12 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
* and write the results back to sourceCatalog[outField_mag].
*
* @param[in] sourceCatalog The source catalog to get instFlux and position from.
* @param[in] instFluxField The instFlux field: Keys of the form "*_instFlux" and "*_instFluxErr" must
* @param[in] instFluxField The instFlux field: Keys of the form "*_flux" and "*_fluxSigma" must
* exist.
* For example: instFluxField = "PsfFlux" -> "PsfFlux_instFlux",
* "PsfFlux_instFluxErr"
* "PsfFlux_fluxSigma"
* @param[in] outField The field to write the magnitudes and magnitude errors to.
* Keys of the form "*_flux", "*_fluxErr", *_mag", and "*_magErr"
* Keys of the form "*_flux", "*_fluxSigma", *_mag", and "*_magErr"
* must exist in the schema.
*
* @warning Not implemented yet: See DM-10155.
Expand Down
28 changes: 24 additions & 4 deletions include/lsst/afw/math/BoundedField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,12 +117,17 @@ class BoundedField : public table::io::PersistableFacade<BoundedField>, public t
* @param[out] image Image to fill.
* @param[in] overlapOnly If true, only modify the region in the intersection of
* image.getBBox(image::PARENT) and this->getBBox().
* @param[in] xStep Distance between grid points in X to evaluate; values
* between grid points will be linearly interpolated.
* @param[in] yStep Distance between grid points in Y to evaluate; values
* between grid points will be linearly interpolated.
*
* @throws pex::exceptions::RuntimeError if the bounding boxes do not overlap
* and overlapOnly=false.
*/
template <typename T>
void fillImage(image::Image<T>& image, bool overlapOnly = false) const;
void fillImage(image::Image<T>& image, bool overlapOnly = false,
int xStep = 1, int yStep = 1) const;

/**
* Add the field or a constant multiple of it to an image in-place
Expand All @@ -131,38 +136,53 @@ class BoundedField : public table::io::PersistableFacade<BoundedField>, public t
* @param[in] scaleBy Multiply the field by this before adding it to the image.
* @param[in] overlapOnly If true, only modify the region in the intersection of
* image.getBBox(image::PARENT) and this->getBBox().
* @param[in] xStep Distance between grid points in X to evaluate; values
* between grid points will be linearly interpolated.
* @param[in] yStep Distance between grid points in Y to evaluate; values
* between grid points will be linearly interpolated.
*
* @throws pex::exceptions::RuntimeError if the bounding boxes do not overlap
* and overlapOnly=false.
*/
template <typename T>
void addToImage(image::Image<T>& image, double scaleBy = 1.0, bool overlapOnly = false) const;
void addToImage(image::Image<T>& image, double scaleBy = 1.0, bool overlapOnly = false,
int xStep = 1, int yStep = 1) const;

/**
* Multiply an image by the field in-place.
*
* @param[out] image Image to fill.
* @param[in] overlapOnly If true, only modify the region in the intersection of
* image.getBBox(image::PARENT) and this->getBBox().
* @param[in] xStep Distance between grid points in X to evaluate; values
* between grid points will be linearly interpolated.
* @param[in] yStep Distance between grid points in Y to evaluate; values
* between grid points will be linearly interpolated.
*
* @throws pex::exceptions::RuntimeError if the bounding boxes do not overlap
* and overlapOnly=false.
*/
template <typename T>
void multiplyImage(image::Image<T>& image, bool overlapOnly = false) const;
void multiplyImage(image::Image<T>& image, bool overlapOnly = false,
int xStep = 1, int yStep = 1) const;

/**
* Divide an image by the field in-place.
*
* @param[out] image Image to fill.
* @param[in] overlapOnly If true, only modify the region in the intersection of
* image.getBBox(image::PARENT) and this->getBBox().
* @param[in] xStep Distance between grid points in X to evaluate; values
* between grid points will be linearly interpolated.
* @param[in] yStep Distance between grid points in Y to evaluate; values
* between grid points will be linearly interpolated.
*
* @throws pex::exceptions::RuntimeError if the bounding boxes do not overlap
* and overlapOnly=false.
*/
template <typename T>
void divideImage(image::Image<T>& image, bool overlapOnly = false) const;
void divideImage(image::Image<T>& image, bool overlapOnly = false,
int xStep = 1, int yStep = 1) const;

/**
* Return a scaled BoundedField
Expand Down
2 changes: 1 addition & 1 deletion python/lsst/afw/image/slicing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ def _getBBoxFromSliceTuple(img, imageSlice):
if not (isinstance(imageSlice, tuple) and len(imageSlice) == 2 and
sum([isinstance(_, (slice, type(Ellipsis), int)) for _ in imageSlice]) == 2):
raise IndexError(
"Images may only be indexed as a 2-D slice not %s", imageSlice)
"Images may only be indexed as a 2-D slice not %s" % imageSlice)

imageSlice, _imageSlice = [], imageSlice
for s, wh in zip(_imageSlice, img.getDimensions()):
Expand Down
25 changes: 16 additions & 9 deletions python/lsst/afw/math/boundedField.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
#include "lsst/afw/table/io/Persistable.h"
#include "lsst/afw/geom/Point.h"
#include "lsst/afw/math/BoundedField.h"
#include "lsst/afw/table/io/python.h"

namespace py = pybind11;

Expand All @@ -43,15 +44,20 @@ namespace afw {
namespace math {
namespace {

using PyClass = py::class_<BoundedField, std::shared_ptr<BoundedField>, lsst::afw::table::io::Persistable>;
using PyClass = py::class_<BoundedField, std::shared_ptr<BoundedField>,
afw::table::io::PersistableFacade<BoundedField>,
afw::table::io::Persistable>;

template <typename PixelT>
void declareTemplates(PyClass &cls) {
cls.def("fillImage", &BoundedField::fillImage<PixelT>, "image"_a, "overlapOnly"_a = false);
cls.def("fillImage", &BoundedField::fillImage<PixelT>,
"image"_a, "overlapOnly"_a = false, "xStep"_a = 1, "yStep"_a = 1);
cls.def("addToImage", &BoundedField::addToImage<PixelT>, "image"_a, "scaleBy"_a = 1.0,
"overlapOnly"_a = false);
cls.def("multiplyImage", &BoundedField::multiplyImage<PixelT>, "image"_a, "overlapOnly"_a = false);
cls.def("divideImage", &BoundedField::divideImage<PixelT>, "image"_a, "overlapOnly"_a = false);
"overlapOnly"_a = false, "xStep"_a = 1, "yStep"_a = 1);
cls.def("multiplyImage", &BoundedField::multiplyImage<PixelT>,
"image"_a, "overlapOnly"_a = false, "xStep"_a = 1, "yStep"_a = 1);
cls.def("divideImage", &BoundedField::divideImage<PixelT>,
"image"_a, "overlapOnly"_a = false, "xStep"_a = 1, "yStep"_a = 1);
}

PYBIND11_PLUGIN(_boundedField) {
Expand All @@ -62,13 +68,14 @@ PYBIND11_PLUGIN(_boundedField) {
return nullptr;
};

afw::table::io::python::declarePersistableFacade<BoundedField>(mod, "BoundedField");
PyClass cls(mod, "BoundedField");

cls.def("__rmul__", [](BoundedField &bf, double const scale) { return bf * scale; }, py::is_operator());
cls.def("__mul__", &BoundedField::operator*);
cls.def("__truediv__", &BoundedField::operator/);
cls.def("__eq__", &BoundedField::operator==);
cls.def("__ne__", &BoundedField::operator!=);
cls.def("__mul__", &BoundedField::operator*, py::is_operator());
cls.def("__truediv__", &BoundedField::operator/, py::is_operator());
cls.def("__eq__", &BoundedField::operator==, py::is_operator());
cls.def("__ne__", &BoundedField::operator!=, py::is_operator());

cls.def("evaluate", (double (BoundedField::*)(double, double) const) & BoundedField::evaluate);
cls.def("evaluate",
Expand Down
11 changes: 4 additions & 7 deletions python/lsst/afw/math/chebyshevBoundedField.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ namespace py = pybind11;

using namespace lsst::afw::math;

using ClsField = py::class_<ChebyshevBoundedField, std::shared_ptr<ChebyshevBoundedField>, BoundedField,
lsst::afw::table::io::PersistableFacade<ChebyshevBoundedField>>;
using ClsField = py::class_<ChebyshevBoundedField, std::shared_ptr<ChebyshevBoundedField>,
BoundedField>;

template <typename PixelT>
void declareTemplates(ClsField & cls) {
Expand All @@ -58,8 +58,6 @@ PYBIND11_PLUGIN(_chebyshevBoundedField) {
};

/* Module level */
lsst::afw::table::io::python::declarePersistableFacade<ChebyshevBoundedField>(mod,
"ChebyshevBoundedField");
py::class_<ChebyshevBoundedFieldControl> clsChebyshevBoundedFieldControl(mod,
"ChebyshevBoundedFieldControl");
ClsField clsChebyshevBoundedField(mod, "ChebyshevBoundedField");
Expand All @@ -72,9 +70,8 @@ PYBIND11_PLUGIN(_chebyshevBoundedField) {
clsChebyshevBoundedField.def(
py::init<lsst::afw::geom::Box2I const &, ndarray::Array<double const, 2, 2> const &>());

/* Operators */
clsChebyshevBoundedField.def("__mul__", &ChebyshevBoundedField::operator*, py::is_operator());
clsChebyshevBoundedField.def("__eq__", &ChebyshevBoundedField::operator==, py::is_operator());
/* Operators are defined only in the BoundedField base class;
we let Python inheritance provide them here. */

/* Members */
LSST_DECLARE_CONTROL_FIELD(clsChebyshevBoundedFieldControl, ChebyshevBoundedFieldControl, orderX);
Expand Down
17 changes: 8 additions & 9 deletions src/image/Calib.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,13 +158,13 @@ inline double convertToFluxErr(double fluxMag0InvSNR, double flux, double magErr
return flux * std::abs(a) * std::sqrt(1 + std::pow(b / a, 2));
}
inline double convertToMag(double fluxMag0, double flux) { return -2.5 * ::log10(flux / fluxMag0); }
inline void convertToMagWithErr(double* mag, double* magErr, double fluxMag0, double fluxMag0InvSNR,
double flux, double fluxErr) {
double const rat = flux / fluxMag0;
double const ratErr = ::sqrt((::pow(fluxErr, 2) + ::pow(flux * fluxMag0InvSNR, 2)) / ::pow(fluxMag0, 2));

*mag = -2.5 * ::log10(rat);
*magErr = 2.5 / ::log(10.0) * ratErr / rat;
inline void convertToMagWithErr(double* mag, double* magErr, double fluxMag0, double fluxMag0Err,
double flux, double fluxErr) {
*mag = -2.5*std::log10(flux/fluxMag0);
double const x = fluxErr/flux;
double const y = fluxMag0Err/fluxMag0;
*magErr = (2.5/std::log(10.0))*std::sqrt(x*x + y*y);
}

} // anonymous namespace
Expand Down Expand Up @@ -233,7 +233,7 @@ std::pair<double, double> Calib::getMagnitude(double const flux, double const fl
}

double mag, magErr;
convertToMagWithErr(&mag, &magErr, _fluxMag0, _fluxMag0Sigma / _fluxMag0, flux, fluxErr);
convertToMagWithErr(&mag, &magErr, _fluxMag0, _fluxMag0Sigma, flux, fluxErr);
return std::make_pair(mag, magErr);
}

Expand Down Expand Up @@ -275,7 +275,6 @@ std::pair<ndarray::Array<double, 1>, ndarray::Array<double, 1>> Calib::getMagnit
ndarray::Array<double, 1>::Iterator magIter = mag.begin();
ndarray::Array<double, 1>::Iterator magErrIter = magErr.begin();
int nonPositive = 0;
double fluxMag0InvSNR = _fluxMag0Sigma / _fluxMag0;
for (; fluxIter != flux.end(); ++fluxIter, ++fluxErrIter, ++magIter, ++magErrIter) {
if (isNegativeFlux(*fluxIter, false)) {
++nonPositive;
Expand All @@ -285,7 +284,7 @@ std::pair<ndarray::Array<double, 1>, ndarray::Array<double, 1>> Calib::getMagnit
continue;
}
double f, df;
convertToMagWithErr(&f, &df, _fluxMag0, fluxMag0InvSNR, *fluxIter, *fluxErrIter);
convertToMagWithErr(&f, &df, _fluxMag0, _fluxMag0Sigma, *fluxIter, *fluxErrIter);
*magIter = f;
*magErrIter = df;
}
Expand Down
16 changes: 8 additions & 8 deletions src/image/PhotoCalib.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,8 @@ Measurement PhotoCalib::instFluxToMaggies(double instFlux, double instFluxErr) c
Measurement PhotoCalib::instFluxToMaggies(afw::table::SourceRecord const &sourceRecord,
std::string const &instFluxField) const {
auto position = sourceRecord.getCentroid();
auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_flux").key;
auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_fluxSigma").key;
return instFluxToMaggies(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
}
ndarray::Array<double, 2, 2> PhotoCalib::instFluxToMaggies(afw::table::SourceCatalog const &sourceCatalog,
Expand Down Expand Up @@ -138,8 +138,8 @@ Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr)
Measurement PhotoCalib::instFluxToMagnitude(afw::table::SourceRecord const &sourceRecord,
std::string const &instFluxField) const {
auto position = sourceRecord.getCentroid();
auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_flux").key;
auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_fluxSigma").key;
return instFluxToMagnitude(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
}

Expand Down Expand Up @@ -264,8 +264,8 @@ void PhotoCalib::instFluxToMaggiesArray(afw::table::SourceCatalog const &sourceC
std::string const &instFluxField,
ndarray::Array<double, 2, 2> result) const {
double instFlux, instFluxErr, maggies, instFluxMag0;
auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_flux").key;
auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_fluxSigma").key;
auto iter = result.begin();
for (auto const &rec : sourceCatalog) {
instFlux = rec.get(instFluxKey);
Expand All @@ -285,8 +285,8 @@ void PhotoCalib::instFluxToMagnitudeArray(afw::table::SourceCatalog const &sourc
std::string const &instFluxField,
ndarray::Array<double, 2, 2> result) const {
double instFlux, instFluxErr, instFluxMag0;
auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_flux").key;
auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_fluxSigma").key;
auto iter = result.begin();
for (auto const &rec : sourceCatalog) {
instFlux = rec.get(instFluxKey);
Expand Down

0 comments on commit 036366c

Please sign in to comment.