LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
Statistics.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 
25 /*
26  * Support statistical operations on images
27  */
28 #include <cassert>
29 #include <cmath>
30 #include <tuple>
31 #include <type_traits>
32 #include <string>
33 #include <vector>
34 #include <utility>
35 
36 #include "lsst/pex/exceptions.h"
37 #include "lsst/afw/image/Image.h"
39 #include "lsst/geom/Angle.h"
40 
41 using namespace std;
43 
44 namespace lsst {
45 namespace afw {
46 namespace math {
47 
48 namespace {
49 double const NaN = std::numeric_limits<double>::quiet_NaN();
50 double const MAX_DOUBLE = std::numeric_limits<double>::max();
51 double const IQ_TO_STDEV = 0.741301109252802; // 1 sigma in units of iqrange (assume Gaussian)
52 
54 class AlwaysTrue {
55 public:
56  template <typename T>
57  bool operator()(T) const {
58  return true;
59  }
60  template <typename Ta, typename Tb>
61  bool operator()(Ta, Tb) const {
62  return true;
63  }
64  template <typename Ta, typename Tb, typename Tc>
65  bool operator()(Ta, Tb, Tc) const {
66  return true;
67  }
68 };
69 
71 class AlwaysFalse {
72 public:
73  template <typename T>
74  bool operator()(T) const {
75  return false;
76  }
77  template <typename Ta, typename Tb>
78  bool operator()(Ta, Tb) const {
79  return false;
80  }
81  template <typename Ta, typename Tb, typename Tc>
82  bool operator()(Ta, Tb, Tc) const {
83  return false;
84  }
85 };
86 
88 class CheckFinite {
89 public:
90  template <typename T>
91  bool operator()(T val) const {
92  return std::isfinite(static_cast<float>(val));
93  }
94 };
95 
97 class CheckValueLtMin {
98 public:
99  template <typename Tval, typename Tmin>
100  bool operator()(Tval val, Tmin min) const {
101  return (static_cast<Tmin>(val) < min);
102  }
103 };
104 
106 class CheckValueGtMax {
107 public:
108  template <typename Tval, typename Tmax>
109  bool operator()(Tval val, Tmax max) const {
110  return (static_cast<Tmax>(val) > max);
111  }
112 };
113 
115 class CheckClipRange {
116 public:
117  template <typename Tval, typename Tcen, typename Tmax>
118  bool operator()(Tval val, Tcen center, Tmax cliplimit) const {
119  Tmax tmp = fabs(val - center);
120  return (tmp <= cliplimit);
121  }
122 };
123 
124 // define some abbreviated typenames for the test templates
125 using ChkFin = CheckFinite;
126 using ChkMin = CheckValueLtMin;
127 using ChkMax = CheckValueGtMax;
128 using ChkClip = CheckClipRange;
129 using AlwaysT = AlwaysTrue;
130 using AlwaysF = AlwaysFalse;
131 
135 inline double varianceError(double const variance, int const n) {
136  return 2 * (n - 1) * variance * variance / static_cast<double>(n * n);
137 }
138 
141 
142 /*
143  * Functions which convert the booleans into calls to the proper templated types, one type per
144  * recursion level
145  */
155 template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
156  bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
157 StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
158  WeightT const &weights, int const, int const nCrude, int const stride,
159  double const meanCrude, double const cliplimit,
160  bool const weightsAreMultiplicative, int const andMask,
161  bool const calcErrorFromInputVariance,
162  std::vector<double> const &maskPropagationThresholds) {
163  int n = 0;
164  double sumw = 0.0; // sum(weight) (N.b. weight will be 1.0 if !useWeights)
165  double sumw2 = 0.0; // sum(weight^2)
166  double sumx = 0; // sum(data*weight)
167  double sumx2 = 0; // sum(data*weight^2)
168 #if 1
169  double sumvw2 = 0.0; // sum(variance*weight^2)
170 #endif
171  double min = (nCrude) ? meanCrude : MAX_DOUBLE;
172  double max = (nCrude) ? meanCrude : -MAX_DOUBLE;
173 
174  image::MaskPixel allPixelOrMask = 0x0;
175 
176  std::vector<double> rejectedWeightsByBit(maskPropagationThresholds.size(), 0.0);
177 
178  for (int iY = 0; iY < img.getHeight(); iY += stride) {
179  typename MaskT::x_iterator mptr = msk.row_begin(iY);
180  typename VarianceT::x_iterator vptr = var.row_begin(iY);
181  typename WeightT::x_iterator wptr = weights.row_begin(iY);
182 
183  for (typename ImageT::x_iterator ptr = img.row_begin(iY), end = ptr + img.getWidth(); ptr != end;
184  ++ptr, ++mptr, ++vptr, ++wptr) {
185  if (IsFinite()(*ptr) && !(*mptr & andMask) &&
186  InClipRange()(*ptr, meanCrude, cliplimit)) { // clip
187 
188  double const delta = (*ptr - meanCrude);
189 
190  if (useWeights) {
191  double weight = *wptr;
192  if (weightsAreMultiplicative) {
193  ;
194  } else {
195  if (*wptr <= 0) {
196  continue;
197  }
198  weight = 1 / weight;
199  }
200 
201  sumw += weight;
202  sumw2 += weight * weight;
203  sumx += weight * delta;
204  sumx2 += weight * delta * delta;
205 
206  if (calcErrorFromInputVariance) {
207  double const var = *vptr;
208  sumvw2 += var * weight * weight;
209  }
210  } else {
211  sumx += delta;
212  sumx2 += delta * delta;
213 
214  if (calcErrorFromInputVariance) {
215  double const var = *vptr;
216  sumvw2 += var;
217  }
218  }
219 
220  allPixelOrMask |= *mptr;
221 
222  if (HasValueLtMin()(*ptr, min)) {
223  min = *ptr;
224  }
225  if (HasValueGtMax()(*ptr, max)) {
226  max = *ptr;
227  }
228  n++;
229  } else { // pixel has been clipped, rejected, etc.
230  for (int bit = 0, nBits = maskPropagationThresholds.size(); bit < nBits; ++bit) {
231  image::MaskPixel mask = 1 << bit;
232  if (*mptr & mask) {
233  double weight = 1.0;
234  if (useWeights) {
235  weight = *wptr;
236  if (!weightsAreMultiplicative) {
237  if (*wptr <= 0) {
238  continue;
239  }
240  weight = 1.0 / weight;
241  }
242  }
243  rejectedWeightsByBit[bit] += weight;
244  }
245  }
246  }
247  }
248  }
249  if (n == 0) {
250  min = NaN;
251  max = NaN;
252  }
253 
254  // estimate of population mean and variance.
255  double mean, variance;
256  if (!useWeights) {
257  sumw = sumw2 = n;
258  }
259 
260  for (int bit = 0, nBits = maskPropagationThresholds.size(); bit < nBits; ++bit) {
261  double hypotheticalTotalWeight = sumw + rejectedWeightsByBit[bit];
262  rejectedWeightsByBit[bit] /= hypotheticalTotalWeight;
263  if (rejectedWeightsByBit[bit] > maskPropagationThresholds[bit]) {
264  allPixelOrMask |= (1 << bit);
265  }
266  }
267 
268  // N.b. if sumw == 0 or sumw*sumw == sumw2 (e.g. n == 1) we'll get NaNs
269  // N.b. the estimator of the variance assumes that the sample points all have the same variance;
270  // otherwise, what is it that we're estimating?
271  mean = sumx / sumw;
272  variance = sumx2 / sumw - ::pow(mean, 2); // biased estimator
273  variance *= sumw * sumw / (sumw * sumw - sumw2); // debias
274 
275  double meanVar; // (standard error of mean)^2
276  if (calcErrorFromInputVariance) {
277  meanVar = sumvw2 / (sumw * sumw);
278  } else {
279  meanVar = variance * sumw2 / (sumw * sumw);
280  }
281 
282  double varVar = varianceError(variance, n); // error in variance; incorrect if useWeights is true
283 
284  sumx += sumw * meanCrude;
285  mean += meanCrude;
286 
287  return StandardReturn(n, sumx, Statistics::Value(mean, meanVar), Statistics::Value(variance, varVar), min,
288  max, allPixelOrMask);
289 }
290 
291 template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
292  bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
293 StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
294  WeightT const &weights, int const flags, int const nCrude, int const stride,
295  double const meanCrude, double const cliplimit,
296  bool const weightsAreMultiplicative, int const andMask,
297  bool const calcErrorFromInputVariance, bool doGetWeighted,
298  std::vector<double> const &maskPropagationThresholds) {
299  if (doGetWeighted) {
300  return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, true>(
301  img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
302  andMask, calcErrorFromInputVariance, maskPropagationThresholds);
303  } else {
304  return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, false>(
305  img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
306  andMask, calcErrorFromInputVariance, maskPropagationThresholds);
307  }
308 }
309 
310 template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
311  bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
312 StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
313  WeightT const &weights, int const flags, int const nCrude, int const stride,
314  double const meanCrude, double const cliplimit,
315  bool const weightsAreMultiplicative, int const andMask,
316  bool const calcErrorFromInputVariance, bool doCheckFinite, bool doGetWeighted,
317  std::vector<double> const &maskPropagationThresholds) {
318  if (doCheckFinite) {
319  return processPixels<CheckFinite, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
320  img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
321  andMask, calcErrorFromInputVariance, doGetWeighted, maskPropagationThresholds);
322  } else {
323  return processPixels<AlwaysTrue, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
324  img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
325  andMask, calcErrorFromInputVariance, doGetWeighted, maskPropagationThresholds);
326  }
327 }
328 
346 template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
347 StandardReturn getStandard(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
348  int const flags, bool const weightsAreMultiplicative, int const andMask,
349  bool const calcErrorFromInputVariance, bool doCheckFinite, bool doGetWeighted,
350  std::vector<double> const &maskPropagationThresholds) {
351  // =====================================================
352  // a crude estimate of the mean, used for numerical stability of variance
353  int nCrude = 0;
354  double meanCrude = 0.0;
355 
356  // for small numbers of values, use a small stride
357  int const nPix = img.getWidth() * img.getHeight();
358  int strideCrude;
359  if (nPix < 100) {
360  strideCrude = 2;
361  } else {
362  strideCrude = 10;
363  }
364 
365  double cliplimit = -1; // unused
366  StandardReturn values = processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
367  img, msk, var, weights, flags, nCrude, strideCrude, meanCrude, cliplimit,
368  weightsAreMultiplicative, andMask, calcErrorFromInputVariance, doCheckFinite, doGetWeighted,
369  maskPropagationThresholds);
370  nCrude = std::get<0>(values);
371  double sumCrude = std::get<1>(values);
372 
373  meanCrude = 0.0;
374  if (nCrude > 0) {
375  meanCrude = sumCrude / nCrude;
376  }
377 
378  // =======================================================
379  // Estimate the full precision variance using that crude mean
380  // - get the min and max as well
381 
382  if (flags & (MIN | MAX)) {
383  return processPixels<ChkFin, ChkMin, ChkMax, AlwaysT, true>(
384  img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
385  andMask, calcErrorFromInputVariance, true, doGetWeighted, maskPropagationThresholds);
386  } else {
387  return processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
388  img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
389  andMask, calcErrorFromInputVariance, doCheckFinite, doGetWeighted, maskPropagationThresholds);
390  }
391 }
392 
411 template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
412 StandardReturn getStandard(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
413  int const flags, std::pair<double, double> const clipinfo,
414 
415  bool const weightsAreMultiplicative, int const andMask,
416  bool const calcErrorFromInputVariance, bool doCheckFinite, bool doGetWeighted,
417  std::vector<double> const &maskPropagationThresholds) {
418  double const center = clipinfo.first;
419  double const cliplimit = clipinfo.second;
420 
421  if (std::isnan(center) || std::isnan(cliplimit)) {
422  return StandardReturn(0, NaN, Statistics::Value(NaN, NaN), Statistics::Value(NaN, NaN), NaN, NaN,
423  ~0x0);
424  }
425 
426  // =======================================================
427  // Estimate the full precision variance using that crude mean
428  int const stride = 1;
429  int nCrude = 0;
430 
431  if (flags & (MIN | MAX)) {
432  return processPixels<ChkFin, ChkMin, ChkMax, ChkClip, true>(
433  img, msk, var, weights, flags, nCrude, stride, center, cliplimit, weightsAreMultiplicative,
434  andMask, calcErrorFromInputVariance, true, doGetWeighted, maskPropagationThresholds);
435  } else { // fast loop ... just the mean & variance
436  return processPixels<ChkFin, AlwaysF, AlwaysF, ChkClip, true>(
437  img, msk, var, weights, flags, nCrude, stride, center, cliplimit, weightsAreMultiplicative,
438  andMask, calcErrorFromInputVariance, doCheckFinite, doGetWeighted, maskPropagationThresholds);
439  }
440 }
441 
450 template <typename Pixel>
452 percentile(std::vector<Pixel> &img, double const fraction) {
453  assert(fraction >= 0.0 && fraction <= 1.0);
454 
455  int const n = img.size();
456 
457  if (n > 1) {
458  double const idx = fraction * (n - 1);
459 
460  // interpolate linearly between the adjacent values
461  // For efficiency:
462  // - if we're asked for a fraction > 0.5,
463  // we'll do the second partial sort on shorter (upper) portion
464  // - otherwise, the shorter portion will be the lower one, we'll partial-sort that.
465 
466  int const q1 = static_cast<int>(idx);
467  int const q2 = q1 + 1;
468 
469  auto mid1 = img.begin() + q1;
470  auto mid2 = img.begin() + q2;
471  if (fraction > 0.5) {
472  std::nth_element(img.begin(), mid1, img.end());
473  std::nth_element(mid1, mid2, img.end());
474  } else {
475  std::nth_element(img.begin(), mid2, img.end());
476  std::nth_element(img.begin(), mid1, mid2);
477  }
478 
479  double val1 = static_cast<double>(*mid1);
480  double val2 = static_cast<double>(*mid2);
481  double w1 = (static_cast<double>(q2) - idx);
482  double w2 = (idx - static_cast<double>(q1));
483  return w1 * val1 + w2 * val2;
484 
485  } else if (n == 1) {
486  return img[0];
487  } else {
488  return NaN;
489  }
490 }
491 
492 namespace {
493  //
494  // Helper function to estimate a floating-point quantile from integer data
495  //
496  // The data has been partially sorted, using nth_element
497  //
498  // begin: iterator to a point which we know to be below our median
499  // end: iterator to a point which we know to be beyond our median
500  // naive: the integer value of the desired quantile
501  // target: the number of points that should be to the left of the quantile.
502  // N.b. if begin isn't the start of the data, this may not be the
503  // desired number of points. Caveat Callor
504  template<typename T>
505  double
506  computeQuantile(typename std::vector<T>::const_iterator begin,
508  T const naive,
509  double const target
510  )
511  {
512  // investigate the cumulative histogram near naive
513  std::size_t left = 0; // number of values less than naive
514  std::size_t middle = 0; // number of values equal to naive
515 
516  for (auto ptr = begin; ptr != end; ++ptr) {
517  auto const val = *ptr;
518  if (val < naive) {
519  ++left;
520  } else if (val == naive) {
521  ++middle;
522  }
523  }
524 
525  return naive - 0.5 + (target - left)/middle;
526  }
527 
536  template <typename Pixel>
538  percentile(std::vector<Pixel> &img, double const fraction) {
539  assert(fraction >= 0.0 && fraction <= 1.0);
540 
541  auto const n = img.size();
542 
543  if (n == 0) {
544  return NaN;
545  } else if (n == 1) {
546  return img[0];
547  } else {
548  // We need to handle ties. The proper way to do this is to analyse the cumulative curve after
549  // building the histograms (which is faster than a generic partitioning algorithm), but it's a
550  // nuisance as we don't know the range of values
551  //
552  // This code looks clean enough, but actually the call to nth_element is expensive
553  // and we *still* have to go through the array a second time
554 
555  double const idx = fraction*(n - 1);
556 
557  auto midP = img.begin() + static_cast<int>(idx);
558  std::nth_element(img.begin(), midP, img.end());
559  auto const naiveP = *midP; // value of desired element
560 
561  return computeQuantile(img.begin(), img.end(), naiveP, fraction*n);
562  }
563  }
564 }
565 
566 using MedianQuartileReturn = std::tuple<double, double, double>;
567 namespace {
575  template <typename Pixel>
576  typename enable_if<!is_integral<Pixel>::value, MedianQuartileReturn>::type
577  medianAndQuartiles(std::vector<Pixel> &img) {
578  int const n = img.size();
579 
580  if (n > 1) {
581  double const idx50 = 0.50 * (n - 1);
582  double const idx25 = 0.25 * (n - 1);
583  double const idx75 = 0.75 * (n - 1);
584 
585  // For efficiency:
586  // - partition at 50th, then partition the two half further to get 25th and 75th
587  // - to get the adjacent points (for interpolation), partition between 25/50, 50/75, 75/end
588  // these should be much smaller partitions
589 
590  int const q50a = static_cast<int>(idx50);
591  int const q50b = q50a + 1;
592  int const q25a = static_cast<int>(idx25);
593  int const q25b = q25a + 1;
594  int const q75a = static_cast<int>(idx75);
595  int const q75b = q75a + 1;
596 
597  auto mid50a = img.begin() + q50a;
598  auto mid50b = img.begin() + q50b;
599  auto mid25a = img.begin() + q25a;
600  auto mid25b = img.begin() + q25b;
601  auto mid75a = img.begin() + q75a;
602  auto mid75b = img.begin() + q75b;
603 
604  // get the 50th percentile, then get the 25th and 75th on the smaller partitions
605  std::nth_element(img.begin(), mid50a, img.end());
606  std::nth_element(mid50a, mid75a, img.end());
607  std::nth_element(img.begin(), mid25a, mid50a);
608 
609  // and the adjacent points for each ... use the smallest segments available.
610  std::nth_element(mid50a, mid50b, mid75a);
611  std::nth_element(mid25a, mid25b, mid50a);
612  std::nth_element(mid75a, mid75b, img.end());
613 
614  // interpolate linearly between the adjacent values
615  double val50a = static_cast<double>(*mid50a);
616  double val50b = static_cast<double>(*mid50b);
617  double w50a = (static_cast<double>(q50b) - idx50);
618  double w50b = (idx50 - static_cast<double>(q50a));
619  double median = w50a * val50a + w50b * val50b;
620 
621  double val25a = static_cast<double>(*mid25a);
622  double val25b = static_cast<double>(*mid25b);
623  double w25a = (static_cast<double>(q25b) - idx25);
624  double w25b = (idx25 - static_cast<double>(q25a));
625  double q1 = w25a * val25a + w25b * val25b;
626 
627  double val75a = static_cast<double>(*mid75a);
628  double val75b = static_cast<double>(*mid75b);
629  double w75a = (static_cast<double>(q75b) - idx75);
630  double w75b = (idx75 - static_cast<double>(q75a));
631  double q3 = w75a * val75a + w75b * val75b;
632 
633  return MedianQuartileReturn(median, q1, q3);
634  } else if (n == 1) {
635  return MedianQuartileReturn(img[0], img[0], img[0]);
636  } else {
637  return MedianQuartileReturn(NaN, NaN, NaN);
638  }
639  }
640 
648  template <typename Pixel>
649  typename enable_if<is_integral<Pixel>::value, MedianQuartileReturn>::type
650  medianAndQuartiles(std::vector<Pixel> &img) {
651  auto const n = img.size();
652 
653  if (n == 0) {
654  return MedianQuartileReturn(NaN, NaN, NaN);
655  } else if (n == 1) {
656  return MedianQuartileReturn(img[0], img[0], img[0]);
657  } else {
658  // We need to handle ties. The proper way to do this is to analyse the cumulative curve after
659  // building the histograms (which is faster than a generic partitioning algorithm), but it's a
660  // nuisance as we don't know the range of values
661  //
662  // This code looks clean enough, but actually the call to nth_element is expensive
663  // and we *still* have to go through the array a second time.
664 
665  // For efficiency:
666  // - partition at 50th, then partition the two halves further to get 25th and 75th
667 
668  auto mid25 = img.begin() + static_cast<int>(0.25*(n - 1));
669  auto mid50 = img.begin() + static_cast<int>(0.50*(n - 1));
670  auto mid75 = img.begin() + static_cast<int>(0.75*(n - 1));
671 
672  // get the 50th percentile, then get the 25th and 75th on the smaller partitions
673  std::nth_element(img.begin(), mid50, img.end());
674  std::nth_element(img.begin(), mid25, mid50);
675  std::nth_element(mid50, mid75, img.end());
676 
677  double const q1 = computeQuantile(img.begin(), mid50, *mid25,
678  0.25*n);
679  double const median = computeQuantile(mid25, mid75, *mid50,
680  0.50*n - (mid25 - img.begin()));
681  double const q3 = computeQuantile(mid50, img.end(), *mid75,
682  0.75*n - (mid50 - img.begin()));
683 
684  return MedianQuartileReturn(median, q1, q3);
685  }
686  }
687 }
688 
696 template <typename IsFinite, typename ImageT, typename MaskT, typename VarianceT>
697 std::shared_ptr<std::vector<typename ImageT::Pixel> > makeVectorCopy(ImageT const &img, MaskT const &msk,
698  VarianceT const &, int const andMask) {
699  // Note need to keep track of allPixelOrMask here ... processPixels() does that
700  // and it always gets called
702 
703  for (int i_y = 0; i_y < img.getHeight(); ++i_y) {
704  typename MaskT::x_iterator mptr = msk.row_begin(i_y);
705  for (typename ImageT::x_iterator ptr = img.row_begin(i_y), end = img.row_end(i_y); ptr != end;
706  ++ptr) {
707  if (IsFinite()(*ptr) && !(*mptr & andMask)) {
708  imgcp->push_back(*ptr);
709  }
710  ++mptr;
711  }
712  }
713 
714  return imgcp;
715 }
716 } // namespace
717 
718 double StatisticsControl::getMaskPropagationThreshold(int bit) const {
719  int oldSize = _maskPropagationThresholds.size();
720  if (oldSize <= bit) {
721  return 1.0;
722  }
723  return _maskPropagationThresholds[bit];
724 }
725 
726 void StatisticsControl::setMaskPropagationThreshold(int bit, double threshold) {
727  int oldSize = _maskPropagationThresholds.size();
728  if (oldSize <= bit) {
729  int newSize = bit + 1;
730  _maskPropagationThresholds.resize(newSize);
731  for (int i = oldSize; i < bit; ++i) {
732  _maskPropagationThresholds[i] = 1.0;
733  }
734  }
735  _maskPropagationThresholds[bit] = threshold;
736 }
737 
739  static std::map<std::string, Property> statisticsProperty;
740  if (statisticsProperty.size() == 0) {
741  statisticsProperty["NOTHING"] = NOTHING;
742  statisticsProperty["ERRORS"] = ERRORS;
743  statisticsProperty["NPOINT"] = NPOINT;
744  statisticsProperty["MEAN"] = MEAN;
745  statisticsProperty["STDEV"] = STDEV;
746  statisticsProperty["VARIANCE"] = VARIANCE;
747  statisticsProperty["MEDIAN"] = MEDIAN;
748  statisticsProperty["IQRANGE"] = IQRANGE;
749  statisticsProperty["MEANCLIP"] = MEANCLIP;
750  statisticsProperty["STDEVCLIP"] = STDEVCLIP;
751  statisticsProperty["VARIANCECLIP"] = VARIANCECLIP;
752  statisticsProperty["MIN"] = MIN;
753  statisticsProperty["MAX"] = MAX;
754  statisticsProperty["SUM"] = SUM;
755  statisticsProperty["MEANSQUARE"] = MEANSQUARE;
756  statisticsProperty["ORMASK"] = ORMASK;
757  statisticsProperty["NCLIPPED"] = NCLIPPED;
758  statisticsProperty["NMASKED"] = NMASKED;
759  }
760  return statisticsProperty[property];
761 }
762 
763 template <typename ImageT, typename MaskT, typename VarianceT>
764 Statistics::Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags,
765  StatisticsControl const &sctrl)
766  : _flags(flags),
767  _mean(NaN, NaN),
768  _variance(NaN, NaN),
769  _min(NaN),
770  _max(NaN),
771  _sum(NaN),
772  _meanclip(NaN, NaN),
773  _varianceclip(NaN, NaN),
774  _median(NaN, NaN),
775  _nClipped(0),
776  _nMasked(0),
777  _iqrange(NaN),
778  _sctrl(sctrl),
779  _weightsAreMultiplicative(false) {
780  doStatistics(img, msk, var, var, _flags, _sctrl);
781 }
782 
783 namespace {
784 template <typename T>
785 bool isEmpty(T const &t) {
786  return t.empty();
787 }
788 
789 template <typename T>
790 bool isEmpty(image::Image<T> const &im) {
791  return (im.getWidth() == 0 && im.getHeight() == 0);
792 }
793 
794 // Asserts that image dimensions are equal
795 template <typename ImageT1, typename ImageT2>
796 void checkDimensions(ImageT1 const &image1, ImageT2 const &image2) {
797  if (image1.getDimensions() != image2.getDimensions()) {
799  (boost::format("Image sizes don't match: %s vs %s") % image1.getDimensions() %
800  image2.getDimensions())
801  .str());
802  }
803 }
804 
805 // Overloads for MaskImposter (which doesn't have a size)
806 template <typename ImageT, typename PixelT>
807 void checkDimensions(ImageT const &image1, MaskImposter<PixelT> const &image2) {}
808 template <typename ImageT, typename PixelT>
809 void checkDimensions(MaskImposter<PixelT> const &image1, ImageT const &image2) {}
810 } // namespace
811 
812 template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
813 Statistics::Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
814  int const flags, StatisticsControl const &sctrl)
815  : _flags(flags),
816  _mean(NaN, NaN),
817  _variance(NaN, NaN),
818  _min(NaN),
819  _max(NaN),
820  _sum(NaN),
821  _meanclip(NaN, NaN),
822  _varianceclip(NaN, NaN),
823  _median(NaN, NaN),
824  _nClipped(0),
825  _nMasked(0),
826  _iqrange(NaN),
827  _sctrl(sctrl),
828  _weightsAreMultiplicative(true) {
829  if (!isEmpty(weights)) {
830  if (_sctrl.getWeightedIsSet() && !_sctrl.getWeighted()) {
832  "You must use the weights if you provide them");
833  }
834 
835  _sctrl.setWeighted(true);
836  }
837  doStatistics(img, msk, var, weights, _flags, _sctrl);
838 }
839 
840 template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
841 void Statistics::doStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var,
842  WeightT const &weights, int const flags, StatisticsControl const &sctrl) {
843  int const num = img.getWidth() * img.getHeight();
844  _n = num;
845  if (_n == 0) {
846  throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
847  }
848  checkDimensions(img, msk);
849  checkDimensions(img, var);
850  if (sctrl.getWeighted()) {
851  checkDimensions(img, weights);
852  }
853 
854  // Check that an int's large enough to hold the number of pixels
855  assert(img.getWidth() * static_cast<double>(img.getHeight()) < std::numeric_limits<int>::max());
856 
857  // get the standard statistics
858  StandardReturn standard =
859  getStandard(img, msk, var, weights, flags, _weightsAreMultiplicative, _sctrl.getAndMask(),
860  _sctrl.getCalcErrorFromInputVariance(), _sctrl.getNanSafe(), _sctrl.getWeighted(),
861  _sctrl._maskPropagationThresholds);
862 
863  _n = std::get<0>(standard);
864  _sum = std::get<1>(standard);
865  _mean = std::get<2>(standard);
866  _variance = std::get<3>(standard);
867  _min = std::get<4>(standard);
868  _max = std::get<5>(standard);
869  _allPixelOrMask = std::get<6>(standard);
870 
871  // ==========================================================
872  // now only calculate it if it's specifically requested - these all cost more!
873 
874  if (flags & NMASKED) {
875  _nMasked = num - _n;
876  }
877 
878  // copy the image for any routines that will use median or quantiles
879  if (flags & (MEDIAN | IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
880  // make a vector copy of the image to get the median and quartiles (will move values)
882  if (_sctrl.getNanSafe()) {
883  imgcp = makeVectorCopy<ChkFin>(img, msk, var, _sctrl.getAndMask());
884  } else {
885  imgcp = makeVectorCopy<AlwaysT>(img, msk, var, _sctrl.getAndMask());
886  }
887 
888  // if we *only* want the median, just use percentile(), otherwise use medianAndQuartiles()
889  if ((flags & (MEDIAN)) && !(flags & (IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP))) {
890  _median = Value(percentile(*imgcp, 0.5), NaN);
891  } else {
892  MedianQuartileReturn mq = medianAndQuartiles(*imgcp);
893  _median = Value(std::get<0>(mq), NaN);
894  _iqrange = std::get<2>(mq) - std::get<1>(mq);
895  }
896 
897  if (flags & (MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
898  for (int i_i = 0; i_i < _sctrl.getNumIter(); ++i_i) {
899  double const center = ((i_i > 0) ? _meanclip : _median).first;
900  double const hwidth = (i_i > 0 && _n > 1)
901  ? _sctrl.getNumSigmaClip() * std::sqrt(_varianceclip.first)
902  : _sctrl.getNumSigmaClip() * IQ_TO_STDEV * _iqrange;
903  std::pair<double, double> const clipinfo(center, hwidth);
904 
905  StandardReturn clipped = getStandard(
906  img, msk, var, weights, flags, clipinfo, _weightsAreMultiplicative,
907  _sctrl.getAndMask(), _sctrl.getCalcErrorFromInputVariance(), _sctrl.getNanSafe(),
908  _sctrl.getWeighted(), _sctrl._maskPropagationThresholds);
909 
910  int const nClip = std::get<0>(clipped); // number after clipping
911  _nClipped = _n - nClip; // number clipped
912  _meanclip = std::get<2>(clipped); // clipped mean
913  double const varClip = std::get<3>(clipped).first; // clipped variance
914 
915  _varianceclip = Value(varClip, varianceError(varClip, nClip));
916  // ... ignore other values
917  }
918  }
919  }
920 }
921 
923  // if iProp == NOTHING try to return their heart's delight, as specified in the constructor
924  Property const prop = static_cast<Property>(((iProp == NOTHING) ? _flags : iProp) & ~ERRORS);
925 
926  if (!(prop & _flags)) { // we didn't calculate it
928  (boost::format("You didn't ask me to calculate %d") % prop).str());
929  }
930 
931  Value ret(NaN, NaN);
932  switch (prop) {
933  case NPOINT:
934  ret.first = static_cast<double>(_n);
935  if (_flags & ERRORS) {
936  ret.second = 0;
937  }
938  break;
939 
940  case NCLIPPED:
941  ret.first = static_cast<double>(_nClipped);
942  if (_flags & ERRORS) {
943  ret.second = 0;
944  }
945  break;
946 
947  case NMASKED:
948  ret.first = static_cast<double>(_nMasked);
949  if (_flags & ERRORS) {
950  ret.second = 0;
951  }
952  break;
953 
954  case SUM:
955  ret.first = static_cast<double>(_sum);
956  if (_flags & ERRORS) {
957  ret.second = 0;
958  }
959  break;
960 
961  // == means ==
962  case MEAN:
963  ret.first = _mean.first;
964  if (_flags & ERRORS) {
965  ret.second = ::sqrt(_mean.second);
966  }
967  break;
968  case MEANCLIP:
969  ret.first = _meanclip.first;
970  if (_flags & ERRORS) {
971  ret.second = ::sqrt(_meanclip.second);
972  }
973  break;
974 
975  // == stdevs & variances ==
976  case VARIANCE:
977  ret.first = _variance.first;
978  if (_flags & ERRORS) {
979  ret.second = ::sqrt(_variance.second);
980  }
981  break;
982  case STDEV:
983  ret.first = sqrt(_variance.first);
984  if (_flags & ERRORS) {
985  ret.second = 0.5 * ::sqrt(_variance.second) / ret.first;
986  }
987  break;
988  case VARIANCECLIP:
989  ret.first = _varianceclip.first;
990  if (_flags & ERRORS) {
991  ret.second = ret.second;
992  }
993  break;
994  case STDEVCLIP:
995  ret.first = sqrt(_varianceclip.first);
996  if (_flags & ERRORS) {
997  ret.second = 0.5 * ::sqrt(_varianceclip.second) / ret.first;
998  }
999  break;
1000 
1001  case MEANSQUARE:
1002  ret.first = (_n - 1) / static_cast<double>(_n) * _variance.first + ::pow(_mean.first, 2);
1003  if (_flags & ERRORS) {
1004  ret.second = ::sqrt(2 * ::pow(ret.first / _n, 2)); // assumes Gaussian
1005  }
1006  break;
1007 
1008  // == other stats ==
1009  case MIN:
1010  ret.first = _min;
1011  if (_flags & ERRORS) {
1012  ret.second = 0;
1013  }
1014  break;
1015  case MAX:
1016  ret.first = _max;
1017  if (_flags & ERRORS) {
1018  ret.second = 0;
1019  }
1020  break;
1021  case MEDIAN:
1022  ret.first = _median.first;
1023  if (_flags & ERRORS) {
1024  ret.second = sqrt(geom::HALFPI * _variance.first / _n); // assumes Gaussian
1025  }
1026  break;
1027  case IQRANGE:
1028  ret.first = _iqrange;
1029  if (_flags & ERRORS) {
1030  ret.second = NaN; // we're not estimating this properly
1031  }
1032  break;
1033 
1034  // no-op to satisfy the compiler
1035  case ERRORS:
1036  break;
1037  // default: redundant as 'ret' is initialized to NaN, NaN
1038  default: // we must have set prop to _flags
1039  assert(iProp == 0);
1041  "getValue() may only be called without a parameter"
1042  " if you asked for only one statistic");
1043  }
1044  return ret;
1045 }
1046 
1047 double Statistics::getValue(Property const prop) const { return getResult(prop).first; }
1048 
1049 double Statistics::getError(Property const prop) const { return getResult(prop).second; }
1050 
1060 template <>
1062  image::Mask<image::MaskPixel> const &var, int const flags,
1063  StatisticsControl const &sctrl)
1064  : _flags(flags),
1065  _mean(NaN, NaN),
1066  _variance(NaN, NaN),
1067  _min(NaN),
1068  _max(NaN),
1069  _meanclip(NaN, NaN),
1070  _varianceclip(NaN, NaN),
1071  _median(NaN, NaN),
1072  _nClipped(0),
1073  _iqrange(NaN),
1074  _sctrl(sctrl) {
1075  if ((flags & ~(NPOINT | SUM)) != 0x0) {
1077  "Statistics<Mask> only supports NPOINT and SUM");
1078  }
1079 
1080  using Mask = image::Mask<>;
1081 
1082  _n = msk.getWidth() * msk.getHeight();
1083  if (_n == 0) {
1084  throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
1085  }
1086 
1087  // Check that an int's large enough to hold the number of pixels
1088  assert(msk.getWidth() * static_cast<double>(msk.getHeight()) < std::numeric_limits<int>::max());
1089 
1090  image::MaskPixel sum = 0x0;
1091  for (int y = 0; y != msk.getHeight(); ++y) {
1092  for (Mask::x_iterator ptr = msk.row_begin(y), end = msk.row_end(y); ptr != end; ++ptr) {
1093  sum |= (*ptr)[0];
1094  }
1095  }
1096  _sum = sum;
1097 }
1098 
1099 /*
1100  * Although short, the definition can't be in the header as it must
1101  * follow the specialization definition
1102  * (g++ complained when this was in the header.)
1103  */
1105  StatisticsControl const &sctrl) {
1106  return Statistics(msk, msk, msk, flags, sctrl);
1107 }
1108 
1109 /*
1110  * Explicit instantiations
1111  *
1112  * explicit Statistics(MaskedImage const& img, int const flags,
1113  * StatisticsControl const& sctrl=StatisticsControl());
1114  */
1116 //
1117 #define STAT Statistics
1118 
1119 using VPixel = image::VariancePixel;
1120 
1121 #define INSTANTIATE_MASKEDIMAGE_STATISTICS(TYPE) \
1122  template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1123  image::Image<VPixel> const &var, int const flags, \
1124  StatisticsControl const &sctrl); \
1125  template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1126  image::Image<VPixel> const &var, image::Image<VPixel> const &weights, \
1127  int const flags, StatisticsControl const &sctrl); \
1128  template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1129  image::Image<VPixel> const &var, ImageImposter<VPixel> const &weights, \
1130  int const flags, StatisticsControl const &sctrl)
1131 
1132 #define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(TYPE) \
1133  template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1134  image::Image<VPixel> const &var, int const flags, \
1135  StatisticsControl const &sctrl); \
1136  template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1137  image::Image<VPixel> const &var, image::Image<VPixel> const &weights, \
1138  int const flags, StatisticsControl const &sctrl)
1139 
1140 #define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(TYPE) \
1141  template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1142  MaskImposter<VPixel> const &var, int const flags, \
1143  StatisticsControl const &sctrl); \
1144  template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1145  MaskImposter<VPixel> const &var, image::Image<VPixel> const &weights, \
1146  int const flags, StatisticsControl const &sctrl); \
1147  template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1148  MaskImposter<VPixel> const &var, ImageImposter<VPixel> const &weights, \
1149  int const flags, StatisticsControl const &sctrl)
1150 
1151 #define INSTANTIATE_REGULARIMAGE_STATISTICS(TYPE) \
1152  template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1153  MaskImposter<VPixel> const &var, int const flags, \
1154  StatisticsControl const &sctrl)
1155 
1156 #define INSTANTIATE_VECTOR_STATISTICS(TYPE) \
1157  template STAT::Statistics(ImageImposter<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1158  MaskImposter<VPixel> const &var, int const flags, \
1159  StatisticsControl const &sctrl); \
1160  template STAT::Statistics(ImageImposter<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1161  MaskImposter<VPixel> const &var, ImageImposter<VPixel> const &weights, \
1162  int const flags, StatisticsControl const &sctrl)
1163 
1164 #define INSTANTIATE_IMAGE_STATISTICS(T) \
1165  INSTANTIATE_MASKEDIMAGE_STATISTICS(T); \
1166  INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(T); \
1167  INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(T); \
1168  INSTANTIATE_REGULARIMAGE_STATISTICS(T); \
1169  INSTANTIATE_VECTOR_STATISTICS(T)
1170 
1171 INSTANTIATE_IMAGE_STATISTICS(double);
1172 INSTANTIATE_IMAGE_STATISTICS(float);
1173 INSTANTIATE_IMAGE_STATISTICS(int);
1174 INSTANTIATE_IMAGE_STATISTICS(std::uint16_t);
1175 INSTANTIATE_IMAGE_STATISTICS(std::uint64_t);
1176 
1178 } // namespace math
1179 } // namespace afw
1180 } // namespace lsst
Key< Flag > const & target
int end
table::Key< int > type
Definition: Detector.cc:163
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< MaskPixelT > > mask
afw::table::Key< afw::table::Array< VariancePixelT > > variance
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:48
T begin(T... args)
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:294
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:296
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition: ImageBase.h:401
x_iterator row_end(int y) const
Return an x_iterator to the end of the y'th row.
Definition: ImageBase.h:404
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
Pass parameters to a Statistics object.
Definition: Statistics.h:92
bool getCalcErrorFromInputVariance() const noexcept
Definition: Statistics.h:139
int getAndMask() const noexcept
Definition: Statistics.h:134
bool getWeightedIsSet() const noexcept
Definition: Statistics.h:138
void setWeighted(bool useWeights) noexcept
Definition: Statistics.h:158
double getNumSigmaClip() const noexcept
Definition: Statistics.h:132
bool getWeighted() const noexcept
Definition: Statistics.h:137
int getNumIter() const noexcept
Definition: Statistics.h:133
bool getNanSafe() const noexcept
Definition: Statistics.h:136
A class to evaluate image statistics.
Definition: Statistics.h:220
Value getResult(Property const prop=NOTHING) const
Return the value and error in the specified statistic (e.g.
Definition: Statistics.cc:922
double getError(Property const prop=NOTHING) const
Return the error in the desired property (if specified in the constructor)
Definition: Statistics.cc:1049
Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Constructor for Statistics object.
Definition: Statistics.cc:764
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
std::pair< double, double > Value
The type used to report (value, error) for desired statistics.
Definition: Statistics.h:223
Reports invalid arguments.
Definition: Runtime.h:66
T end(T... args)
T fabs(T... args)
T isfinite(T... args)
T isnan(T... args)
T left(T... args)
T max(T... args)
T min(T... args)
float VariancePixel
default type for MaskedImage variance images
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:359
Property
control what is calculated
Definition: Statistics.h:62
@ ORMASK
get the or-mask of all pixels used.
Definition: Statistics.h:79
@ ERRORS
Include errors of requested quantities.
Definition: Statistics.h:64
@ VARIANCECLIP
estimate sample N-sigma clipped variance (N set in StatisticsControl, default=3)
Definition: Statistics.h:73
@ MEANSQUARE
find mean value of square of pixel values
Definition: Statistics.h:78
@ MIN
estimate sample minimum
Definition: Statistics.h:75
@ NCLIPPED
number of clipped points
Definition: Statistics.h:80
@ NOTHING
We don't want anything.
Definition: Statistics.h:63
@ STDEV
estimate sample standard deviation
Definition: Statistics.h:67
@ NMASKED
number of masked points
Definition: Statistics.h:81
@ STDEVCLIP
estimate sample N-sigma clipped stdev (N set in StatisticsControl, default=3)
Definition: Statistics.h:72
@ VARIANCE
estimate sample variance
Definition: Statistics.h:68
@ MEDIAN
estimate sample median
Definition: Statistics.h:69
@ MAX
estimate sample maximum
Definition: Statistics.h:76
@ SUM
find sum of pixels in the image
Definition: Statistics.h:77
@ IQRANGE
estimate sample inter-quartile range
Definition: Statistics.h:70
@ MEAN
estimate sample mean
Definition: Statistics.h:66
@ MEANCLIP
estimate sample N-sigma clipped mean (N set in StatisticsControl, default=3)
Definition: Statistics.h:71
@ NPOINT
number of sample points
Definition: Statistics.h:65
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:738
constexpr double HALFPI
Definition: Angle.h:41
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
STL namespace.
T nth_element(T... args)
T pow(T... args)
T quiet_NaN(T... args)
T size(T... args)
T sqrt(T... args)
afw::table::Key< double > weight
ImageT val
Definition: CR.cc:146