include/interpolator.h

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (c) 1983-2001 Kevin Rosenberg
   4 **
   5 **  $Id: interpolator.h,v 1.8 2003/01/30 21:53:16 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
  22 class CubicSplineInterpolator {
  23 private:
  24   double *m_pdY2;  // second differential of y data
  25
  26   const double* const m_pdY;
  27   const int m_n;
  28
  29 public:
  30   CubicSplineInterpolator (const double* const y, int n);
  31
  32   ~CubicSplineInterpolator ();
  33
  34   double interpolate (double x);
  35 };
  36
  37
  38 class CubicPolyInterpolator {
  39 private:
  40   const double* const m_pdY;
  41   const int m_n;
  42
  43 public:
  44   CubicPolyInterpolator (const double* const y, int n);
  45
  46   ~CubicPolyInterpolator ();
  47
  48   double interpolate (double x);
  49 };
  50
  51
  52 template<class T>
  53 class BilinearInterpolator {
  54 private:
  55   T** const m_ppMatrix;
  56   const int m_nx;
  57   const int m_ny;
  58
  59 public:
  60   BilinearInterpolator (T** ppMatrix, unsigned int nx, unsigned int ny)
  61   : m_ppMatrix(ppMatrix), m_nx(nx), m_ny(ny)
  62   {}
  63
  64   T interpolate (double dXPos, double dYPos)
  65   {
  66     int iFloorX = static_cast<int>(floor(dXPos));
  67     int iFloorY = static_cast<int>(floor (dYPos));
  68     double dXFrac = dXPos - iFloorX;
  69     double dYFrac = dYPos - iFloorY;
  70
  71     T result = 0;
  72
  73     if (iFloorX < 0 || iFloorY < 0 || iFloorX > m_nx-1 || iFloorY > m_ny-1)
  74       result = 0;
  75     else if (iFloorX == m_nx - 1 && iFloorY == m_ny - 1)
  76       result = static_cast<T>(m_ppMatrix[m_nx-1][m_ny-1]);
  77     else if (iFloorX == m_nx - 1)
  78       result = static_cast<T>(m_ppMatrix[iFloorX][iFloorY] + dYFrac * (m_ppMatrix[iFloorX][iFloorY+1] - m_ppMatrix[iFloorX][iFloorY]));
  79     else if (iFloorY == m_ny - 1)
  80       result = static_cast<T>(m_ppMatrix[iFloorX][iFloorY] + dXFrac * (m_ppMatrix[iFloorX+1][iFloorY] - m_ppMatrix[iFloorX][iFloorY]));
  81     else
  82       result = static_cast<T>((1 - dXFrac) * (1 - dYFrac) * m_ppMatrix[iFloorX][iFloorY] +
  83         dXFrac * (1 - dYFrac) * m_ppMatrix[iFloorX+1][iFloorY] +
  84         dYFrac * (1 - dXFrac) * m_ppMatrix[iFloorX][iFloorY+1] +
  85         dXFrac * dYFrac * m_ppMatrix[iFloorX+1][iFloorY+1]);
  86
  87     return result;
  88   }
  89 };
  90
  91
  92 template<class T>
  93 class BicubicPolyInterpolator {
  94 private:
  95   T** const m_ppMatrix;
  96   const unsigned int m_nx;
  97   const unsigned int m_ny;
  98
  99 public:
 100   BicubicPolyInterpolator (T** ppMatrix, unsigned int nx, unsigned int ny)
 101   : m_ppMatrix(ppMatrix), m_nx(nx), m_ny(ny)
 102   {}
 103
 104   T interpolate (double dXPos, double dYPos)
 105   {
 106     // int iFloorX = static_cast<int>(floor (dXPos));
 107     // int iFloorY = static_cast<int>(floor (dYPos));
 108     // double dXFrac = dXPos - iFloorX;
 109     // double dYFrac = dYPos - iFloorY;
 110
 111     T result = 0;
 112
 113     // Need to add code
 114
 115     return result;
 116   }
 117 };
 118
 119
 120 template<class T>
 121 class LinearInterpolator {
 122 private:
 123   T* const m_pX;
 124   T* const m_pY;
 125   const int m_n;
 126   const bool m_bZeroOutsideRange;
 127
 128 public:
 129   LinearInterpolator (T* pY, unsigned int n, bool bZeroOutside = true)
 130   : m_pX(0), m_pY(pY), m_n(n), m_bZeroOutsideRange(bZeroOutside)
 131   {}
 132
 133   LinearInterpolator (T* pX, T* pY, unsigned int n, bool bZeroOutside = true)
 134   : m_pX(pX), m_pY(pY), m_n(n), m_bZeroOutsideRange(bZeroOutside)
 135   {}
 136
 137   double interpolate (double dX, int* piLastFloor = NULL)
 138   {
 139     double result = 0;
 140
 141     if (! m_pX) {
 142       if (dX == 0)
 143         result = m_pY[0];
 144       else if (dX < 0) {
 145         if (m_bZeroOutsideRange)
 146           result = 0;
 147         else
 148           result = m_pY[0];
 149       } else if (dX == m_n - 1)
 150         result = m_pY[m_n-1];
 151       else if (dX > m_n - 1) {
 152         if (m_bZeroOutsideRange)
 153           result = 0;
 154         else
 155           result = m_pY[m_n - 1];
 156       } else {
 157        int iFloor = static_cast<int>(floor(dX));
 158        result = m_pY[iFloor] + (m_pY[iFloor+1] - m_pY[iFloor]) * (dX - iFloor);
 159       }
 160     } else {
 161       int iLower = -1;
 162       int iUpper = m_n;
 163       if (piLastFloor && *piLastFloor >= 0 && m_pX[*piLastFloor] < dX)
 164         iLower = *piLastFloor;
 165
 166       while (iUpper - iLower > 1) {
 167         int iMiddle = (iUpper + iLower) >> 1;
 168         if (dX >= m_pX[iMiddle])
 169           iLower = iMiddle;
 170         else
 171           iUpper = iMiddle;
 172       }
 173       if (dX == m_pX[0])
 174         result = m_pY[0];
 175       else if (dX < m_pX[0]) {
 176         if (m_bZeroOutsideRange)
 177           result = 0;
 178         else
 179           result = m_pY[0];
 180       } else if (dX == m_pX[m_n-1])
 181         result = m_pY[m_n-1];
 182       else if (dX > m_pX[m_n - 1]) {
 183         if (m_bZeroOutsideRange)
 184           result = 0;
 185         else
 186           result = m_pY[m_n - 1];
 187       } else {
 188         if (iLower < 0 || iLower >= m_n) {
 189           sys_error (ERR_SEVERE, "Coordinate out of range [linearInterpolate]");
 190           return 0;
 191         }
 192
 193        if (piLastFloor)
 194          *piLastFloor = iLower;
 195        result = m_pY[iLower] + (m_pY[iUpper] - m_pY[iLower]) * ((dX - m_pX[iLower]) / (m_pX[iUpper] - m_pX[iLower]));
 196       }
 197     }
 198
 199     return result;
 200   }
 201 };
 202