LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
39#include "lsst/geom/Angle.h"
40
41using namespace std;
43
44namespace lsst {
45namespace afw {
46namespace math {
47
48namespace {
50double const MAX_DOUBLE = std::numeric_limits<double>::max();
51double const IQ_TO_STDEV = 0.741301109252802; // 1 sigma in units of iqrange (assume Gaussian)
52
54class AlwaysTrue {
55public:
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
71class AlwaysFalse {
72public:
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
88class CheckFinite {
89public:
90 template <typename T>
91 bool operator()(T val) const {
92 return std::isfinite(static_cast<float>(val));
93 }
94};
95
97class CheckValueLtMin {
98public:
99 template <typename Tval, typename Tmin>
100 bool operator()(Tval val, Tmin min) const {
101 return (static_cast<Tmin>(val) < min);
102 }
103};
104
106class CheckValueGtMax {
107public:
108 template <typename Tval, typename Tmax>
109 bool operator()(Tval val, Tmax max) const {
110 return (static_cast<Tmax>(val) > max);
111 }
112};
113
115class CheckClipRange {
116public:
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
125using ChkFin = CheckFinite;
126using ChkMin = CheckValueLtMin;
127using ChkMax = CheckValueGtMax;
128using ChkClip = CheckClipRange;
129using AlwaysT = AlwaysTrue;
130using AlwaysF = AlwaysFalse;
131
135inline 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 */
155template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
156 bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
157StandardReturn 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 bool const calcErrorMosaicMode,
163 std::vector<double> const &maskPropagationThresholds) {
164 int n = 0;
165 double sumw = 0.0; // sum(weight) (N.b. weight will be 1.0 if !useWeights)
166 double sumw2 = 0.0; // sum(weight^2)
167 double sumx = 0; // sum(data*weight)
168 double sumx2 = 0; // sum(data*weight^2)
169#if 1
170 double sumvw2 = 0.0; // sum(variance*weight^2)
171 double sumvw = 0.0; // sum(variance*weight)
172#endif
173 double min = (nCrude) ? meanCrude : MAX_DOUBLE;
174 double max = (nCrude) ? meanCrude : -MAX_DOUBLE;
175
176 image::MaskPixel allPixelOrMask = 0x0;
177
178 std::vector<double> rejectedWeightsByBit(maskPropagationThresholds.size(), 0.0);
179
180 for (int iY = 0; iY < img.getHeight(); iY += stride) {
181 typename MaskT::x_iterator mptr = msk.row_begin(iY);
182 typename VarianceT::x_iterator vptr = var.row_begin(iY);
183 typename WeightT::x_iterator wptr = weights.row_begin(iY);
184
185 for (typename ImageT::x_iterator ptr = img.row_begin(iY), end = ptr + img.getWidth(); ptr != end;
186 ++ptr, ++mptr, ++vptr, ++wptr) {
187 if (IsFinite()(*ptr) && !(*mptr & andMask) &&
188 InClipRange()(*ptr, meanCrude, cliplimit)) { // clip
189
190 double const delta = (*ptr - meanCrude);
191
192 if (useWeights) {
193 double weight = *wptr;
194 if (weightsAreMultiplicative) {
195 ;
196 } else {
197 if (*wptr <= 0) {
198 continue;
199 }
200 weight = 1 / weight;
201 }
202
203 sumw += weight;
204 sumw2 += weight * weight;
205 sumx += weight * delta;
206 sumx2 += weight * delta * delta;
207
208 if (calcErrorFromInputVariance) {
209 double const var = *vptr;
210 sumvw2 += var * weight * weight;
211 }
212 if (calcErrorMosaicMode) {
213 double const var = *vptr;
214 sumvw += var * weight;
215 }
216 } else {
217 sumx += delta;
218 sumx2 += delta * delta;
219
220 if (calcErrorFromInputVariance) {
221 double const var = *vptr;
222 sumvw2 += var;
223 }
224 if (calcErrorMosaicMode) {
225 double const var = *vptr;
226 sumvw += var;
227 }
228 }
229
230 allPixelOrMask |= *mptr;
231
232 if (HasValueLtMin()(*ptr, min)) {
233 min = *ptr;
234 }
235 if (HasValueGtMax()(*ptr, max)) {
236 max = *ptr;
237 }
238 n++;
239 } else { // pixel has been clipped, rejected, etc.
240 for (int bit = 0, nBits = maskPropagationThresholds.size(); bit < nBits; ++bit) {
241 image::MaskPixel mask = 1 << bit;
242 if (*mptr & mask) {
243 double weight = 1.0;
244 if (useWeights) {
245 weight = *wptr;
246 if (!weightsAreMultiplicative) {
247 if (*wptr <= 0) {
248 continue;
249 }
250 weight = 1.0 / weight;
251 }
252 }
253 rejectedWeightsByBit[bit] += weight;
254 }
255 }
256 }
257 }
258 }
259 if (n == 0) {
260 min = NaN;
261 max = NaN;
262 }
263
264 // estimate of population mean and variance.
265 double mean, variance;
266 if (!useWeights) {
267 sumw = sumw2 = n;
268 }
269
270 for (int bit = 0, nBits = maskPropagationThresholds.size(); bit < nBits; ++bit) {
271 double hypotheticalTotalWeight = sumw + rejectedWeightsByBit[bit];
272 rejectedWeightsByBit[bit] /= hypotheticalTotalWeight;
273 if (rejectedWeightsByBit[bit] > maskPropagationThresholds[bit]) {
274 allPixelOrMask |= (1 << bit);
275 }
276 }
277
278 // N.b. if sumw == 0 or sumw*sumw == sumw2 (e.g. n == 1) we'll get NaNs
279 // N.b. the estimator of the variance assumes that the sample points all have the same variance;
280 // otherwise, what is it that we're estimating?
281 mean = sumx / sumw;
282 variance = sumx2 / sumw - ::pow(mean, 2); // biased estimator
283 variance *= sumw * sumw / (sumw * sumw - sumw2); // debias
284
285 double meanVar; // (standard error of mean)^2
286 if (calcErrorFromInputVariance) {
287 meanVar = sumvw2 / (sumw * sumw);
288 } else if (calcErrorMosaicMode){
289 meanVar = sumvw / sumw;
290 } else {
291 meanVar = variance * sumw2 / (sumw * sumw);
292 }
293
294 double varVar = varianceError(variance, n); // error in variance; incorrect if useWeights is true
295
296 sumx += sumw * meanCrude;
297 mean += meanCrude;
298
299 return StandardReturn(n, sumx, Statistics::Value(mean, meanVar), Statistics::Value(variance, varVar), min,
300 max, allPixelOrMask);
301}
302
303template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
304 bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
305StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
306 WeightT const &weights, int const flags, int const nCrude, int const stride,
307 double const meanCrude, double const cliplimit,
308 bool const weightsAreMultiplicative, int const andMask,
309 bool const calcErrorFromInputVariance, bool const calcErrorMosaicMode,
310 bool doGetWeighted, std::vector<double> const &maskPropagationThresholds) {
311 if (doGetWeighted) {
312 return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, true>(
313 img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
314 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, maskPropagationThresholds);
315 } else {
316 return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, false>(
317 img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
318 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, maskPropagationThresholds);
319 }
320}
321
322template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
323 bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
324StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
325 WeightT const &weights, int const flags, int const nCrude, int const stride,
326 double const meanCrude, double const cliplimit,
327 bool const weightsAreMultiplicative, int const andMask,
328 bool const calcErrorFromInputVariance,
329 bool const calcErrorMosaicMode, bool doCheckFinite, bool doGetWeighted,
330 std::vector<double> const &maskPropagationThresholds) {
331 if (doCheckFinite) {
332 return processPixels<CheckFinite, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
333 img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
334 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, doGetWeighted, maskPropagationThresholds);
335 } else {
336 return processPixels<AlwaysTrue, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
337 img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
338 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, doGetWeighted, maskPropagationThresholds);
339 }
340}
341
362template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
363StandardReturn getStandard(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
364 int const flags, bool const weightsAreMultiplicative, int const andMask,
365 bool const calcErrorFromInputVariance, bool const calcErrorMosaicMode,
366 bool doCheckFinite, bool doGetWeighted,
367 std::vector<double> const &maskPropagationThresholds) {
368 // =====================================================
369 // a crude estimate of the mean, used for numerical stability of variance
370 int nCrude = 0;
371 double meanCrude = 0.0;
372
373 // for small numbers of values, use a small stride
374 int const nPix = img.getWidth() * img.getHeight();
375 int strideCrude;
376 if (nPix < 100) {
377 strideCrude = 2;
378 } else {
379 strideCrude = 10;
380 }
381
382 double cliplimit = -1; // unused
383 StandardReturn values = processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
384 img, msk, var, weights, flags, nCrude, strideCrude, meanCrude, cliplimit,
385 weightsAreMultiplicative, andMask, calcErrorFromInputVariance, calcErrorMosaicMode,
386 doCheckFinite, doGetWeighted, maskPropagationThresholds);
387 nCrude = std::get<0>(values);
388 double sumCrude = std::get<1>(values);
389
390 meanCrude = 0.0;
391 if (nCrude > 0) {
392 meanCrude = sumCrude / nCrude;
393 }
394
395 // =======================================================
396 // Estimate the full precision variance using that crude mean
397 // - get the min and max as well
398
399 if (flags & (MIN | MAX)) {
400 return processPixels<ChkFin, ChkMin, ChkMax, AlwaysT, true>(
401 img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
402 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, true, doGetWeighted,
403 maskPropagationThresholds);
404 } else {
405 return processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
406 img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
407 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, doCheckFinite, doGetWeighted,
408 maskPropagationThresholds);
409 }
410}
411
433template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
434StandardReturn getStandard(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
435 int const flags, std::pair<double, double> const clipinfo,
436 bool const weightsAreMultiplicative, int const andMask,
437 bool const calcErrorFromInputVariance, bool const calcErrorMosaicMode,
438 bool doCheckFinite, bool doGetWeighted,
439 std::vector<double> const &maskPropagationThresholds) {
440 double const center = clipinfo.first;
441 double const cliplimit = clipinfo.second;
442
443 if (std::isnan(center) || std::isnan(cliplimit)) {
444 return StandardReturn(0, NaN, Statistics::Value(NaN, NaN), Statistics::Value(NaN, NaN), NaN, NaN,
445 ~0x0);
446 }
447
448 // =======================================================
449 // Estimate the full precision variance using that crude mean
450 int const stride = 1;
451 int nCrude = 0;
452
453 if (flags & (MIN | MAX)) {
454 return processPixels<ChkFin, ChkMin, ChkMax, ChkClip, true>(
455 img, msk, var, weights, flags, nCrude, stride, center, cliplimit, weightsAreMultiplicative,
456 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, true, doGetWeighted,
457 maskPropagationThresholds);
458 } else { // fast loop ... just the mean & variance
459 return processPixels<ChkFin, AlwaysF, AlwaysF, ChkClip, true>(
460 img, msk, var, weights, flags, nCrude, stride, center, cliplimit, weightsAreMultiplicative,
461 andMask, calcErrorFromInputVariance, calcErrorMosaicMode, doCheckFinite, doGetWeighted,
462 maskPropagationThresholds);
463 }
464}
465
474template <typename Pixel>
475typename enable_if<!is_integral<Pixel>::value, double>::type
476percentile(std::vector<Pixel> &img, double const fraction) {
477 assert(fraction >= 0.0 && fraction <= 1.0);
478
479 int const n = img.size();
480
481 if (n > 1) {
482 double const idx = fraction * (n - 1);
483
484 // interpolate linearly between the adjacent values
485 // For efficiency:
486 // - if we're asked for a fraction > 0.5,
487 // we'll do the second partial sort on shorter (upper) portion
488 // - otherwise, the shorter portion will be the lower one, we'll partial-sort that.
489
490 int const q1 = static_cast<int>(idx);
491 int const q2 = q1 + 1;
492
493 auto mid1 = img.begin() + q1;
494 auto mid2 = img.begin() + q2;
495 if (fraction > 0.5) {
496 std::nth_element(img.begin(), mid1, img.end());
497 std::nth_element(mid1, mid2, img.end());
498 } else {
499 std::nth_element(img.begin(), mid2, img.end());
500 std::nth_element(img.begin(), mid1, mid2);
501 }
502
503 double val1 = static_cast<double>(*mid1);
504 double val2 = static_cast<double>(*mid2);
505 double w1 = (static_cast<double>(q2) - idx);
506 double w2 = (idx - static_cast<double>(q1));
507 return w1 * val1 + w2 * val2;
508
509 } else if (n == 1) {
510 return img[0];
511 } else {
512 return NaN;
513 }
514}
515
516namespace {
517 //
518 // Helper function to estimate a floating-point quantile from integer data
519 //
520 // The data has been partially sorted, using nth_element
521 //
522 // begin: iterator to a point which we know to be below our median
523 // end: iterator to a point which we know to be beyond our median
524 // naive: the integer value of the desired quantile
525 // target: the number of points that should be to the left of the quantile.
526 // N.b. if begin isn't the start of the data, this may not be the
527 // desired number of points. Caveat Callor
528 template<typename T>
529 double
530 computeQuantile(typename std::vector<T>::const_iterator begin,
531 typename std::vector<T>::const_iterator end,
532 T const naive,
533 double const target
534 )
535 {
536 // investigate the cumulative histogram near naive
537 std::size_t left = 0; // number of values less than naive
538 std::size_t middle = 0; // number of values equal to naive
539
540 for (auto ptr = begin; ptr != end; ++ptr) {
541 auto const val = *ptr;
542 if (val < naive) {
543 ++left;
544 } else if (val == naive) {
545 ++middle;
546 }
547 }
548
549 return naive - 0.5 + (target - left)/middle;
550 }
551
560 template <typename Pixel>
561 typename enable_if<is_integral<Pixel>::value, double>::type
562 percentile(std::vector<Pixel> &img, double const fraction) {
563 assert(fraction >= 0.0 && fraction <= 1.0);
564
565 auto const n = img.size();
566
567 if (n == 0) {
568 return NaN;
569 } else if (n == 1) {
570 return img[0];
571 } else {
572 // We need to handle ties. The proper way to do this is to analyse the cumulative curve after
573 // building the histograms (which is faster than a generic partitioning algorithm), but it's a
574 // nuisance as we don't know the range of values
575 //
576 // This code looks clean enough, but actually the call to nth_element is expensive
577 // and we *still* have to go through the array a second time
578
579 double const idx = fraction*(n - 1);
580
581 auto midP = img.begin() + static_cast<int>(idx);
582 std::nth_element(img.begin(), midP, img.end());
583 auto const naiveP = *midP; // value of desired element
584
585 return computeQuantile(img.begin(), img.end(), naiveP, fraction*n);
586 }
587 }
588}
589
590using MedianQuartileReturn = std::tuple<double, double, double>;
591namespace {
599 template <typename Pixel>
600 typename enable_if<!is_integral<Pixel>::value, MedianQuartileReturn>::type
601 medianAndQuartiles(std::vector<Pixel> &img) {
602 int const n = img.size();
603
604 if (n > 1) {
605 double const idx50 = 0.50 * (n - 1);
606 double const idx25 = 0.25 * (n - 1);
607 double const idx75 = 0.75 * (n - 1);
608
609 // For efficiency:
610 // - partition at 50th, then partition the two half further to get 25th and 75th
611 // - to get the adjacent points (for interpolation), partition between 25/50, 50/75, 75/end
612 // these should be much smaller partitions
613
614 int const q50a = static_cast<int>(idx50);
615 int const q50b = q50a + 1;
616 int const q25a = static_cast<int>(idx25);
617 int const q25b = q25a + 1;
618 int const q75a = static_cast<int>(idx75);
619 int const q75b = q75a + 1;
620
621 auto mid50a = img.begin() + q50a;
622 auto mid50b = img.begin() + q50b;
623 auto mid25a = img.begin() + q25a;
624 auto mid25b = img.begin() + q25b;
625 auto mid75a = img.begin() + q75a;
626 auto mid75b = img.begin() + q75b;
627
628 // get the 50th percentile, then get the 25th and 75th on the smaller partitions
629 std::nth_element(img.begin(), mid50a, img.end());
630 std::nth_element(mid50a, mid75a, img.end());
631 std::nth_element(img.begin(), mid25a, mid50a);
632
633 // and the adjacent points for each ... use the smallest segments available.
634 std::nth_element(mid50a, mid50b, mid75a);
635 std::nth_element(mid25a, mid25b, mid50a);
636 std::nth_element(mid75a, mid75b, img.end());
637
638 // interpolate linearly between the adjacent values
639 double val50a = static_cast<double>(*mid50a);
640 double val50b = static_cast<double>(*mid50b);
641 double w50a = (static_cast<double>(q50b) - idx50);
642 double w50b = (idx50 - static_cast<double>(q50a));
643 double median = w50a * val50a + w50b * val50b;
644
645 double val25a = static_cast<double>(*mid25a);
646 double val25b = static_cast<double>(*mid25b);
647 double w25a = (static_cast<double>(q25b) - idx25);
648 double w25b = (idx25 - static_cast<double>(q25a));
649 double q1 = w25a * val25a + w25b * val25b;
650
651 double val75a = static_cast<double>(*mid75a);
652 double val75b = static_cast<double>(*mid75b);
653 double w75a = (static_cast<double>(q75b) - idx75);
654 double w75b = (idx75 - static_cast<double>(q75a));
655 double q3 = w75a * val75a + w75b * val75b;
656
657 return MedianQuartileReturn(median, q1, q3);
658 } else if (n == 1) {
659 return MedianQuartileReturn(img[0], img[0], img[0]);
660 } else {
661 return MedianQuartileReturn(NaN, NaN, NaN);
662 }
663 }
664
672 template <typename Pixel>
673 typename enable_if<is_integral<Pixel>::value, MedianQuartileReturn>::type
674 medianAndQuartiles(std::vector<Pixel> &img) {
675 auto const n = img.size();
676
677 if (n == 0) {
678 return MedianQuartileReturn(NaN, NaN, NaN);
679 } else if (n == 1) {
680 return MedianQuartileReturn(img[0], img[0], img[0]);
681 } else {
682 // We need to handle ties. The proper way to do this is to analyse the cumulative curve after
683 // building the histograms (which is faster than a generic partitioning algorithm), but it's a
684 // nuisance as we don't know the range of values
685 //
686 // This code looks clean enough, but actually the call to nth_element is expensive
687 // and we *still* have to go through the array a second time.
688
689 // For efficiency:
690 // - partition at 50th, then partition the two halves further to get 25th and 75th
691
692 auto mid25 = img.begin() + static_cast<int>(0.25*(n - 1));
693 auto mid50 = img.begin() + static_cast<int>(0.50*(n - 1));
694 auto mid75 = img.begin() + static_cast<int>(0.75*(n - 1));
695
696 // get the 50th percentile, then get the 25th and 75th on the smaller partitions
697 std::nth_element(img.begin(), mid50, img.end());
698 std::nth_element(img.begin(), mid25, mid50);
699 std::nth_element(mid50, mid75, img.end());
700
701 double const q1 = computeQuantile(img.begin(), mid50, *mid25,
702 0.25*n);
703 double const median = computeQuantile(mid25, mid75, *mid50,
704 0.50*n - (mid25 - img.begin()));
705 double const q3 = computeQuantile(mid50, img.end(), *mid75,
706 0.75*n - (mid50 - img.begin()));
707
708 return MedianQuartileReturn(median, q1, q3);
709 }
710 }
711}
712
720template <typename IsFinite, typename ImageT, typename MaskT, typename VarianceT>
721std::shared_ptr<std::vector<typename ImageT::Pixel> > makeVectorCopy(ImageT const &img, MaskT const &msk,
722 VarianceT const &, int const andMask) {
723 // Note need to keep track of allPixelOrMask here ... processPixels() does that
724 // and it always gets called
726
727 for (int i_y = 0; i_y < img.getHeight(); ++i_y) {
728 typename MaskT::x_iterator mptr = msk.row_begin(i_y);
729 for (typename ImageT::x_iterator ptr = img.row_begin(i_y), end = img.row_end(i_y); ptr != end;
730 ++ptr) {
731 if (IsFinite()(*ptr) && !(*mptr & andMask)) {
732 imgcp->push_back(*ptr);
733 }
734 ++mptr;
735 }
736 }
737
738 return imgcp;
739}
740} // namespace
741
743 int oldSize = _maskPropagationThresholds.size();
744 if (oldSize <= bit) {
745 return 1.0;
746 }
747 return _maskPropagationThresholds[bit];
748}
749
750void StatisticsControl::setMaskPropagationThreshold(int bit, double threshold) {
751 int oldSize = _maskPropagationThresholds.size();
752 if (oldSize <= bit) {
753 int newSize = bit + 1;
754 _maskPropagationThresholds.resize(newSize);
755 for (int i = oldSize; i < bit; ++i) {
756 _maskPropagationThresholds[i] = 1.0;
757 }
758 }
759 _maskPropagationThresholds[bit] = threshold;
760}
761
763 static std::map<std::string, Property> statisticsProperty;
764 if (statisticsProperty.size() == 0) {
765 statisticsProperty["NOTHING"] = NOTHING;
766 statisticsProperty["ERRORS"] = ERRORS;
767 statisticsProperty["NPOINT"] = NPOINT;
768 statisticsProperty["MEAN"] = MEAN;
769 statisticsProperty["STDEV"] = STDEV;
770 statisticsProperty["VARIANCE"] = VARIANCE;
771 statisticsProperty["MEDIAN"] = MEDIAN;
772 statisticsProperty["IQRANGE"] = IQRANGE;
773 statisticsProperty["MEANCLIP"] = MEANCLIP;
774 statisticsProperty["STDEVCLIP"] = STDEVCLIP;
775 statisticsProperty["VARIANCECLIP"] = VARIANCECLIP;
776 statisticsProperty["MIN"] = MIN;
777 statisticsProperty["MAX"] = MAX;
778 statisticsProperty["SUM"] = SUM;
779 statisticsProperty["MEANSQUARE"] = MEANSQUARE;
780 statisticsProperty["ORMASK"] = ORMASK;
781 statisticsProperty["NCLIPPED"] = NCLIPPED;
782 statisticsProperty["NMASKED"] = NMASKED;
783 }
784 return statisticsProperty[property];
785}
786
787template <typename ImageT, typename MaskT, typename VarianceT>
788Statistics::Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags,
789 StatisticsControl const &sctrl)
790 : _flags(flags),
791 _mean(NaN, NaN),
792 _variance(NaN, NaN),
793 _min(NaN),
794 _max(NaN),
795 _sum(NaN),
796 _meanclip(NaN, NaN),
797 _varianceclip(NaN, NaN),
798 _median(NaN, NaN),
799 _nClipped(0),
800 _nMasked(0),
801 _iqrange(NaN),
802 _sctrl(sctrl),
803 _weightsAreMultiplicative(false) {
804 doStatistics(img, msk, var, var, _flags, _sctrl);
805}
806
807namespace {
808template <typename T>
809bool isEmpty(T const &t) {
810 return t.empty();
811}
812
813template <typename T>
814bool isEmpty(image::Image<T> const &im) {
815 return (im.getWidth() == 0 && im.getHeight() == 0);
816}
817
818// Asserts that image dimensions are equal
819template <typename ImageT1, typename ImageT2>
820void checkDimensions(ImageT1 const &image1, ImageT2 const &image2) {
821 if (image1.getDimensions() != image2.getDimensions()) {
823 (boost::format("Image sizes don't match: %s vs %s") % image1.getDimensions() %
824 image2.getDimensions())
825 .str());
826 }
827}
828
829// Overloads for MaskImposter (which doesn't have a size)
830template <typename ImageT, typename PixelT>
831void checkDimensions(ImageT const &image1, MaskImposter<PixelT> const &image2) {}
832template <typename ImageT, typename PixelT>
833void checkDimensions(MaskImposter<PixelT> const &image1, ImageT const &image2) {}
834} // namespace
835
836template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
837Statistics::Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
838 int const flags, StatisticsControl const &sctrl)
839 : _flags(flags),
840 _mean(NaN, NaN),
841 _variance(NaN, NaN),
842 _min(NaN),
843 _max(NaN),
844 _sum(NaN),
845 _meanclip(NaN, NaN),
846 _varianceclip(NaN, NaN),
847 _median(NaN, NaN),
848 _nClipped(0),
849 _nMasked(0),
850 _iqrange(NaN),
851 _sctrl(sctrl),
852 _weightsAreMultiplicative(true) {
853 if (!isEmpty(weights)) {
854 if (_sctrl.getWeightedIsSet() && !_sctrl.getWeighted()) {
856 "You must use the weights if you provide them");
857 }
858
859 _sctrl.setWeighted(true);
860 }
861 doStatistics(img, msk, var, weights, _flags, _sctrl);
862}
863
864template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
865void Statistics::doStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var,
866 WeightT const &weights, int const flags, StatisticsControl const &sctrl) {
867
868 if (_sctrl.getCalcErrorFromInputVariance() && _sctrl.getCalcErrorMosaicMode()) {
870 "Both calcErrorFromInputVariance and calcErrorMosaicMode are True");
871 }
872
873 int const num = img.getWidth() * img.getHeight();
874 _n = num;
875 if (_n == 0) {
876 throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
877 }
878 checkDimensions(img, msk);
879 checkDimensions(img, var);
880 if (sctrl.getWeighted()) {
881 checkDimensions(img, weights);
882 }
883
884 // Check that an int's large enough to hold the number of pixels
885 assert(img.getWidth() * static_cast<double>(img.getHeight()) < std::numeric_limits<int>::max());
886
887 // get the standard statistics
888 StandardReturn standard =
889 getStandard(img, msk, var, weights, flags, _weightsAreMultiplicative, _sctrl.getAndMask(),
891 _sctrl.getNanSafe(), _sctrl.getWeighted(),
892 _sctrl._maskPropagationThresholds);
893
894 _n = std::get<0>(standard);
895 _sum = std::get<1>(standard);
896 _mean = std::get<2>(standard);
897 _variance = std::get<3>(standard);
898 _min = std::get<4>(standard);
899 _max = std::get<5>(standard);
900 _allPixelOrMask = std::get<6>(standard);
901
902 // ==========================================================
903 // now only calculate it if it's specifically requested - these all cost more!
904
905 if (flags & NMASKED) {
906 _nMasked = num - _n;
907 }
908
909 // copy the image for any routines that will use median or quantiles
910 if (flags & (MEDIAN | IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
911 // make a vector copy of the image to get the median and quartiles (will move values)
913 if (_sctrl.getNanSafe()) {
914 imgcp = makeVectorCopy<ChkFin>(img, msk, var, _sctrl.getAndMask());
915 } else {
916 imgcp = makeVectorCopy<AlwaysT>(img, msk, var, _sctrl.getAndMask());
917 }
918
919 // if we *only* want the median, just use percentile(), otherwise use medianAndQuartiles()
920 if ((flags & (MEDIAN)) && !(flags & (IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP))) {
921 _median = Value(percentile(*imgcp, 0.5), NaN);
922 } else {
923 MedianQuartileReturn mq = medianAndQuartiles(*imgcp);
924 _median = Value(std::get<0>(mq), NaN);
925 _iqrange = std::get<2>(mq) - std::get<1>(mq);
926 }
927
928 if (flags & (MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
929 for (int i_i = 0; i_i < _sctrl.getNumIter(); ++i_i) {
930 double const center = ((i_i > 0) ? _meanclip : _median).first;
931 double const hwidth = (i_i > 0 && _n > 1)
932 ? _sctrl.getNumSigmaClip() * std::sqrt(_varianceclip.first)
933 : _sctrl.getNumSigmaClip() * IQ_TO_STDEV * _iqrange;
934 std::pair<double, double> const clipinfo(center, hwidth);
935
936 StandardReturn clipped = getStandard(
937 img, msk, var, weights, flags, clipinfo, _weightsAreMultiplicative,
938 _sctrl.getAndMask(), _sctrl.getCalcErrorFromInputVariance(),
939 _sctrl.getCalcErrorMosaicMode(), _sctrl.getNanSafe(),
940 _sctrl.getWeighted(), _sctrl._maskPropagationThresholds);
941
942 int const nClip = std::get<0>(clipped); // number after clipping
943 _nClipped = _n - nClip; // number clipped
944 _meanclip = std::get<2>(clipped); // clipped mean
945 double const varClip = std::get<3>(clipped).first; // clipped variance
946
947 _varianceclip = Value(varClip, varianceError(varClip, nClip));
948 // ... ignore other values
949 }
950 }
951 }
952}
953
955 // if iProp == NOTHING try to return their heart's delight, as specified in the constructor
956 Property const prop = static_cast<Property>(((iProp == NOTHING) ? _flags : iProp) & ~ERRORS);
957
958 if (!(prop & _flags)) { // we didn't calculate it
960 (boost::format("You didn't ask me to calculate %d") % prop).str());
961 }
962
963 Value ret(NaN, NaN);
964 switch (prop) {
965 case NPOINT:
966 ret.first = static_cast<double>(_n);
967 if (_flags & ERRORS) {
968 ret.second = 0;
969 }
970 break;
971
972 case NCLIPPED:
973 ret.first = static_cast<double>(_nClipped);
974 if (_flags & ERRORS) {
975 ret.second = 0;
976 }
977 break;
978
979 case NMASKED:
980 ret.first = static_cast<double>(_nMasked);
981 if (_flags & ERRORS) {
982 ret.second = 0;
983 }
984 break;
985
986 case SUM:
987 ret.first = static_cast<double>(_sum);
988 if (_flags & ERRORS) {
989 ret.second = 0;
990 }
991 break;
992
993 // == means ==
994 case MEAN:
995 ret.first = _mean.first;
996 if (_flags & ERRORS) {
997 ret.second = ::sqrt(_mean.second);
998 }
999 break;
1000 case MEANCLIP:
1001 ret.first = _meanclip.first;
1002 if (_flags & ERRORS) {
1003 ret.second = ::sqrt(_meanclip.second);
1004 }
1005 break;
1006
1007 // == stdevs & variances ==
1008 case VARIANCE:
1009 ret.first = _variance.first;
1010 if (_flags & ERRORS) {
1011 ret.second = ::sqrt(_variance.second);
1012 }
1013 break;
1014 case STDEV:
1015 ret.first = sqrt(_variance.first);
1016 if (_flags & ERRORS) {
1017 ret.second = 0.5 * ::sqrt(_variance.second) / ret.first;
1018 }
1019 break;
1020 case VARIANCECLIP:
1021 ret.first = _varianceclip.first;
1022 if (_flags & ERRORS) {
1023 ret.second = ret.second;
1024 }
1025 break;
1026 case STDEVCLIP:
1027 ret.first = sqrt(_varianceclip.first);
1028 if (_flags & ERRORS) {
1029 ret.second = 0.5 * ::sqrt(_varianceclip.second) / ret.first;
1030 }
1031 break;
1032
1033 case MEANSQUARE:
1034 ret.first = (_n - 1) / static_cast<double>(_n) * _variance.first + ::pow(_mean.first, 2);
1035 if (_flags & ERRORS) {
1036 ret.second = ::sqrt(2 * ::pow(ret.first / _n, 2)); // assumes Gaussian
1037 }
1038 break;
1039
1040 // == other stats ==
1041 case MIN:
1042 ret.first = _min;
1043 if (_flags & ERRORS) {
1044 ret.second = 0;
1045 }
1046 break;
1047 case MAX:
1048 ret.first = _max;
1049 if (_flags & ERRORS) {
1050 ret.second = 0;
1051 }
1052 break;
1053 case MEDIAN:
1054 ret.first = _median.first;
1055 if (_flags & ERRORS) {
1056 ret.second = sqrt(geom::HALFPI * _variance.first / _n); // assumes Gaussian
1057 }
1058 break;
1059 case IQRANGE:
1060 ret.first = _iqrange;
1061 if (_flags & ERRORS) {
1062 ret.second = NaN; // we're not estimating this properly
1063 }
1064 break;
1065
1066 // no-op to satisfy the compiler
1067 case ERRORS:
1068 break;
1069 // default: redundant as 'ret' is initialized to NaN, NaN
1070 default: // we must have set prop to _flags
1071 assert(iProp == 0);
1073 "getValue() may only be called without a parameter"
1074 " if you asked for only one statistic");
1075 }
1076 return ret;
1077}
1078
1079double Statistics::getValue(Property const prop) const { return getResult(prop).first; }
1080
1081double Statistics::getError(Property const prop) const { return getResult(prop).second; }
1082
1092template <>
1094 image::Mask<image::MaskPixel> const &var, int const flags,
1095 StatisticsControl const &sctrl)
1096 : _flags(flags),
1097 _mean(NaN, NaN),
1098 _variance(NaN, NaN),
1099 _min(NaN),
1100 _max(NaN),
1101 _meanclip(NaN, NaN),
1102 _varianceclip(NaN, NaN),
1103 _median(NaN, NaN),
1104 _nClipped(0),
1105 _iqrange(NaN),
1106 _sctrl(sctrl) {
1107 if ((flags & ~(NPOINT | SUM)) != 0x0) {
1109 "Statistics<Mask> only supports NPOINT and SUM");
1110 }
1111
1112 using Mask = image::Mask<>;
1113
1114 _n = msk.getWidth() * msk.getHeight();
1115 if (_n == 0) {
1116 throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
1117 }
1118
1119 // Check that an int's large enough to hold the number of pixels
1120 assert(msk.getWidth() * static_cast<double>(msk.getHeight()) < std::numeric_limits<int>::max());
1121
1122 image::MaskPixel sum = 0x0;
1123 for (int y = 0; y != msk.getHeight(); ++y) {
1124 for (Mask::x_iterator ptr = msk.row_begin(y), end = msk.row_end(y); ptr != end; ++ptr) {
1125 sum |= (*ptr)[0];
1126 }
1127 }
1128 _sum = sum;
1129}
1130
1131/*
1132 * Although short, the definition can't be in the header as it must
1133 * follow the specialization definition
1134 * (g++ complained when this was in the header.)
1135 */
1137 StatisticsControl const &sctrl) {
1138 return Statistics(msk, msk, msk, flags, sctrl);
1139}
1140
1141/*
1142 * Explicit instantiations
1143 *
1144 * explicit Statistics(MaskedImage const& img, int const flags,
1145 * StatisticsControl const& sctrl=StatisticsControl());
1146 */
1148//
1149#define STAT Statistics
1150
1151using VPixel = image::VariancePixel;
1152
1153#define INSTANTIATE_MASKEDIMAGE_STATISTICS(TYPE) \
1154 template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1155 image::Image<VPixel> const &var, int const flags, \
1156 StatisticsControl const &sctrl); \
1157 template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1158 image::Image<VPixel> const &var, image::Image<VPixel> const &weights, \
1159 int const flags, StatisticsControl const &sctrl); \
1160 template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1161 image::Image<VPixel> const &var, ImageImposter<VPixel> const &weights, \
1162 int const flags, StatisticsControl const &sctrl)
1163
1164#define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(TYPE) \
1165 template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1166 image::Image<VPixel> const &var, int const flags, \
1167 StatisticsControl const &sctrl); \
1168 template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1169 image::Image<VPixel> const &var, image::Image<VPixel> const &weights, \
1170 int const flags, StatisticsControl const &sctrl)
1171
1172#define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(TYPE) \
1173 template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1174 MaskImposter<VPixel> const &var, int const flags, \
1175 StatisticsControl const &sctrl); \
1176 template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1177 MaskImposter<VPixel> const &var, image::Image<VPixel> const &weights, \
1178 int const flags, StatisticsControl const &sctrl); \
1179 template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
1180 MaskImposter<VPixel> const &var, ImageImposter<VPixel> const &weights, \
1181 int const flags, StatisticsControl const &sctrl)
1182
1183#define INSTANTIATE_REGULARIMAGE_STATISTICS(TYPE) \
1184 template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1185 MaskImposter<VPixel> const &var, int const flags, \
1186 StatisticsControl const &sctrl)
1187
1188#define INSTANTIATE_VECTOR_STATISTICS(TYPE) \
1189 template STAT::Statistics(ImageImposter<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1190 MaskImposter<VPixel> const &var, int const flags, \
1191 StatisticsControl const &sctrl); \
1192 template STAT::Statistics(ImageImposter<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
1193 MaskImposter<VPixel> const &var, ImageImposter<VPixel> const &weights, \
1194 int const flags, StatisticsControl const &sctrl)
1195
1196#define INSTANTIATE_IMAGE_STATISTICS(T) \
1197 INSTANTIATE_MASKEDIMAGE_STATISTICS(T); \
1198 INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(T); \
1199 INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(T); \
1200 INSTANTIATE_REGULARIMAGE_STATISTICS(T); \
1201 INSTANTIATE_VECTOR_STATISTICS(T)
1202
1203INSTANTIATE_IMAGE_STATISTICS(double);
1204INSTANTIATE_IMAGE_STATISTICS(float);
1205INSTANTIATE_IMAGE_STATISTICS(int);
1206INSTANTIATE_IMAGE_STATISTICS(std::uint16_t);
1207INSTANTIATE_IMAGE_STATISTICS(std::uint64_t);
1208
1210} // namespace math
1211} // namespace afw
1212} // namespace lsst
Key< Flag > const & target
int min
int end
int max
#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:95
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:83
double getMaskPropagationThreshold(int bit) const
When pixels with the given bit are rejected, we count what fraction the rejected pixels would have co...
bool getCalcErrorFromInputVariance() const noexcept
Definition Statistics.h:131
int getAndMask() const noexcept
Definition Statistics.h:126
bool getWeightedIsSet() const noexcept
Definition Statistics.h:130
void setMaskPropagationThreshold(int bit, double threshold)
void setWeighted(bool useWeights) noexcept
Definition Statistics.h:151
double getNumSigmaClip() const noexcept
Definition Statistics.h:124
bool getWeighted() const noexcept
Definition Statistics.h:129
int getNumIter() const noexcept
Definition Statistics.h:125
bool getCalcErrorMosaicMode() const noexcept
Definition Statistics.h:132
bool getNanSafe() const noexcept
Definition Statistics.h:128
A class to evaluate image statistics.
Definition Statistics.h:222
Value getResult(Property const prop=NOTHING) const
Return the value and error in the specified statistic (e.g.
double getError(Property const prop=NOTHING) const
Return the error in the desired property (if specified in the constructor)
Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Constructor for Statistics object.
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
std::pair< double, double > Value
The type used to report (value, error) for desired statistics.
Definition Statistics.h:225
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)
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:361
Property
control what is calculated
Definition Statistics.h:53
@ ORMASK
get the or-mask of all pixels used.
Definition Statistics.h:70
@ ERRORS
Include errors of requested quantities.
Definition Statistics.h:55
@ VARIANCECLIP
estimate sample N-sigma clipped variance (N set in StatisticsControl, default=3)
Definition Statistics.h:64
@ MEANSQUARE
find mean value of square of pixel values
Definition Statistics.h:69
@ MIN
estimate sample minimum
Definition Statistics.h:66
@ NCLIPPED
number of clipped points
Definition Statistics.h:71
@ NOTHING
We don't want anything.
Definition Statistics.h:54
@ STDEV
estimate sample standard deviation
Definition Statistics.h:58
@ NMASKED
number of masked points
Definition Statistics.h:72
@ STDEVCLIP
estimate sample N-sigma clipped stdev (N set in StatisticsControl, default=3)
Definition Statistics.h:63
@ VARIANCE
estimate sample variance
Definition Statistics.h:59
@ MEDIAN
estimate sample median
Definition Statistics.h:60
@ MAX
estimate sample maximum
Definition Statistics.h:67
@ SUM
find sum of pixels in the image
Definition Statistics.h:68
@ IQRANGE
estimate sample inter-quartile range
Definition Statistics.h:61
@ MEAN
estimate sample mean
Definition Statistics.h:57
@ MEANCLIP
estimate sample N-sigma clipped mean (N set in StatisticsControl, default=3)
Definition Statistics.h:62
@ NPOINT
number of sample points
Definition Statistics.h:56
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
double constexpr HALFPI
Definition Angle.h:42
STL namespace.
T nth_element(T... args)
T pow(T... args)
T quiet_NaN(T... args)
T resize(T... args)
T size(T... args)
T sqrt(T... args)
afw::table::Key< double > weight
ImageT val
Definition CR.cc:146