Skip to content

Commit b913b68

Browse files
pasabanovpasabanov
authored
Added new quadratic spline types and unified calculations (#30)
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 b913b68

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)