23 #include "ndarray/eigen.h" 32 namespace lsst {
namespace meas {
namespace base {
34 FlagDefinitionList flagDefinitions;
41 return flagDefinitions;
45 #define USE_WEIGHT 0 // zweight is only set, not used. It isn't even set if this is false 55 typedef Eigen::Matrix<double, FittedModel::NPARAM, FittedModel::NPARAM> Matrix;
56 typedef Eigen::Matrix<double, FittedModel::NPARAM, 1> Vector;
58 template<
typename PixelT>
59 explicit Fit2d(afw::image::Image<PixelT>
const& im) :
wide(32),
rast(
wide*
wide) {
60 ncols = im.getWidth();
61 nrows = im.getHeight();
118 static void dcalc2(Fit2d *fit,
119 Fit2d::Vector
const &el,
120 Fit2d::Matrix *
alpha,
135 }
else if (s <= 0.0) {
150 for (
int i = 0, nrast = fit->rast.size(); i != nrast; i++) {
151 double const dx = (fit->rast[i].x - x0)/s;
152 double const dy = (fit->rast[i].y - y0)/s;
153 double const d = hypot(dx, dy);
155 if (d >= fit->dc2zmin && d <= fit->
dc2zmax) {
156 double const arg = exp(-d*d/2.0);
157 double const funct = a*arg +
b;
160 df << arg, 1.0, a*arg*dx/
s, a*arg*dy/
s, a*arg*(dx*dx/s + dy*dy/
s);
162 double const r = fit->rast[i].zo - funct;
165 *alpha += df*df.transpose();
175 fit->stnew = sqrt(fit->chnew/nu);
182 static void curf2(Fit2d *fit) {
187 if (fit->iter == 0) {
188 dcalc2(fit, fit->elnew, &fit->alnew, &fit->benew);
189 if (fit->status != 0) {
192 fit->nref = fit->nin;
198 fit->chold = fit->chnew;
199 fit->stold = fit->stnew;
200 fit->elold = fit->elnew;
201 fit->sgold = fit->sgnew;
202 fit->beold = fit->benew;
203 fit->alold = fit->alnew;
208 for (
int poor = 0; poor != NPLIM; ++poor) {
210 if (fit->alnew(j, j) == 0.0) {
215 fit->alpha(j, k) = fit->alnew(j, k)/sqrt(fit->alnew(j, j)*fit->alnew(k, k));
217 fit->alpha(j, j) = 1.0 + fit->flamd;
223 fit->alpha.computeInverse();
225 Eigen::FullPivLU<Fit2d::Matrix> alphaLU(fit->alpha);
226 if (!alphaLU.isInvertible()) {
230 fit->alpha = alphaLU.inverse();
238 fit->elnew[j] += fit->benew[k]*fit->alpha(j,k)/sqrt(fit->alnew(j, j)*fit->alnew(k, k));
244 dcalc2(fit, fit->elnew, &fit->alnew, &fit->benew);
251 fit->sgnew[j] = fit->stnew*sqrt(fabs(fit->alpha(j, j)/fit->alnew(j, j)));
256 if (fit->status == 0.0 && fit->stnew <= fit->stold) {
257 fit->flamd = fit->flamd/fit->xlamd;
263 fit->flamd = 3.0*fit->flamd;
264 fit->chnew = fit->chold;
265 fit->stnew = fit->stold;
266 fit->elnew = fit->elold;
267 fit->sgnew = fit->sgold;
268 fit->benew = fit->beold;
269 fit->alnew = fit->alold;
277 static void cond2(Fit2d *fit) {
281 if (fit->iter <= 3) {
285 if (fit->flamd < fit->flamdok) {
292 if (fit->status < 0) {
299 if (fit->chnew <= 0.0) {
315 if (ex > fit->lost || ey > fit->lost) {
322 double const ratio = (fit->stnew - fit->stold)/fit->stnew;
323 if (ratio > fit->ratiomin) {
325 }
else if (fit->iter > fit->nitmax) {
336 template<
typename PixelT>
337 static void fg2(afw::image::Image<PixelT>
const& im,
338 double x0,
double y0,
344 int ix0 =
static_cast<int>(x0 + 0.5);
345 int iy0 =
static_cast<int>(y0 + 0.5);
346 if(ix0 < fit->
tooclose || im.getWidth() - ix0 < fit->tooclose ||
347 iy0 < fit->tooclose || im.getHeight() - iy0 < fit->tooclose) {
353 double peak = im(jx0, jy0);
357 double const w = fit->wide/2;
358 int xl =
static_cast<int>(ix0 - w + 0.5);
362 int xh =
static_cast<int>(xl + 2*w + 0.5);
363 if (xh > im.getWidth()) {
366 int yl =
static_cast<int>(iy0 - w + 0.5);
370 int yh =
static_cast<int>(yl + 2*w + 0.5);
371 if (yh > im.getHeight()) {
375 double sky = im(xl, yl);
377 for (
int y = yl;
y != yh; ++
y) {
378 for (
int x = xl;
x != xh; ++
x) {
379 fit->rast[put].x =
x;
380 fit->rast[put].y =
y;
381 fit->rast[put].zo = im(
x,
y);
383 fit->rast[put].zweight = 1.0;
385 if (im(
x,
y) < sky) {
388 double const ex =
x - ix0;
389 double const ey =
y - iy0;
390 double const er = hypot(ex, ey);
391 if (er <= fit->
lost) {
392 if (im(
x,
y) > peak) {
401 fit->rast.resize(put);
408 double xmax = xh - 1;
410 double ymax = yh - 1;
411 double const test = 0.5*(im(ix0, iy0) + sky);
412 for (
int x = ix0;
x < xh - 1; ++
x) {
413 if (im(
x + 1, iy0) <= test) {
414 double a = test - im(
x, iy0);
415 double b = im(
x + 1, iy0) - im(
x, iy0);
417 xmax =
x + ((b == 0.0) ? 0.5 : a/
b);
421 for (
int x = ix0;
x > 0; --
x) {
422 if (im(
x - 1, iy0) <= test) {
423 double a = test - im(
x, iy0);
424 double b = im(
x - 1, iy0) - im(
x, iy0);
426 xmin =
x - ((b == 0.0) ? 0.5 : a/
b);
430 for (
int y = iy0;
y < yh - 1; ++
y) {
431 if (im(ix0,
y + 1) <= test) {
432 double a = test - im(ix0,
y);
433 double b = im(ix0,
y + 1) - im(ix0,
y);
435 ymax =
y + ((b == 0.0) ? 0.5 : a/
b);
439 for (
int y = iy0;
y > 0; --
y) {
440 if (im(ix0,
y - 1) <= test) {
441 double a = test - im(ix0,
y);
442 double b = im(ix0,
y - 1) - im(ix0,
y);
444 ymin =
y - ((b == 0.0) ? 0.5 : a/
b);
459 fit->elnew = fit->param;
468 template<
typename PixelT>
469 FittedModel twodg(afw::image::Image<PixelT>
const& im,
480 fg2(im, x0, y0, &fit);
481 if (fit.status != 0) {
483 std::copy(&fit.elnew[0], &fit.elnew[0] + fit.elnew.size(), params.
begin());
485 return FittedModel(fit.status, params);
490 for (fit.iter = 0; fit.status == 0; ++fit.iter) {
496 std::copy(&fit.elnew[0], &fit.elnew[0] + fit.elnew.size(), params.
begin());
498 return FittedModel(fit.status, params, fit.iter, fit.flamd, fit.chnew);
503 #define MAKE_TWODG(IMAGE_T) \ 504 template FittedModel twodg(IMAGE_T const& im, double x0, double y0) 529 template<
typename PixelT>
548 result.
x = center.getX();
549 result.
y = center.getY();
554 int x =
static_cast<int>(center.getX() + 0.5);
555 int y =
static_cast<int>(center.getY() + 0.5);
583 #define MAKE_FIT_CENTROID(IMAGE_T) \ 584 template afw::geom::Point2D GaussianCentroidAlgorithm::fitCentroid(IMAGE_T const &, double, double) 598 if (flag == GaussianCentroidAlgorithm::FAILURE)
continue;
599 if (mapper.getInputSchema().getNames().count(
600 mapper.getInputSchema().join(
name, flag.
name)) == 0)
continue;
603 mapper.addMapping(key);
std::size_t size() const
return the current size (number of defined elements) of the collection
Defines the fields and offsets for a table.
GaussianCentroidAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
std::vector< Raster > rast
double indexToPosition(double ind)
Convert image index to image position.
table::Key< std::string > name
afw::table::Schema schema
A mapping between the keys of two Schemas, used to copy data between them.
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...
A reusable struct for centroid measurements.
static afw::geom::Point2D fitCentroid(afw::image::Image< PixelT > const &im, double x0, double y0)
Compute centroids with 2-D Gaussian fitter.
CentroidElement x
x (column) coordinate of the measured position
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=nullptr) const
Handle an exception thrown by the current algorithm by setting flags in the given record...
static FlagDefinition const FAILURE
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Utility class for handling flag fields that indicate the failure modes of an algorithm.
A base class for image defects.
MaskedImageT getMaskedImage()
Return the MaskedImage.
std::vector< double > params
Point< double, 2 > Point2D
A C++ control class to handle GaussianCentroidAlgorithm's configuration.
CentroidChecker _centroidChecker
#define LSST_EXCEPT(type,...)
Create an exception with a given type and message and optionally other arguments (dependent on the ty...
#define MAKE_TWODG(IMAGE_T)
Algorithm provides no uncertainy information at all.
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
static FlagDefinitionList const & getFlagDefinitions()
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
A FunctorKey for CentroidResult.
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
CentroidElement y
y (row) coordinate of the measured position
#define MAKE_FIT_CENTROID(IMAGE_T)
Record class that contains measurements made on a single exposure.
static FlagDefinition const NO_PEAK
vector-type utility class to build a collection of FlagDefinitions
A class to represent a 2-dimensional array of pixels.
CentroidResultKey _centroidKey
SafeCentroidExtractor _centroidExtractor