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

DM-46720: Add uncalibrateImage method to afw.image.PhotoCalib #748

Merged
merged 6 commits into from
Oct 31, 2024
Merged

DM-46720: Add uncalibrateImage method to afw.image.PhotoCalib #748

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5d93fea
Add uncalibrateImage to PhotoCalib
parejkoj Oct 16, 2024
8a367cb
clang-format pass
parejkoj Oct 16, 2024
5273137
Fix incorrect subimage test
parejkoj Oct 17, 2024
adb6e90
Refactor to make flux/instFlux ratio more explicit
parejkoj Oct 30, 2024
e1fa458
Update PhotoCalib test values
parejkoj Oct 30, 2024
f753d4d
Use more useful/appropriate image and variance in test
parejkoj Oct 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions include/lsst/afw/image/PhotoCalib.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -350,6 +350,22 @@ class PhotoCalib final : public table::io::PersistableFacade<PhotoCalib>, public
MaskedImage<float> calibrateImage(MaskedImage<float> const &maskedImage,
bool includeScaleUncertainty = true) const;

/**
* Return a un-calibrated image, with pixel values in ADU (or whatever the original input to
* this photoCalib was).
*
* Mask pixels are propagated directly from the input image.
*
* @param maskedImage The masked image with pixel units of nJy to uncalibrate.
* @param includeScaleUncertainty Remove uncertainty from this calibration that was previously propagated
* into the variance plane. Must have the same value as the parameter of the same name when
* calibrateImage was called on this image.
*
* @return The uncalibrated masked image.
*/
MaskedImage<float> uncalibrateImage(MaskedImage<float> const &maskedImage,
bool includeScaleUncertainty = true) const;

/**
* Return a flux calibrated catalog, with new `_flux`, `_fluxErr`, `_mag`, and `_magErr` fields.
*
Expand Down
18 changes: 10 additions & 8 deletions python/lsst/afw/image/_photoCalib.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,11 @@ void declarePhotoCalib(lsst::cpputils::python::WrapperCollection &wrappers) {

/* Members - nanojansky */
cls.def("instFluxToNanojansky",
(double (PhotoCalib::*)(double, lsst::geom::Point<double, 2> const &) const) &
(double(PhotoCalib::*)(double, lsst::geom::Point<double, 2> const &) const) &
PhotoCalib::instFluxToNanojansky,
"instFlux"_a, "point"_a);
cls.def("instFluxToNanojansky",
(double (PhotoCalib::*)(double) const) & PhotoCalib::instFluxToNanojansky,
(double(PhotoCalib::*)(double) const) & PhotoCalib::instFluxToNanojansky,
"instFlux"_a);

