Skip to content

Instantly share code, notes, and snippets.

@parejkoj

parejkoj/PhotoCalib.cpp

Last active Feb 8, 2017
Embed
What would you like to do?
Header for new PhotoCalib object.
/**
* @brief The photometric calibration of an exposure.
*
* A PhotoCalib is a BoundedField (a function with a specified domain) that converts between calibrated
* counts-on-chip (ADU) to flux and magnitude. It is defined in terms of "maggies", which are a linear
* unit defined in SDSS: http://www.sdss.org/dr12/algorithms/magnitudes/#nmgy
*
* PhotoCalib is immutable.
*
* The spatially varying flux/magnitude zero point is defined such that,
* at a position (x,y) in the domain of the boundedField zeroPoint
* and for a given measured source counts:
* zeroPoint(x,y) * counts = flux (in maggies)
* while the errors (constant on the domain) are defined as:
* sqrt(countsSigma^2 + zeroPointSigma^2) = fluxSigma (in maggies)
*/
class PhotoCalib {
public:
PhotoCalib(PhotoCalib const &) = delete;
PhotoCalib(PhotoCalib &&) = delete;
PhotoCalib & operator=(PhotoCalib const &) = delete;
PhotoCalib & operator=(PhotoCalib &&) = delete;
/**
* @brief Create a non-spatially-varying calibration.
*
* @param[in] fluxMag0 The constant flux/magnitude zero point (counts at magnitude 0).
* @param[in] fluxMag0Sigma The error on the zero point.
*/
PhotoCalib(double fluxMag0, double fluxMag0Sigma=0) const;
/**
* @brief Create a spatially-varying calibration.
*
* @param[in] zeroPoint The spatially varying photometric zero point.
* @param[in] fluxMag0Sigma The error on the zero point.
*/
PhotoCalib(std::shared_ptr<lsst::afw::math::BoundedField> zeroPoint, double fluxMag0Sigma=0) const;
/**
* @brief Convert counts in ADU to maggies.
*
* If passed point, use the exact calculation at that point, otherwise, use the mean scaling factor.
*
* @param[in] counts The source flux in ADU.
* @param[in] point The point that flux is measured at (must be within the domain of the
* BoundedField of this PhotoCalib).
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The flux in maggies.
*/
double countsToMaggies(double counts,
const lsst::afw::geom::Point &point=nullptr,
bool throwOnNegativeFlux=false) const;
/**
* @brief Convert counts and error in counts (ADU) to maggies and maggies error.
*
* If passed point, use the exact calculation at that point, otherwise, use the mean scaling factor.
*
* @param[in] counts The source flux in ADU.
* @param[in] countsSigma The counts error (sigma).
* @param[in] point The point that flux is measured at (must be within the domain of the
* BoundedField of this PhotoCalib).
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The flux in maggies and error (sigma).
*/
std::pair<double, double> countsToMaggies(double counts, double countsSigma,
const lsst::afw::geom::Point &point=nullptr,
bool throwOnNegativeFlux=false) const;
/**
* @brief Convert sourceRecord[fluxField] (ADU) at location
* (sourceRecord.get('x'), sourceRecord.get('y')) (pixels) to maggies and maggie error.
*
* @param[in] sourceRecord The source record to get flux and position from.
* @param[in] fluxField The flux field: Keys of the form "*_flux" and "*_fluxSigma" must exist.
* For example: fluxField = "PsfFlux" -> "PsfFlux_flux", "PsfFlux_fluxSigma"
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The flux in maggies and error (sigma) for this source.
*/
std::pair<double, double> countsToMaggies(const lsst::afw::table::SourceRecord &sourceRecord,
const std::string &fluxField,
bool throwOnNegativeFlux=false) const;
//@{
/**
* @brief Convert sourceCatalog[fluxField] (ADU) at locations
* (sourceCatalog.get('x'), sourceCatalog.get('y')) (pixels) to maggies.
*
* @param[in] sourceCatalog The source catalog to get flux and position from.
* @param[in] fluxField The flux field: Keys of the form "*_flux" and "*_fluxSigma" must exist.
* For example: fluxField = "PsfFlux" -> "PsfFlux_flux", "PsfFlux_fluxSigma"
* @param[out] maggies, maggiesSigma Pass pre-allocated arrays to put the output in, instead of
* returning it.
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The flux in maggies and error (sigma) for this source.
*/
std::pair< ndarray::Array<double,1>, ndarray::Array<double,1> >
countsToMaggies(const lsst::afw::table::SourceCatalog &sourceCatalog,
const std::string &fluxField,
bool throwOnNegativeFlux=false) const;
void countsToMaggies(const lsst::afw::table::SourceCatalog &sourceCatalog,
const std::string &fluxField,
ndarray::Array<double,1> maggies,
ndarray::Array<double,1> maggiesSigma,
bool throwOnNegativeFlux=false) const;
//@}
/**
* @brief Convert sourceCatalog[fluxField_flux] (ADU) at locations
* (sourceCatalog.get('x'), sourceCatalog.get('x')) (pixels) to maggies
* and write the results back to sourceCatalog[outField_mag].
*
* @param[in] sourceCatalog The source catalog to get flux and position from.
* @param[in] fluxField The flux field: Keys of the form "*_flux" and "*_fluxSigma" must exist.
* For example: fluxField = "PsfFlux" -> "PsfFlux_flux", "PsfFlux_fluxSigma"
* @param[in] outField The field to write the maggies and maggie errors to.
* Keys of the form "*_flux" and "*_fluxSigma" must exist in the schema.
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*/
void countsToMaggies(lsst::afw::table::sourceCatalog &sourceCatalog,
const std::string &fluxField,
const std::string &outField,
bool throwOnNegativeFlux=false) const;
/**
* @brief Convert counts in ADU to AB magnitude.
*
* If passed point, use the exact calculation at that point, otherwise, use the mean scaling factor.
*
* @param[in] counts The source flux in ADU.
* @param[in] point The point that flux is measured at (must be within the domain of the
* BoundedField of this PhotoCalib).
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The AB magnitude.
*/
double countsToMagnitude(double counts,
const lsst::afw::geom::Point &point=nullptr,
bool throwOnNegativeFlux=false) const;
/**
* @brief Convert counts and error in counts (ADU) to AB magnitude and magnitude error.
*
* If passed point, use the exact calculation at that point, otherwise, use the mean scaling factor.
*
* @param[in] counts The source flux in ADU.
* @param[in] countsSigma The counts error (standard deviation).
* @param[in] point The point that flux is measured at (must be within the domain of the
* BoundedField of this PhotoCalib).
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The AB magnitude and error (sigma).
*/
std::pair<double, double> countsToMagnitude(double counts, double countsSigma,
const lsst::afw::geom::Point &point=nullptr,
bool throwOnNegativeFlux=false) const;
/**
* @brief Convert sourceRecord[fluxField] (ADU) at location
* (sourceRecord.get('x'), sourceRecord.get('y')) (pixels) to AB magnitude.
*
* @param[in] sourceRecord The source record to get flux and position from.
* @param[in] fluxField The flux field: Keys of the form "*_flux" and "*_fluxSigma" must exist.
* For example: fluxField = "PsfFlux" -> "PsfFlux_flux", "PsfFlux_fluxSigma"
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The magnitude and magnitude error for this source.
*/
std::pair<double, double> countsToMagnitude(const lsst::afw::table::SourceRecord &sourceRecord,
const std::string &fluxField,
bool throwOnNegativeFlux=false) const;
//@{
/**
* @brief Convert sourceCatalog[fluxField] (ADU) at locations
* (sourceCatalog.get('x'), sourceCatalog.get('y')) (pixels) to AB magnitudes.
*
* @param[in] sourceCatalog The source catalog to get flux and position from.
* @param[in] fluxField The flux field: Keys of the form "*_flux" and "*_fluxSigma" must exist.
* For example: fluxField = "PsfFlux" -> "PsfFlux_flux", "PsfFlux_fluxSigma"
* @param[out] mag, magSigma Pass pre-allocated arrays to put the output in, instead of
* returning it.
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*
* @return The magnitudes and magnitude errors for the sources.
*/
std::pair< ndarray::Array<double,1>, ndarray::Array<double,1> >
countsToMagnitude(const lsst::afw::table::SourceCatalog &sourceCatalog,
const std::string &fluxField,
bool throwOnNegativeFlux=false) const;
void countsToMagnitude(const lsst::afw::table::SourceCatalog &sourceCatalog,
const std::string &fluxField,
ndarray::Array<double,1>, mag,
ndarray::Array<double,1> magSigma,
bool throwOnNegativeFlux=false) const;
//@}
/**
* @brief Convert sourceCatalog[fluxField_flux] (ADU) at locations
* (sourceCatalog.get('x'), sourceCatalog.get('x')) (pixels) to AB magnitudes
* and write the results back to sourceCatalog[outField_mag].
*
* @param[in] sourceCatalog The source catalog to get flux and position from.
* @param[in] fluxField The flux field: Keys of the form "*_flux" and "*_fluxSigma" must exist.
* For example: fluxField = "PsfFlux" -> "PsfFlux_flux", "PsfFlux_fluxSigma"
* @param[in] outField The field to write the magnitudes and magnitude errors to.
* Keys of the form "*_mag" and "*_magSigma" must exist in the schema.
* @param[in] throwOnNegativeFlux Raise an exception when passed negative counts.
*/
void countsToMagnitude(lsst::afw::table::sourceCatalog &sourceCatalog,
const std::string &fluxField,
const std::string &outField,
bool throwOnNegativeFlux=false) const;
/**
* @brief Convert AB magnitude to counts (ADU), using the mean flux/magnitude scaling factor.
*
* @param[in] magnitude The AB magnitude to convert.
*
* @return Source counts in ADU.
*/
double magnitudeToCounts(double magnitude) const;
/**
* @brief Get the mean flux/magnitude zero point.
*
* This value is defined, for counts at (x,y), such that:
* getFluxMag0() * counts * computeScaledZeroPoint()(x,y) = countsToMaggies(counts, (x,y))
*
* @see PhotoCalib::computeScaledZeroPoint()
*
* @return The flux magnitude zero point.
*/
double getFluxMag0() const;
/**
* @brief Calculates the spatially-variable zero point, normalized by the mean in the valid domain.
*
* This value is defined, for counts at (x,y), such that:
* getFluxMag0() * counts * computeScaledZeroPoint()(x,y) = countsToMaggies(counts, (x,y))
*
* @see PhotoCalib::getFluxMag0()
*
* @return The normalized spatially-variable zero point.
*/
std::shared_ptr<lsst::afw::math::BoundedField> computeScaledZeroPoint() const;
/**
* @brief Calculates the scaling between this PhotoCalib and another PhotoCalib.
*
* With:
* c = counts at position (x,y)
* this = this PhotoCalib
* other = other PhotoCalib
* return = BoundedField returned by this method
* the return value from this method is defined as:
* this.countsToMaggies(c, (x,y)) * return(x, y) = other.countsToMaggies(c, (x,y))
*
* @param[in] other The PhotoCalib to scale to.
*
* @return The BoundedField as defined above.
*/
std::shared_ptr<lsst::afw::math::BoundedField> computeScalingTo(std::shared_ptr<PhotoCalib> other) const;
/**
* @brief Compare two PhotoCalibs for equality.
*
* @param rhs The PhotoCalib to compare to.
*
* @return True if both PhotoCalibs have identical representations.
*/
bool operator==(PhotoCalib const& rhs) const;
/**
* @brief Compare two PhotoCalibs for non-equality.
*
* @param rhs The PhotoCalib to compare to.
*
* @return True if both PhotoCalibs do not have identical representations.
*/
bool operator!=(PhotoCalib const& rhs) const { return !(*this == rhs); }
private:
std::shared_ptr<lsst::afw::math::BoundedField> _zeroPoint;
double _fluxMag0Sigma;
// The "mean" zero point, defined as the geometric mean of _zeroPoint evaluated over _zeroPoint's bbox.
// Computed on instantiation as a convinience.
// Also, the actual zeroPoint for a spatially-constant calibration.
double _fluxMag0;
}
@TallJimbo

