libctsupport/mathfuncs.cpp

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