cls.def("instFluxToNanojansky",
Expand All @@ -109,18 +109,18 @@ void declarePhotoCalib(lsst::cpputils::python::WrapperCollection &wrappers) {
"sourceCatalog"_a, "instFluxField"_a);

cls.def("instFluxToNanojansky",
(void (PhotoCalib::*)(afw::table::SourceCatalog &, std::string const &,
std::string const &) const) &
(void(PhotoCalib::*)(afw::table::SourceCatalog &, std::string const &,
std::string const &) const) &
PhotoCalib::instFluxToNanojansky,
"sourceCatalog"_a, "instFluxField"_a, "outField"_a);

/* Members - magnitudes */
cls.def("instFluxToMagnitude",
(double (PhotoCalib::*)(double, lsst::geom::Point<double, 2> const &) const) &
(double(PhotoCalib::*)(double, lsst::geom::Point<double, 2> const &) const) &
PhotoCalib::instFluxToMagnitude,
"instFlux"_a, "point"_a);
cls.def("instFluxToMagnitude",
(double (PhotoCalib::*)(double) const) & PhotoCalib::instFluxToMagnitude,
(double(PhotoCalib::*)(double) const) & PhotoCalib::instFluxToMagnitude,
"instFlux"_a);

cls.def("instFluxToMagnitude",
Expand All @@ -145,8 +145,8 @@ void declarePhotoCalib(lsst::cpputils::python::WrapperCollection &wrappers) {
"sourceCatalog"_a, "instFluxField"_a);

cls.def("instFluxToMagnitude",
(void (PhotoCalib::*)(afw::table::SourceCatalog &, std::string const &,
std::string const &) const) &
(void(PhotoCalib::*)(afw::table::SourceCatalog &, std::string const &,
std::string const &) const) &
PhotoCalib::instFluxToMagnitude,
"sourceCatalog"_a, "instFluxField"_a, "outField"_a);

Expand All @@ -171,6 +171,8 @@ void declarePhotoCalib(lsst::cpputils::python::WrapperCollection &wrappers) {

cls.def("calibrateImage", &PhotoCalib::calibrateImage, "maskedImage"_a,
"includeScaleUncertainty"_a = true);
cls.def("uncalibrateImage", &PhotoCalib::uncalibrateImage, "maskedImage"_a,
"includeScaleUncertainty"_a = true);

cls.def("calibrateCatalog",
py::overload_cast<afw::table::SourceCatalog const &,
Expand Down
65 changes: 57 additions & 8 deletions src/image/PhotoCalib.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,9 @@ int const SERIALIZATION_VERSION = 1;

double toNanojansky(double instFlux, double scale) { return instFlux * scale; }

double toMagnitude(double instFlux, double scale) { return cpputils::nanojanskyToABMagnitude(instFlux * scale); }
double toMagnitude(double instFlux, double scale) {
return cpputils::nanojanskyToABMagnitude(instFlux * scale);
}

double toInstFluxFromMagnitude(double magnitude, double scale) {
// Note: flux[nJy] / scale = instFlux[counts]
Expand Down Expand Up @@ -90,8 +92,33 @@ void toNanojanskyVariance(ndarray::Array<float const, 2, 1> const &instFlux,
auto eigenInstFluxVar = ndarray::asEigen<Eigen::ArrayXpr>(instFluxVar);
auto eigenInstFlux = ndarray::asEigen<Eigen::ArrayXpr>(instFlux);
auto eigenOut = ndarray::asEigen<Eigen::ArrayXpr>(out);
eigenOut = eigenFlux.square() *
(eigenInstFluxVar / eigenInstFlux.square() + (scaleErr / eigenFlux * eigenInstFlux).square());
eigenOut = eigenFlux.square() * (eigenInstFluxVar / eigenInstFlux.square() +
(scaleErr / (eigenFlux / eigenInstFlux)).square());
}

/**
* Compute the variance of an array of fluxes, for calculations on nJy calibrated MaskedImages.
*
* MaskedImage stores the variance instead of the standard deviation, so we can skip a sqrt().
* Usage in calibrateImage() is to compute instFlux (ADU) directly, so that the calibration scale is
* implictly `flux/instFlux` (and thus not passed as an argument).
*
* @param flux[in] The flux in nJy.
* @param fluxVar[in] The variance of the fluxes.
* @param scaleErr[in] The error on the calibration scale.
* @param instFlux[in] The instrumental fluxes calculated from the fluxes.
* @param out[out] The output array to fill with the instrumental variance values.
*/
void fromNanojanskyVariance(ndarray::Array<float const, 2, 1> const &flux,
ndarray::Array<float const, 2, 1> const &fluxVar, float scaleErr,
ndarray::Array<float const, 2, 1> const &instFlux,
ndarray::Array<float, 2, 1> out) {
auto eigenFlux = ndarray::asEigen<Eigen::ArrayXpr>(flux);
auto eigenFluxVar = ndarray::asEigen<Eigen::ArrayXpr>(fluxVar);
auto eigenInstFlux = ndarray::asEigen<Eigen::ArrayXpr>(instFlux);
auto eigenOut = ndarray::asEigen<Eigen::ArrayXpr>(out);
eigenOut = eigenInstFlux.square() *
(eigenFluxVar / eigenFlux.square() - (scaleErr / (eigenFlux / eigenInstFlux)).square());
}

double toMagnitudeErr(double instFlux, double instFluxErr, double scale, double scaleErr) {
Expand Down Expand Up @@ -280,6 +307,28 @@ MaskedImage<float> PhotoCalib::calibrateImage(MaskedImage<float> const &maskedIm
return result;
}

MaskedImage<float> PhotoCalib::uncalibrateImage(MaskedImage<float> const &maskedImage,
bool includeScaleUncertainty) const {
// Deep copy construct, as we're mutiplying in-place.
auto result = MaskedImage<float>(maskedImage, true);

if (_isConstant) {
*(result.getImage()) /= _calibrationMean;
} else {
_calibration->divideImage(*(result.getImage()), true); // only in the overlap region
}
if (includeScaleUncertainty) {
fromNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(),
_calibrationErr, result.getImage()->getArray(),
result.getVariance()->getArray());
} else {
fromNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(), 0,
result.getImage()->getArray(), result.getVariance()->getArray());
}

return result;
}

afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog,
std::vector<std::string> const &instFluxFields) const {
auto const &inSchema = catalog.getSchema();
Expand Down Expand Up @@ -399,8 +448,8 @@ class PhotoCalibSchema {

class PhotoCalibFactory : public table::io::PersistableFactory {
public:
std::shared_ptr<table::io::Persistable>
read(InputArchive const &archive, CatalogVector const &catalogs) const override {
std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
CatalogVector const &catalogs) const override {
table::BaseRecord const &record = catalogs.front().front();
PhotoCalibSchema const &keys = PhotoCalibSchema::get();
int version = getVersion(record);
Expand Down Expand Up @@ -482,7 +531,7 @@ class CalibFactory : public table::io::PersistableFactory {
int tableVersion = 1;
try {
catalogs.front().getSchema().find<double>(EXPTIME_FIELD_NAME);
} catch (pex::exceptions::NotFoundError const&) {
} catch (pex::exceptions::NotFoundError const &) {
tableVersion = CALIB_TABLE_CURRENT_VERSION;
}

Expand All @@ -493,8 +542,8 @@ class CalibFactory : public table::io::PersistableFactory {
table::BaseRecord const &record = catalogs.front().front();

double calibration = cpputils::referenceFlux / record.get(keys.fluxMag0);
double calibrationErr =
cpputils::referenceFlux * record.get(keys.fluxMag0Err) / std::pow(record.get(keys.fluxMag0), 2);
double calibrationErr = cpputils::referenceFlux * record.get(keys.fluxMag0Err) /
std::pow(record.get(keys.fluxMag0), 2);
return std::make_shared<PhotoCalib>(calibration, calibrationErr);
}

Expand Down
45 changes: 30 additions & 15 deletions tests/test_photoCalib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,17 +84,17 @@ def setUp(self):
self.pointYShift = lsst.geom.Point2D(0, -10)
self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(-100, -100), lsst.geom.Point2I(100, 100))

# calibration and instFlux1 are selected to produce calibrated flux of 1.
self.calibration = 1e-3
self.calibrationErr = 1e-4
self.instFlux1 = 1000.
self.instFluxErr1 = 10.
self.flux1 = 1.0
self.calibration = 1000.0
# A 1% error on the calibration.
self.calibrationErr = 10.0
self.instFlux1 = 1.0
self.instFluxErr1 = 0.1
self.flux1 = 1000.0 # nJy
self.mag1 = (self.flux1*u.nJy).to_value(u.ABmag)

# useful reference points: 575.44 nJy ~= 24.5 mag, 3630.78 * 10^9 nJy ~= 0 mag
self.flux2 = 575.44
self.instFlux2 = self.instFlux1*self.flux2
self.flux2 = 575.44 * self.flux1
self.instFlux2 = self.flux2/self.calibration
self.mag2 = (self.flux2*u.nJy).to_value(u.ABmag)

self.schema = lsst.afw.table.SourceTable.makeMinimalSchema()
Expand Down Expand Up @@ -185,7 +185,7 @@ def _testPhotoCalibCenter(self, photoCalib, calibrationErr):
self.calibration,
self.flux1)
result = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1)
self.assertEqual(1, result.value)
self.assertEqual(self.flux1, result.value)
self.assertFloatsAlmostEqual(errFlux1, result.error)
result = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1, self.point0)
self.assertFloatsAlmostEqual(self.flux1, result.value)
Expand Down Expand Up @@ -458,14 +458,15 @@ def testLinearXBoundedField(self):
def testComputeScaledCalibration(self):
photoCalib = lsst.afw.image.PhotoCalib(self.calibration, bbox=self.bbox)
scaledCalib = lsst.afw.image.PhotoCalib(photoCalib.computeScaledCalibration())
self.assertEqual(1, scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean())
self.assertEqual(self.flux1,
scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean())
self.assertEqual(photoCalib.instFluxToNanojansky(self.instFlux1),
scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean())

photoCalib = lsst.afw.image.PhotoCalib(self.constantCalibration)
scaledCalib = lsst.afw.image.PhotoCalib(photoCalib.computeScaledCalibration())

self.assertEqual(1, scaledCalib.instFluxToNanojansky(self.instFlux1*self.calibration))
self.assertEqual(self.flux1, scaledCalib.instFluxToNanojansky(self.instFlux1*self.calibration))
self.assertEqual(photoCalib.instFluxToNanojansky(self.instFlux1),
scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean())

Expand Down Expand Up @@ -577,9 +578,11 @@ def setupImage(self):
dim = (5, 6)
npDim = (dim[1], dim[0]) # numpy and afw have a different x/y order
sigma = 10.0
image = np.random.normal(loc=1000.0, scale=sigma, size=npDim).astype(np.float32)
# An image with increasing pixel values, to more easily check the
# values in specific calculations.
image = np.arange(dim[0]*dim[1]).astype(np.float32).reshape(npDim) + 1000
mask = np.zeros(npDim, dtype=np.int32)
variance = (np.random.normal(loc=0.0, scale=sigma, size=npDim).astype(np.float32))**2
variance = np.ones(npDim, dtype=np.float32)*sigma
maskedImage = lsst.afw.image.makeMaskedImageFromArrays(image, mask, variance)
maskedImage.mask[0, 0] = True # set one mask bit to check propagation of mask bits.

Expand All @@ -594,11 +597,15 @@ def testCalibrateImageConstant(self):
photoCalib = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr)
result = photoCalib.calibrateImage(maskedImage)
self.assertMaskedImagesAlmostEqual(expect, result)
uncalibrated = photoCalib.uncalibrateImage(result)
self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated)

# same test, but without using the calibration error
expect = makeCalibratedMaskedImageNoCalibrationError(image, mask, variance, outImage)
result = photoCalib.calibrateImage(maskedImage, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(expect, result)
uncalibrated = photoCalib.uncalibrateImage(result, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated)

def testCalibrateImageNonConstant(self):
"""Test a spatially-varying calibration."""
Expand All @@ -612,11 +619,15 @@ def testCalibrateImageNonConstant(self):
photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr)
result = photoCalib.calibrateImage(maskedImage)
self.assertMaskedImagesAlmostEqual(expect, result)
uncalibrated = photoCalib.uncalibrateImage(result)
self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated)

# same test, but without using the calibration error
expect = makeCalibratedMaskedImageNoCalibrationError(image, mask, variance, outImage)
result = photoCalib.calibrateImage(maskedImage, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(expect, result)
uncalibrated = photoCalib.uncalibrateImage(result, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated)

def testCalibrateImageNonConstantSubimage(self):
"""Test a non-constant calibration on a sub-image, to ensure we're
Expand All @@ -634,11 +645,15 @@ def testCalibrateImageNonConstantSubimage(self):
photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr)
result = photoCalib.calibrateImage(subImage)
self.assertMaskedImagesAlmostEqual(expect.subset(subBox), result)
uncalibrated = photoCalib.uncalibrateImage(result)
self.assertMaskedImagesAlmostEqual(subImage, uncalibrated)

# same test, but without using the calibration error
expect = makeCalibratedMaskedImageNoCalibrationError(image, mask, variance, outImage)
result = photoCalib.calibrateImage(maskedImage, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(expect, result)
result = photoCalib.calibrateImage(subImage, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(expect.subset(subBox), result)
uncalibrated = photoCalib.uncalibrateImage(result, includeScaleUncertainty=False)
self.assertMaskedImagesAlmostEqual(subImage, uncalibrated)

def testNonPositiveMeans(self):
# no negative calibrations
Expand Down
Loading