This comment has been minimized.

Copy link

@TallJimbo TallJimbo commented Feb 2, 2017

A few minor stylistic things:

  • Pass string and SourceCatalog by const reference.
  • I think all of these methods should be const.

We should take this opportunity to make throwOnNegativeFlux a non-static data member. I think it's mostly used via a Python context manager, which should still be able to behave the same way.

I think the big question is how to handle Calibs with conversions that are not spatially variable. I think we'll try to measure on images that don't have spatially-varying photometric calibration whenever possible, by re-flat-fielding the images with the spatially-varying part before we measure; I'm hoping we'll be able to construct all of our coadds this way. So there will still be a big role for non-spatially varying Calibs to play, and I'd like to avoid having to pass points to those when it's not needed. The cleanest solution might inheritance of some kind, but that's potentially tricky: FixedCalib is a SpatiallyVaryingCalib if and only if both are immutable (which we may or may not want anyway), and having them inherit from a common base class implies the base class doesn't have a usable interface by itself. I suppose we could do SpatiallyVaryingCalib is a FixedCalib if we pretend that the converters that don't take points actually do the same thing, which is a lie (but maybe an acceptable one).

@parejkoj

This comment has been minimized.

Copy link
Owner Author

@parejkoj parejkoj commented Feb 2, 2017

