Engauge Digitizer 2
Loading...
Searching...
No Matches
Public Member Functions | Friends | List of all members
Spline Class Reference

Cubic interpolation given independent and dependent value vectors. More...

#include <Spline.h>

Collaboration diagram for Spline:
Collaboration graph

Public Member Functions

 Spline (const std::vector< double > &t, const std::vector< SplinePair > &xy, SplineTCheck splineTCheck=SPLINE_ENABLE_T_CHECK)
 Initialize spline with independent (t) and dependent (x and y) value vectors.
 
virtual ~Spline ()
 
void computeUntranslatedCoefficients (double aTranslated, double bTranslated, double cTranslated, double dTranslated, double tI, double &aUntranslated, double &bUntranslated, double &cUntranslated, double &dUntranslated) const
 From coefficients in xy=d*(t-ti)^3+c*(t-ti)^2+b*(t-ti)+a we compute and return the coefficients in xy=d* t ^3+c* t ^2+b* t +a.
 
SplinePair findSplinePairForFunctionX (double x, int numIterations) const
 Use bisection algorithm to iteratively find the SplinePair interpolated to best match the specified x value.
 
SplinePair interpolateCoeff (double t) const
 Return interpolated y for specified x.
 
SplinePair interpolateControlPoints (double t) const
 Return interpolated y for specified x, for testing.
 
SplinePair p1 (unsigned int i) const
 Bezier p1 control point for specified interval. P0 is m_xy[i] and P3 is m_xy[i+1].
 
SplinePair p2 (unsigned int i) const
 Bezier p2 control point for specified interval. P0 is m_xy[i] and P3 is m_xy[i+1].
 

Friends

class TestSpline
 For unit testing.
 

Detailed Description

Cubic interpolation given independent and dependent value vectors.

X is handled as a dependent variable based on the unitless independent parameter t so curves are not restricted to x(i)!=x(i+1).

This class has two modes that can be run side by side:

  1. Coefficient mode, where each interval has coefficients a,b,c,d in xy=ai+bi*(t-ti)+ci*(t-ti)^2+di*(t-ti)^3
  2. Bezier mode, where each interval has 2 endpoints and 2 control points in P=(1-s)^3*P0+3*(1-s)^2*s*P1+3*(1-s)*s^2*P2+s^3*P3. Although interpolation can be performed in bezier mode, this interpolation was just to verify consistent operation between the two modes. The real purpose of this mode is to produce reliable control points, p1 and p2, for each interval. The control points can be used by external code that relies on control points to perform its own interpolation

Definition at line 29 of file Spline.h.

Constructor & Destructor Documentation

◆ Spline()

Spline::Spline ( const std::vector< double > & t,
const std::vector< SplinePair > & xy,
SplineTCheck splineTCheck = SPLINE_ENABLE_T_CHECK )

Initialize spline with independent (t) and dependent (x and y) value vectors.

Besides initializing the a,b,c,d coefficients for each interval, this constructor initializes bezier points (P1 and P2) for each interval, where P0 and P3 are the start and end points for each interval.

Definition at line 13 of file Spline.cpp.

16{
17 ENGAUGE_ASSERT (t.size() == xy.size());
18 ENGAUGE_ASSERT (xy.size() > 0); // Need at least one point for this class to not fail with a crash
19
21 // In normal production this check should be performed
22 checkTIncrements (t);
23 }
24 computeCoefficientsForIntervals (t, xy);
25 computeControlPointsForIntervals ();
26}
const int INNER_RADIUS_MIN
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT if defined(QT_NO_DEBUG) && !defined(QT_FORCE_ASSERTS) define ENGAUGE...
@ SPLINE_ENABLE_T_CHECK
Definition Spline.h:15

◆ ~Spline()

Spline::~Spline ( )
virtual

Definition at line 28 of file Spline.cpp.

29{
30}

Member Function Documentation

◆ computeUntranslatedCoefficients()

void Spline::computeUntranslatedCoefficients ( double aTranslated,
double bTranslated,
double cTranslated,
double dTranslated,
double tI,
double & aUntranslated,
double & bUntranslated,
double & cUntranslated,
double & dUntranslated ) const

From coefficients in xy=d*(t-ti)^3+c*(t-ti)^2+b*(t-ti)+a we compute and return the coefficients in xy=d* t ^3+c* t ^2+b* t +a.

Definition at line 145 of file Spline.cpp.

154{
155 // Expanding the equation using t translated as (t-ti) we get the equation using untranslated (t) as follows
156 // xy=d*t^3+
157 // (c-3*d*ti)*t^2+
158 // (b-2*c*ti+3*d*ti^2)*t+
159 // (a-b*ti+c*ti^2-d*ti^3)
160 // which matches up with
161 // xy=d*t^3+
162 // c*t^2+
163 // b*t+
164 // a
165 // Both equations are given in the header file documentation for this method
167 bUntranslated = bTranslated - 2.0 * cTranslated * tI + 3.0 * dTranslated * tI * tI;
170}

◆ findSplinePairForFunctionX()

SplinePair Spline::findSplinePairForFunctionX ( double x,
int numIterations ) const

Use bisection algorithm to iteratively find the SplinePair interpolated to best match the specified x value.

This assumes the curve is a function since otherwise there is the potential for multiple solutions

Definition at line 172 of file Spline.cpp.

