LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
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 
37 #include <algorithm>
38 #include <cassert>
39 #include <limits>
40 #include "boost/iterator/iterator_adaptor.hpp"
41 #include "boost/tuple/tuple.hpp"
42 #include <memory>
45 
46 namespace lsst {
47 namespace afw {
48 namespace image {
49 template <typename>
50 class Image;
51 template <typename, typename, typename>
52 class MaskedImage;
53 } // namespace image
54 namespace math {
55 template <typename>
56 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  NCLIPPED = 0x8000,
82  NMASKED = 0x10000
83 };
86 
94 public:
95  enum WeightsBoolean { WEIGHTS_FALSE = 0, WEIGHTS_TRUE = 1, WEIGHTS_NONE }; // initial state is NONE
96 
97  StatisticsControl(double numSigmaClip = 3.0,
98  int numIter = 3,
100  0x0,
101  bool isNanSafe = true,
102  WeightsBoolean useWeights =
103  WEIGHTS_NONE
104  )
105  : _numSigmaClip(numSigmaClip),
106  _numIter(numIter),
107  _andMask(andMask),
108  _noGoodPixelsMask(0x0),
109  _isNanSafe(isNanSafe),
110  _useWeights(useWeights),
111  _calcErrorFromInputVariance(false),
112  _maskPropagationThresholds() {
113  try {
114  _noGoodPixelsMask = lsst::afw::image::Mask<>::getPlaneBitMask("NO_DATA");
116  ; // Mask has no NO_DATA plane defined
117  }
118 
119  assert(_numSigmaClip > 0);
120  assert(_numIter > 0);
121  }
122 
124 
129  double getMaskPropagationThreshold(int bit) const;
130  void setMaskPropagationThreshold(int bit, double threshold);
132 
133  double getNumSigmaClip() const noexcept { return _numSigmaClip; }
134  int getNumIter() const noexcept { return _numIter; }
135  int getAndMask() const noexcept { return _andMask; }
136  int getNoGoodPixelsMask() const noexcept { return _noGoodPixelsMask; }
137  bool getNanSafe() const noexcept { return _isNanSafe; }
138  bool getWeighted() const noexcept { return _useWeights == WEIGHTS_TRUE ? true : false; }
139  bool getWeightedIsSet() const noexcept { return _useWeights != WEIGHTS_NONE ? true : false; }
140  bool getCalcErrorFromInputVariance() const noexcept { return _calcErrorFromInputVariance; }
141 
142  void setNumSigmaClip(double numSigmaClip) {
143  assert(numSigmaClip > 0);
144  _numSigmaClip = numSigmaClip;
145  }
146  void setNumIter(int numIter) {
147  assert(numIter > 0);
148  _numIter = numIter;
149  }
150  void setAndMask(int andMask) { _andMask = andMask; }
151  void setNoGoodPixelsMask(int noGoodPixelsMask) { _noGoodPixelsMask = noGoodPixelsMask; }
152  void setNanSafe(bool isNanSafe) noexcept { _isNanSafe = isNanSafe; }
153  void setWeighted(bool useWeights) noexcept { _useWeights = useWeights ? WEIGHTS_TRUE : WEIGHTS_FALSE; }
154  void setCalcErrorFromInputVariance(bool calcErrorFromInputVariance) noexcept {
155  _calcErrorFromInputVariance = calcErrorFromInputVariance;
156  }
157 
158 private:
159  friend class Statistics;
160 
161  double _numSigmaClip; // Number of standard deviations to clip at
162  int _numIter; // Number of iterations
163  int _andMask; // and-Mask to specify which mask planes to ignore
164  int _noGoodPixelsMask; // mask to set if no values are acceptable
165  bool _isNanSafe; // Check for NaNs & Infs before running (slower)
166  WeightsBoolean _useWeights; // Calculate weighted statistics (enum because of 3-valued logic)
167  bool _calcErrorFromInputVariance; // Calculate errors from the input variances, if available
168  std::vector<double> _maskPropagationThresholds; // Thresholds for when to propagate mask bits,
169  // treated like a dict (unset bits are set to 1.0)
170 };
171 
215 class Statistics final {
216 public:
219 
232  template <typename ImageT, typename MaskT, typename VarianceT>
233  explicit Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags,
234  StatisticsControl const &sctrl = StatisticsControl());
235 
244  template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
245  explicit Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
246  int const flags, StatisticsControl const &sctrl = StatisticsControl());
247 
248  Statistics(Statistics const &) = default;
249  Statistics(Statistics &&) = default;
250  Statistics &operator=(Statistics const &) = default;
251  Statistics &operator=(Statistics &&) = default;
252  ~Statistics() noexcept = default;
253 
267  Value getResult(Property const prop = NOTHING) const;
268 
276  double getError(Property const prop = NOTHING) const;
282  double getValue(Property const prop = NOTHING) const;
283  lsst::afw::image::MaskPixel getOrMask() const noexcept { return _allPixelOrMask; }
284 
285 private:
286  long _flags; // The desired calculation
287 
288  int _n; // number of pixels in the image
289  Value _mean; // the image's mean
290  Value _variance; // the image's variance
291  double _min; // the image's minimum
292  double _max; // the image's maximum
293  double _sum; // the sum of all the image's pixels
294  Value _meanclip; // the image's N-sigma clipped mean
295  Value _varianceclip; // the image's N-sigma clipped variance
296  Value _median; // the image's median
297  int _nClipped; // number of pixels clipped
298  int _nMasked; // number of pixels masked
299  double _iqrange; // the image's interquartile range
300  lsst::afw::image::MaskPixel _allPixelOrMask; // the 'or' of all masked pixels
301 
302  StatisticsControl _sctrl; // the control structure
303  bool _weightsAreMultiplicative; // Multiply by weights rather than dividing by them
304 
313  template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
314  void doStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
315  int const flags, StatisticsControl const &sctrl);
316 };
317 
318 /* ************************************ The factory functions ********************************* */
323 template <typename ValueT>
324 class infinite_iterator : public boost::iterator_adaptor<infinite_iterator<ValueT>, const ValueT *,
325  const ValueT, boost::forward_traversal_tag> {
326 public:
327  infinite_iterator() : infinite_iterator::iterator_adaptor_(0) {}
328  explicit infinite_iterator(const ValueT *p) : infinite_iterator::iterator_adaptor_(p) {}
329 
330 private:
331  friend class boost::iterator_core_access;
332  void increment() noexcept { ; } // never actually advance the iterator
333 };
338 template <typename ValueT>
340 public:
342  explicit MaskImposter(ValueT val = 0) noexcept { _val[0] = val; }
343  x_iterator row_begin(int) const noexcept { return x_iterator(_val); }
344 
345 private:
346  ValueT _val[1];
347 };
348 
353 template <typename Pixel>
355  lsst::afw::image::Mask<image::MaskPixel> const &msk, int const flags,
356  StatisticsControl const &sctrl = StatisticsControl()) {
358  return Statistics(img, msk, var, flags, sctrl);
359 }
360 
365 template <typename ImageT, typename MaskT, typename VarianceT>
366 Statistics makeStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags,
367  StatisticsControl const &sctrl = StatisticsControl()) {
368  return Statistics(img, msk, var, flags, sctrl);
369 }
370 
375 template <typename Pixel>
377  StatisticsControl const &sctrl = StatisticsControl()) {
378  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
379  return Statistics(*mimg.getImage(), *mimg.getMask(), *mimg.getVariance(), flags, sctrl);
380  } else {
382  return Statistics(*mimg.getImage(), *mimg.getMask(), var, flags, sctrl);
383  }
384 }
385 
390 template <typename Pixel>
392  lsst::afw::image::Image<WeightPixel> const &weights, int const flags,
393  StatisticsControl const &sctrl = StatisticsControl()) {
394  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance() ||
395  (!sctrl.getWeightedIsSet() && (weights.getWidth() != 0 && weights.getHeight() != 0))) {
396  return Statistics(*mimg.getImage(), *mimg.getMask(), *mimg.getVariance(), weights, flags, sctrl);
397  } else {
399  return Statistics(*mimg.getImage(), *mimg.getMask(), var, weights, flags, sctrl);
400  }
401 }
402 
413  StatisticsControl const &sctrl = StatisticsControl());
414 
419 template <typename Pixel>
421  lsst::afw::image::Image<Pixel> const &img,
422  int const flags,
423  StatisticsControl const &sctrl = StatisticsControl()
424 ) {
425  // make a phony mask that will be compiled out
427  MaskImposter<WeightPixel> const var;
428  return Statistics(img, msk, var, flags, sctrl);
429 }
430 
435 template <typename ValueT>
436 class ImageImposter final {
437 public:
438  // types we'll use in Statistics
441  typedef ValueT Pixel;
442 
443  // constructors for std::vector<>, and copy constructor
444  // These are both shallow! ... no actual copying of values
445  explicit ImageImposter(std::vector<ValueT> const &v) : _v(v) {}
446  explicit ImageImposter(ImageImposter<ValueT> const &img) : _v(img._getVector()) {}
447 
448  // The methods we'll use in Statistics
449  x_iterator row_begin(int) const noexcept { return _v.begin(); }
450  x_iterator row_end(int) const noexcept { return _v.end(); }
451  int getWidth() const noexcept { return _v.size(); }
452  int getHeight() const noexcept { return 1; }
454  return lsst::geom::Extent2I(getWidth(), getHeight());
455  }
456 
457  bool empty() const noexcept { return _v.empty(); }
458 
459 private:
460  std::vector<ValueT> const &_v; // a private reference to the data
461  std::vector<ValueT> const &_getVector() const { return _v; } // get the ref for the copyCon
462 };
463 
468 template <typename EntryT>
470  int const flags,
471  StatisticsControl const &sctrl = StatisticsControl()
472 ) {
473  ImageImposter<EntryT> img(v); // wrap the vector in a fake image
474  MaskImposter<lsst::afw::image::MaskPixel> msk; // instantiate a fake mask that will be compiled out.
476  return Statistics(img, msk, var, flags, sctrl);
477 }
478 
483 template <typename EntryT>
485  std::vector<WeightPixel> const &vweights,
486  int const flags,
487  StatisticsControl const &sctrl = StatisticsControl()
488 ) {
489  ImageImposter<EntryT> img(v); // wrap the vector in a fake image
490  MaskImposter<lsst::afw::image::MaskPixel> msk; // instantiate a fake mask that will be compiled out.
492 
493  ImageImposter<WeightPixel> weights(vweights);
494 
495  return Statistics(img, msk, var, weights, flags, sctrl);
496 }
497 
502 template <typename EntryT>
504  int const flags,
505  StatisticsControl const &sctrl = StatisticsControl()
506 ) {
507  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
508  return Statistics(*mv.getImage(), *mv.getMask(), *mv.getVariance(), flags, sctrl);
509  } else {
511  return Statistics(*mv.getImage(), *mv.getMask(), var, flags, sctrl);
512  }
513 }
514 
519 template <typename EntryT>
521  std::vector<WeightPixel> const &vweights,
522  int const flags,
523  StatisticsControl const &sctrl = StatisticsControl()
524 ) {
525  ImageImposter<WeightPixel> weights(vweights);
526 
527  if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
528  return Statistics(*mv.getImage(), *mv.getMask(), *mv.getVariance(), weights, flags, sctrl);
529  } else {
531  return Statistics(*mv.getImage(), *mv.getMask(), var, weights, flags, sctrl);
532  }
533 }
534 } // namespace math
535 } // namespace afw
536 } // namespace lsst
537 
538 #endif
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:1088
double getNumSigmaClip() const noexcept
Definition: Statistics.h:133
MaskImposter(ValueT val=0) noexcept
Definition: Statistics.h:342
ImageImposter(ImageImposter< ValueT > const &img)
Definition: Statistics.h:446
T empty(T... args)
std::vector< ValueT >::const_iterator x_iterator
Definition: Statistics.h:439
bool getNanSafe() const noexcept
Definition: Statistics.h:137
x_iterator row_begin(int) const noexcept
Definition: Statistics.h:343
lsst::afw::image::VariancePixel WeightPixel
Definition: Statistics.h:56
get the or-mask of all pixels used.
Definition: Statistics.h:80
void setNumSigmaClip(double numSigmaClip)
Definition: Statistics.h:142
estimate sample minimum
Definition: Statistics.h:76
bool getCalcErrorFromInputVariance() const noexcept
Definition: Statistics.h:140
bool empty() const noexcept
Definition: Statistics.h:457
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:315
number of clipped points
Definition: Statistics.h:81
estimate sample standard deviation
Definition: Statistics.h:68
int getHeight() const noexcept
Definition: Statistics.h:452
find mean value of square of pixel values
Definition: Statistics.h:79
ImageT val
Definition: CR.cc:146
int getNumIter() const noexcept
Definition: Statistics.h:134
estimate sample maximum
Definition: Statistics.h:77
T end(T... args)
lsst::geom::Extent2I getDimensions() const noexcept
Definition: Statistics.h:453
void setNoGoodPixelsMask(int noGoodPixelsMask)
Definition: Statistics.h:151
We don&#39;t want anything.
Definition: Statistics.h:64
lsst::afw::image::MaskedImage< EntryT >::ImagePtr getImage() const
Definition: MaskedVector.h:88
Include errors of requested quantities.
Definition: Statistics.h:65
x_iterator row_end(int) const noexcept
Definition: Statistics.h:450
A class to evaluate image statistics.
Definition: Statistics.h:215
void setNanSafe(bool isNanSafe) noexcept
Definition: Statistics.h:152
number of masked points
Definition: Statistics.h:82
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:1055
STL class.
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<>
Definition: Statistics.h:520
std::vector< ValueT >::const_iterator fast_iterator
Definition: Statistics.h:440
estimate sample N-sigma clipped stdev (N set in StatisticsControl, default=3)
Definition: Statistics.h:73
std::pair< double, double > Value
The type used to report (value, error) for desired statistics.
Definition: Statistics.h:218
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:747
infinite_iterator< ValueT > x_iterator
Definition: Statistics.h:341
StatisticsControl(double numSigmaClip=3.0, int numIter=3, lsst::afw::image::MaskPixel andMask=0x0, bool isNanSafe=true, WeightsBoolean useWeights=WEIGHTS_NONE)
Definition: Statistics.h:97
estimate sample median
Definition: Statistics.h:70
Pass parameters to a Statistics object.
Definition: Statistics.h:93
estimate sample N-sigma clipped mean (N set in StatisticsControl, default=3)
Definition: Statistics.h:72
A base class for image defects.
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
estimate sample N-sigma clipped variance (N set in StatisticsControl, default=3)
Definition: Statistics.h:74
A Mask wrapper to provide an infinite_iterator for Mask::row_begin().
Definition: Statistics.h:339
int getWidth() const noexcept
Definition: Statistics.h:451
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
number of sample points
Definition: Statistics.h:66
x_iterator row_begin(int) const noexcept
Definition: Statistics.h:449
bool getWeightedIsSet() const noexcept
Definition: Statistics.h:139
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:379
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
int getAndMask() const noexcept
Definition: Statistics.h:135
lsst::afw::image::MaskedImage< EntryT >::VariancePtr getVariance() const
Definition: MaskedVector.h:94
This iterator will never increment.
Definition: Statistics.h:324
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:1067
lsst::afw::image::MaskedImage< EntryT >::MaskPtr getMask() const
Definition: MaskedVector.h:91
T size(T... args)
Extent< int, 2 > Extent2I
Definition: Extent.h:397
infinite_iterator(const ValueT *p)
Definition: Statistics.h:328
estimate sample inter-quartile range
Definition: Statistics.h:71
T begin(T... args)
estimate sample mean
Definition: Statistics.h:67
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:313
estimate sample variance
Definition: Statistics.h:69
void setCalcErrorFromInputVariance(bool calcErrorFromInputVariance) noexcept
Definition: Statistics.h:154
Reports invalid arguments.
Definition: Runtime.h:66
A vector wrapper to provide a vector with the necessary methods and typedefs to be processed by Stati...
Definition: Statistics.h:436
void setWeighted(bool useWeights) noexcept
Definition: Statistics.h:153
find sum of pixels in the image
Definition: Statistics.h:78
bool getWeighted() const noexcept
Definition: Statistics.h:138
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
int getNoGoodPixelsMask() const noexcept
Definition: Statistics.h:136
float VariancePixel
default type for MaskedImage variance images
Property
control what is calculated
Definition: Statistics.h:63
ImageImposter(std::vector< ValueT > const &v)
Definition: Statistics.h:445