Thanks for the comments. I hadn't thought at all about constness/references, this is just a sketch of the API.

I've been thinking about the question of how to make what you're calling FixedCalib easier to work with. I think we can do it with this same object, by adding a setFluxMag0-type of method that just makes a constant BoundedField. The interface above has countsToMagnitude(double counts), which should behave identically to the current Calib in the case of a constant BoundedField.

@r-owen

This comment has been minimized.

Copy link

@r-owen r-owen commented Feb 2, 2017

Do you want boundedField(x,y) * counts * fluxMag0 = flux (in maggies) (as you have now) or fluxMag0(x,y) * counts = flux (in maggies) where fluxMag0 is a BoundedField? In your design boundedField(x,y) is a dimensionless multiplicative correction on the constant fluxMag0. I think I prefer the latter. The meaning of boundedField and fluxMag0 are both difficult to interpret, especially if the former is not very nearly 1 at most points. If you can come up with a better name for boundedField that might ameliorate the problem. Nonetheless, fluxMag0 as the bounded field is simplest to understand.

As to Jim's points...rather than a separate class, I'd be tempted to have a boundedField that is a constant and have users able to ask if the bounded field model is intrinsically constant ( e.g. ConstantBoundedField, not a spatially varying model whose spatial coefficients happen to be zero). Then if code really cared it could ask. Though in most cases I imagine the overhead of providing a position would be small.

