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
Statistics.h
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 
25 #if !defined(LSST_AFW_MATH_STATISTICS_H)
26 #define LSST_AFW_MATH_STATISTICS_H
27 
40 #include <algorithm>
41 #include <cassert>
42 #include <limits>
43 #include "boost/iterator/iterator_adaptor.hpp"
44 #include "boost/tuple/tuple.hpp"
45 #include "boost/shared_ptr.hpp"
48 
49 namespace lsst {
50 namespace afw {
51  namespace image {
52  template<typename> class Image;
53  template<typename, typename, typename> class MaskedImage;
54  }
55 namespace math {
56  template<typename> class MaskedVector; // forward declaration
57 
58  typedef lsst::afw::image::VariancePixel WeightPixel; // Type used for weights
59 
63 enum Property {
64  NOTHING = 0x0,
65  ERRORS = 0x1,
66  NPOINT = 0x2,
67  MEAN = 0x4,
68  STDEV = 0x8,
69  VARIANCE = 0x10,
70  MEDIAN = 0x20,
71  IQRANGE = 0x40,
72  MEANCLIP = 0x80,
73  STDEVCLIP = 0x100,
74  VARIANCECLIP = 0x200,
75  MIN = 0x400,
77  MAX = 0x800,
78  SUM = 0x1000,
79  MEANSQUARE = 0x2000,
80  ORMASK = 0x4000
81 };
82 Property stringToStatisticsProperty(std::string const property);
83 
84 
93  typedef enum { WEIGHTS_FALSE=0, WEIGHTS_TRUE=1, WEIGHTS_NONE } WeightsBoolean; // initial state is NONE
94 public:
95 
96  typedef boost::shared_ptr<StatisticsControl> Ptr;
97  typedef boost::shared_ptr<StatisticsControl> const ConstPtr;
98 
100  double numSigmaClip = 3.0,
101  int numIter = 3,
102  lsst::afw::image::MaskPixel andMask = 0x0,
103  bool isNanSafe = true,
104  int useWeights = WEIGHTS_NONE
105  ) :
106  _numSigmaClip(numSigmaClip),
107  _numIter(numIter),
108  _andMask(andMask),
110  _isNanSafe(isNanSafe),
111  _useWeights(useWeights == 0 ? WEIGHTS_FALSE : (useWeights == 1) ? WEIGHTS_TRUE : WEIGHTS_NONE),
114  {
115  try {
117  } catch(lsst::pex::exceptions::InvalidParameterError) {
118  ; // Mask has no NO_DATA plane defined
119  }
120 
121  assert(_numSigmaClip > 0);
122  assert(_numIter > 0);
123  }
124 
126 
131  double getMaskPropagationThreshold(int bit) const;
132  void setMaskPropagationThreshold(int bit, double threshold);
134 
135  double getNumSigmaClip() const { return _numSigmaClip; }
136  int getNumIter() const { return _numIter; }
137  int getAndMask() const { return _andMask; }
138  int getNoGoodPixelsMask() const { return _noGoodPixelsMask; }
139  bool getNanSafe() const { return _isNanSafe; }
140  bool getWeighted() const { return _useWeights == WEIGHTS_TRUE ? true : false; }
141  bool getWeightedIsSet() const { return _useWeights != WEIGHTS_NONE ? true : false; }
143 
144  void setNumSigmaClip(double numSigmaClip) { assert(numSigmaClip > 0); _numSigmaClip = numSigmaClip; }
145  void setNumIter(int numIter) { assert(numIter > 0); _numIter = numIter; }
146  void setAndMask(int andMask) { _andMask = andMask; }
147  void setNoGoodPixelsMask(int noGoodPixelsMask) { _noGoodPixelsMask = noGoodPixelsMask; }
148  void setNanSafe(bool isNanSafe) { _isNanSafe = isNanSafe; }
149  void setWeighted(bool useWeights) { _useWeights = useWeights ? WEIGHTS_TRUE : WEIGHTS_FALSE; }
150  void setCalcErrorFromInputVariance(bool calcErrorFromInputVariance) {
151  _calcErrorFromInputVariance = calcErrorFromInputVariance;
152  }
153 
154 private:
155 
156  friend class Statistics;
157 
158  double _numSigmaClip; // Number of standard deviations to clip at
159  int _numIter; // Number of iterations
160  int _andMask; // and-Mask to specify which mask planes to ignore
161  int _noGoodPixelsMask; // mask to set if no values are acceptable
162  bool _isNanSafe; // Check for NaNs & Infs before running (slower)
163  WeightsBoolean _useWeights; // Calculate weighted statistics (enum because of 3-valued logic)
164  bool _calcErrorFromInputVariance; // Calculate errors from the input variances, if available
165  std::vector<double> _maskPropagationThresholds; // Thresholds for when to propagate mask bits,
166  // treated like a dict (unset bits are set to 1.0)
167 };
168 
212 class Statistics {
213 public:
215  typedef std::pair<double, double> Value;
216 
217  template<typename ImageT, typename MaskT, typename VarianceT>
218  explicit Statistics(ImageT const &img,
219  MaskT const &msk,
220  VarianceT const &var,
221  int const flags,
222  StatisticsControl const& sctrl = StatisticsControl());
223 
224  template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
225  explicit Statistics(ImageT const &img,
226  MaskT const &msk,
227  VarianceT const &var,
228  WeightT const &weights,
229  int const flags,
230  StatisticsControl const& sctrl = StatisticsControl());
231 
232  Value getResult(Property const prop = NOTHING) const;
233 
234  double getError(Property const prop = NOTHING) const;
235  double getValue(Property const prop = NOTHING) const;
237  return _allPixelOrMask;
238  }
239 
240 private:
241  long _flags; // The desired calculation
242 
243  int _n; // number of pixels in the image
244  Value _mean; // the image's mean
245  Value _variance; // the image's variance
246  double _min; // the image's minimum
247  double _max; // the image's maximum
248  double _sum; // the sum of all the image's pixels
249  Value _meanclip; // the image's N-sigma clipped mean
250  Value _varianceclip; // the image's N-sigma clipped variance
251  Value _median; // the image's median
252  double _iqrange; // the image's interquartile range
253  lsst::afw::image::MaskPixel _allPixelOrMask; // the 'or' of all masked pixels
254 
255  StatisticsControl _sctrl; // the control structure
256  bool _weightsAreMultiplicative; // Multiply by weights rather than dividing by them
257 
258  template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
259  void doStatistics(ImageT const &img,
260  MaskT const &msk,
261  VarianceT const &var,
262  WeightT const &weights,
263  int const flags,
264  StatisticsControl const& sctrl);
265 };
266 
267 /************************************* The factory functions **********************************/
272 template <typename ValueT>
274  : public boost::iterator_adaptor<infinite_iterator<ValueT>,
275  const ValueT*, const ValueT,
276  boost::forward_traversal_tag> {
277 public:
278  infinite_iterator() : infinite_iterator::iterator_adaptor_(0) {}
279  explicit infinite_iterator(const ValueT* p) : infinite_iterator::iterator_adaptor_(p) {}
280 private:
282  void increment() { ; } // never actually advance the iterator
283 };
288 template<typename ValueT>
290 public:
292  explicit MaskImposter(ValueT val = 0) { _val[0] = val; }
293  x_iterator row_begin(int) const { return x_iterator(_val); }
294 private:
295  ValueT _val[1];
296 };
297 
298 
303 template<typename Pixel>
306  int const flags,
307  StatisticsControl const& sctrl = StatisticsControl()
308  ) {
310  return Statistics(img, msk, var, flags, sctrl);
311 }
312 
313 
318 template<typename ImageT, typename MaskT, typename VarianceT>
319 Statistics makeStatistics(ImageT const &img,
320  MaskT const &msk,
321  VarianceT const &var,
322  int const flags,
323  StatisticsControl const& sctrl = StatisticsControl()
324  ) {
325  return Statistics(img, msk, var, flags, sctrl);
326 }
327 
332 template<typename Pixel>
335  int const flags,
336  StatisticsControl const& sctrl = StatisticsControl()
337  )
338 {
339  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
340  return Statistics(*mimg.getImage(), *mimg.getMask(), *mimg.getVariance(), flags, sctrl);
341  } else {
343  return Statistics(*mimg.getImage(), *mimg.getMask(), var, flags, sctrl);
344  }
345 }
346 
351 template<typename Pixel>
355  int const flags,
356  StatisticsControl const& sctrl = StatisticsControl()
357  )
358 {
359  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance() ||
360  (!sctrl.getWeightedIsSet() && (weights.getWidth() != 0 && weights.getHeight() != 0))) {
361  return Statistics(*mimg.getImage(), *mimg.getMask(), *mimg.getVariance(), weights, flags, sctrl);
362  } else {
364  return Statistics(*mimg.getImage(), *mimg.getMask(), var, weights, flags, sctrl);
365  }
366 }
367 
374  int const flags,
375  StatisticsControl const& sctrl = StatisticsControl());
376 
377 
378 
383 template<typename Pixel>
385  lsst::afw::image::Image<Pixel> const &img,
386  int const flags,
387  StatisticsControl const& sctrl = StatisticsControl()
388 ) {
389  // make a phony mask that will be compiled out
391  MaskImposter<WeightPixel> const var;
392  return Statistics(img, msk, var, flags, sctrl);
393 }
394 
395 
400 template<typename ValueT>
402 public:
403 
404  // types we'll use in Statistics
405  typedef typename std::vector<ValueT>::const_iterator x_iterator;
406  typedef typename std::vector<ValueT>::const_iterator fast_iterator;
407  typedef ValueT Pixel;
408 
409  // constructors for std::vector<>, and copy constructor
410  // These are both shallow! ... no actual copying of values
411  explicit ImageImposter(std::vector<ValueT> const &v) : _v(v) { }
412  explicit ImageImposter(ImageImposter<ValueT> const &img) : _v(img._getVector()) {}
413 
414  // The methods we'll use in Statistics
415  x_iterator row_begin(int) const { return _v.begin(); }
416  x_iterator row_end(int) const { return _v.end(); }
417  int getWidth() const { return _v.size(); }
418  int getHeight() const { return 1; }
419 
420  bool empty() const { return _v.empty(); }
421 private:
422  std::vector<ValueT> const &_v; // a private reference to the data
423  std::vector<ValueT> const &_getVector() const { return _v; } // get the ref for the copyCon
424 };
425 
430 template<typename EntryT>
431 Statistics makeStatistics(std::vector<EntryT> const &v,
432  int const flags,
433  StatisticsControl const& sctrl = StatisticsControl()
434  ) {
435  ImageImposter<EntryT> img(v); // wrap the vector in a fake image
436  MaskImposter<lsst::afw::image::MaskPixel> msk; // instantiate a fake mask that will be compiled out.
438  return Statistics(img, msk, var, flags, sctrl);
439 }
440 
445 template<typename EntryT>
446 Statistics makeStatistics(std::vector<EntryT> const &v,
447  std::vector<WeightPixel> const &vweights,
448  int const flags,
449  StatisticsControl const& sctrl = StatisticsControl()
450  ) {
451  ImageImposter<EntryT> img(v); // wrap the vector in a fake image
452  MaskImposter<lsst::afw::image::MaskPixel> msk; // instantiate a fake mask that will be compiled out.
454 
455  ImageImposter<WeightPixel> weights(vweights);
456 
457  return Statistics(img, msk, var, weights, flags, sctrl);
458 }
459 
464 template<typename EntryT>
466  int const flags,
467  StatisticsControl const& sctrl = StatisticsControl()
468  ) {
469  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
470  return Statistics(*mv.getImage(), *mv.getMask(), *mv.getVariance(), flags, sctrl);
471  } else {
473  return Statistics(*mv.getImage(), *mv.getMask(), var, flags, sctrl);
474  }
475 }
476 
481 template<typename EntryT>
483  std::vector<WeightPixel> const &vweights,
484  int const flags,
485  StatisticsControl const& sctrl = StatisticsControl()
486  ) {
487  ImageImposter<WeightPixel> weights(vweights);
488 
489  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
490  return Statistics(*mv.getImage(), *mv.getMask(), *mv.getVariance(), weights, flags, sctrl);
491  } else {
493  return Statistics(*mv.getImage(), *mv.getMask(), var, weights, flags, sctrl);
494  }
495 }
496 
497 }}}
498 
499 #endif
ImageImposter(ImageImposter< ValueT > const &img)
Definition: Statistics.h:412
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1009
boost::uint16_t MaskPixel
double getMaskPropagationThreshold(int bit) const
Definition: Statistics.cc:691
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
Statistics makeStatistics(std::vector< EntryT > const &v, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle std::vector&lt;&gt;
Definition: Statistics.h:431
boost::shared_ptr< StatisticsControl > Ptr
Definition: Statistics.h:96
lsst::afw::image::MaskedImage< EntryT >::VariancePtr getVariance() const
Definition: MaskedVector.h:89
lsst::afw::image::MaskedImage< EntryT >::MaskPtr getMask() const
Definition: MaskedVector.h:86
estimate sample mean
Definition: Statistics.h:67
std::pair< double, double > Value
The type used to report (value, error) for desired statistics.
Definition: Statistics.h:215
std::vector< ValueT > const & _getVector() const
Definition: Statistics.h:423
estimate sample N-sigma clipped mean (N set in StatisticsControl, default=3)
Definition: Statistics.h:72
void setNumIter(int numIter)
Definition: Statistics.h:145
StatisticsControl(double numSigmaClip=3.0, int numIter=3, lsst::afw::image::MaskPixel andMask=0x0, bool isNanSafe=true, int useWeights=WEIGHTS_NONE)
Definition: Statistics.h:99
infinite_iterator< ValueT > x_iterator
Definition: Statistics.h:291
void setAndMask(int andMask)
Definition: Statistics.h:146
x_iterator row_begin(int) const
Definition: Statistics.h:293
MaskImposter(ValueT val=0)
Definition: Statistics.h:292
Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Constructor for Statistics object.
Definition: Statistics.cc:746
estimate sample standard deviation
Definition: Statistics.h:68
estimate sample variance
Definition: Statistics.h:69
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle regular (non-masked) Images.
Definition: Statistics.h:384
find sum of pixels in the image
Definition: Statistics.h:78
Value getResult(Property const prop=NOTHING) const
Return the value and error in the specified statistic (e.g. MEAN)
Definition: Statistics.cc:887
ImageImposter(std::vector< ValueT > const &v)
Definition: Statistics.h:411
lsst::afw::image::VariancePixel WeightPixel
Definition: Statistics.h:56
boost::shared_ptr< StatisticsControl > const ConstPtr
Definition: Statistics.h:97
Statistics makeStatistics(lsst::afw::math::MaskedVector< EntryT > const &mv, std::vector< WeightPixel > const &vweights, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle lsst::afw::math::MaskedVector&lt;&gt;
Definition: Statistics.h:482
estimate sample N-sigma clipped stdev (N set in StatisticsControl, default=3)
Definition: Statistics.h:73
Statistics makeStatistics(lsst::afw::math::MaskedVector< EntryT > const &mv, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle lsst::afw::math::MaskedVector&lt;&gt;
Definition: Statistics.h:465
int const x0
Definition: saturated.cc:45
find mean value of square of pixel values
Definition: Statistics.h:79
lsst::afw::image::MaskedImage< EntryT >::ImagePtr getImage() const
Definition: MaskedVector.h:83
infinite_iterator(const ValueT *p)
Definition: Statistics.h:279
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Include errors of requested quantities.
Definition: Statistics.h:65
We don&#39;t want anything.
Definition: Statistics.h:64
VariancePtr getVariance(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:890
lsst::afw::image::MaskPixel getOrMask() const
Definition: Statistics.h:236
Statistics makeStatistics(lsst::afw::image::MaskedImage< Pixel > const &mimg, lsst::afw::image::Image< WeightPixel > const &weights, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle MaskedImages, just pass the getImage() and getMask() values right on through.
Definition: Statistics.h:352
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
Statistics makeStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a straight front-end to the constructor.
Definition: Statistics.h:319
Pass parameters to a Statistics objectA class to pass parameters which control how the stats are calc...
Definition: Statistics.h:92
void doStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights, int const flags, StatisticsControl const &sctrl)
Definition: Statistics.cc:793
static MaskPixelT getPlaneBitMask(const std::vector< std::string > &names)
Return the bitmask corresponding to a vector of plane names OR&#39;d together.
Definition: Mask.cc:860
void setNumSigmaClip(double numSigmaClip)
Definition: Statistics.h:144
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:93
void setNanSafe(bool isNanSafe)
Definition: Statistics.h:148
std::vector< ValueT >::const_iterator fast_iterator
Definition: Statistics.h:406
A Mask wrapper to provide an infinite_iterator for Mask::row_begin(). This allows a fake Mask to be p...
Definition: Statistics.h:289
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
Statistics makeStatistics(std::vector< EntryT > const &v, std::vector< WeightPixel > const &vweights, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle std::vector&lt;&gt;
Definition: Statistics.h:446
StatisticsControl _sctrl
Definition: Statistics.h:255
void setCalcErrorFromInputVariance(bool calcErrorFromInputVariance)
Definition: Statistics.h:150
float VariancePixel
! default type for Masks and MaskedImage Masks
estimate sample maximum
Definition: Statistics.h:77
number of sample points
Definition: Statistics.h:66
double getError(Property const prop=NOTHING) const
Return the error in the desired property (if specified in the constructor)
Definition: Statistics.cc:1020
std::vector< ValueT > const & _v
Definition: Statistics.h:422
estimate sample minimum
Definition: Statistics.h:76
This iterator will never increment. It is returned by row_begin() in the MaskImposter class (below) t...
Definition: Statistics.h:273
estimate sample median
Definition: Statistics.h:70
std::vector< ValueT >::const_iterator x_iterator
Definition: Statistics.h:405
bool getCalcErrorFromInputVariance() const
Definition: Statistics.h:142
MaskPtr getMask(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:879
void setMaskPropagationThreshold(int bit, double threshold)
Definition: Statistics.cc:699
void setWeighted(bool useWeights)
Definition: Statistics.h:149
void setNoGoodPixelsMask(int noGoodPixelsMask)
Definition: Statistics.h:147
estimate sample inter-quartile range
Definition: Statistics.h:71
lsst::afw::image::MaskPixel _allPixelOrMask
Definition: Statistics.h:253
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
A vector wrapper to provide a vector with the necessary methods and typedefs to be processed by Stati...
Definition: Statistics.h:401
get the or-mask of all pixels used.
Definition: Statistics.h:80
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:715
Implementation of the Class MaskedImage.
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
x_iterator row_end(int) const
Definition: Statistics.h:416
Property
control what is calculated
Definition: Statistics.h:63
Statistics makeStatistics(lsst::afw::image::MaskedImage< Pixel > const &mimg, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle MaskedImages, just pass the getImage() and getMask() values right on through.
Definition: Statistics.h:333
ImageT val
Definition: CR.cc:154
friend class boost::iterator_core_access
Definition: Statistics.h:281
std::vector< double > _maskPropagationThresholds
Definition: Statistics.h:165
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:304
x_iterator row_begin(int) const
Definition: Statistics.h:415