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
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
|
/************************************************************************/
/* */
/* 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 verify.cc
\brief verify bspline interpolation against polynomial
A b-spline is a piecewise polynomial function. Therefore, it should
model a polynomial of the same degree precisely. This program tests
this assumption.
While the test should hold throughout, we have to limit it to
'reasonable' degrees, because we have to create the spline over
a range sufficiently large to make margin errors disappear,
so even if we want to, say, only look at the spline's values
between 0 and 1, we have to have a few ten or even 100 values
to the left and right of this interval, because the polynomial
does not exhibit convenient features like periodicity or
mirroring on the bounds. But since a polynomial, outside
[-1,1], grows with the power of it's degree, the values around
the test interval become very large very quickly with high
degrees. We can't reasonable expect to calculate a meaningful
spline over such data. The test shows how the measured fidelity
degrades with higher degrees due to this effect.
still , with 'reasonable' degrees, we see that the spline fits the
signal very well indeed, demonstrating that the spline can
faithfully represent a polynomial of equal degree.
*/
#include <iostream>
#include <typeinfo>
#include <random>
#include <vigra/multi_array.hxx>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>
#include <vspline/vspline.h>
template < class dtype >
struct random_polynomial
{
int degree ;
std::vector < dtype > coefficient ;
// set up a polynomial with random coefficients in [0,1]
random_polynomial ( int _degree )
: degree ( _degree ) ,
coefficient ( _degree + 1 )
{
std::random_device rd ;
std::mt19937 gen ( rd() ) ;
std::uniform_real_distribution<> dis ( -1.0 , 1.0 ) ;
for ( auto & e : coefficient )
e = dis ( gen ) ;
}
// evaluate the polynomial at x
dtype operator() ( dtype x )
{
dtype power = 1 ;
dtype result = 0 ;
for ( auto e : coefficient )
{
result += e * power ;
power *= x ;
}
return result ;
}
} ;
template < class dtype >
void polynominal_test ( int degree , const char * type_name )
{
// this is the function we want to model:
random_polynomial < long double > rp ( degree ) ;
// we evaluate this function in the range [-200,200[
// note that for high degrees, the signal will grow very
// large outside [-1,1], 'spoiling' the test
vigra::MultiArray < 1 , dtype > data ( vigra::Shape1 ( 400 ) ) ;
for ( int i = 0 ; i < 400 ; i++ )
{
dtype x = ( i - 200 ) ;
data[i] = rp ( x ) ;
}
// we create a b-spline over these data
vspline::bspline < dtype , 1 > bspl ( 400 , degree , vspline::NATURAL , 0.0 ) ;
bspl.prefilter ( data ) ;
auto ev = vspline::make_evaluator < decltype(bspl), dtype > ( bspl ) ;
// to test the spline against the polynomial, we generate random
// arguments in [-2,2]
std::random_device rd ;
std::mt19937 gen ( rd() ) ;
std::uniform_real_distribution<long double> dis ( -2.0 , 2.0 ) ;
long double signal = 0 ;
long double spline = 0 ;
long double noise = 0 ;
long double error ;
long double max_error = 0 ;
// now we evaluate the spline and the polynomial at equal arguments
// and do the statistics
for ( int i = 0 ; i < 10000 ; i++ )
{
long double x = dis ( gen ) ;
long double p = rp ( dtype ( x ) ) ;
// note how we have to translate x to spline coordinates here
long double s = ev ( dtype ( x + 200 ) ) ;
error = std::fabs ( p - s ) ;
if ( error > max_error )
max_error = error ;
signal += std::fabs ( p ) ;
spline += std::fabs ( s ) ;
noise += error ;
}
// note: with optimized code, this does not work:
if ( std::isnan ( noise ) || std::isnan ( noise ) )
{
std::cout << type_name << " aborted due to numeric overflow" << std::endl ;
return ;
}
long double mean_error = noise / 10000.0L ;
// finally we echo the results of the test
std::cout << type_name
<< " Mean error: "
<< mean_error
<< " Maximum error: "
<< max_error
<< " SNR "
<< int ( 20 * std::log10 ( signal / noise ) )
<< "dB"
<< std::endl ;
}
int main ( int argc , char * argv[] )
{
std::cout << std::fixed << std::showpos << std::showpoint
<< std::setprecision(18) ;
for ( int degree = 1 ; degree < 15 ; degree++ )
{
std::cout << "testing spline against polynomial, degree " << degree << std::endl ;
polynominal_test < float > ( degree , "using float........" ) ;
polynominal_test < double > ( degree , "using double......." ) ;
polynominal_test < long double > ( degree , "using long double.." ) ;
}
}
|