I think it partly comes down to this: is it useful for a spatially varying Calib to offer an interface that does not accept points? I hope so. Doing so implies that we usually expect the spatial variation to be small enough that positional correction is only required for extremely high accuracy work.

@r-owen

This comment has been minimized.

Copy link

@r-owen r-owen commented Feb 2, 2017

I'm uneasy about the catalog interface. There are two flavors of countsToMagnitude: one returns a vector of <mag, mag err> and the other sets a field. I assume you did it this way because it can be a headache adding fields to a catalog in order to set them, and sometimes it suffices to just get a vector of data. Also, if the catalog is not contiguous, the vector is more efficient to work with in Python.

If we can get away with it I would drop the vector version, because a vector is difficult to keep in synch with a catalog. That said, there would be a lot less suffering if it was easier to add fields, or the python wrapper this worked with astropy tables. Jim, what do you think? This sort of thing comes up a lot, e.g. in star selectors.

John: are you intending to fill out the magnitudesToCounts methods to match the countsToMagnitudes, or is intentional that there is just one spatially non-varying version?

@parejkoj

This comment has been minimized.

Copy link
Owner Author

@parejkoj parejkoj commented Feb 3, 2017

I've updated the interface with constructors, countsToMaggies methods, the countsTo* "output array" versions and the exception context manager. I think most of the necessary methods are now described. More comments welcome.

Responding to Russell: it was Jim's suggestion to break it up into a spatial term, and a constant term (though I'd been thinking along those lines). However, we could hide that aspect of it as an implementation detail or just compute and cache the constant zero point for when we need it, while keeping the double fluxMag0 constructor (which builds a constant BoundedField internally). I could be convinced about this, but I feel like explicitly separating the spatial and constant parts is conceptually useful. It could be useful in persistence too, because we can then save the constant part to a FITS header (as we do now), and then do with the BoundedField as we like.

... is it useful for a spatially varying Calib to offer an interface that does not accept points? I hope so. Doing so implies that we usually expect the spatial variation to be small enough that positional correction is only required for extremely high accuracy work.

Actually, this suggests that we should be able to separate the spatial and constant terms in a sensible manner. I'd be surprised if that caused a problem in practice, worries about someone doing something non-sensible aside.

I think there would be only the one magnitudeToCounts: Jim says it's really only used for reference catalogs. If we think we need a magnitudesToCounts(counts, point) etc., I can certainly add them.

