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