X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=include%2Finterpolator.h;h=c1f454288c055d5f1d3db56384fc77c45b49897d;hp=cea4961809cda3300db2c8b227d7639c3d1d734e;hb=8a7697ce57b56cdc43698cd1241ad98d49f9b5ac;hpb=3ea498d51ce4597e9649cd21f155b51175ea0bea diff --git a/include/interpolator.h b/include/interpolator.h index cea4961..c1f4542 100644 --- a/include/interpolator.h +++ b/include/interpolator.h @@ -2,7 +2,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: interpolator.h,v 1.2 2001/03/21 21:45:31 kevin Exp $ +** $Id$ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -34,6 +34,7 @@ public: double interpolate (double x); }; + class CubicPolyInterpolator { private: const double* const m_pdY; @@ -50,38 +51,134 @@ public: template class BilinearInterpolator { +private: + T** const m_ppMatrix; + const int m_nx; + const int m_ny; + +public: + BilinearInterpolator (T** ppMatrix, unsigned int nx, unsigned int ny) + : m_ppMatrix(ppMatrix), m_nx(nx), m_ny(ny) + {} + + T interpolate (double dXPos, double dYPos) +{ + int iFloorX = static_cast(floor(dXPos)); + int iFloorY = static_cast(floor (dYPos)); + double dXFrac = dXPos - iFloorX; + double dYFrac = dYPos - iFloorY; + + T result = 0; + + if (iFloorX < 0 || iFloorY < 0 || iFloorX > m_nx-1 || iFloorY > m_ny-1) + result = 0; + else if (iFloorX == m_nx - 1 && iFloorY == m_ny - 1) + result = static_cast(m_ppMatrix[m_nx-1][m_ny-1]); + else if (iFloorX == m_nx - 1) + result = static_cast(m_ppMatrix[iFloorX][iFloorY] + dYFrac * (m_ppMatrix[iFloorX][iFloorY+1] - m_ppMatrix[iFloorX][iFloorY])); + else if (iFloorY == m_ny - 1) + result = static_cast(m_ppMatrix[iFloorX][iFloorY] + dXFrac * (m_ppMatrix[iFloorX+1][iFloorY] - m_ppMatrix[iFloorX][iFloorY])); + else + result = static_cast + ((1 - dXFrac) * (1 - dYFrac) * m_ppMatrix[iFloorX][iFloorY] + + dXFrac * (1 - dYFrac) * m_ppMatrix[iFloorX+1][iFloorY] + + dYFrac * (1 - dXFrac) * m_ppMatrix[iFloorX][iFloorY+1] + + dXFrac * dYFrac * m_ppMatrix[iFloorX+1][iFloorY+1]); + + return result; +} + }; + + +template +class BilinearPolarInterpolator { +private: + T** const m_ppMatrix; + const int m_nAngle; + const int m_nPos; + int m_nCenterPos; + +public: + BilinearPolarInterpolator (T** ppMatrix, unsigned int nAngle, + unsigned int nPos) + : m_ppMatrix(ppMatrix), m_nAngle(nAngle), m_nPos(nPos) + { + if (m_nPos %2) + m_nCenterPos = (m_nPos - 1) / 2; + else + m_nCenterPos = m_nPos / 2; + } + + T interpolate (double dAngle, double dPos) +{ + int iFloorAngle = static_cast(floor(dAngle)); + int iFloorPos = static_cast(floor (dPos)); + double dAngleFrac = dAngle - iFloorAngle; + double dPosFrac = dPos - iFloorPos; + + T result = 0; + + if (iFloorAngle < -1 || iFloorPos < 0 || iFloorAngle > m_nAngle-1 || iFloorPos > m_nPos-1) + result = 0; + else if (iFloorAngle == -1 && iFloorPos == m_nPos-1) + result = static_cast(m_ppMatrix[0][m_nPos-1] + dAngleFrac * (m_ppMatrix[m_nAngle-1][iFloorPos] - m_ppMatrix[0][iFloorPos])); + else if (iFloorAngle == m_nAngle - 1 && iFloorPos == m_nPos-1) + result = static_cast(m_ppMatrix[m_nAngle-1][m_nPos-1] + dAngleFrac * (m_ppMatrix[0][iFloorPos] - m_ppMatrix[m_nAngle-1][iFloorPos])); + else if (iFloorPos == m_nPos - 1) + result = static_cast(m_ppMatrix[iFloorAngle][iFloorPos] + dAngleFrac * (m_ppMatrix[iFloorAngle+1][iFloorPos] - m_ppMatrix[iFloorAngle][iFloorPos])); + else { + if (iFloorAngle == m_nAngle-1) { + int iUpperAngle = 0; + int iLowerPos = (m_nPos-1) - iFloorPos; + int iUpperPos = (m_nPos-1) - (iFloorPos+1); + result = static_cast + ((1-dAngleFrac) * (1-dPosFrac) * m_ppMatrix[iFloorAngle][iFloorPos] + + dAngleFrac * (1-dPosFrac) * m_ppMatrix[iUpperAngle][iLowerPos] + + dPosFrac * (1-dAngleFrac) * m_ppMatrix[iFloorAngle][iFloorPos+1] + + dAngleFrac * dPosFrac * m_ppMatrix[iUpperAngle][iUpperPos]); + } else if (iFloorAngle == -1) { + int iLowerAngle = m_nAngle - 1; + int iLowerPos = (m_nPos-1) - iFloorPos; + int iUpperPos = (m_nPos-1) - (iFloorPos+1); + result = static_cast + ((1-dAngleFrac) * (1-dPosFrac) * m_ppMatrix[iLowerAngle][iLowerPos] + + dAngleFrac * (1-dPosFrac) * m_ppMatrix[iFloorAngle+1][iFloorPos] + + dPosFrac * (1-dAngleFrac) * m_ppMatrix[iLowerAngle][iUpperPos] + + dAngleFrac * dPosFrac * m_ppMatrix[iFloorAngle+1][iFloorPos+1]); + } else + result = static_cast + ((1-dAngleFrac) * (1-dPosFrac) * m_ppMatrix[iFloorAngle][iFloorPos] + + dAngleFrac * (1-dPosFrac) * m_ppMatrix[iFloorAngle+1][iFloorPos] + + dPosFrac * (1-dAngleFrac) * m_ppMatrix[iFloorAngle][iFloorPos+1] + + dAngleFrac * dPosFrac * m_ppMatrix[iFloorAngle+1][iFloorPos+1]); + } + return result; +} +}; + + +template +class BicubicPolyInterpolator { private: T** const m_ppMatrix; const unsigned int m_nx; const unsigned int m_ny; public: - BilinearInterpolator (T** ppMatrix, unsigned int nx, unsigned int ny) + BicubicPolyInterpolator (T** ppMatrix, unsigned int nx, unsigned int ny) : m_ppMatrix(ppMatrix), m_nx(nx), m_ny(ny) {} T interpolate (double dXPos, double dYPos) { - int iFloorX = floor (dXPos); - int iFloorY = floor (dYPos); - double dXFrac = dXPos - iFloorX; - double dYFrac = dYPos - iFloorY; + // int iFloorX = static_cast(floor (dXPos)); + // int iFloorY = static_cast(floor (dYPos)); + // double dXFrac = dXPos - iFloorX; + // double dYFrac = dYPos - iFloorY; T result = 0; - if (iFloorX < 0 || iFloorY < 0 || iFloorX > m_nx-1 || iFloorY > m_ny-1) - result = 0; - else if (iFloorX == m_nx - 1 && iFloorY == m_ny - 1) - result = m_ppMatrix[m_nx-1][m_ny-1]; - else if (iFloorX == m_nx - 1) - result = m_ppMatrix[iFloorX][iFloorY] + dYFrac * (m_ppMatrix[iFloorX][iFloorY+1] - m_ppMatrix[iFloorX][iFloorY]); - else if (iFloorY == m_ny - 1) - result = m_ppMatrix[iFloorX][iFloorY] + dXFrac * (m_ppMatrix[iFloorX+1][iFloorY] - m_ppMatrix[iFloorX][iFloorY]); - else - result = (1 - dXFrac) * (1 - dYFrac) * m_ppMatrix[iFloorX][iFloorY] + - dXFrac * (1 - dYFrac) * m_ppMatrix[iFloorX+1][iFloorY] + - dYFrac * (1 - dXFrac) * m_ppMatrix[iFloorX][iFloorY+1] + - dXFrac * dYFrac * m_ppMatrix[iFloorX+1][iFloorY+1]; + // Need to add code return result; } @@ -91,17 +188,18 @@ public: template class LinearInterpolator { private: - T* const m_pY; T* const m_pX; - const unsigned int m_n; + T* const m_pY; + const int m_n; + const bool m_bZeroOutsideRange; public: - LinearInterpolator (T* pY, unsigned int n) - : m_pY(pY), m_pX(0), m_n(n) + LinearInterpolator (T* pY, unsigned int n, bool bZeroOutside = true) + : m_pX(0), m_pY(pY), m_n(n), m_bZeroOutsideRange(bZeroOutside) {} - LinearInterpolator (T* pY, T* pX, unsigned int n) - : m_pY(pY), m_pX(pX), m_n(n) + LinearInterpolator (T* pX, T* pY, unsigned int n, bool bZeroOutside = true) + : m_pX(pX), m_pY(pY), m_n(n), m_bZeroOutsideRange(bZeroOutside) {} double interpolate (double dX, int* piLastFloor = NULL) @@ -111,14 +209,20 @@ public: if (! m_pX) { if (dX == 0) result = m_pY[0]; - else if (dX < 0) - result = 0; - else if (dX == m_n - 1) + else if (dX < 0) { + if (m_bZeroOutsideRange) + result = 0; + else + result = m_pY[0]; + } else if (dX == m_n - 1) result = m_pY[m_n-1]; - else if (dX > m_n - 1) - result = 0; - else { - int iFloor = floor(dX); + else if (dX > m_n - 1) { + if (m_bZeroOutsideRange) + result = 0; + else + result = m_pY[m_n - 1]; + } else { + int iFloor = static_cast(floor(dX)); result = m_pY[iFloor] + (m_pY[iFloor+1] - m_pY[iFloor]) * (dX - iFloor); } } else { @@ -136,13 +240,19 @@ public: } if (dX == m_pX[0]) result = m_pY[0]; - else if (dX < m_pX[0]) - result = 0; - else if (dX == m_pX[m_n-1]) + else if (dX < m_pX[0]) { + if (m_bZeroOutsideRange) + result = 0; + else + result = m_pY[0]; + } else if (dX == m_pX[m_n-1]) result = m_pY[m_n-1]; - else if (dX > m_pX[m_n-1]) - result = 0; - else { + else if (dX > m_pX[m_n - 1]) { + if (m_bZeroOutsideRange) + result = 0; + else + result = m_pY[m_n - 1]; + } else { if (iLower < 0 || iLower >= m_n) { sys_error (ERR_SEVERE, "Coordinate out of range [linearInterpolate]"); return 0;