LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Statistics.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 
33 #include <iostream>
34 #include <limits>
35 #include <cmath>
36 #include <cassert>
37 #include "boost/shared_ptr.hpp"
38 #include "lsst/pex/exceptions.h"
39 #include "lsst/afw/image/Image.h"
41 #include "lsst/utils/ieee.h"
42 #include "lsst/afw/geom/Angle.h"
43 
44 using namespace std;
45 namespace afwImage = lsst::afw::image;
46 namespace afwMath = lsst::afw::math;
47 namespace afwGeom = lsst::afw::geom;
48 namespace pexExceptions = lsst::pex::exceptions;
49 
50 namespace {
51  double const NaN = std::numeric_limits<double>::quiet_NaN();
52  double const MAX_DOUBLE = std::numeric_limits<double>::max();
53  double const IQ_TO_STDEV = 0.741301109252802; // 1 sigma in units of iqrange (assume Gaussian)
54 
58  class AlwaysTrue {
59  public:
60  template<typename T>
61  bool operator()(T) const {
62  return true;
63  }
64  template<typename Ta, typename Tb>
65  bool operator()(Ta, Tb) const {
66  return true;
67  }
68  template<typename Ta, typename Tb, typename Tc>
69  bool operator()(Ta, Tb, Tc) const {
70  return true;
71  }
72  };
73 
77  class AlwaysFalse {
78  public:
79  template<typename T>
80  bool operator()(T) const {
81  return false;
82  }
83  template<typename Ta, typename Tb>
84  bool operator()(Ta, Tb) const {
85  return false;
86  }
87  template<typename Ta, typename Tb, typename Tc>
88  bool operator()(Ta, Tb, Tc) const {
89  return false;
90  }
91  };
92 
96  class CheckFinite {
97  public:
98  template<typename T>
99  bool operator()(T val) const {
100  return lsst::utils::isfinite(static_cast<float>(val));
101  }
102  };
103 
107  class CheckValueLtMin {
108  public:
109  template<typename Tval, typename Tmin>
110  bool operator()(Tval val, Tmin min) const {
111  return (static_cast<Tmin>(val) < min);
112  }
113  };
114 
118  class CheckValueGtMax {
119  public:
120  template<typename Tval, typename Tmax>
121  bool operator()(Tval val, Tmax max) const {
122  return (static_cast<Tmax>(val) > max);
123  }
124  };
125 
129  class CheckClipRange {
130  public:
131  template<typename Tval, typename Tcen, typename Tmax>
132  bool operator()(Tval val, Tcen center, Tmax cliplimit) const {
133  Tmax tmp = fabs(val - center);
134  return (tmp <= cliplimit);
135  }
136  };
137 
138  // define some abbreviated typenames for the test templates
139  typedef CheckFinite ChkFin;
140  typedef CheckValueLtMin ChkMin;
141  typedef CheckValueGtMax ChkMax;
142  typedef CheckClipRange ChkClip;
143  typedef AlwaysTrue AlwaysT;
144  typedef AlwaysFalse AlwaysF;
145 
146  // Return the variance of a variance, assuming a Gaussian
147  // There is apparently an attempt to correct for bias in the factor (n - 1)/n. RHL
148  inline double varianceError(double const variance, int const n)
149  {
150  return 2*(n - 1)*variance*variance/static_cast<double>(n*n);
151  }
152 
153  /*********************************************************************************************************/
154  // return type for processPixels
155  typedef boost::tuple<int, // n
156  double, // sum
158  afwMath::Statistics::Value, // variance
159  double, // min
160  double, // max
161  lsst::afw::image::MaskPixel // allPixelOrMask
162  > StandardReturn;
163 
164  /*
165  * Functions which convert the booleans into calls to the proper templated types, one type per
166  * recursion level
167  */
168  /*
169  * This function handles the inner summation loop, with tests templated
170  *
171  * The idea here is to allow different conditionals in the inner loop, but avoid repeating code.
172  * Each test is actually a functor which is handled through a template. If the
173  * user requests a test (eg check for NaNs), the function is instantiated with the appropriate functor.
174  * Otherwise, an 'AlwaysTrue' or 'AlwaysFalse' object is passed in. The compiler then compiles-out
175  * a test which is always false, or removes the conditional for a test which is always true.
176  */
177  template<typename IsFinite,
178  typename HasValueLtMin,
179  typename HasValueGtMax,
180  typename InClipRange,
181  bool useWeights,
182  typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
183  StandardReturn processPixels(ImageT const &img,
184  MaskT const &msk,
185  VarianceT const &var,
186  WeightT const &weights,
187  int const,
188  int const nCrude,
189  int const stride,
190  double const meanCrude,
191  double const cliplimit,
192  bool const weightsAreMultiplicative,
193  int const andMask,
194  bool const calcErrorFromInputVariance,
195  std::vector<double> const & maskPropagationThresholds
196  )
197  {
198  int n = 0;
199  double sumw = 0.0; // sum(weight) (N.b. weight will be 1.0 if !useWeights)
200  double sumw2 = 0.0; // sum(weight^2)
201  double sumx = 0; // sum(data*weight)
202  double sumx2 = 0; // sum(data*weight^2)
203 #if 1
204  double sumvw2 = 0.0; // sum(variance*weight^2)
205 #endif
206  double min = (nCrude) ? meanCrude : MAX_DOUBLE;
207  double max = (nCrude) ? meanCrude : -MAX_DOUBLE;
208 
209  afwImage::MaskPixel allPixelOrMask = 0x0;
210 
211  std::vector<double> rejectedWeightsByBit(maskPropagationThresholds.size(), 0.0);
212 
213  for (int iY = 0; iY < img.getHeight(); iY += stride) {
214 
215  typename MaskT::x_iterator mptr = msk.row_begin(iY);
216  typename VarianceT::x_iterator vptr = var.row_begin(iY);
217  typename WeightT::x_iterator wptr = weights.row_begin(iY);
218 
219  for (typename ImageT::x_iterator ptr = img.row_begin(iY), end = ptr + img.getWidth();
220  ptr != end; ++ptr, ++mptr, ++vptr, ++wptr) {
221 
222  if (IsFinite()(*ptr) && !(*mptr & andMask) &&
223  InClipRange()(*ptr, meanCrude, cliplimit) ) { // clip
224 
225  double const delta = (*ptr - meanCrude);
226 
227  if (useWeights) {
228  double weight = *wptr;
229  if (weightsAreMultiplicative) {
230  ;
231  } else {
232  if (*wptr <= 0) {
233  continue;
234  }
235  weight = 1/weight;
236  }
237 
238  sumw += weight;
239  sumw2 += weight*weight;
240  sumx += weight*delta;
241  sumx2 += weight*delta*delta;
242 
243  if (calcErrorFromInputVariance) {
244  double const var = *vptr;
245  sumvw2 += var*weight*weight;
246  }
247  } else {
248  sumx += delta;
249  sumx2 += delta*delta;
250 
251  if (calcErrorFromInputVariance) {
252  double const var = *vptr;
253  sumvw2 += var;
254  }
255  }
256 
257  allPixelOrMask |= *mptr;
258 
259  if (HasValueLtMin()(*ptr, min)) { min = *ptr; }
260  if (HasValueGtMax()(*ptr, max)) { max = *ptr; }
261  n++;
262  } else { // pixel has been clipped, rejected, etc.
263  for (int bit = 0, nBits=maskPropagationThresholds.size(); bit < nBits; ++bit) {
264  lsst::afw::image::MaskPixel mask = 1 << bit;
265  if (*mptr & mask) {
266  double weight = 1.0;
267  if (useWeights) {
268  weight = *wptr;
269  if (!weightsAreMultiplicative) {
270  if (*wptr <= 0) {
271  continue;
272  }
273  weight = 1.0 / weight;
274  }
275  }
276  rejectedWeightsByBit[bit] += weight;
277  }
278  }
279  }
280  }
281  }
282  if (n == 0) {
283  min = NaN;
284  max = NaN;
285  }
286 
287  // estimate of population mean and variance.
288  double mean, variance;
289  if (!useWeights) {
290  sumw = sumw2 = n;
291  }
292 
293  for (int bit = 0, nBits=maskPropagationThresholds.size(); bit < nBits; ++bit) {
294  double hypotheticalTotalWeight = sumw + rejectedWeightsByBit[bit];
295  rejectedWeightsByBit[bit] /= hypotheticalTotalWeight;
296  if (rejectedWeightsByBit[bit] > maskPropagationThresholds[bit]) {
297  allPixelOrMask |= (1 << bit);
298  }
299  }
300 
301 
302  // N.b. if sumw == 0 or sumw*sumw == sumw2 (e.g. n == 1) we'll get NaNs
303  // N.b. the estimator of the variance assumes that the sample points all have the same variance;
304  // otherwise, what is it that we're estimating?
305  mean = sumx/sumw;
306  variance = sumx2/sumw - ::pow(mean, 2); // biased estimator
307  variance *= sumw*sumw/(sumw*sumw - sumw2); // debias
308 
309  double meanVar; // (standard error of mean)^2
310  if (calcErrorFromInputVariance) {
311  meanVar = sumvw2/(sumw*sumw);
312  } else {
313  meanVar = variance*sumw2/(sumw*sumw);
314  }
315 
316  double varVar = varianceError(variance, n); // error in variance; incorrect if useWeights is true
317 
318  sumx += sumw*meanCrude;
319  mean += meanCrude;
320 
321  return StandardReturn(n, sumx,
322  afwMath::Statistics::Value(mean, meanVar),
323  afwMath::Statistics::Value(variance, varVar), min, max, allPixelOrMask);
324  }
325 
326  template<typename IsFinite,
327  typename HasValueLtMin,
328  typename HasValueGtMax,
329  typename InClipRange,
330  bool useWeights,
331  typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
332  StandardReturn processPixels(ImageT const &img,
333  MaskT const &msk,
334  VarianceT const &var,
335  WeightT const &weights,
336  int const flags,
337  int const nCrude,
338  int const stride,
339  double const meanCrude,
340  double const cliplimit,
341  bool const weightsAreMultiplicative,
342  int const andMask,
343  bool const calcErrorFromInputVariance,
344  bool doGetWeighted,
345  std::vector<double> const & maskPropagationThresholds
346  )
347  {
348  if (doGetWeighted) {
349  return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, true>(
350  img, msk, var, weights,
351  flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative, andMask,
352  calcErrorFromInputVariance, maskPropagationThresholds);
353  } else {
354  return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, false>(
355  img, msk, var, weights,
356  flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative, andMask,
357  calcErrorFromInputVariance, maskPropagationThresholds);
358  }
359  }
360 
361  template<typename IsFinite,
362  typename HasValueLtMin,
363  typename HasValueGtMax,
364  typename InClipRange,
365  bool useWeights,
366  typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
367  StandardReturn processPixels(ImageT const &img,
368  MaskT const &msk,
369  VarianceT const &var,
370  WeightT const &weights,
371  int const flags,
372  int const nCrude,
373  int const stride,
374  double const meanCrude,
375  double const cliplimit,
376  bool const weightsAreMultiplicative,
377  int const andMask,
378  bool const calcErrorFromInputVariance,
379  bool doCheckFinite,
380  bool doGetWeighted,
381  std::vector<double> const & maskPropagationThresholds
382  )
383  {
384  if (doCheckFinite) {
385  return processPixels<CheckFinite, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
386  img, msk, var, weights,
387  flags, nCrude, 1, meanCrude, cliplimit,
388  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
389  doGetWeighted, maskPropagationThresholds);
390  } else {
391  return processPixels<AlwaysTrue, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
392  img, msk, var, weights,
393  flags, nCrude, 1, meanCrude, cliplimit,
394  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
395  doGetWeighted, maskPropagationThresholds);
396  }
397  }
398 
407  template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
408  StandardReturn getStandard(ImageT const &img, // image
409  MaskT const &msk, // mask
410  VarianceT const &var, // variance
411  WeightT const &weights, // weights to apply to each pixel
412  int const flags, // what to measure
413  bool const weightsAreMultiplicative, // weights are multiplicative (not inverse)
414  int const andMask, // mask of bad pixels
415  bool const calcErrorFromInputVariance, // estimate errors from variance
416  bool doCheckFinite, // check for NaN/Inf
417  bool doGetWeighted, // use the weights
418  std::vector<double> const & maskPropagationThresholds
419  )
420  {
421  // =====================================================
422  // a crude estimate of the mean, used for numerical stability of variance
423  int nCrude = 0;
424  double meanCrude = 0.0;
425 
426  // for small numbers of values, use a small stride
427  int const nPix = img.getWidth()*img.getHeight();
428  int strideCrude;
429  if (nPix < 100) {
430  strideCrude = 2;
431  } else {
432  strideCrude = 10;
433  }
434 
435  double cliplimit = -1; // unused
436  StandardReturn values = processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
437  img, msk, var, weights,
438  flags, nCrude, strideCrude, meanCrude,
439  cliplimit,
440  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
441  doCheckFinite, doGetWeighted,
442  maskPropagationThresholds);
443  nCrude = values.get<0>();
444  double sumCrude = values.get<1>();
445 
446  meanCrude = 0.0;
447  if (nCrude > 0) {
448  meanCrude = sumCrude/nCrude;
449  }
450 
451  // =======================================================
452  // Estimate the full precision variance using that crude mean
453  // - get the min and max as well
454 
455  if (flags & (afwMath::MIN | afwMath::MAX)) {
456  return processPixels<ChkFin, ChkMin, ChkMax, AlwaysT, true>(
457  img, msk, var, weights,
458  flags, nCrude, 1, meanCrude,
459  cliplimit,
460  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
461  true, doGetWeighted, maskPropagationThresholds);
462  } else {
463  return processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT,true>(
464  img, msk, var, weights,
465  flags, nCrude, 1, meanCrude,
466  cliplimit,
467  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
468  doCheckFinite, doGetWeighted, maskPropagationThresholds);
469  }
470  }
471 
477  template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
478  StandardReturn getStandard(ImageT const &img, // image
479  MaskT const &msk, // mask
480  VarianceT const &var, // variance
481  WeightT const &weights, // weights to apply to each pixel
482  int const flags, // what to measure
483  std::pair<double, double> const clipinfo, // the center and cliplimit for the
484  // first clip iteration
485  bool const weightsAreMultiplicative, // weights are multiplicative (not inverse)
486  int const andMask, // mask of bad pixels
487  bool const calcErrorFromInputVariance, // estimate errors from variance
488  bool doCheckFinite, // check for NaN/Inf
489  bool doGetWeighted, // use the weights,
490  std::vector<double> const & maskPropagationThresholds
491  )
492  {
493  double const center = clipinfo.first;
494  double const cliplimit = clipinfo.second;
495 
496  if (lsst::utils::isnan(center) || lsst::utils::isnan(cliplimit)) {
497  return StandardReturn(0, NaN,
498  afwMath::Statistics::Value(NaN, NaN),
499  afwMath::Statistics::Value(NaN, NaN), NaN, NaN, ~0x0);
500  }
501 
502  // =======================================================
503  // Estimate the full precision variance using that crude mean
504  int const stride = 1;
505  int nCrude = 0;
506 
507  if (flags & (afwMath::MIN | afwMath::MAX)) {
508  return processPixels<ChkFin, ChkMin, ChkMax, ChkClip, true>(
509  img, msk, var, weights,
510  flags, nCrude, stride,
511  center, cliplimit,
512  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
513  true, doGetWeighted, maskPropagationThresholds);
514  } else { // fast loop ... just the mean & variance
515  return processPixels<ChkFin, AlwaysF, AlwaysF, ChkClip, true>(
516  img, msk, var, weights,
517  flags, nCrude, stride,
518  center, cliplimit,
519  weightsAreMultiplicative, andMask, calcErrorFromInputVariance,
520  doCheckFinite, doGetWeighted, maskPropagationThresholds);
521  }
522  }
523 
532  template<typename Pixel>
533  double percentile(std::vector<Pixel> &img,
534  double const fraction)
535  {
536  assert (fraction >= 0.0 && fraction <= 1.0);
537 
538  int const n = img.size();
539 
540  if (n > 1) {
541 
542  double const idx = fraction*(n - 1);
543 
544  // interpolate linearly between the adjacent values
545  // For efficiency:
546  // - if we're asked for a fraction > 0.5,
547  // we'll do the second partial sort on shorter (upper) portion
548  // - otherwise, the shorter portion will be the lower one, we'll partial-sort that.
549 
550  int const q1 = static_cast<int>(idx);
551  int const q2 = q1 + 1;
552 
553  typename std::vector<Pixel>::iterator mid1 = img.begin() + q1;
554  typename std::vector<Pixel>::iterator mid2 = img.begin() + q2;
555  if (fraction > 0.5) {
556  std::nth_element(img.begin(), mid1, img.end());
557  std::nth_element(mid1, mid2, img.end());
558  } else {
559  std::nth_element(img.begin(), mid2, img.end());
560  std::nth_element(img.begin(), mid1, mid2);
561  }
562 
563  double val1 = static_cast<double>(*mid1);
564  double val2 = static_cast<double>(*mid2);
565  double w1 = (static_cast<double>(q2) - idx);
566  double w2 = (idx - static_cast<double>(q1));
567  return w1*val1 + w2*val2;
568 
569  } else if (n == 1) {
570  return img[0];
571  } else {
572  return NaN;
573  }
574  }
575 
576 
585  typedef boost::tuple<double, double, double> MedianQuartileReturn;
586 
587  template<typename Pixel>
588  MedianQuartileReturn medianAndQuartiles(std::vector<Pixel> &img)
589  {
590  int const n = img.size();
591 
592  if (n > 1) {
593 
594  double const idx50 = 0.50*(n - 1);
595  double const idx25 = 0.25*(n - 1);
596  double const idx75 = 0.75*(n - 1);
597 
598  // For efficiency:
599  // - partition at 50th, then partition the two half further to get 25th and 75th
600  // - to get the adjacent points (for interpolation), partition between 25/50, 50/75, 75/end
601  // these should be much smaller partitions
602 
603  int const q50a = static_cast<int>(idx50);
604  int const q50b = q50a + 1;
605  int const q25a = static_cast<int>(idx25);
606  int const q25b = q25a + 1;
607  int const q75a = static_cast<int>(idx75);
608  int const q75b = q75a + 1;
609 
610  typename std::vector<Pixel>::iterator mid50a = img.begin() + q50a;
611  typename std::vector<Pixel>::iterator mid50b = img.begin() + q50b;
612  typename std::vector<Pixel>::iterator mid25a = img.begin() + q25a;
613  typename std::vector<Pixel>::iterator mid25b = img.begin() + q25b;
614  typename std::vector<Pixel>::iterator mid75a = img.begin() + q75a;
615  typename std::vector<Pixel>::iterator mid75b = img.begin() + q75b;
616 
617  // get the 50th percentile, then get the 25th and 75th on the smaller partitions
618  std::nth_element(img.begin(), mid50a, img.end());
619  std::nth_element(mid50a, mid75a, img.end());
620  std::nth_element(img.begin(), mid25a, mid50a);
621 
622  // and the adjacent points for each ... use the smallest segments available.
623  std::nth_element(mid50a, mid50b, mid75a);
624  std::nth_element(mid25a, mid25b, mid50a);
625  std::nth_element(mid75a, mid75b, img.end());
626 
627  // interpolate linearly between the adjacent values
628  double val50a = static_cast<double>(*mid50a);
629  double val50b = static_cast<double>(*mid50b);
630  double w50a = (static_cast<double>(q50b) - idx50);
631  double w50b = (idx50 - static_cast<double>(q50a));
632  double median = w50a*val50a + w50b*val50b;
633 
634  double val25a = static_cast<double>(*mid25a);
635  double val25b = static_cast<double>(*mid25b);
636  double w25a = (static_cast<double>(q25b) - idx25);
637  double w25b = (idx25 - static_cast<double>(q25a));
638  double q1 = w25a*val25a + w25b*val25b;
639 
640  double val75a = static_cast<double>(*mid75a);
641  double val75b = static_cast<double>(*mid75b);
642  double w75a = (static_cast<double>(q75b) - idx75);
643  double w75b = (idx75 - static_cast<double>(q75a));
644  double q3 = w75a*val75a + w75b*val75b;
645 
646  return MedianQuartileReturn(median, q1, q3);
647  } else if (n == 1) {
648  return MedianQuartileReturn(img[0], img[0], img[0]);
649  } else {
650  return MedianQuartileReturn(NaN, NaN, NaN);
651  }
652  }
653 
654  /*********************************************************************************************************/
662  template<typename IsFinite, typename ImageT, typename MaskT, typename VarianceT>
663  boost::shared_ptr<std::vector<typename ImageT::Pixel> > makeVectorCopy(ImageT const &img,
664  MaskT const &msk,
665  VarianceT const &,
666  int const andMask
667  )
668  {
669  // Note need to keep track of allPixelOrMask here ... processPixels() does that
670  // and it always gets called
671  boost::shared_ptr<std::vector<typename ImageT::Pixel> >
672  imgcp(new std::vector<typename ImageT::Pixel>(0));
673 
674  for (int i_y = 0; i_y < img.getHeight(); ++i_y) {
675  typename MaskT::x_iterator mptr = msk.row_begin(i_y);
676  for (typename ImageT::x_iterator ptr = img.row_begin(i_y), end = img.row_end(i_y);
677  ptr != end; ++ptr) {
678  if (IsFinite()(*ptr) && !(*mptr & andMask)) {
679  imgcp->push_back(*ptr);
680  }
681  ++mptr;
682  }
683  }
684 
685  return imgcp;
686  }
687 }
688 
689 
690 
692  int oldSize = _maskPropagationThresholds.size();
693  if (oldSize < bit) {
694  return 1.0;
695  }
696  return _maskPropagationThresholds[bit];
697 }
698 
700  int oldSize = _maskPropagationThresholds.size();
701  if (oldSize <= bit) {
702  int newSize = bit + 1;
703  _maskPropagationThresholds.resize(newSize);
704  for (int i = oldSize; i < bit; ++i) {
705  _maskPropagationThresholds[i] = 1.0;
706  }
707  }
708  _maskPropagationThresholds[bit] = threshold;
709 }
710 
711 
716  static std::map<std::string, Property> statisticsProperty;
717  if (statisticsProperty.size() == 0) {
718  statisticsProperty["NOTHING"] = NOTHING;
719  statisticsProperty["ERRORS"] = ERRORS;
720  statisticsProperty["NPOINT"] = NPOINT;
721  statisticsProperty["MEAN"] = MEAN;
722  statisticsProperty["STDEV"] = STDEV;
723  statisticsProperty["VARIANCE"] = VARIANCE;
724  statisticsProperty["MEDIAN"] = MEDIAN;
725  statisticsProperty["IQRANGE"] = IQRANGE;
726  statisticsProperty["MEANCLIP"] = MEANCLIP;
727  statisticsProperty["STDEVCLIP"] = STDEVCLIP;
728  statisticsProperty["VARIANCECLIP"] = VARIANCECLIP;
729  statisticsProperty["MIN"] = MIN;
730  statisticsProperty["MAX"] = MAX;
731  statisticsProperty["SUM"] = SUM;
732  statisticsProperty["MEANSQUARE"] = MEANSQUARE;
733  statisticsProperty["ORMASK"] = ORMASK;
734  }
735  return statisticsProperty[property];
736 }
737 
745 template<typename ImageT, typename MaskT, typename VarianceT>
747  ImageT const &img,
748  MaskT const &msk,
749  VarianceT const &var,
750  int const flags,
751  afwMath::StatisticsControl const& sctrl
752  ) :
753  _flags(flags), _mean(NaN, NaN), _variance(NaN, NaN), _min(NaN), _max(NaN), _sum(NaN),
754  _meanclip(NaN, NaN), _varianceclip(NaN, NaN), _median(NaN, NaN), _iqrange(NaN),
755  _sctrl(sctrl), _weightsAreMultiplicative(false)
756 {
757  doStatistics(img, msk, var, var, _flags, _sctrl);
758 }
759 
760 namespace {
761  template<typename T>
762  bool isEmpty(T const& t) { return t.empty(); }
763 
764  template<typename T>
765  bool isEmpty(afwImage::Image<T> const& im) { return (im.getWidth() == 0 && im.getHeight() == 0); }
766 }
767 
768 template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
770  ImageT const &img,
771  MaskT const &msk,
772  VarianceT const &var,
773  WeightT const &weights,
774  int const flags,
775  afwMath::StatisticsControl const& sctrl
776  ) :
777  _flags(flags), _mean(NaN, NaN), _variance(NaN, NaN), _min(NaN), _max(NaN), _sum(NaN),
778  _meanclip(NaN, NaN), _varianceclip(NaN, NaN), _median(NaN, NaN), _iqrange(NaN),
779  _sctrl(sctrl), _weightsAreMultiplicative(true)
780 {
781  if (!isEmpty(weights)) {
783  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
784  "You must use the weights if you provide them");
785  }
786 
787  _sctrl.setWeighted(true);
788  }
789  doStatistics(img, msk, var, weights, _flags, _sctrl);
790 }
791 
792 template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
794  ImageT const &img,
795  MaskT const &msk,
796  VarianceT const &var,
797  WeightT const &weights,
798  int const flags,
799  afwMath::StatisticsControl const& sctrl
800  )
801 {
802  _n = img.getWidth()*img.getHeight();
803  if (_n == 0) {
804  throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
805  }
806 
807  // Check that an int's large enough to hold the number of pixels
808  assert(img.getWidth()*static_cast<double>(img.getHeight()) < std::numeric_limits<int>::max());
809 
810  // get the standard statistics
811  StandardReturn standard = getStandard(img, msk, var, weights, flags,
812  _weightsAreMultiplicative,
813  _sctrl.getAndMask(),
814  _sctrl.getCalcErrorFromInputVariance(),
815  _sctrl.getNanSafe(), _sctrl.getWeighted(),
816  _sctrl._maskPropagationThresholds);
817 
818  _n = standard.get<0>();
819  _sum = standard.get<1>();
820  _mean = standard.get<2>();
821  _variance = standard.get<3>();
822  _min = standard.get<4>();
823  _max = standard.get<5>();
824  _allPixelOrMask = standard.get<6>();
825 
826  // ==========================================================
827  // now only calculate it if it's specifically requested - these all cost more!
828 
829  // copy the image for any routines that will use median or quantiles
830  if (flags & (MEDIAN | IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
831 
832  // make a vector copy of the image to get the median and quartiles (will move values)
833  boost::shared_ptr<std::vector<typename ImageT::Pixel> > imgcp;
834  if (_sctrl.getNanSafe()) {
835  imgcp = makeVectorCopy<ChkFin>(img, msk, var, _sctrl.getAndMask());
836  } else {
837  imgcp = makeVectorCopy<AlwaysT>(img, msk, var, _sctrl.getAndMask());
838  }
839 
840  // if we *only* want the median, just use percentile(), otherwise use medianAndQuartiles()
841  if ( (flags & (MEDIAN)) && !(flags & (IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP)) ) {
842  _median = Value(percentile(*imgcp, 0.5), NaN);
843  } else {
844  MedianQuartileReturn mq = medianAndQuartiles(*imgcp);
845  _median = Value(mq.get<0>(), NaN);
846  _iqrange = mq.get<2>() - mq.get<1>();
847  }
848 
849 
850  if (flags & (MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
851  for (int i_i = 0; i_i < _sctrl.getNumIter(); ++i_i) {
852  double const center = ((i_i > 0) ? _meanclip : _median).first;
853  double const hwidth = (i_i > 0 && _n > 1) ?
854  _sctrl.getNumSigmaClip()*std::sqrt(_varianceclip.first) :
855  _sctrl.getNumSigmaClip()*IQ_TO_STDEV*_iqrange;
856  std::pair<double, double> const clipinfo(center, hwidth);
857 
858  StandardReturn clipped = getStandard(img, msk, var, weights, flags, clipinfo,
859  _weightsAreMultiplicative,
860  _sctrl.getAndMask(),
861  _sctrl.getCalcErrorFromInputVariance(),
862  _sctrl.getNanSafe(), _sctrl.getWeighted(),
863  _sctrl._maskPropagationThresholds);
864 
865  int const nClip = clipped.get<0>();
866  _meanclip = clipped.get<2>(); // clipped mean
867  double const varClip = clipped.get<3>().first; // clipped variance
868 
869  _varianceclip = Value(varClip, varianceError(varClip, nClip));
870  // ... ignore other values
871  }
872  }
873  }
874 }
875 
876 /************************************************************************************************************/
887 std::pair<double, double> afwMath::Statistics::getResult(
888  afwMath::Property const iProp
889  ) const {
893 
894  // if iProp == NOTHING try to return their heart's delight, as specified in the constructor
895  afwMath::Property const prop =
896  static_cast<afwMath::Property>(((iProp == NOTHING) ? _flags : iProp) & ~ERRORS);
897 
898  if (!(prop & _flags)) { // we didn't calculate it
899  throw LSST_EXCEPT(pexExceptions::InvalidParameterError,
900  (boost::format("You didn't ask me to calculate %d") % prop).str());
901  }
902 
903 
904  Value ret(NaN, NaN);
905  switch (prop) {
906 
907  case NPOINT:
908  ret.first = static_cast<double>(_n);
909  if (_flags & ERRORS) {
910  ret.second = 0;
911  }
912  break;
913 
914  case SUM:
915  ret.first = static_cast<double>(_sum);
916  if (_flags & ERRORS) {
917  ret.second = 0;
918  }
919  break;
920 
921  // == means ==
922  case MEAN:
923  ret.first = _mean.first;
924  if (_flags & ERRORS) {
925  ret.second = ::sqrt(_mean.second);
926  }
927  break;
928  case MEANCLIP:
929  ret.first = _meanclip.first;
930  if ( _flags & ERRORS ) {
931  ret.second = ::sqrt(_meanclip.second);
932  }
933  break;
934 
935  // == stdevs & variances ==
936  case VARIANCE:
937  ret.first = _variance.first;
938  if (_flags & ERRORS) {
939  ret.second = ::sqrt(_variance.second);
940  }
941  break;
942  case STDEV:
943  ret.first = sqrt(_variance.first);
944  if (_flags & ERRORS) {
945  ret.second = 0.5*::sqrt(_variance.second)/ret.first;
946  }
947  break;
948  case VARIANCECLIP:
949  ret.first = _varianceclip.first;
950  if (_flags & ERRORS) {
951  ret.second = ret.second;
952  }
953  break;
954  case STDEVCLIP:
955  ret.first = sqrt(_varianceclip.first);
956  if (_flags & ERRORS) {
957  ret.second = 0.5*::sqrt(_varianceclip.second)/ret.first;
958  }
959  break;
960 
961  case MEANSQUARE:
962  ret.first = (_n - 1)/static_cast<double>(_n)*_variance.first + ::pow(_mean.first, 2);
963  if (_flags & ERRORS) {
964  ret.second = ::sqrt(2*::pow(ret.first/_n, 2)); // assumes Gaussian
965  }
966  break;
967 
968  // == other stats ==
969  case MIN:
970  ret.first = _min;
971  if ( _flags & ERRORS ) {
972  ret.second = 0;
973  }
974  break;
975  case MAX:
976  ret.first = _max;
977  if ( _flags & ERRORS ) {
978  ret.second = 0;
979  }
980  break;
981  case MEDIAN:
982  ret.first = _median.first;
983  if ( _flags & ERRORS ) {
984  ret.second = sqrt(afwGeom::HALFPI*_variance.first/_n); // assumes Gaussian
985  }
986  break;
987  case IQRANGE:
988  ret.first = _iqrange;
989  if ( _flags & ERRORS ) {
990  ret.second = NaN; // we're not estimating this properly
991  }
992  break;
993 
994  // no-op to satisfy the compiler
995  case ERRORS:
996  break;
997  // default: redundant as 'ret' is initialized to NaN, NaN
998  default: // we must have set prop to _flags
999  assert (iProp == 0);
1000  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
1001  "getValue() may only be called without a parameter"
1002  " if you asked for only one statistic");
1003  }
1004  return ret;
1005 }
1006 
1010  afwMath::Property const prop
1011  ) const {
1014  return getResult(prop).first;
1015 }
1016 
1017 
1021  afwMath::Property const prop
1022  ) const {
1026  return getResult(prop).second;
1027 }
1028 
1029 
1030 /************************************************************************************************/
1035 namespace lsst {
1036 namespace afw {
1037 namespace math {
1038 
1039 template<>
1044  int const flags,
1045  StatisticsControl const& sctrl
1046  ) :
1047  _flags(flags),
1048  _mean(NaN, NaN), _variance(NaN, NaN), _min(NaN), _max(NaN),
1049  _meanclip(NaN, NaN), _varianceclip(NaN, NaN), _median(NaN, NaN), _iqrange(NaN),
1050  _sctrl(sctrl) {
1051 
1052  if ((flags & ~(NPOINT | SUM)) != 0x0) {
1053  throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Statistics<Mask> only supports NPOINT and SUM");
1054  }
1055 
1057 
1058  _n = msk.getWidth()*msk.getHeight();
1059  if (_n == 0) {
1060  throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
1061  }
1062 
1063  // Check that an int's large enough to hold the number of pixels
1064  assert(msk.getWidth()*static_cast<double>(msk.getHeight()) < std::numeric_limits<int>::max());
1065 
1066  afwImage::MaskPixel sum = 0x0;
1067  for (int y = 0; y != msk.getHeight(); ++y) {
1068  for (Mask::x_iterator ptr = msk.row_begin(y), end = msk.row_end(y); ptr != end; ++ptr) {
1069  sum |= (*ptr)[0];
1070  }
1071  }
1072  _sum = sum;
1073 }
1074 
1084  int const flags,
1085  StatisticsControl const& sctrl
1086  ) {
1087  return Statistics(msk, msk, msk, flags, sctrl);
1088 }
1089 
1090 }}}
1091 
1092 /****************************************************************************************************/
1093 /*
1094  * Explicit instantiations
1095  *
1096  * explicit Statistics(MaskedImage const& img, int const flags,
1097  * StatisticsControl const& sctrl=StatisticsControl());
1098  */
1100 //
1101 #define STAT afwMath::Statistics
1102 
1103 typedef afwImage::VariancePixel VPixel;
1104 
1105 #define INSTANTIATE_MASKEDIMAGE_STATISTICS(TYPE) \
1106  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1107  afwImage::Mask<afwImage::MaskPixel> const &msk, \
1108  afwImage::Image<VPixel> const &var, \
1109  int const flags, StatisticsControl const& sctrl); \
1110  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1111  afwImage::Mask<afwImage::MaskPixel> const &msk, \
1112  afwImage::Image<VPixel> const &var, \
1113  afwImage::Image<VPixel> const &weights, \
1114  int const flags, StatisticsControl const& sctrl); \
1115  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1116  afwImage::Mask<afwImage::MaskPixel> const &msk, \
1117  afwImage::Image<VPixel> const &var, \
1118  afwMath::ImageImposter<VPixel> const &weights, \
1119  int const flags, StatisticsControl const& sctrl)
1120 
1121 #define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(TYPE) \
1122  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1123  afwMath::MaskImposter<afwImage::MaskPixel> const &msk, \
1124  afwImage::Image<VPixel> const &var, \
1125  int const flags, StatisticsControl const& sctrl); \
1126  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1127  afwMath::MaskImposter<afwImage::MaskPixel> const &msk, \
1128  afwImage::Image<VPixel> const &var, \
1129  afwImage::Image<VPixel> const &weights, \
1130  int const flags, StatisticsControl const& sctrl)
1131 
1132 #define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(TYPE) \
1133  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1134  afwImage::Mask<afwImage::MaskPixel> const &msk, \
1135  afwMath::MaskImposter<VPixel> const &var, \
1136  int const flags, StatisticsControl const& sctrl); \
1137  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1138  afwImage::Mask<afwImage::MaskPixel> const &msk, \
1139  afwMath::MaskImposter<VPixel> const &var, \
1140  afwImage::Image<VPixel> const &weights, \
1141  int const flags, StatisticsControl const& sctrl); \
1142  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1143  afwImage::Mask<afwImage::MaskPixel> const &msk, \
1144  afwMath::MaskImposter<VPixel> const &var, \
1145  afwMath::ImageImposter<VPixel> const &weights, \
1146  int const flags, StatisticsControl const& sctrl)
1147 
1148 #define INSTANTIATE_REGULARIMAGE_STATISTICS(TYPE) \
1149  template STAT::Statistics(afwImage::Image<TYPE> const &img, \
1150  afwMath::MaskImposter<afwImage::MaskPixel> const &msk, \
1151  afwMath::MaskImposter<VPixel> const &var, \
1152  int const flags, StatisticsControl const& sctrl)
1153 
1154 #define INSTANTIATE_VECTOR_STATISTICS(TYPE) \
1155  template STAT::Statistics(afwMath::ImageImposter<TYPE> const &img, \
1156  afwMath::MaskImposter<afwImage::MaskPixel> const &msk, \
1157  afwMath::MaskImposter<VPixel> const &var, \
1158  int const flags, StatisticsControl const& sctrl); \
1159  template STAT::Statistics(afwMath::ImageImposter<TYPE> const &img, \
1160  afwMath::MaskImposter<afwImage::MaskPixel> const &msk, \
1161  afwMath::MaskImposter<VPixel> const &var, \
1162  afwMath::ImageImposter<VPixel> const &weights, \
1163  int const flags, StatisticsControl const& sctrl)
1164 
1165 #define INSTANTIATE_IMAGE_STATISTICS(T) \
1166  INSTANTIATE_MASKEDIMAGE_STATISTICS(T); \
1167  INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(T); \
1168  INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(T); \
1169  INSTANTIATE_REGULARIMAGE_STATISTICS(T); \
1170  INSTANTIATE_VECTOR_STATISTICS(T)
1171 
1172 INSTANTIATE_IMAGE_STATISTICS(double);
1173 INSTANTIATE_IMAGE_STATISTICS(float);
1174 INSTANTIATE_IMAGE_STATISTICS(int);
1175 INSTANTIATE_IMAGE_STATISTICS(boost::uint16_t);
1176 INSTANTIATE_IMAGE_STATISTICS(boost::uint64_t);
1177 
int y
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1009
boost::uint16_t MaskPixel
double getMaskPropagationThreshold(int bit) const
Definition: Statistics.cc:691
estimate sample mean
Definition: Statistics.h:67
std::pair< double, double > Value
The type used to report (value, error) for desired statistics.
Definition: Statistics.h:215
x_iterator row_begin(int y) const
Definition: Image.h:319
estimate sample N-sigma clipped mean (N set in StatisticsControl, default=3)
Definition: Statistics.h:72
tbl::Key< double > weight
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Constructor for Statistics object.
Definition: Statistics.cc:746
estimate sample standard deviation
Definition: Statistics.h:68
estimate sample variance
Definition: Statistics.h:69
find sum of pixels in the image
Definition: Statistics.h:78
Value getResult(Property const prop=NOTHING) const
Return the value and error in the specified statistic (e.g. MEAN)
Definition: Statistics.cc:887
int isnan(T t)
Definition: ieee.h:110
estimate sample N-sigma clipped stdev (N set in StatisticsControl, default=3)
Definition: Statistics.h:73
int const x0
Definition: saturated.cc:45
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
find mean value of square of pixel values
Definition: Statistics.h:79
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Include errors of requested quantities.
Definition: Statistics.h:65
We don&#39;t want anything.
Definition: Statistics.h:64
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
Pass parameters to a Statistics objectA class to pass parameters which control how the stats are calc...
Definition: Statistics.h:92
ImageT::Pixel _sum
Definition: CR.cc:313
void doStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights, int const flags, StatisticsControl const &sctrl)
Definition: Statistics.cc:793
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:93
Support for 2-D images.
StatisticsControl _sctrl
Definition: Statistics.h:255
float VariancePixel
! default type for Masks and MaskedImage Masks
estimate sample maximum
Definition: Statistics.h:77
number of sample points
Definition: Statistics.h:66
double getError(Property const prop=NOTHING) const
Return the error in the desired property (if specified in the constructor)
Definition: Statistics.cc:1020
int isfinite(T t)
Definition: ieee.h:100
estimate sample minimum
Definition: Statistics.h:76
estimate sample median
Definition: Statistics.h:70
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
void setMaskPropagationThreshold(int bit, double threshold)
Definition: Statistics.cc:699
void setWeighted(bool useWeights)
Definition: Statistics.h:149
double const HALFPI
Definition: Angle.h:20
int _n
estimate sample inter-quartile range
Definition: Statistics.h:71
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
Compute Image Statistics.
get the or-mask of all pixels used.
Definition: Statistics.h:80
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:715
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
Property
control what is calculated
Definition: Statistics.h:63
ImageT val
Definition: CR.cc:154
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
Include files required for standard LSST Exception handling.