Skip to content

Commit e00dad9

Browse files
pasabanovpasabanov
committed
Added new quadratic spline types and unified calculations
1. Added "ClampedStart", "ClampedEnd", "SemiClamped", "FixedSecondStart", "FixedSecondEnd", "SemiFixedSecond", "Clamped", "FixedSecond", and "NotAKnot" quadratic spline types. 2. Added minimal documentation for these types. 3. Unified calculations in `QuadraticSpline::compute_coeffs`.
1 parent 61a30e3 commit e00dad9

File tree

3 files changed

+199
-94
lines changed

3 files changed

+199
-94
lines changed

‎ALFI/ALFI/spline/cubic.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ namespace alfi::spline {
7070
*/
7171
struct NotAKnotEnd final {};
7272
/**
73-
The arithmetic mean_array between NotAKnotStart and NotAKnotEnd.
73+
The arithmetic mean of NotAKnotStart and NotAKnotEnd.
7474
*/
7575
struct SemiNotAKnot final {};
7676
using Default = Natural;

‎ALFI/ALFI/spline/quadratic.h‎

Lines changed: 187 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ namespace alfi::spline {
2828
*/
2929
struct NotAKnotEnd final {};
3030
/**
31-
The arithmetic mean_array between NotAKnotStart and NotAKnotEnd.
31+
The arithmetic mean of NotAKnotStart and NotAKnotEnd.
3232
*/
3333
struct SemiNotAKnot final {};
3434
/**
@@ -40,13 +40,80 @@ namespace alfi::spline {
4040
*/
4141
struct NaturalEnd final {};
4242
/**
43-
The arithmetic mean_array between NaturalStart and NaturalEnd.
43+
The arithmetic mean of NaturalStart and NaturalEnd.
4444
*/
4545
struct SemiNatural final {};
4646
/**
47-
The arithmetic mean_array between SemiNotAKnot and SemiNatural.
47+
The arithmetic mean of SemiNotAKnot and SemiNatural.
4848
*/
4949
struct SemiSemi final {};
50+
/**
51+
The first derivative equals `df` at the first point.
52+
*/
53+
struct ClampedStart final {
54+
explicit ClampedStart(Number df) : df(df) {}
55+
Number df;
56+
};
57+
/**
58+
The first derivative equals `df` at the last point.
59+
*/
60+
struct ClampedEnd final {
61+
explicit ClampedEnd(Number df) : df(df) {}
62+
Number df;
63+
};
64+
/**
65+
The arithmetic mean of ClampedStart and ClampedEnd.
66+
*/
67+
struct SemiClamped final {
68+
SemiClamped(Number df_1, Number df_n) : df_1(df_1), df_n(df_n) {}
69+
Number df_1, df_n;
70+
};
71+
/**
72+
The second derivative equals `d2f` on the first segment.
73+
*/
74+
struct FixedSecondStart final {
75+
explicit FixedSecondStart(Number d2f) : d2f(d2f) {}
76+
Number d2f;
77+
};
78+
/**
79+
The second derivative equals `d2f` on the last segment.
80+
*/
81+
struct FixedSecondEnd final {
82+
explicit FixedSecondEnd(Number d2f) : d2f(d2f) {}
83+
Number d2f;
84+
};
85+
/**
86+
The arithmetic mean of FixedSecondStart and FixedSecondEnd.
87+
*/
88+
struct SemiFixedSecond final {
89+
SemiFixedSecond(Number d2f_1, Number d2f_n) : d2f_1(d2f_1), d2f_n(d2f_n) {}
90+
Number d2f_1, d2f_n;
91+
};
92+
/**
93+
The first derivative equals `df` at the point with index `point_idx`.
94+
*/
95+
struct Clamped final {
96+
Clamped(SizeT point_idx, Number df) : point_idx(point_idx), df(df) {}
97+
SizeT point_idx;
98+
Number df;
99+
};
100+
/**
101+
The second derivative equals `d2f` on the segment with index `segment_idx`.
102+
*/
103+
struct FixedSecond final {
104+
FixedSecond(SizeT segment_idx, Number d2f) : segment_idx(segment_idx), d2f(d2f) {}
105+
SizeT segment_idx;
106+
Number d2f;
107+
};
108+
/**
109+
The second derivative is continuous at the point with index `point_idx`.\n
110+
This also means, that the second derivative is the same on the segments with indices `point_idx-1` and `point_idx`.\n
111+
If number of points is less than or equal to 2, then `point_idx` is ignored.
112+
*/
113+
struct NotAKnot final {
114+
explicit NotAKnot(SizeT point_idx) : point_idx(point_idx) {}
115+
SizeT point_idx;
116+
};
50117
using Default = SemiNotAKnot;
51118
};
52119

@@ -56,7 +123,16 @@ namespace alfi::spline {
56123
typename Types::NaturalStart,
57124
typename Types::NaturalEnd,
58125
typename Types::SemiNatural,
59-
typename Types::SemiSemi>;
126+
typename Types::SemiSemi,
127+
typename Types::ClampedStart,
128+
typename Types::ClampedEnd,
129+
typename Types::SemiClamped,
130+
typename Types::FixedSecondStart,
131+
typename Types::FixedSecondEnd,
132+
typename Types::SemiFixedSecond,
133+
typename Types::Clamped,
134+
typename Types::FixedSecond,
135+
typename Types::NotAKnot>;
60136

61137
static Container<Number> compute_coeffs(const auto& X, const auto& Y, Type type = typename Types::Default{}) {
62138
if (X.size() != Y.size()) {
@@ -73,97 +149,115 @@ namespace alfi::spline {
73149
return util::spline::simple_spline<Number,Container>(X, Y, 2);
74150
}
75151

152+
if (std::holds_alternative<typename Types::NotAKnotStart>(type)) {
153+
return compute_coeffs(X, Y, typename Types::NotAKnot(1));
154+
} else if (std::holds_alternative<typename Types::NotAKnotEnd>(type)) {
155+
return compute_coeffs(X, Y, typename Types::NotAKnot(n - 2));
156+
} else if (std::holds_alternative<typename Types::SemiNotAKnot>(type)) {
157+
return util::arrays::mean(compute_coeffs(X, Y, typename Types::NotAKnot{1}), compute_coeffs(X, Y, typename Types::NotAKnot{n - 2}));
158+
} else if (std::holds_alternative<typename Types::NaturalStart>(type)) {
159+
return compute_coeffs(X, Y, typename Types::FixedSecond(0, 0));
160+
} else if (std::holds_alternative<typename Types::NaturalEnd>(type)) {
161+
return compute_coeffs(X, Y, typename Types::FixedSecond(n - 2, 0));
162+
} else if (std::holds_alternative<typename Types::SemiNatural>(type)) {
163+
return util::arrays::mean(compute_coeffs(X, Y, typename Types::NaturalStart{}), compute_coeffs(X, Y, typename Types::NaturalEnd{}));
164+
} else if (std::holds_alternative<typename Types::SemiSemi>(type)) {
165+
return util::arrays::mean(compute_coeffs(X, Y, typename Types::SemiNotAKnot{}), compute_coeffs(X, Y, typename Types::SemiNatural{}));
166+
} else if (const auto* cs = std::get_if<typename Types::ClampedStart>(&type)) {
167+
return compute_coeffs(X, Y, typename Types::Clamped(0, cs->df));
168+
} else if (const auto* ce = std::get_if<typename Types::ClampedEnd>(&type)) {
169+
return compute_coeffs(X, Y, typename Types::Clamped(n - 1, ce->df));
170+
} else if (const auto* sc = std::get_if<typename Types::SemiClamped>(&type)) {
171+
return util::arrays::mean(compute_coeffs(X, Y, typename Types::Clamped(0, sc->df_1)), compute_coeffs(X, Y, typename Types::Clamped(n - 1, sc->df_n)));
172+
} else if (const auto* fss = std::get_if<typename Types::FixedSecondStart>(&type)) {
173+
return compute_coeffs(X, Y, typename Types::FixedSecond(0, fss->d2f));
174+
} else if (const auto* fse = std::get_if<typename Types::FixedSecondEnd>(&type)) {
175+
return compute_coeffs(X, Y, typename Types::FixedSecond(n - 2, fse->d2f));
176+
} else if (const auto* sfs = std::get_if<typename Types::SemiFixedSecond>(&type)) {
177+
return util::arrays::mean(compute_coeffs(X, Y, typename Types::FixedSecond(0, sfs->d2f_1)), compute_coeffs(X, Y, typename Types::FixedSecond(n - 2, sfs->d2f_n)));
178+
}
179+
76180
Container<Number> coeffs(3 * (n - 1));
77181

78-
std::visit(util::misc::overload{
79-
[&](const typename Types::NotAKnotStart&) {
80-
if (n <= 3) {
81-
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 2);
82-
return;
83-
}
84-
{ // first three points
85-
const Number h1 = X[1] - X[0], h2 = X[2] - X[1];
86-
const Number d1 = (Y[1] - Y[0]) / h1, d2 = (Y[2] - Y[1]) / h2;
87-
coeffs[0] = (d2 - d1) / (h1 + h2); // a1
88-
coeffs[1] = d1 - h1 * coeffs[0]; // b1
89-
coeffs[2] = Y[0]; // c1
90-
coeffs[3] = coeffs[0]; // a2
91-
coeffs[4] = d1 + h1 * coeffs[0]; // b2
92-
coeffs[5] = Y[1]; // c2
93-
}
94-
for (SizeT i = 2; i < n - 1; ++i) {
95-
const Number hi = X[i+1] - X[i], hi_1 = X[i] - X[i-1];
96-
const Number d = 2 * coeffs[3*(i-1)] * hi_1 + coeffs[3*(i-1)+1];
97-
coeffs[3*i] = ((Y[i+1] - Y[i])/hi - d) / hi;
98-
coeffs[3*i+1] = d;
99-
coeffs[3*i+2] = Y[i];
100-
}
101-
},
102-
[&](const typename Types::NotAKnotEnd&) {
103-
if (n <= 3) {
104-
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 2);
105-
return;
106-
}
107-
{ // last three points
108-
const Number hn_2 = X[n-2] - X[n-3], hn_1 = X[n-1] - X[n-2];
109-
const Number dn_2 = (Y[n-2] - Y[n-3]) / hn_2, dn_1 = (Y[n-1] - Y[n-2]) / hn_1;
110-
coeffs[3*(n-3)] = (dn_1 - dn_2) / (hn_2 + hn_1); // an-2
111-
coeffs[3*(n-3)+1] = dn_2 - hn_2 * coeffs[3*(n-3)]; // bn-2
112-
coeffs[3*(n-3)+2] = Y[n-3]; // cn-2
113-
coeffs[3*(n-2)] = coeffs[3*(n-3)]; // an-1
114-
coeffs[3*(n-2)+1] = dn_2 + hn_2 * coeffs[3*(n-3)]; // bn-1
115-
coeffs[3*(n-2)+2] = Y[n-2]; // cn-1
116-
}
117-
for (SizeT iter = 0; iter <= n - 4; ++iter) {
118-
const SizeT i = n - 4 - iter;
119-
const Number hi = X[i+1] - X[i];
120-
coeffs[3*i] = (coeffs[3*(i+1)+1] - (Y[i+1] - Y[i])/hi) / hi;
121-
coeffs[3*i+1] = coeffs[3*(i+1)+1] - 2*coeffs[3*i]*hi;
122-
// another way: coeffs[3*i+1] = (Y[i+1] - Y[i])/hi - coeffs[3*i]*hi;
123-
coeffs[3*i+2] = Y[i];
124-
}
125-
},
126-
[&](const typename Types::SemiNotAKnot&) {
127-
coeffs = util::arrays::mean(compute_coeffs(X, Y, typename Types::NotAKnotStart{}), compute_coeffs(X, Y, typename Types::NotAKnotEnd{}));
128-
},
129-
[&](const typename Types::NaturalStart&) {
130-
if (n <= 2) {
131-
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 2);
132-
return;
133-
}
134-
coeffs[0] = 0; // a1 = 0
135-
coeffs[1] = (Y[1] - Y[0]) / (X[1] - X[0]); // b1
136-
coeffs[2] = Y[0]; // c1
137-
for (SizeT i = 1; i < n - 1; ++i) {
138-
coeffs[3*i+1] = 2 * (Y[i] - Y[i-1]) / (X[i] - X[i-1]) - coeffs[3*(i-1)+1];
139-
coeffs[3*i] = ((Y[i+1] - Y[i]) / (X[i+1] - X[i]) - coeffs[3*i+1]) / (X[i+1] - X[i]);
140-
// another way: coeffs[3*i] = (coeffs[3*(i+1)+1] - coeffs[3*i+1]) / 2 / (X[i+1] - X[i]);
141-
coeffs[3*i+2] = Y[i];
142-
}
143-
},
144-
[&](const typename Types::NaturalEnd&) {
145-
if (n <= 2) {
146-
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 2);
147-
return;
148-
}
149-
coeffs[3*(n-2)] = 0; // an-1 = 0
150-
coeffs[3*(n-2)+1] = (Y[n-1] - Y[n-2]) / (X[n-1] - X[n-2]); // bn-1
151-
coeffs[3*(n-2)+2] = Y[n-2]; // cn-1
152-
for (SizeT iter = 0; iter <= n - 3; ++iter) {
153-
const SizeT i = n - 3 - iter;
154-
coeffs[3*i+1] = 2 * (Y[i+1] - Y[i]) / (X[i+1] - X[i]) - coeffs[3*(i+1)+1];
155-
// another way: coeffs[3*i] = ((Y[i+1] - Y[i]) / (X[i+1] - X[i]) - coeffs[3*i+1]) / (X[i+1] - X[i]);
156-
coeffs[3*i] = (coeffs[3*(i+1)+1] - coeffs[3*i+1]) / 2 / (X[i+1] - X[i]);
157-
coeffs[3*i+2] = Y[i];
158-
}
159-
},
160-
[&](const typename Types::SemiNatural&) {
161-
coeffs = util::arrays::mean(compute_coeffs(X, Y, typename Types::NaturalStart{}), compute_coeffs(X, Y, typename Types::NaturalEnd{}));
162-
},
163-
[&](const typename Types::SemiSemi&) {
164-
coeffs = util::arrays::mean(compute_coeffs(X, Y, typename Types::SemiNotAKnot{}), compute_coeffs(X, Y, typename Types::SemiNatural{}));
165-
},
166-
}, type);
182+
SizeT i1 = 0, i2 = n - 1;
183+
184+
if (const auto* clamped = std::get_if<typename Types::Clamped>(&type)) {
185+
const auto i = clamped->point_idx;
186+
if (i < 0 || i >= n) {
187+
std::cerr << "Error in function " << __FUNCTION__
188+
<< ": point index for type 'Clamped' is out of bounds:"
189+
<< "expected to be non-negative and less than " << n << ", but got " << i
190+
<< ". Returning an empty array..." << std::endl;
191+
return {};
192+
}
193+
if (i > 0) {
194+
const auto dx1 = X[i] - X[i-1], dy1 = Y[i] - Y[i-1];
195+
coeffs[3*(i-1)+2] = Y[i-1];
196+
coeffs[3*(i-1)+1] = 2*dy1/dx1 - clamped->df;
197+
coeffs[3*(i-1)] = dy1/dx1/dx1 - coeffs[3*(i-1)+1]/dx1;
198+
i1 = i - 1;
199+
}
200+
if (i < n - 1) {
201+
const auto dx = X[i+1] - X[i], dy = Y[i+1] - Y[i];
202+
coeffs[3*i+2] = Y[i];
203+
coeffs[3*i+1] = clamped->df;
204+
coeffs[3*i] = dy/dx/dx - coeffs[3*i+1]/dx;
205+
i2 = i;
206+
}
207+
} else if (const auto* fixed_second = std::get_if<typename Types::FixedSecond>(&type)) {
208+
const auto i = fixed_second->segment_idx;
209+
if (i < 0 || i >= n - 1) {
210+
std::cerr << "Error in function " << __FUNCTION__
211+
<< ": point index for type 'FixedSecond' is out of bounds:"
212+
<< "expected to be non-negative and less than " << n - 1 << ", but got " << i
213+
<< ". Returning an empty array..." << std::endl;
214+
return {};
215+
}
216+
coeffs[3*i] = fixed_second->d2f / 2;
217+
coeffs[3*i+1] = (Y[i+1]-Y[i])/(X[i+1]-X[i]) - coeffs[3*i]*(X[i+1]-X[i]);
218+
coeffs[3*i+2] = Y[i];
219+
i1 = i2 = i;
220+
} else if (const auto* not_a_knot = std::get_if<typename Types::NotAKnot>(&type)) {
221+
if (n <= 2) {
222+
return util::spline::simple_spline<Number,Container>(X, Y, 2);
223+
}
224+
const auto i = not_a_knot->point_idx;
225+
if (i < 1 || i >= n - 1) {
226+
std::cerr << "Error in function " << __FUNCTION__
227+
<< ": point index for type 'NotAKnot' is out of bounds:"
228+
<< "expected to be positive and less than " << n - 1 << ", but got " << i
229+
<< ". Returning an empty array..." << std::endl;
230+
return {};
231+
}
232+
const auto dx1 = X[i] - X[i-1], dy1 = Y[i] - Y[i-1];
233+
const auto dx = X[i+1] - X[i], dy = Y[i+1] - Y[i];
234+
coeffs[3*(i-1)+2] = Y[i-1];
235+
coeffs[3*i+2] = Y[i];
236+
coeffs[3*(i-1)] = coeffs[3*i] = (dy/dx - dy1/dx1) / (dx1 + dx);
237+
coeffs[3*(i-1)+1] = dy1/dx1 - coeffs[3*(i-1)]*dx1;
238+
coeffs[3*i+1] = dy/dx - coeffs[3*i]*dx;
239+
i1 = i - 1;
240+
i2 = i;
241+
} else {
242+
std::cerr << "Error in function " << __FUNCTION__ << ": Unknown type. This should not have happened."
243+
<< " Please report this to the developers. Returning an empty array..." << std::endl;
244+
return {};
245+
}
246+
247+
for (SizeT iter = 0; iter < i1; ++iter) {
248+
const auto i = i1 - 1 - iter;
249+
const auto dx = X[i+1] - X[i], dy = Y[i+1] - Y[i];
250+
coeffs[3*i] = coeffs[3*(i+1)+1]/dx - dy/dx/dx;
251+
coeffs[3*i+1] = 2*dy/dx - coeffs[3*(i+1)+1];
252+
coeffs[3*i+2] = Y[i];
253+
}
254+
255+
for (SizeT i = i2 + 1; i < n - 1; ++i) {
256+
const auto dx = X[i+1] - X[i], dy = Y[i+1] - Y[i];
257+
coeffs[3*i+2] = Y[i];
258+
coeffs[3*i+1] = coeffs[3*(i-1)+1] + 2*coeffs[3*(i-1)]*(X[i]-X[i-1]);
259+
coeffs[3*i] = dy/dx/dx - coeffs[3*i+1]/dx;
260+
}
167261

