1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
/************************************************************************/
/* */
/* vspline - a set of generic tools for creation and evaluation */
/* of uniform b-splines */
/* */
/* Copyright 2015 - 2023 by Kay F. Jahnke */
/* */
/* Permission is hereby granted, free of charge, to any person */
/* obtaining a copy of this software and associated documentation */
/* files (the "Software"), to deal in the Software without */
/* restriction, including without limitation the rights to use, */
/* copy, modify, merge, publish, distribute, sublicense, and/or */
/* sell copies of the Software, and to permit persons to whom the */
/* Software is furnished to do so, subject to the following */
/* conditions: */
/* */
/* The above copyright notice and this permission notice shall be */
/* included in all copies or substantial portions of the */
/* Software. */
/* */
/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
/* OTHER DEALINGS IN THE SOFTWARE. */
/* */
/************************************************************************/
/// \file splinus.cc
///
/// \brief compare a periodic b-spline with a sine
///
/// This is a simple example using a periodic b-spline
/// over just two values: 1 and -1. This spline is used to approximate
/// a sine function. You pass the spline's desired degree on the command
/// line. Next you enter a number (interpreted as degrees) and the program
/// will output the sine and the 'splinus' of the given angle.
/// As you can see when playing with higher degrees, the higher the spline's
/// degree, the closer the match with the sine. So apart from serving as a
/// very simple demonstration of using a 1D periodic b-spline, it teaches us
/// that a periodic b-spline can approximate a sine.
/// To show off, we use long double as the spline's data type.
/// This program also shows a bit of functional programming magic, putting
/// together the 'splinus' functor from several of vspline's functional
/// bulding blocks.
///
/// compile with: clang++ -pthread -O3 -std=c++11 splinus.cc -o splinus
#include <iostream>
#include <assert.h>
#include <vspline/vspline.h>
int main ( int argc , char * argv[] )
{
if ( argc < 2 )
{
std::cerr << "please pass the spline degree on the command line"
<< std::endl ;
exit ( 1 ) ;
}
int degree = std::atoi ( argv[1] ) ;
if ( degree < 4 )
{
std::cout << "rising degree to 4, the minimum for this program" << std::endl ;
degree = 4 ;
}
assert ( degree >= 4 && degree <= vspline_constants::max_degree ) ;
// create the bspline object
typedef vspline::bspline < long double , 1 > spline_type ;
spline_type bsp ( 2 , // two values
degree , // degree as per command line
vspline::PERIODIC , // periodic boundary conditions
0.0 ) ; // no tolerance
// the bspline object's 'core' is a MultiArrayView to the knot point
// data, which we set one by one for this simple example:
bsp.core[0] = 1.0L ;
bsp.core[1] = -1.0L ;
// now we prefilter the data
bsp.prefilter() ;
// we build 'splinus' as a functional construct. Inside the brace,
// we 'chain' several vspline::unary_functors:
// - a 'domain' which scales and shifts input to the spline's range.
// with the 'incoming' range of [ 90 , 450 ] and the spline's
// range of [ 0 , 2 ], the translation of incoming coordinates is:
// x' = 0 + 2 * ( x - 90 ) / ( 450 - 90 )
// - a 'safe' evaluator for the spline. since the spline has been
// built with PERIODIC boundary conditions, this evaluator will
// map incoming coordinates into the first period, [0,2].
auto splinus =
( vspline::domain ( 90.0L , 450.0L , bsp )
+ vspline::make_safe_evaluator < spline_type , long double > ( bsp ) ) ;
// alternatively we can use this construct. This will work just about
// the same, but has a potential flaw: If arithmetic imprecision should
// land the output of the periodic gate just slightly below 90.0, the
// domain may produce a value just below 0.0, resulting in a slightly
// out-of-bounds access. So the construct above is preferable.
// Just to demonstrate that vspline::grok also produces an object that
// can be used with function call syntax, we use vspline::grok here.
// Note how 'grokking' the chain of functors produces a simply typed
// object, rather than the complexly typed result of the chaining
// operation inside the brace:
vspline::grok_type < long double , long double > splinus2
= vspline::grok
( vspline::periodic ( 90.0L , 450.0L )
+ vspline::domain ( 90.0L , 450.0L , bsp )
+ vspline::evaluator < long double , long double > ( bsp ) ) ;
// we throw derivatives in the mix. If our spline models a sine,
// it's derivatives should model the sine's derivatives, cos etc.
// note how this is obscured by the higher 'steepness' of the spline,
// which is over [ 0 , 2 ] , not [ 0 , 2 * pi ]. Hence the derivatives
// come out amplified with a power of pi, which we compensate for.
vigra::TinyVector < int , 1 > derivative ;
derivative[0] = 1 ;
vspline::evaluator < long double , long double > evd1 ( bsp , derivative ) ;
derivative[0] = 2 ;
vspline::evaluator < long double , long double > evd2 ( bsp , derivative ) ;
derivative[0] = 3 ;
vspline::evaluator < long double , long double > evd3 ( bsp , derivative ) ;
auto splinus_d1
= ( vspline::domain ( 90.0L , 450.0L , bsp )
+ vspline::periodic ( 0.0L , 2.0L )
+ evd1 ) ;
auto splinus_d2
= ( vspline::domain ( 90.0L , 450.0L , bsp )
+ vspline::periodic ( 0.0L , 2.0L )
+ evd2 ) ;
auto splinus_d3
= ( vspline::domain ( 90.0L , 450.0L , bsp )
+ vspline::periodic ( 0.0L , 2.0L )
+ evd3 ) ;
// now we let the user enter arbitrary angles, calculate the sine
// and the 'splinus' and print the result and difference:
while ( std::cin.good() )
{
std::cout << " angle in degrees > " ;
long double x ;
std::cin >> x ; // get an angle
if ( std::cin.eof() )
{
std::cout << std::endl ;
break ;
}
long double xs = x * M_PI / 180.0L ; // note: sin() uses radians
// finally we can produce both results. Note how we can use periodic_ev,
// the combination of gate and evaluator, like an ordinary function.
std::cout << "sin(" << x << ") = "
<< sin ( xs ) << std::endl
<< "cos(" << x << ") = "
<< cos ( xs ) << std::endl
<< "splinus(" << x << ") = "
<< splinus ( x ) << std::endl
<< "splinus2(" << x << ") = "
<< splinus2 ( x ) << std::endl
<< "difference sin/splinus: "
<< sin ( xs ) - splinus ( x ) << std::endl
<< "difference sin/splinus2: "
<< sin ( xs ) - splinus2 ( x ) << std::endl
<< "difference splinus/splinus2: "
<< splinus2 ( x ) - splinus ( x ) << std::endl
<< "derivatives: " << splinus_d1 ( x ) / M_PI
<< " " << splinus_d2 ( x ) / ( M_PI * M_PI )
<< " " << splinus_d3 ( x ) / ( M_PI * M_PI * M_PI ) << std::endl ;
}
}
|
