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.

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-35265: Reduce usage of MeasurementError #213

Merged
merged 6 commits into from
Jun 23, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 6 additions & 2 deletions src/LocalBackground.cc
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,9 @@ void LocalBackgroundAlgorithm::measure(afw::table::SourceRecord& measRecord,
// Define pixels in annulus
auto const psf = exposure.getPsf();
if (!psf) {
throw LSST_EXCEPT(MeasurementError, NO_PSF.doc, NO_PSF.number);
_flagHandler.setValue(measRecord, NO_PSF.number, true);
_flagHandler.setValue(measRecord, LocalBackgroundAlgorithm::FAILURE.number, true);
return;
}
float const psfSigma = psf->computeShape(center).getDeterminantRadius();

Expand Down Expand Up @@ -95,7 +97,9 @@ void LocalBackgroundAlgorithm::measure(afw::table::SourceRecord& measRecord,
}

if (values.size() == 0) {
throw LSST_EXCEPT(MeasurementError, NO_GOOD_PIXELS.doc, NO_GOOD_PIXELS.number);
_flagHandler.setValue(measRecord, NO_GOOD_PIXELS.number, true);
_flagHandler.setValue(measRecord, LocalBackgroundAlgorithm::FAILURE.number, true);
return;
}

// Measure the background
Expand Down
8 changes: 6 additions & 2 deletions src/NaiveCentroid.cc
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,9 @@ void NaiveCentroidAlgorithm::measure(afw::table::SourceRecord& measRecord,
y -= image.getY0();

if (x < 1 || x >= image.getWidth() - 1 || y < 1 || y >= image.getHeight() - 1) {
throw LSST_EXCEPT(MeasurementError, EDGE.doc, EDGE.number);
_flagHandler.setValue(measRecord, EDGE.number, true);
_flagHandler.setValue(measRecord, NaiveCentroidAlgorithm::FAILURE.number, true);
return;
}

ImageT::xy_locator im = image.xy_at(x, y);
Expand All @@ -76,7 +78,9 @@ void NaiveCentroidAlgorithm::measure(afw::table::SourceRecord& measRecord,
9 * _ctrl.background;

if (sum == 0.0) {
throw LSST_EXCEPT(MeasurementError, NO_COUNTS.doc, NO_COUNTS.number);
_flagHandler.setValue(measRecord, NO_COUNTS.number, true);
_flagHandler.setValue(measRecord, NaiveCentroidAlgorithm::FAILURE.number, true);
return;
}

double const sum_x = -im(-1, 1) + im(1, 1) + -im(-1, 0) + im(1, 0) + -im(-1, -1) + im(1, -1);
Expand Down
4 changes: 3 additions & 1 deletion src/PsfFlux.cc
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,9 @@ void PsfFluxAlgorithm::measure(afw::table::SourceRecord& measRecord,
->clippedTo(exposure.getMaskedImage().getMask()->getBBox()));
}
if (fitRegion.getArea() == 0) {
throw LSST_EXCEPT(MeasurementError, NO_GOOD_PIXELS.doc, NO_GOOD_PIXELS.number);
_flagHandler.setValue(measRecord, NO_GOOD_PIXELS.number, true);
_flagHandler.setValue(measRecord, FAILURE.number, true);
return;
timj marked this conversation as resolved.
Show resolved Hide resolved
}
typedef afw::detection::Psf::Pixel PsfPixel;
// SpanSet::flatten returns a new ndarray::Array, which must stay in scope
Expand Down
109 changes: 53 additions & 56 deletions src/SdssCentroid.cc
Original file line number Diff line number Diff line change
Expand Up @@ -128,15 +128,15 @@ float astrom_errors(float skyVar, // variance of pixels at the sky level
*/

template <typename ImageXy_locatorT, typename VarImageXy_locatorT>
void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
double *dxc, // output; error in xCenter
double *yCenter, // output; y-position of object
double *dyc, // output; error in yCenter
double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
ImageXy_locatorT im, // Locator for the pixel values
VarImageXy_locatorT vim, // Locator for the image containing the variance
double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
FlagHandler flagHandler) {
int doMeasureCentroidImpl(double *xCenter, // output; x-position of object
double *dxc, // output; error in xCenter
double *yCenter, // output; y-position of object
double *dyc, // output; error in yCenter
double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
ImageXy_locatorT im, // Locator for the pixel values
VarImageXy_locatorT vim, // Locator for the image containing the variance
double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
FlagHandler flagHandler) {
/*
* find a first quadratic estimate
*/
Expand All @@ -146,25 +146,17 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-positio
double const sy = 0.5 * (im(0, 1) - im(0, -1));

if (d2x == 0.0 || d2y == 0.0) {
throw LSST_EXCEPT(MeasurementError, SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.doc,
SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number);
return SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number;
}
if (d2x < 0.0 || d2y < 0.0) {
throw LSST_EXCEPT(MeasurementError,
SdssCentroidAlgorithm::NOT_AT_MAXIMUM.doc +
(boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
SdssCentroidAlgorithm::NOT_AT_MAXIMUM.number);
return SdssCentroidAlgorithm::NOT_AT_MAXIMUM.number;
}

double const dx0 = sx / d2x;
double const dy0 = sy / d2y; // first guess

if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
throw LSST_EXCEPT(
MeasurementError,
SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE.doc +
(boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE.number);
return SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number;
}

double vpk = im(0, 0) + 0.5 * (sx * dx0 + sy * dy0); // height of peak in image
Expand Down Expand Up @@ -238,18 +230,19 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-positio

*sizeX2 = tauX2; // return the estimates of the (object size)^2
*sizeY2 = tauY2;
return 0;
}

template <typename MaskedImageXy_locatorT>
void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
double *dxc, // output; error in xCenter
double *yCenter, // output; y-position of object
double *dyc, // output; error in yCenter
double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
double *peakVal, // output; peak of object
MaskedImageXy_locatorT mim, // Locator for the pixel values
double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
bool negative, FlagHandler flagHandler) {
int doMeasureCentroidImpl(double *xCenter, // output; x-position of object
double *dxc, // output; error in xCenter
double *yCenter, // output; y-position of object
double *dyc, // output; error in yCenter
double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
double *peakVal, // output; peak of object
MaskedImageXy_locatorT mim, // Locator for the pixel values
double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
bool negative, FlagHandler flagHandler) {
/*
* find a first quadratic estimate
*/
Expand All @@ -259,25 +252,17 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-positio
double const sy = 0.5 * (mim.image(0, 1) - mim.image(0, -1));

if (d2x == 0.0 || d2y == 0.0) {
throw LSST_EXCEPT(MeasurementError, SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.doc,
SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number);
return SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number;
}
if ((!negative && (d2x < 0.0 || d2y < 0.0)) || (negative && (d2x > 0.0 || d2y > 0.0))) {
throw LSST_EXCEPT(MeasurementError,
SdssCentroidAlgorithm::NOT_AT_MAXIMUM.doc +
(boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
SdssCentroidAlgorithm::NOT_AT_MAXIMUM.number);
return SdssCentroidAlgorithm::NOT_AT_MAXIMUM.number;
}

double const dx0 = sx / d2x;
double const dy0 = sy / d2y; // first guess

if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
throw LSST_EXCEPT(
MeasurementError,
SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.doc +
(boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE.number);
return SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number;
}

double vpk = mim.image(0, 0) + 0.5 * (sx * dx0 + sy * dy0); // height of peak in image
Expand Down Expand Up @@ -354,12 +339,13 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-positio
*sizeY2 = tauY2;

*peakVal = vpk;
return 0;
}

template <typename MaskedImageT>
std::pair<MaskedImageT, double> smoothAndBinImage(std::shared_ptr<afw::detection::Psf const> psf, int const x,
const int y, MaskedImageT const &mimage, int binX, int binY,
FlagHandler _flagHandler) {
std::tuple<MaskedImageT, double, int> smoothAndBinImage(std::shared_ptr<afw::detection::Psf const> psf, int const x,
const int y, MaskedImageT const &mimage, int binX, int binY,
FlagHandler _flagHandler) {
geom::Point2D const center(x + mimage.getX0(), y + mimage.getY0());
afw::geom::ellipses::Quadrupole const &shape = psf->computeShape(center);
double const smoothingSigma = shape.getDeterminantRadius();
Expand All @@ -378,12 +364,10 @@ std::pair<MaskedImageT, double> smoothAndBinImage(std::shared_ptr<afw::detection

// image to smooth, a shallow copy
std::shared_ptr<MaskedImageT> subImage;
try {
subImage.reset(new MaskedImageT(mimage, bbox, afw::image::LOCAL));
} catch (pex::exceptions::LengthError &err) {
throw LSST_EXCEPT(MeasurementError, SdssCentroidAlgorithm::EDGE.doc,
SdssCentroidAlgorithm::EDGE.number);
if (!mimage.getBBox(afw::image::LOCAL).contains(bbox)) {
return std::make_tuple(mimage, 0, SdssCentroidAlgorithm::EDGE.number);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want to remove this try/catch as well, I think the replacement is:

if (!mimage.getBBox(afw::image::LOCAL).contains(bbox)) {
    return std::make_tuple(mimage, 0, SdssCentroidAlgorithm::EDGE.number);
}
subImage.reset(new MaskedImageT(mimage, bbox, afw::image::LOCAL));

}
subImage.reset(new MaskedImageT(mimage, bbox, afw::image::LOCAL));
std::shared_ptr<MaskedImageT> binnedImage = afw::math::binImage(*subImage, binX, binY, afw::math::MEAN);
binnedImage->setXY0(subImage->getXY0());
// image to smooth into, a deep copy.
Expand All @@ -397,7 +381,7 @@ std::pair<MaskedImageT, double> smoothAndBinImage(std::shared_ptr<afw::detection
*smoothedImage.getVariance() *= binX * binY * nEffective; // We want the per-pixel variance, so undo the
// effects of binning and smoothing

return std::make_pair(smoothedImage, smoothingSigma);
return std::make_tuple(smoothedImage, smoothingSigma, 0);
}

} // end anonymous namespace
Expand Down Expand Up @@ -436,7 +420,9 @@ void SdssCentroidAlgorithm::measure(afw::table::SourceRecord &measRecord,
int const y = image.positionToIndex(center.getY(), afw::image::Y).first;

if (!image.getBBox().contains(geom::Extent2I(x, y) + image.getXY0())) {
throw LSST_EXCEPT(meas::base::MeasurementError, EDGE.doc, EDGE.number);
_flagHandler.setValue(measRecord, EDGE.number, true);
_flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
return;
}

// Algorithm uses a least-squares fit (implemented via a convolution) to a symmetrized PSF model.
Expand All @@ -449,19 +435,30 @@ void SdssCentroidAlgorithm::measure(afw::table::SourceRecord &measRecord,
int binY = 1;
double xc = 0., yc = 0., dxc = 0., dyc = 0.; // estimated centre and error therein
for (int binsize = 1; binsize <= _ctrl.binmax; binsize *= 2) {
std::pair<MaskedImageT, double> result =
smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
MaskedImageT const smoothedImage = result.first;
double const smoothingSigma = result.second;
std::tuple<MaskedImageT, double, int> smoothResult =
smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
int errorFlag = std::get<2>(smoothResult);
if (errorFlag > 0) {
_flagHandler.setValue(measRecord, errorFlag, true);
_flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
return;
}
MaskedImageT const smoothedImage = std::get<0>(smoothResult);
double const smoothingSigma = std::get<1>(smoothResult);

MaskedImageT::xy_locator mim =
smoothedImage.xy_at(smoothedImage.getWidth() / 2, smoothedImage.getHeight() / 2);

double sizeX2, sizeY2; // object widths^2 in x and y directions
double peakVal; // peak intensity in image

doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim, smoothingSigma, negative,
_flagHandler);
errorFlag = doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim, smoothingSigma, negative,
_flagHandler);
if (errorFlag > 0) {
_flagHandler.setValue(measRecord, errorFlag, true);
_flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
return;
}

if (binsize > 1) {
// dilate from the lower left corner of central pixel
Expand Down