168262
return coeffs;
169263
}

‎examples/interpolation/interpolation.cpp‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,17 @@ class PlotWindow final : public QWidget {
5454
{"Natural end", alfi::spline::QuadraticSpline<>::Types::NaturalEnd{}},
5555
{"Semi-natural", alfi::spline::QuadraticSpline<>::Types::SemiNatural{}},
5656
{"Semi-semi", alfi::spline::QuadraticSpline<>::Types::SemiSemi{}},
57+
{"Clamped-start(10)", alfi::spline::QuadraticSpline<>::Types::ClampedStart{10}},
58+
{"Clamped-start(-10)", alfi::spline::QuadraticSpline<>::Types::ClampedStart{-10}},
59+
{"Clamped-end(10)", alfi::spline::QuadraticSpline<>::Types::ClampedEnd{10}},
60+
{"Clamped-end(-10)", alfi::spline::QuadraticSpline<>::Types::ClampedEnd{-10}},
61+
{"Semi-clamped(10, 10)", alfi::spline::QuadraticSpline<>::Types::SemiClamped{10, 10}},
62+
{"Fixed-second-start(10)", alfi::spline::QuadraticSpline<>::Types::FixedSecondStart{10}},
63+
{"Fixed-second-end(10)", alfi::spline::QuadraticSpline<>::Types::FixedSecondEnd{10}},
64+
{"Semi-fixed-second(10, 10)", alfi::spline::QuadraticSpline<>::Types::SemiFixedSecond{10, 10}},
65+
{"Clamped(10, 10)", alfi::spline::QuadraticSpline<>::Types::Clamped{10, 10}},
66+
{"Fixed-second(10, 10)", alfi::spline::QuadraticSpline<>::Types::FixedSecond{10, 10}},
67+
{"Not-a-knot(10)", alfi::spline::QuadraticSpline<>::Types::NotAKnot{10}},
5768
};
5869
static const inline size_t quadratic_spline_default_type_index = 2;
5970

0 commit comments

Comments
 (0)