Skip to content

Commit

Permalink
Set specific MeasurementErrors for SdssCentroidAlgorithm
Browse files Browse the repository at this point in the history 
These new flags correspond to the errors in the original algorithm:

not_at_maximum when the centroid is at a point lower than surround
no_2nd_derivative when the 2nd derivative around the point is <= 0
almost_not_at_2nd_derivate when it is too small around the point
  • Loading branch information
@pgee2000
pgee authored and pgee2000 committed Mar 29, 2015
1 parent 671d8eb commit 6f11e7d
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 56 deletions.
4 changes: 3 additions & 1 deletion include/lsst/meas/base/SdssCentroid.h
View file
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
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
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

0 comments on commit 6f11e7d

Please sign in to comment.