LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
GaussianCentroid.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2013 LSST Corporation.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 #include "ndarray/eigen.h"
24 #include "lsst/afw/table/Source.h"
26 
27 #include <algorithm>
28 #include "Eigen/Core"
29 #include "Eigen/LU"
30 
31 
32 namespace lsst { namespace meas { namespace base {
33 namespace {
34 #define USE_WEIGHT 0 // zweight is only set, not used. It isn't even set if this is false
35 struct Raster {
36  int x, y; // x, y index of rastered pixel
37  double zo; // pixel's value
38 #if USE_WEIGHT
39  double zweight; // weight for the pixel
40 #endif
41 };
42 
43 struct Fit2d {
44  typedef Eigen::Matrix<double, FittedModel::NPARAM, FittedModel::NPARAM> Matrix;
45  typedef Eigen::Matrix<double, FittedModel::NPARAM, 1> Vector;
46 
47  template<typename PixelT>
48  explicit Fit2d(afw::image::Image<PixelT> const& im) : wide(32), rast(*new std::vector<Raster>(wide*wide)) {
49  ncols = im.getWidth();
50  nrows = im.getHeight();
51 
52  dc2zmin = 0.3;
53  dc2zmax = 3.5;
54  tooclose = 3;
55  toosmall = 1.5;
56  nitmax = 15;
57  flamdok = 1.0e-7;
58  ratiomin = -1.0e-7;
59  lost = 3.5;
60  xlamd = 5.0;
61 
62  iter = 0;
63  flamd = 1.0;
64  }
65  ~Fit2d() {
66  delete &rast;
67  }
68 
69  int wide; // patch of image being fitted is wide*wide
70  std::vector<Raster> &rast; // The pixels being fitted, thought of as a raster scan
71  int nin;
72  int nref;
73  int iter;
74  int tooclose;
75  Vector param;
76  Matrix alpha;
77  Vector beta;
78  Vector elold;
79  Vector elnew;
80  Vector sgold;
81  Vector sgnew;
82  Matrix alold;
83  Matrix alnew;
84  Vector beold;
85  Vector benew;
86  double chisq;
87  double chold;
88  double stold;
89  double chnew;
90  double stnew;
91  double flamd;
92  double flamdok;
93  double xlamd;
94  double ratiomin;
95  int nitmax;
96  double dc2zmin;
97  double dc2zmax;
98  double lost;
99  double toosmall; // minimum acceptable width
100  int ncols; // number of columns in image
101  int nrows; // number of rows in image
102  int status; // status of processing
103 };
104 
105 /************************************************************************************************************/
110 static void dcalc2(Fit2d *fit, // the struct carrying all the fit information
111  Fit2d::Vector const &el, // current parameters of the fit
112  Fit2d::Matrix *alpha, // desired alpha matrix
113  Fit2d::Vector *beta // desired beta vector
114  ) {
115 /*
116  * Pass arguments and initialize
117  */
118  double const a = el[FittedModel::PEAK];
119  double const b = el[FittedModel::SKY];
120  double const x0 = el[FittedModel::X0];
121  double const y0 = el[FittedModel::Y0];
122  double const s = el[FittedModel::SIGMA];
123 
124  if (a <= 0.0) {
125  fit->status = FittedModel::BAD_A;
126  return;
127  } else if (s <= 0.0) {
128  fit->status = FittedModel::BAD_WIDTH;
129  return;
130  } else {
131  fit->status = 0;
132  }
133 
134  beta->setZero();
135  alpha->setZero();
136 
137  fit->nin = 0;
138  fit->chnew = 0.0;
139 /*
140  * Examine all pixels in raster
141  */
142  for (int i = 0, nrast = fit->rast.size(); i != nrast; i++) {
143  double const dx = (fit->rast[i].x - x0)/s;
144  double const dy = (fit->rast[i].y - y0)/s;
145  double const d = hypot(dx, dy);
146 
147  if (d >= fit->dc2zmin && d <= fit->dc2zmax) {
148  double const arg = exp(-d*d/2.0);
149  double const funct = a*arg + b;
150 
151  Fit2d::Vector df;
152  df << arg, 1.0, a*arg*dx/s, a*arg*dy/s, a*arg*(dx*dx/s + dy*dy/s);
153 
154  double const r = fit->rast[i].zo - funct; // residual
155  fit->chnew += r*r;
156  *beta += r*df;
157  *alpha += df*df.transpose(); // outer product
158 
159  fit->nin++;
160  }
161  }
162  int nu = fit->nin - (FittedModel::NPARAM + 1); // number of degrees of freedom
163  if (nu <= 0) {
164  fit->status = FittedModel::TOO_FEW;
165  return;
166  }
167  fit->stnew = sqrt(fit->chnew/nu);
168 }
169 
170 /************************************************************************************************************/
171 /*
172  * From Bevington's CURFIT
173  */
174 static void curf2(Fit2d *fit) {
175 /*
176  * Initialization
177  */
178  fit->status = 0;
179  if (fit->iter == 0) {
180  dcalc2(fit, fit->elnew, &fit->alnew, &fit->benew);
181  if (fit->status != 0) {
182  return;
183  }
184  fit->nref = fit->nin;
185  fit->stnew = sqrt(fit->chnew/(fit->nref - (FittedModel::NPARAM + 1)));
186  }
187 /*
188  * Save Current Parameters
189  */
190  fit->chold = fit->chnew;
191  fit->stold = fit->stnew;
192  fit->elold = fit->elnew;
193  fit->sgold = fit->sgnew;
194  fit->beold = fit->benew;
195  fit->alold = fit->alnew;
196 /*
197  * Previous Call To DCALC Used To Fill Current
198  */
199  int const NPLIM = 4;
200  for (int poor = 0; poor != NPLIM; ++poor) {
201  for (int j = 0; j != FittedModel::NPARAM; ++j) {
202  if (fit->alnew(j, j) == 0.0) {
203  fit->status = FittedModel::DIAGONAL;
204  return;
205  }
206  for (int k = 0; k != FittedModel::NPARAM; ++k) {
207  fit->alpha(j, k) = fit->alnew(j, k)/sqrt(fit->alnew(j, j)*fit->alnew(k, k));
208  }
209  fit->alpha(j, j) = 1.0 + fit->flamd;
210  }
211 /*
212  * Inversion.
213  */
214 #if 0
215  fit->alpha.computeInverse(); // has no way to check if inverse succeeded
216 #else
217  Eigen::FullPivLU<Fit2d::Matrix> alphaLU(fit->alpha);
218  if (!alphaLU.isInvertible()) {
219  fit->status = -1;
220  return;
221  }
222  fit->alpha = alphaLU.inverse();
223 #endif
224 
225 /*
226  * New Elements
227  */
228  for (int j = 0; j != FittedModel::NPARAM; ++j) {
229  for (int k = 0; k != FittedModel::NPARAM; ++k) {
230  fit->elnew[j] += fit->benew[k]*fit->alpha(j,k)/sqrt(fit->alnew(j, j)*fit->alnew(k, k));
231  }
232  }
233 /*
234  * Compute Current Chi-Squared And Next ALPHA+BETA
235  */
236  dcalc2(fit, fit->elnew, &fit->alnew, &fit->benew);
237  if (fit->status == FittedModel::TOO_FEW) {
238  return;
239  }
240  fit->stnew = sqrt(fit->chnew/(fit->nref - (FittedModel::NPARAM + 1)));
241 
242  for (int j = 0; j != FittedModel::NPARAM; ++j) {
243  fit->sgnew[j] = fit->stnew*sqrt(fabs(fit->alpha(j, j)/fit->alnew(j, j)));
244  }
245 /*
246  * Quick Return If Solution Got Better
247  */
248  if (fit->status == 0.0 && fit->stnew <= fit->stold) {
249  fit->flamd = fit->flamd/fit->xlamd;
250  return;
251  }
252 /*
253  * Undo Solution And Try Again
254  */
255  fit->flamd = 3.0*fit->flamd;
256  fit->chnew = fit->chold;
257  fit->stnew = fit->stold;
258  fit->elnew = fit->elold;
259  fit->sgnew = fit->sgold;
260  fit->benew = fit->beold;
261  fit->alnew = fit->alold;
262  }
263 }
264 
265 /************************************************************************************************************/
266 /*
267  * Test of convergence or fatal errors
268  */
269 static void cond2(Fit2d *fit) {
270 /*
271  * Ignore CURF errors for these conditions
272  */
273  if (fit->iter <= 3) {
274  fit->status = 0;
275  return;
276  }
277  if (fit->flamd < fit->flamdok) {
278  fit->status = FittedModel::ALMOST;
279  return;
280  }
281 /*
282  * Catch Fatal Errors
283  */
284  if (fit->status < 0) {
285  return;
286  }
287  if (fit->nin < FittedModel::NPARAM + 1) {
288  fit->status = FittedModel::TOO_FEW;
289  return;
290  }
291  if (fit->chnew <= 0.0) {
292  fit->status = FittedModel::CHI_SQUARED;
293  return;
294  }
295  if (fit->elnew[FittedModel::X0] < 0.0 || fit->elnew[FittedModel::X0] > fit->ncols ||
296  fit->elnew[FittedModel::Y0] < 0.0 || fit->elnew[FittedModel::Y0] > fit->nrows) {
297  fit->status = FittedModel::RANGE;
298  return;
299  }
300  if (fit->elnew[FittedModel::SIGMA] < 0.0) {
301  fit->status = FittedModel::BAD_WIDTH;
302  return;
303  }
304 
305  double const ex = fabs(fit->param[FittedModel::X0] - fit->elnew[FittedModel::X0]);
306  double const ey = fabs(fit->param[FittedModel::Y0] - fit->elnew[FittedModel::Y0]);
307  if (ex > fit->lost || ey > fit->lost) {
308  fit->status = FittedModel::LOST;
309  return;
310  }
311 /*
312  * Test for convergence
313  */
314  double const ratio = (fit->stnew - fit->stold)/fit->stnew;
315  if (ratio > fit->ratiomin) {
316  fit->status = (fit->flamd < 0.001) ? FittedModel::CONVERGE : FittedModel::POOR;
317  } else if (fit->iter > fit->nitmax) {
318  fit->status = FittedModel::ITERATE;
319  } else {
320  fit->status = 0;
321  }
322 }
323 
324 /************************************************************************************************************/
325 /*
326  * First guess for 2-D Gaussian
327  */
328 template<typename PixelT>
329 static void fg2(afw::image::Image<PixelT> const& im,
330  double x0, double y0,
331  Fit2d *fit
332  ) {
333 /*
334  * Sanity Check
335  */
336  int ix0 = static_cast<int>(x0 + 0.5);
337  int iy0 = static_cast<int>(y0 + 0.5);
338  if(ix0 < fit->tooclose || im.getWidth() - ix0 < fit->tooclose ||
339  iy0 < fit->tooclose || im.getHeight() - iy0 < fit->tooclose) {
340  fit->status = FittedModel::BAD_GUESS;
341  return;
342  }
343  int jx0 = ix0;
344  int jy0 = iy0;
345  double peak = im(jx0, jy0);
346 /*
347  * Extract interesting portion
348  */
349  double const w = fit->wide/2;
350  int xl = static_cast<int>(ix0 - w + 0.5);
351  if (xl < 0) {
352  xl = 0;
353  }
354  int xh = static_cast<int>(xl + 2*w + 0.5);
355  if (xh > im.getWidth()) {
356  xh = im.getWidth();
357  }
358  int yl = static_cast<int>(iy0 - w + 0.5);
359  if (yl < 0) {
360  yl = 0;
361  }
362  int yh = static_cast<int>(yl + 2*w + 0.5);
363  if (yh > im.getHeight()) {
364  yh = im.getHeight();
365  }
366 
367  double sky = im(xl, yl);
368  int put = 0;
369  for (int y = yl; y != yh; ++y) {
370  for (int x = xl; x != xh; ++x) {
371  fit->rast[put].x = x;
372  fit->rast[put].y = y;
373  fit->rast[put].zo = im(x, y);
374 #if USE_WEIGHT
375  fit->rast[put].zweight = 1.0;
376 #endif
377  if (im(x, y) < sky) {
378  sky = im(x, y);
379  }
380  double const ex = x - ix0;
381  double const ey = y - iy0;
382  double const er = hypot(ex, ey);
383  if (er <= fit->lost) {
384  if (im(x, y) > peak) {
385  jx0 = x;
386  jy0 = y;
387  peak = im(x, y);
388  }
389  }
390  put++;
391  }
392  }
393  fit->rast.resize(put);
394  ix0 = jx0;
395  iy0 = jy0;
396 /*
397  * Look For FWHM
398  */
399  double xmin = xl;
400  double xmax = xh - 1;
401  double ymin = yl;
402  double ymax = yh - 1;
403  double const test = 0.5*(im(ix0, iy0) + sky); // pixel value of half-maximum above sky
404  for (int x = ix0; x < xh - 1; ++x) {
405  if (im(x + 1, iy0) <= test) {
406  double a = test - im(x, iy0);
407  double b = im(x + 1, iy0) - im(x, iy0);
408 
409  xmax = x + ((b == 0.0) ? 0.5 : a/b);
410  break;
411  }
412  }
413  for (int x = ix0; x > 0; --x) {
414  if (im(x - 1, iy0) <= test) {
415  double a = test - im(x, iy0);
416  double b = im(x - 1, iy0) - im(x, iy0);
417 
418  xmin = x - ((b == 0.0) ? 0.5 : a/b);
419  break;
420  }
421  }
422  for (int y = iy0; y < yh - 1; ++y) {
423  if (im(ix0, y + 1) <= test) {
424  double a = test - im(ix0, y);
425  double b = im(ix0, y + 1) - im(ix0, y);
426 
427  ymax = y + ((b == 0.0) ? 0.5 : a/b);
428  break;
429  }
430  }
431  for (int y = iy0; y > 0; --y) {
432  if (im(ix0, y - 1) <= test) {
433  double a = test - im(ix0, y);
434  double b = im(ix0, y - 1) - im(ix0, y);
435 
436  ymin = y - ((b == 0.0) ? 0.5 : a/b);
437  break;
438  }
439  }
440 /*
441  * Assemble the fitting vector
442  */
443  fit->param[FittedModel::PEAK] = im(ix0, iy0) - sky;
444  fit->param[FittedModel::SKY] = sky;
445  fit->param[FittedModel::X0] = x0;
446  fit->param[FittedModel::Y0] = y0;
447  fit->param[FittedModel::SIGMA] = 0.5*((xmax - xmin)+(ymax - ymin))/2.354;
448  if (fit->param[FittedModel::SIGMA] < fit->toosmall) {
449  fit->param[FittedModel::SIGMA] = fit->toosmall;
450  }
451  fit->elnew = fit->param;
452 
453  fit->status = 0;
454 }
455 
456 /************************************************************************************************************/
460 template<typename PixelT>
461 FittedModel twodg(afw::image::Image<PixelT> const& im,
462  double x0,
463  double y0
464  ) {
465 /*
466  * Initialize the fitting structure
467  */
468  Fit2d fit(im);
469 /*
470  * Initialization
471  */
472  fg2(im, x0, y0, &fit);
473  if (fit.status != 0) {
474  std::vector<double> params(FittedModel::NPARAM);
475  std::copy(&fit.elnew[0], &fit.elnew[0] + fit.elnew.size(), params.begin());
476 
477  return FittedModel(fit.status, params);
478  }
479 /*
480  * Find best-fit model
481  */
482  for (fit.iter = 0; fit.status == 0; ++fit.iter) {
483  curf2(&fit);
484  cond2(&fit);
485  }
486 
487  std::vector<double> params(FittedModel::NPARAM);
488  std::copy(&fit.elnew[0], &fit.elnew[0] + fit.elnew.size(), params.begin());
489 
490  return FittedModel(fit.status, params, fit.iter, fit.flamd, fit.chnew);
491 }
492 //
493 // Explicit instantiations
494 //
495 #define MAKE_TWODG(IMAGE_T) \
496  template FittedModel twodg(IMAGE_T const& im, double x0, double y0)
497 
498 MAKE_TWODG(afw::image::Image<float>);
499 MAKE_TWODG(afw::image::Image<double>);
500 MAKE_TWODG(afw::image::Image<int>);
501 
502 boost::array<FlagDefinition,GaussianCentroidAlgorithm::N_FLAGS> const & getFlagDefinitions() {
503  static boost::array<FlagDefinition,GaussianCentroidAlgorithm::N_FLAGS> const flagDefs = {{
504  {"flag", "general failure flag, set if anything went wrong"},
505  {"flag_noPeak", "Fitted Centroid has a negative peak"}
506  }};
507  return flagDefs;
508 }
509 
510 } // end anonymous namespace
511 
513  Control const & ctrl,
514  std::string const & name,
516 ) : _ctrl(ctrl),
517  _centroidKey(
518  CentroidResultKey::addFields(schema, name, "centroid from Gaussian Centroid algorithm", NO_UNCERTAINTY)
519  ),
520  _centroidExtractor(schema, name, true)
521 {
522  _flagHandler = FlagHandler::addFields(schema, name,
523  getFlagDefinitions().begin(), getFlagDefinitions().end());
524 }
525 
526 
527 template<typename PixelT>
529  double x,
530  double y
531 ) {
532  base::FittedModel fit = base::twodg(image, x, y);
533  double const xCen = fit.params[base::FittedModel::X0];
534  double const yCen = fit.params[base::FittedModel::Y0];
535  return afw::geom::Point2D(xCen, yCen);
536 }
537 
538 
540  afw::table::SourceRecord & measRecord,
541  afw::image::Exposure<float> const & exposure
542 ) const {
543  // get our current best guess about the centroid: either a centroider measurement or peak.
544  afw::geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
545  CentroidResult result;
546  result.x = center.getX();
547  result.y = center.getY();
548  measRecord.set(_centroidKey, result); // better than NaN
549  typedef afw::image::Image<float> ImageT;
550  ImageT const& image = *exposure.getMaskedImage().getImage();
551 
552  int x = static_cast<int>(center.getX() + 0.5);
553  int y = static_cast<int>(center.getY() + 0.5);
554 
555  x -= image.getX0(); // work in image Pixel coordinates
556  y -= image.getY0();
557 
558  FittedModel fit = twodg(image, x, y); // here's the fitter
559  if (fit.params[FittedModel::PEAK] <= 0) {
560  throw LSST_EXCEPT(
563  NO_PEAK
564  );
565  }
566  result.x = lsst::afw::image::indexToPosition(image.getX0()) + fit.params[FittedModel::X0];
567  result.y = lsst::afw::image::indexToPosition(image.getY0()) + fit.params[FittedModel::Y0];
568  _flagHandler.setValue(measRecord, FAILURE, false); // if we had a suspect flag, we'd set that instead
569  measRecord.set(_centroidKey, result);
570 }
571 
572 
574  _flagHandler.handleFailure(measRecord, error);
575 }
576 
577 //
578 // Explicit instantiations
579 //
580 #define MAKE_FIT_CENTROID(IMAGE_T) \
581  template afw::geom::Point2D GaussianCentroidAlgorithm::fitCentroid(IMAGE_T const &, double, double)
582 
583 MAKE_FIT_CENTROID(afw::image::Image<float>);
584 MAKE_FIT_CENTROID(afw::image::Image<double>);
585 
587  Control const & ctrl,
588  std::string const & name,
589  afw::table::SchemaMapper & mapper
590 ) :
591  CentroidTransform{name, mapper}
592 {
593  for (auto flag = getFlagDefinitions().begin() + 1; flag < getFlagDefinitions().end(); ++flag) {
594  mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>(
595  mapper.getInputSchema().join(name, flag->name)).key);
596  }
597 }
598 
599 }}} // namespace lsst::meas::base
600 
601 
int y
int iter
Defines the fields and offsets for a table.
Definition: Schema.h:46
int tooclose
GaussianCentroidAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:150
int nref
int ncols
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
double xlamd
Matrix alold
double indexToPosition(double ind)
Convert image index to image position.
Definition: ImageUtils.h:54
table::Key< std::string > name
Definition: ApCorrMap.cc:71
Eigen matrix objects that present a view into an ndarray::Array.
Vector beold
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=NULL) const
std::vector< Raster > & rast
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:19
A reusable struct for centroid measurements.
Vector elold
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
Vector sgold
static afw::geom::Point2D fitCentroid(afw::image::Image< PixelT > const &im, double x0, double y0)
x0, y0 is an initial guess for position, column
CentroidElement x
x (column) coordinate of the measured position
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
Vector benew
int nrows
int const x0
Definition: saturated.cc:45
Matrix alpha
Vector param
Key< T > addMapping(Key< T > const &inputKey, bool doReplace=false)
Add a new field to the output Schema that is a copy of a field in the input Schema.
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Vector sgnew
double stnew
def error
Definition: log.py:108
Matrix alnew
FlagDefinition getDefinition(int i) const
Definition: FlagHandler.h:66
double toosmall
double w
Definition: CoaddPsf.cc:57
std::vector< double > params
double dc2zmin
double flamd
A C++ control class to handle GaussianCentroidAlgorithm&#39;s configuration.
int wide
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
void setValue(afw::table::BaseRecord &record, int i, bool value) const
Definition: FlagHandler.h:72
double dc2zmax
tbl::Schema schema
double chold
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
double stold
int x
int nin
double chisq
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
GaussianCentroidTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
double flamdok
int nitmax
double zo
#define MAKE_TWODG(IMAGE_T)
int status
double chnew
Algorithm provides no uncertainy information at all.
Definition: constants.h:42
Vector elnew
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:136
A FunctorKey for CentroidResult.
afw::table::Key< double > b
CentroidElement y
y (row) coordinate of the measured position
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=NULL) const
Definition: FlagHandler.cc:59
Point< double, 2 > Point2D
Definition: Point.h:286
#define MAKE_FIT_CENTROID(IMAGE_T)
Record class that contains measurements made on a single exposure.
Definition: Source.h:81
double ratiomin
int const y0
Definition: saturated.cc:45
double lost
static FlagHandler addFields(afw::table::Schema &schema, std::string const &prefix, FlagDefinition const *begin, FlagDefinition const *end)
Definition: FlagHandler.cc:28
Vector beta