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
PsfAttributes.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
32 #include <cmath>
33 #include "lsst/base.h"
34 #include "lsst/afw/geom/Point.h"
35 #include "lsst/afw/geom/Angle.h"
39 #include "lsst/afw/detection/Psf.h"
41 #include "lsst/afw/table/Source.h"
44 
45 /************************************************************************************************************/
46 
47 namespace afwDetection = lsst::afw::detection;
48 namespace afwGeom = lsst::afw::geom;
49 namespace afwImage = lsst::afw::image;
50 namespace afwMath = lsst::afw::math;
51 
52 namespace lsst {
53 namespace meas {
54 namespace algorithms {
60  int const iX,
61  int const iY
62  )
63 {
64  // N.b. (iX, iY) are ints so that we know this image is centered in the central pixel of _psfImage
65  _psfImage = psf->computeImage(afwGeom::PointD(iX, iY));
66 }
67 
73  lsst::afw::geom::Point2I const& cen
74  ) :
75  // N.b. cen is a PointI so that we know this image is centered in the central pixel of _psfImage
76  _psfImage(psf->computeImage(afwGeom::PointD(cen)))
77 {
78 }
79 
80 namespace {
81 
82 /*
83  * Return an estimate of <r> == <sqrt(x^2 + y^2)> for an image (i.e. sum(I*r)/sum(I))
84  *
85  * For a Gaussian N(0, alpha^2), <r> = sqrt(pi/2) alpha
86  */
87 template<typename ImageT>
88 double
89 computeFirstMoment(ImageT const& image, // the data to process
90  float const xCen, float const yCen // centre of object
91  )
92 {
93  double sum = 0.0;
94  double norm = 0.0;
95  for (int iY = 0; iY != image->getHeight(); ++iY) {
96  int iX = 0;
97  for (afwImage::Image<double>::x_iterator ptr = image->row_begin(iY),
98  end = image->row_end(iY); ptr != end; ++ptr, ++iX) {
99  double const x = iX - xCen;
100  double const y = iY - yCen;
101  double const r = std::sqrt( x*x + y*y );
102  double const m = (*ptr)*r;
103  norm += *ptr;
104  sum += m;
105  }
106  }
107 
108  std::string errmsg("");
109  if (sum < 0.0) {
110  errmsg = "sum(I*r) is negative. ";
111  }
112  if (norm <= 0.0) {
113  errmsg += "sum(I) is <= 0.";
114  }
115  if (errmsg != "") {
116  throw LSST_EXCEPT(lsst::pex::exceptions::DomainError, errmsg);
117  }
118 
119  return sum/norm;
120 }
121 
122 /*
123  * Return an estimate of <r^2> == <x^2 + y^2> for an image (i.e. sum(I*r^2)/sum(I))
124  *
125  * For a Gaussian N(0, alpha^2), <r^2> = 2 alpha^2
126  */
127 template<typename ImageT>
128 double
129 computeSecondMoment(ImageT const& image, // the data to process
130  float const xCen, float const yCen // centre of object
131  )
132 {
133  double sum = 0.0;
134  double norm = 0.0;
135  for (int iY = 0; iY != image->getHeight(); ++iY) {
136  int iX = 0;
137  for (afwImage::Image<double>::x_iterator ptr = image->row_begin(iY),
138  end = image->row_end(iY); ptr != end; ++ptr, ++iX) {
139  double const x = iX - xCen;
140  double const y = iY - yCen;
141  double const r2 = x*x + y*y;
142  double const m = (*ptr)*r2;
143  norm += *ptr;
144  sum += m;
145  }
146  }
147 
148  std::string errmsg("");
149  if (sum < 0.0) {
150  errmsg = "sum(I*r*r) is negative. ";
151  }
152  if (norm <= 0.0) {
153  errmsg += "sum(I) is <= 0.";
154  }
155  if (errmsg != "") {
156  throw LSST_EXCEPT(lsst::pex::exceptions::DomainError, errmsg);
157  }
158 
159  return sum/norm;
160 }
161 
162 /*****************************************************************************/
163 /*
164  * Calculate weighted moments of an object up to 2nd order
165  */
166 template<typename ImageT>
167 std::pair<bool, double>
168 calcmom(ImageT const& image, // the image data
169  float const xCen, float const yCen, // centre of object
170  double const w11 // weights
171  )
172 {
173  assert(w11 >= 0); /* i.e. it was set */
174  if (fabs(w11) > 1e6) {
175  return std::make_pair(false, std::numeric_limits<double>::quiet_NaN());
176  }
177 
178  double sum = 0, sumrr = 0.0;
179 
180  for (int i = 0; i < image.getHeight(); ++i) {
181  float const y = i - yCen;
182  float const y2 = y*y;
183 
184  typename ImageT::x_iterator ptr = image.row_begin(i);
185  for (int j = 0; j < image.getWidth(); ++j, ++ptr) {
186  float const x = j - xCen;
187  float const x2 = x*x;
188  float const expon = (x2 + y2)*w11;
189 
190  if (expon <= 14.0) {
191  float const weight = exp(-0.5*expon);
192  float const tmod = *ptr;
193  float const ymod = tmod*weight;
194  sum += ymod;
195  sumrr += (x2 + y2)*ymod;
196  }
197  }
198  }
199 
200  if (sum <= 0 || sumrr < 0) {
201  return std::make_pair(false, std::numeric_limits<double>::quiet_NaN());
202  }
203 
204  return std::make_pair(true, 0.5*sumrr/sum); // 0.5: 1-D moment
205 }
206 
207 /*
208  * Return an estimate of <r^2> == <x^2 + y^2> for an image using adaptive moments
209  *
210  * For a Gaussian N(0, alpha^2), <r^2> = 2 alpha^2
211  *
212  * This is basically the SdssShape code simplified for a circularly symmetrical case. I don't want to call
213  * the shape code here as this class may well be moving to afw along with Psf
214  */
215 template<typename ImageT>
216 double
217 computeSecondMomentAdaptive(ImageT const& image, // the data to process
218  float const xCen, float const yCen // centre of object
219  )
220 {
221  int const MAXIT = 100; // \todo from Policy XXX
222  float const TOL = 0.0001;
223  double w11 = 0.5; // current weight for moments
224  float sigma11_ow_old = 1e6; // previous version of sigma11_ow
225 
226  bool unweighted = false; // do we need to use an unweighted moment?
227  int iter = 0; // iteration number
228  for (; iter < MAXIT; ++iter) {
229  std::pair<bool, double> moments = calcmom(*image, xCen, yCen, w11);
230 
231  if (not moments.first) {
232  unweighted = true;
233  break;
234  }
235 /*
236  * Did we converge?
237  */
238  float const sigma11_ow = moments.second; // quadratic moments of weight*object
239 
240  if (iter > 0 && fabs(sigma11_ow/sigma11_ow_old - 1.0) < TOL) {
241  break; /* yes; we converged */
242  }
243 
244  sigma11_ow_old = sigma11_ow;
245 /*
246  * Didn't converge, calculate new values for weighting function
247  *
248  * The product of two Gaussians is a Gaussian, the inverse-variances add
249  *
250  * We know sigma11_ow and sigma11W == 1/w11, the variances of the weighted object
251  * and of the weights themselves. We can estimate the object's variance as
252  * 1/sigma11_ow - 1/sigma11W
253  * and, as we want to find a set of weights with the _same_ covariance as the
254  * object we take this to be the an estimate of our correct weights.
255  *
256  * N.b. This assumes that the object is roughly Gaussian.
257  * Consider the object:
258  * O == delta(x + p) + delta(x - p)
259  * the covariance of the weighted object is equal to that of the unweighted
260  * object, and this prescription fails badly. If we detect this, we set
261  * unweighted, and calculate the UNweighted moments
262  * instead.
263  */
264  w11 = 1/sigma11_ow - w11; // inverse of new sigma11_ow
265 
266  if (w11 <= 0) { // product-of-Gaussians assumption failed
267  unweighted = true;
268  break;
269  }
270  }
271 /*
272  * Problems; try calculating the un-weighted moments
273  */
274  double sigma11W; // quadratic moment of the weighting function
275 
276  if (iter == MAXIT || unweighted) {
277  w11 = 0; // => unweighted
278  std::pair<bool, double> moments = calcmom(*image, xCen, yCen, w11);
279 
280  if (moments.first) {
281  sigma11W = moments.second; // estimate of object moment
282  } else {
283  sigma11W = 1/12.0; // a single pixel
284  }
285  } else {
286  sigma11W = 1/w11;
287  }
288 
289  return 2*sigma11W; // 2: <x^2> + <y^2>
290 }
291 
292 }
293 
299  /*
300  * Estimate the PSF's center. This *really* needs to be rewritten to avoid using MeasureSources;
301  * we shouldn't need to instantiate source objects just to measure an adaptive centroid!
302  */
304  typedef afwImage::Exposure<double> Exposure;
305  Exposure::Ptr exposure = makeExposure(mi);
306  afwDetection::Footprint::Ptr foot = boost::make_shared<afwDetection::Footprint>(exposure->getBBox(
307  afwImage::LOCAL));
308 
309  afwGeom::Point2D center(_psfImage->getX0() + _psfImage->getWidth()/2,
310  _psfImage->getY0() + _psfImage->getHeight()/2);
311  double x(center.getX());
312  double y(center.getY());
314  double const xCen = fittedCenter.getX();
315  double const yCen = fittedCenter.getY();
316  switch (how) {
317  case ADAPTIVE_MOMENT:
318  return ::sqrt(0.5*computeSecondMomentAdaptive(_psfImage, xCen, yCen));
319  case FIRST_MOMENT:
320  return ::sqrt(2.0/afwGeom::PI)*computeFirstMoment(_psfImage, xCen, yCen);
321  case SECOND_MOMENT:
322  return ::sqrt(0.5*computeSecondMoment(_psfImage, xCen, yCen));
323  case NOISE_EQUIVALENT:
324  return ::sqrt(computeEffectiveArea()/(4*afwGeom::PI));
325  case BICKERTON:
326  {
327  double sum = 0.0;
328  double norm = 0.0;
329  for (int iY = 0; iY != _psfImage->getHeight(); ++iY) {
330  int iX = 0;
331  for (afwImage::Image<double>::x_iterator ptr = _psfImage->row_begin(iY),
332  end = _psfImage->row_end(iY); ptr != end;
333  ++ptr, ++iX) {
334  double const x = iX - xCen;
335  double const y = iY - yCen;
336  double const r = std::sqrt( x*x + y*y );
337  double const m = (*ptr)*r;
338  norm += (*ptr)*(*ptr);
339  sum += m*m;
340  }
341  }
342  return sqrt(sum/norm);
343  }
344  default:
345  abort();
346  }
347 }
348 
354 
355  double sum = 0.0;
356  double sumsqr = 0.0;
357  for (int iY = 0; iY != _psfImage->getHeight(); ++iY) {
359  for (afwImage::Image<double>::x_iterator ptr = _psfImage->row_begin(iY); ptr != end; ++ptr) {
360  sum += *ptr;
361  sumsqr += (*ptr)*(*ptr);
362  }
363  }
364  return sum*sum/sumsqr;
365 }
366 
367 }}}
int y
int iter
Calculate width as sqrt(n_eff/(4 pi))
Definition: PSF.h:71
#define MAXIT
Definition: Function2D.cc:646
A coordinate class intended to represent absolute positions.
PsfAttributes(boost::shared_ptr< lsst::afw::detection::Psf const > psf, int const iX, int const iY)
Constructor for PsfAttributes.
Represent a set of pixels of an arbitrary shape and size.
x_iterator row_begin(int y) const
Definition: Image.h:319
tbl::Key< double > weight
A class to contain the data, WCS, and other information needed to describe an image of the sky...
Definition: Exposure.h:48
#define CONST_PTR(...)
Definition: base.h:47
T norm(const T &x)
Definition: Integrate.h:191
static afw::geom::Point2D fitCentroid(afw::image::Image< PixelT > const &im, double x0, double y0)
x0, y0 is an initial guess for position, column
boost::shared_ptr< Footprint > Ptr
Definition: Footprint.h:67
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Calculate width using &lt;r^2&gt;
Definition: PSF.h:70
boost::shared_ptr< lsst::afw::image::Image< double > > _psfImage
boost::shared_ptr< lsst::afw::image::Image< double > > _psfImage
Definition: PSF.h:82
Support for PCA analysis of 2-D images.
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
double computeEffectiveArea() const
Compute the effective area of the psf ( sum(I)^2/sum(I^2) )
Calculate width using &lt;r&gt;
Definition: PSF.h:69
int x
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Weight &lt;r^2&gt; by I^2 to avoid negative fluxes.
Definition: PSF.h:72
tuple m
Definition: lsstimport.py:48
Calculate width using adaptive Gaussian weights.
Definition: PSF.h:68
boost::shared_ptr< Image > computeImage(geom::Point2D position=makeNullPoint(), image::Color color=image::Color(), ImageOwnerEnum owner=COPY) const
Return an Image of the PSF, in a form that can be compared directly with star images.
Class to ensure constraints for spatial modeling.
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
Definition: Exposure.h:308
Point< double, 2 > PointD
Definition: Point.h:285
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:68
double computeGaussianWidth(Method how=ADAPTIVE_MOMENT) const
Compute the &#39;sigma&#39; value for an equivalent gaussian psf.
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43