libctsupport/mathfuncs.cpp

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (c) 1983-2009 Kevin Rosenberg
   4 **
   5 **  This program is free software; you can redistribute it and/or modify
   6 **  it under the terms of the GNU General Public License (version 2) as
   7 **  published by the Free Software Foundation.
   8 **
   9 **  This program is distributed in the hope that it will be useful,
  10 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 **  GNU General Public License for more details.
  13 **
  14 **  You should have received a copy of the GNU General Public License
  15 **  along with this program; if not, write to the Free Software
  16 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  17 ******************************************************************************/
  18
  19 #include "ctsupport.h"
  20
  21
  22 /* NAME
  23 *    integrateSimpson         Integrate array of data by Simpson's rule
  24 *
  25 * SYNOPSIS
  26 *    double integrateSimpson (xmin, xmax, y, np)
  27 *    double xmin, xmax          Extent of integration
  28 *    double y[]         Function values to be integrated
  29 *    int np                     number of data points
  30 *                               (must be an odd number and at least 3)
  31 *
  32 * RETURNS
  33 *    integrand of function
  34 */
  35
  36 double
  37 integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
  38 {
  39   if (np < 2)
  40     return (0.);
  41   else if (np == 2)
  42     return ((xmax - xmin) * (y[0] + y[1]) / 2);
  43
  44   double area = 0;
  45   int nDiv = (np - 1) / 2;  // number of divisions
  46   double width = (xmax - xmin) / (double) (np - 1);     // width of cells
  47
  48   for (int i = 1; i <= nDiv; i++) {
  49     int xr = 2 * i;
  50     int xl = xr - 2;       // 2 * (i - 1) == 2 * i - 2 == xr - 2
  51     int xm = xr - 1;       // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
  52
  53     area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
  54   }
  55
  56   if ((np & 1) == 0)            /* do last trapazoid */
  57     area += width * (y[np-2] + y[np-1]) / 2;
  58
  59   return (area);
  60 }
  61
  62
  63 /* NAME
  64 *    normalizeAngle       Normalize angle to 0 to 2 * PI range
  65 *
  66 * SYNOPSIS
  67 *    t = normalizeAngle (theta)
  68 *    double t          Normalized angle
  69 *    double theta     Input angle
  70 */
  71
  72 double
  73 normalizeAngle (double theta)
  74 {
  75   while (theta < 0.)
  76     theta += TWOPI;
  77   while (theta >= TWOPI)
  78     theta -= TWOPI;
  79
  80   return (theta);
  81 }
  82
  83
  84 void
  85 vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)
  86 {
  87   if (nPoints <= 0)
  88     return;
  89
  90   mean = 0;
  91   min = vec[0];
  92   max = vec[0];
  93   int i;
  94   for (i = 0; i < nPoints; i++) {
  95     double v = vec[i];
  96     if (v > max)
  97       max = v;
  98     if (v < min)
  99       min = v;
 100     mean += v;
 101   }
 102   mean /= nPoints;
 103
 104   static const int nbin = 1024;
 105   int hist[ nbin ] = {0};
 106   double spread = max - min;
 107   mode = 0;
 108   stddev = 0;
 109   for (i = 0; i < nPoints; i++) {
 110     double v = vec[i];
 111     int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
 112     hist[b]++;
 113     double diff = (v - mean);
 114     stddev += diff * diff;
 115   }
 116   stddev = sqrt (stddev / nPoints);
 117
 118   int max_binindex = 0;
 119   int max_bin = -1;
 120   for (int ibin = 0; ibin < nbin; ibin++) {
 121     if (hist[ibin] > max_bin) {
 122       max_bin = hist[ibin];
 123       max_binindex = ibin;
 124     }
 125   }
 126
 127   mode = (max_binindex * spread / (nbin - 1)) + min;
 128
 129   std::sort(vec.begin(), vec.end());
 130
 131   if (nPoints % 2)  // Odd
 132     median = vec[((nPoints - 1) / 2)];
 133   else        // Even
 134     median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
 135 }
 136
 137