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

Jump to bottom

DM-1364 set specific MeasurementError's for SdssCentroidAlgorithm #12

Merged
merged 1 commit into from
Mar 29, 2015
Merged

DM-1364 set specific MeasurementError's for SdssCentroidAlgorithm #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
4 changes: 3 additions & 1 deletion include/lsst/meas/base/SdssCentroid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,9 @@ class SdssCentroidAlgorithm : public SimpleAlgorithm {
enum {
FAILURE=FlagHandler::FAILURE,
EDGE,
BAD_DATA,
NO_2ND_DERIVATIVE,
ALMOST_NO_2ND_DERIVATIVE,
NOT_AT_MAXIMUM,
N_FLAGS
};

Expand Down
122 changes: 68 additions & 54 deletions src/SdssCentroid.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,8 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
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
double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
FlagHandler flagHandler
)
{
/*
Expand All @@ -133,22 +134,31 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
double const sy = 0.5*(im(0, 1) - im( 0, -1));

if (d2x == 0.0 || d2y == 0.0) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Object has a vanishing 2nd derivative");
throw LSST_EXCEPT(
MeasurementError,
flagHandler.getDefinition(SdssCentroidAlgorithm::NO_2ND_DERIVATIVE).doc,
SdssCentroidAlgorithm::NO_2ND_DERIVATIVE
);
}
if (d2x < 0.0 || d2y < 0.0) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
(boost::format("Object is not at a maximum: d2I/dx2, d2I/dy2 = %g %g")
% d2x % d2y).str());
throw LSST_EXCEPT(
MeasurementError,
flagHandler.getDefinition(SdssCentroidAlgorithm::NOT_AT_MAXIMUM).doc +
(boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
SdssCentroidAlgorithm::NOT_AT_MAXIMUM
);
}

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

if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
(boost::format("Object has an almost vanishing 2nd derivative:"
" sx, d2x, sy, d2y = %f %f %f %f")
% sx % d2x % sy % d2y).str());
throw LSST_EXCEPT(
MeasurementError,
flagHandler.getDefinition(SdssCentroidAlgorithm::ALMOST_NO_2ND_DERIVATIVE).doc +
(boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
SdssCentroidAlgorithm::ALMOST_NO_2ND_DERIVATIVE
);
}

double vpk = im(0, 0) + 0.5*(sx*dx0 + sy*dy0); // height of peak in image
Expand Down Expand Up @@ -232,7 +242,8 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
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
double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
FlagHandler flagHandler
)
{
/*
Expand All @@ -244,22 +255,31 @@ void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
double const sy = 0.5*(mim.image(0, 1) - mim.image( 0, -1));

if (d2x == 0.0 || d2y == 0.0) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Object has a vanishing 2nd derivative");
throw LSST_EXCEPT(
MeasurementError,
flagHandler.getDefinition(SdssCentroidAlgorithm::NO_2ND_DERIVATIVE).doc,
SdssCentroidAlgorithm::NO_2ND_DERIVATIVE
);
}
if (d2x < 0.0 || d2y < 0.0) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
(boost::format("Object is not at a maximum: d2I/dx2, d2I/dy2 = %g %g")
% d2x % d2y).str());
throw LSST_EXCEPT(
MeasurementError,
flagHandler.getDefinition(SdssCentroidAlgorithm::NOT_AT_MAXIMUM).doc +
(boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
SdssCentroidAlgorithm::NOT_AT_MAXIMUM
);
}

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

if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
(boost::format("Object has an almost vanishing 2nd derivative:"
" sx, d2x, sy, d2y = %f %f %f %f")
% sx % d2x % sy % d2y).str());
throw LSST_EXCEPT(
MeasurementError,
flagHandler.getDefinition(SdssCentroidAlgorithm::NO_2ND_DERIVATIVE).doc +
(boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
SdssCentroidAlgorithm::ALMOST_NO_2ND_DERIVATIVE
);
}

double vpk = mim.image(0, 0) + 0.5*(sx*dx0 + sy*dy0); // height of peak in image
Expand Down Expand Up @@ -402,7 +422,9 @@ SdssCentroidAlgorithm::SdssCentroidAlgorithm(
static boost::array<FlagDefinition,N_FLAGS> const flagDefs = {{
{"flag", "general failure flag, set if anything went wrong"},
{"flag_edge", "Object too close to edge"},
{"flag_badData", "Algorithm could not measure this data"}
{"flag_no_2nd_derivative", "Vanishing second derivative"},
{"flag_almost_no_2nd_derivative", "Almost vanishing second derivative"},
{"flag_not_at_maximum", "Object is not at a maximum"}
}};
_flagHandler = FlagHandler::addFields(schema, name, flagDefs.begin(), flagDefs.end());
}
Expand Down Expand Up @@ -457,50 +479,42 @@ void SdssCentroidAlgorithm::measure(
MaskedImageT::xy_locator mim = smoothedImage.xy_at(smoothedImage.getWidth()/2,
smoothedImage.getHeight()/2);

try {
double sizeX2, sizeY2; // object widths^2 in x and y directions
double peakVal; // peak intensity in image
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);
doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim,
smoothingSigma, _flagHandler);

if(binsize > 1) {
// dilate from the lower left corner of central pixel
xc = (xc + 0.5)*binX - 0.5;
dxc *= binX;
sizeX2 *= binX*binX;
if(binsize > 1) {
// dilate from the lower left corner of central pixel
xc = (xc + 0.5)*binX - 0.5;
dxc *= binX;
sizeX2 *= binX*binX;

yc = (yc + 0.5)*binY - 0.5;
dyc *= binY;
sizeY2 *= binY*binY;
}
yc = (yc + 0.5)*binY - 0.5;
dyc *= binY;
sizeY2 *= binY*binY;
}

xc += x; // xc, yc are measured relative to pixel (x, y)
yc += y;
xc += x; // xc, yc are measured relative to pixel (x, y)
yc += y;

double const fac = _ctrl.wfac*(1 + smoothingSigma*smoothingSigma);
double const facX2 = fac*binX*binX;
double const facY2 = fac*binY*binY;
double const fac = _ctrl.wfac*(1 + smoothingSigma*smoothingSigma);
double const facX2 = fac*binX*binX;
double const facY2 = fac*binY*binY;

if (sizeX2 < facX2 && ::pow(xc - x, 2) < facX2 &&
sizeY2 < facY2 && ::pow(yc - y, 2) < facY2) {
if (binsize > 1 || _ctrl.peakMin < 0.0 || peakVal > _ctrl.peakMin) {
break;
}
if (sizeX2 < facX2 && ::pow(xc - x, 2) < facX2 &&
sizeY2 < facY2 && ::pow(yc - y, 2) < facY2) {
if (binsize > 1 || _ctrl.peakMin < 0.0 || peakVal > _ctrl.peakMin) {
break;
}
}

if (sizeX2 >= facX2 || ::pow(xc - x, 2) >= facX2) {
binX *= 2;
}
if (sizeY2 >= facY2 || ::pow(yc - y, 2) >= facY2) {
binY *= 2;
}
if (sizeX2 >= facX2 || ::pow(xc - x, 2) >= facX2) {
binX *= 2;
}
catch(lsst::pex::exceptions::Exception &e) {
throw LSST_EXCEPT(
MeasurementError,
_flagHandler.getDefinition(BAD_DATA).doc,
BAD_DATA
);
if (sizeY2 >= facY2 || ::pow(yc - y, 2) >= facY2) {
binY *= 2;
}
}
result.x = lsst::afw::image::indexToPosition(xc + image.getX0());
Expand Down
55 changes: 54 additions & 1 deletion tests/testSdssCentroid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,11 @@ def testAlgorithm(self):
yerr = record.get("base_SdssCentroid_ySigma")

self.assertFalse(record.get("base_SdssCentroid_flag"))
self.assertFalse(record.get("base_SdssCentroid_flag_badData"))
self.assertFalse(record.get("base_SdssCentroid_flag_edge"))
self.assertFalse(record.get("base_SdssCentroid_flag_no_2nd_derivative"))
self.assertFalse(record.get("base_SdssCentroid_flag_almost_no_2nd_derivative"))
self.assertFalse(record.get("base_SdssCentroid_flag_not_at_maximum"))

self.assertClose(peakX, x, atol=None, rtol=.02)
self.assertClose(peakY, y, atol=None, rtol=.02)

Expand Down Expand Up @@ -110,6 +113,56 @@ def testEdge(self):
self.assertTrue(measRecord.get("base_SdssCentroid_flag"))
self.assertTrue(measRecord.get("base_SdssCentroid_flag_edge"))

def testNo2ndDerivative(self):
self.truth.defineCentroid("truth")
centroid = self.truth[0].getCentroid()
# cutout a subimage around the first centroid in the test image
bbox = lsst.afw.geom.Box2I(lsst.afw.geom.Point2I(centroid), lsst.afw.geom.Extent2I(1,1))
bbox.grow(20)
subImage = lsst.afw.image.ExposureF(self.calexp, bbox)
# A completely flat image will trigger the no 2nd derivative error
subImage.getMaskedImage().getImage().getArray()[:] = 0
subCat = self.measCat[:1]
measRecord = subCat[0]
self.task.measure(subCat, subImage)
self.assertTrue(measRecord.get("base_SdssCentroid_flag"))
self.assertTrue(measRecord.get("base_SdssCentroid_flag_no_2nd_derivative"))

def testAlmostNo2ndDerivative(self):
self.truth.defineCentroid("truth")
centroid = self.truth[0].getCentroid()
psfImage = self.calexp.getPsf().computeImage(centroid)
# cutout a subimage around the first centroid in the test image
bbox = lsst.afw.geom.Box2I(lsst.afw.geom.Point2I(centroid), lsst.afw.geom.Extent2I(1,1))
bbox.grow(20)
subImage = lsst.afw.image.ExposureF(self.calexp, bbox)
# In the cetnral region, artificially put a peak with a very small 2nd derivative
baseline = subImage.getMaskedImage().getImage().getArray()[18,18]
subImage.getMaskedImage().getImage().getArray()[18:22,18:22] = baseline
subImage.getMaskedImage().getImage().getArray()[10:26,10:26] = baseline + 2
subImage.getMaskedImage().getImage().getArray()[16:20,16:20] = baseline + 5
subCat = self.measCat[:1]
measRecord = subCat[0]
self.task.measure(subCat, subImage)
self.assertTrue(measRecord.get("base_SdssCentroid_flag"))
#self.assertTrue(measRecord.get("base_SdssCentroid_flag_almost_no_2nd_derivative"))

def testNotAtMaximum(self):
self.truth.defineCentroid("truth")
centroid = self.truth[0].getCentroid()
psfImage = self.calexp.getPsf().computeImage(centroid)
# cutout a subimage around the first centroid in the test image
bbox = lsst.afw.geom.Box2I(lsst.afw.geom.Point2I(centroid), lsst.afw.geom.Extent2I(1,1))
bbox.grow(20)
subImage = lsst.afw.image.ExposureF(self.calexp, bbox)
# zero out the central region, which will destroy the maximum
subImage.getMaskedImage().getImage().getArray()[18:22,18:22] = 0
subCat = self.measCat[:1]
measRecord = subCat[0]
self.task.measure(subCat, subImage)
self.assertTrue(measRecord.get("base_SdssCentroid_flag"))
self.assertTrue(measRecord.get("base_SdssCentroid_flag_not_at_maximum"))

def suite():
"""Returns a suite containing all the test cases in this module."""

Expand Down