Main Lemma Repository
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

CubicSplineInterpolator.h 5.3KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 09/25/2013 08:20:14 AM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@xri-geo.com
  14. * @copyright Copyright (c) 2013, XRI Geophysics, LLC
  15. * @copyright Copyright (c) 2013, Trevor Irons
  16. */
  17. #ifndef CUBICSPLINEINTERPOLATOR_INC
  18. #define CUBICSPLINEINTERPOLATOR_INC
  19. #include "LemmaObject.h"
  20. namespace Lemma {
  21. // Simple struct to hold spline terms
  22. struct SplineSet{
  23. VectorXr a;
  24. VectorXr b;
  25. VectorXr c;
  26. VectorXr d;
  27. VectorXr x;
  28. SplineSet( ) {
  29. }
  30. SplineSet(const int&n) {
  31. a = VectorXr::Zero(n+1);
  32. b = VectorXr::Zero(n);
  33. c = VectorXr::Zero(n+1);
  34. d = VectorXr::Zero(n);
  35. x = VectorXr::Zero(n+1);
  36. }
  37. };
  38. /**
  39. \brief Real 1D Natural cubic spline interpolator.
  40. \details Splines are fit between knots \f$j\f$ according to the forulae
  41. \f[ S_j(x) = a_j + b_j(x - x_j) + c_j(x-x_j)^2 + d_j(x-x_y)^3 \f]
  42. The spline must satisfy the following conditions
  43. \f{eqnarray} {
  44. S_i(x_i) & = & y_i = S_{i-1}(x_i), i = 1,..., n-1 \\
  45. S'_i(x_i) & = & S'_{i-1}(x_i), i = 1,..., n-1 \\
  46. S''_i(x_i) & = & S''_{i-1}(x_i), i = 1,..., n-1 \\
  47. S''_0(x_0) & = & S''_{n-1}(x_n) = 0
  48. \f}
  49. */
  50. class CubicSplineInterpolator : public LemmaObject {
  51. friend std::ostream &operator<<(std::ostream &stream,
  52. const CubicSplineInterpolator &ob);
  53. public:
  54. // ==================== LIFECYCLE =======================
  55. /**
  56. * Factory method for generating concrete class.
  57. * @return a std::shared_ptr of type CubicSplineInterpolator
  58. */
  59. static std::shared_ptr<CubicSplineInterpolator> NewSP();
  60. // ==================== OPERATORS =======================
  61. // ==================== OPERATIONS =======================
  62. /** Sets the knots to use for interpolation.
  63. @param[in] x are the absissa values
  64. @param[in] y are the ordinate values
  65. */
  66. void SetKnots(const VectorXr& x, const VectorXr& y);
  67. /** Resets the knots to use for interpolation, when abscissa values haven't changed.
  68. @param[in] y are the ordinate values
  69. */
  70. void ResetKnotOrdinate( const VectorXr& y );
  71. /** Interpolate a monotonically increasing ordered set.
  72. @param[in] x are the interpolation abscissa points
  73. @return the ordinate values at x
  74. */
  75. VectorXr InterpolateOrderedSet(const VectorXr& x);
  76. /** integrates the spline from x0 to x1. Uses composite Simpson's rule and n is the number of segments
  77. * @param[in] x0 is left argument
  78. * @param[in] x1 is right argument
  79. * @param[in] n is the number of points, must be even
  80. */
  81. Real Integrate(const Real& x0, const Real& x1, const int& n);
  82. /** integrates using cubic spline values. Taken from AMRIRA P223F project code Leroi, which in turn was based on
  83. This is a modification of the FUNCTION PPVALU in the book
  84. "A PRACTICAL GUIDE TO SPLINES" by C. DE BOOR
  85. */
  86. Real Integrate(const Real& x0, const Real& x1);
  87. /** @returns the know abscissa values
  88. */
  89. VectorXr GetKnotAbscissa();
  90. /** @returns the know abscissa values
  91. */
  92. VectorXr GetKnotOrdinate();
  93. /** Interpolation at a single point.
  94. @param[in] x is the interpolation abscissa point
  95. @param[in] i is an optional index to start searching at. Defaults to zero
  96. @return the ordinate value at x
  97. */
  98. Real Interpolate(const Real& x, int& i);
  99. /** Interpolation at a single point.
  100. @param[in] x is the interpolation abscissa point
  101. @return the ordinate value at x
  102. */
  103. Real Interpolate(const Real& x);
  104. // ==================== ACCESS =======================
  105. // ==================== INQUIRY =======================
  106. /** Returns the name of the underlying class, similiar to Python's type */
  107. virtual inline std::string GetName() const {
  108. return CName;
  109. }
  110. protected:
  111. // ==================== LIFECYCLE =======================
  112. /** Default protected constructor, use New */
  113. CubicSplineInterpolator ( );
  114. /** Default protected destructor, smart pointers auto delete */
  115. ~CubicSplineInterpolator ();
  116. /**
  117. * @copybrief LemmaObject::Release()
  118. * @copydetails LemmaObject::Release()
  119. */
  120. void Release();
  121. // ==================== OPERATIONS =======================
  122. /**
  123. Performs spline fitting.
  124. */
  125. //void Spline() {
  126. // }
  127. /** Finds the interval of knots in spline to use for integration.
  128. */
  129. int Interval(const Real& x);
  130. private:
  131. /** ASCII string representation of the class name */
  132. static constexpr auto CName = "CubicSplineInterpolator";
  133. SplineSet Spline;
  134. int ilo;
  135. int mflag;
  136. // ==================== DATA MEMBERS =========================
  137. }; // ----- end of class CubicSplineInterpolator -----
  138. } // ----- end of Lemma name -----
  139. #endif // ----- #ifndef CUBICSPLINEINTERPOLATOR_INC -----