I'll wait to see what Jim has to say about the catalog interface. I can certainly see utility in being able to just get an array out, but it does clutter the interface a bit.

@TallJimbo

This comment has been minimized.

Copy link

@TallJimbo TallJimbo commented Feb 3, 2017

Some C++/immutability things:

  • Maybe instead of throwOnNegativeFlux, we could add a boolean argument to all of these? Now that we can easily do kwargs via pybind11, I think that'll be much cleaner, and it'll let the class be completely immutable, which is fairly important.
  • If the class is fully immutable, I think we can delete the copy and move constructors as well, since there won't be any point to copying or moving, as we'll be passing these by shared_ptr and you can just copy/move the shared_ptr since there's no point to having your own copy.
  • We should pass and hold the BoundedField by shared_ptr. BoundedFields are already immutable, so there'll be no need to copy when we take ownership, and it'll be a polymorphic object so we have to hold it by some kind of pointer.
  • Instead of individual functions that both take output arrays and return new ones, I think we should have a version that takes output arrays and a version that returns a ones for each case. In Python, that will naturally appear like the same function with an optional keyword arg, but it'll be more natural in C++.
  • I'm not sure operator== is going to be useful here, except for in test code (assuming it doesn't do fuzzy comparisons, which would probably belong in another method). It might be useful enough there to leave it, but it could also go in a test utility module.

More meaning comments:

  • I do like being able to just get arrays for magnitudes. Easier field-adding in afw.table isn't coming anytime soon, and the array-returning version is the one I'd expect us to use in analysis tools, even if we're using astropy tables for that (it's pretty straightforward now to make an astropy table from an afw catalog, and then add a pair of numpy array columns to that). I actually imagine the version that sets other fields in a SourceCatalog will be used less, but I still think it's worth keeping for more pipeline-oriented use where we set the schema up well in advance.
  • I assume the _fluxMag0 and _fluxMag0Sigma data members were supposed to be double, not BoundedField?
  • I think it's important that this interface doesn't demand that we separate the calibration into a constant fluxMag0 and spatially-varying correction. But I do think that's the representation that will turn out to be most convenient and probably efficient: I think we need to d exactly a o that decomposition for persistence, so we can be the best constant approximation to the calib in a single number in the header, and I think we'll frequently want to multiply images by the spatially-varying correction, creating an image that has an exactly constant Calib with no spatially-varying part. I think it really is an implementation detail whether we represent a spatially-constant Calib via a null BoundedField pointer or a ConstantBoundedField object (but I suspect the nullptr is simpler and no less efficient).
  • As per my last point, I think we want a shared_ptr<BoundedField> computeScalingTo(shared_ptr<Calib> other) const method that returns a the BoundedField that should be used to multiply an image with Calib == this to create an image with Calib == other. A form that returns a pair of shared_ptr<BoundedField>, shared_ptr<Calib> and takes no arguments would also be useful; I'm thinking this would create a spatially-constant Calib with just fluxMag0.
@parejkoj

This comment has been minimized.

Copy link
Owner Author

@parejkoj parejkoj commented Feb 4, 2017

Thanks for the follow-up, Jim. I also chatted with Russell, and have changed the constructor API to take either a constant, or a boundingField and optional BoundingBox over which it is valid, and over which the mean zero point should be computed.

I've also taken care of I think all the rest of your suggestions, though I may not have correctly understood your last one re:computeScalingTo.

Note that I did define two different countsToMagnitudes(SourceCatalog), one that returning the arrays, the other returning void, I just stuck them together and tweaked the docstring. Maybe there's a better way to do that that doesn't result in more already-over-duplicated docstrings?

@TallJimbo

This comment has been minimized.

Copy link

@TallJimbo TallJimbo commented Feb 7, 2017

The "bounded" part of BoundedField is that it's a field that knows the bounding box over which it is valid, so I don't think it make sense to have an overload that specifies both a BoundedField and a bounding box. It probably does make sense to add a BoundedField method that returns a new one with a smaller bounding box.

There are some Doxygen tricks that let you share text over multiple functions, but it's probably not worth bothering worth now. I would recommend that you make the output array arguments separate arguments rather than a pair, though.

Your math for computeScalingTo matches my intent, though it might best to rename the second overload since there's no argument for the "To" to refer to. Maybe call it "decompose"?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.