174{
176
177 double tLow = m_t[0];
178 double tHigh = m_t[m_xy.size() - 1];
179
180 // This method implicitly assumes that the x values are monotonically increasing - except for the
181 // "export raw points" exception right here
182
183 // Subtle stuff here - we first look for exact matches in m_elements. If a match is found, we exit
184 // immediately. This strategy works for "export raw points" option with functions having overlaps,
185 // in which case users expect interpolation to be used (issue #303)
186 for (unsigned int i = 0; i < m_xy.size(); i++) {
187 const SplinePair &sp = m_xy.at (i);
188 if (sp.x() == x) {
189 return sp;
190 }
191 }
192
193 // Extrapolation that is performed if x is out of bounds. As a starting point, we assume that the t
194 // values and x values behave the same, which is linearly. This assumption works best when user
195 // has set the points so the spline line is linear at the endpoints - which is also preferred since
196 // higher order polynomials are typically unstable and can "explode" off into unwanted directions
197 double x0 = interpolateCoeff (m_t[0]).x();
198 double xNm1 = interpolateCoeff (m_t[m_xy.size() - 1]).x();
199 if (x < x0) {
200
201 // Extrapolate with x < x(0) < x(N-1) which correspond to s, s0 and sNm1
202 double x1 = interpolateCoeff (m_t[1]).x();
203 double tStart = (x - x0) / (x1 - x0); // This will be less than zero. x=x0 for t=0 and x=x1 for t=1
204 tLow = 2.0 * tStart;
205 tHigh = 0.0;
206
207 } else if (xNm1 < x) {
208
209 // Extrapolate with x(0) < x(N-1) < x which correspond to s0, sNm1 and s
210 double xNm2 = interpolateCoeff (m_t[m_xy.size() - 2]).x();
211 double tStart = tHigh + (x - xNm1) / (xNm1 - xNm2); // This is greater than one. x=xNm2 for t=0 and x=xNm1 for t=1
212 tLow = m_xy.size() - 1;
213 tHigh = tHigh + 2.0 * (tStart - tLow);
214
215 }
216
217 // Interpolation using bisection search
218 double tCurrent = (tHigh + tLow) / 2.0;
219 double tDelta = (tHigh - tLow) / 4.0;
220 for (int iteration = 0; iteration < numIterations; iteration++) {
222 if (spCurrent.x() > x) {
223 tCurrent -= tDelta;
224 } else {
225 tCurrent += tDelta;
226 }
227 tDelta /= 2.0;
228 }
229
230 return spCurrent;
231}
Single X/Y pair for cubic spline interpolation initialization and calculations.
Definition SplinePair.h:14
double x() const
Get method for x.
SplinePair interpolateCoeff(double t) const
Return interpolated y for specified x.
Definition Spline.cpp:233

◆ interpolateCoeff()

SplinePair Spline::interpolateCoeff ( double t) const

Return interpolated y for specified x.

The appropriate interval is selected from the entire set of piecewise-defined intervals, then the corresponding a,b,c,d coefficients are applied

Definition at line 233 of file Spline.cpp.

234{
235 ENGAUGE_ASSERT (m_elements.size() != 0);
236
237 vector<SplineCoeff>::const_iterator itr;
238 itr = lower_bound(m_elements.begin(), m_elements.end(), t);
239 if (itr != m_elements.begin()) {
240 itr--;
241 }
242
243 return itr->eval(t);
244}

◆ interpolateControlPoints()

SplinePair Spline::interpolateControlPoints ( double t) const

Return interpolated y for specified x, for testing.

This uses the bezier points. If the t values are not separated by +1 consistently then this algorithm will probably need additional effort to work right

Definition at line 246 of file Spline.cpp.

247{
248 ENGAUGE_ASSERT (m_xy.size() != 0);
249
250 for (unsigned int i = 0; i < unsigned (m_xy.size() - 1); i++) {
251
252 if (m_t[i] <= t && t <= m_t[i+1]) {
253
254 SplinePair s ((t - m_t[i]) / (m_t[i + 1] - m_t[i]));
255 SplinePair onems (SplinePair (1.0) - s);
256 SplinePair xy = onems * onems * onems * m_xy [i] +
257 SplinePair (3.0) * onems * onems * s * m_p1 [i] +
258 SplinePair (3.0) * onems * s * s * m_p2 [i] +
259 s * s * s * m_xy[i + 1];
260 return xy;
261 }
262 }
263
264 // Input argument is out of bounds
265 ENGAUGE_ASSERT (false);
266 return m_t[0];
267}

◆ p1()

SplinePair Spline::p1 ( unsigned int i) const

Bezier p1 control point for specified interval. P0 is m_xy[i] and P3 is m_xy[i+1].

Definition at line 269 of file Spline.cpp.

270{
271 ENGAUGE_ASSERT (i < m_p1.size ());
272
273 return m_p1 [i];
274}

◆ p2()

SplinePair Spline::p2 ( unsigned int i) const

Bezier p2 control point for specified interval. P0 is m_xy[i] and P3 is m_xy[i+1].

Definition at line 276 of file Spline.cpp.

277{
278 ENGAUGE_ASSERT (i < m_p2.size ());
279
280 return m_p2 [i];
281}

Friends And Related Symbol Documentation

◆ TestSpline

For unit testing.

Definition at line 32 of file Spline.h.


The documentation for this class was generated from the following files: