/usr/include/vspline/brace.h is in vspline-dev 0.3.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 | /************************************************************************/
/* */
/* vspline - a set of generic tools for creation and evaluation */
/* of uniform b-splines */
/* */
/* Copyright 2015 - 2017 by Kay F. Jahnke */
/* */
/* The git repository for this software is at */
/* */
/* https://bitbucket.org/kfj/vspline */
/* */
/* Please direct questions, bug reports, and contributions to */
/* */
/* kfjahnke+vspline@gmail.com */
/* */
/* 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 brace.h
\brief This file provides code for 'bracing' the spline's coefficient array.
Note that this isn't really user code, it's code used by class vspline::bspline.
Inspired by libeinspline, I wrote code to 'brace' the spline coefficients. The concept is
this: while the IIR filter used to calculate the coefficients has infinite support (though
arithmetic precision limits this in real-world applications), the evaluation of the spline
at a specific location only looks at a small window of coefficients (compact, finite support).
This fact can be exploited by taking note of how large the support area is and providing
a few more coefficients in a frame around the 'core' coefficients to allow the evaluation
to proceed without having to check for boundary conditions. While the difference is not
excessive (the main computational cost is the actual evaluation itself), it's still
nice to be able to code the evaluation without boundary checking, which makes the code
very straightforward and legible.
There is another aspect to bracing: In my implementation of vectorized evaluation,
the window into the coefficient array used to pick out coefficients to evaluate at
a specific location is coded as a set of offsets from it's 'low' corner. This way,
several such windows can be processed in parallel. This mechanism can only function
efficiently in a braced coefficient array, since it would otherwise have to give up
if any of the windows accessed by the vector of coordinates had members outside the
(unbraced) coefficient array and submit the coordinate vector to individual processing.
I consider the logic to code this and the loss in performance too much of a bother
to go down this path; all my evaluation code uses braced coefficient arrays. Of course
the user is free to omit bracing, but then they have to use their own evaluation
code.
What's in the brace? Of course this depends on the boundary conditions chosen.
In vspline, I offer code for several boundary conditions, but most have something
in common: the original, finite sequence is extrapolated into an infinite periodic
signal. With straight PERIODIC boundary conditions, the initial sequence is
immediately followed and preceded by copies of itself. The other boundary conditions
mirror the signal in some way and then repeat the mirrored signal periodically.
Using boundary conditions like these, both the extrapolated signal and the
coefficients share the same periodicity and mirroring.
There are two ways of arriving at a braced coeffcient array: We can start from the
extrapolated signal, pick a section large enough to make margin effects vanish
(due to limited arithmetic precision), prefilter it and pick out a subsection containing
the 'core' coefficients and their support. Alternatively, we can work only on the core
coefficients, calculate suitable initial causal and anticausal coeffcients (where the
calculation considers the extrapolated signal, which remains implicit), apply the filter
and *then* surround the core coefficient array with more coeffcients (the brace) following
the same extrapolation pattern as we imposed on the signal, but now on the coefficients
rather than on the initial knot point values.
The bracing can be performed without any solver-related maths by simply copying
(possibly trivially modified) slices of the core coefficients to the margin area.
Following the 'implicit' scheme, my default modus operandi braces after the
prefiltering. Doing so, it is posible to calculate the inital causal and anticausal
coefficient for the prefilter exactly. But this exactness is still, eventually,
subject to discretization and can only be represented after quantization. If,
instead, we prefilter a suitably extrapolated signal, now with arbitrary boundary
conditions, the margin effects will vanish towards the center (due to the characteristics
of the filter), and the 'core' coefficients will end up the same as in the first
approach. So we might as well extrapolate the signal 'far enough', pick any boundary
conditions we like (even zero padding), prefilter, and discard the margin outside the
area which is unaffected by margin effects. The result is, within arithmetic precision,
the same. Both approaches have advantages and disadvantages:
Implicit extrapolation needs less memory - we only need to provide storage for the
core coeffcients, which is just as much as we need for the original signal, so we can
operate in-place. The disadvantage of the implicit scheme is that we have to capture
the implicit extrapolation in code to calculate the initial causal/anticausal coefficients,
which is non-trivial and requires separate routines for each case, as can be seen in my
prefiltering code. And if, after prefiltering, we want to brace the core coeffcients
for efficient evaluation, we still need additional memory, which, if it hasn't been
allocated around the core before prefiltering, even requires us to copy the data out
into a larger memory area.
Explicit extrapolation needs more memory. A typical scheme would be to anticipate the
space needed for the explicit extrapolation, allocate enough memory for the extrapolated
signal, place the same into the center of the allocated memory, perform the extrapolation
and then prefilter. The advantage is that we can run the prefilter with arbitrary initial
causal/anticausal coefficients. No matter what the extrapolation looks like, we can always
use the same code. And we can extrapolate in any way we see fit, without having to produce
code to deal with our choice. If we pick the frame of extrapolated values large enough,
we can even pick out the 'braced' coefficient array from the result of the filter.
Obviously, there is no one 'right' way of doing this. Offered several choices
of implicit extrapolation, the user can choose between the implicit and explicit scheme.
The code in this file is useful for both choices: for the implicit scheme, bracing is
applied after prefiltering to enable evaluation with vspline. For the explicit scheme,
bracing may be used on the original data before prefiltering with arbitrary boundary
conditions, if the user's extrapolation scheme is covered by the code given here.
When using the higher-level access methods (via bspline objects), using the explicit or
implicit scheme becomes a matter of passing in the right flag, so at this level, a deep
understanding of the extrapolation mechanism isn't needed at all. I use the implicit scheme
as the default.
Since the bracing mainly requires copying data or trivial maths we can do the operations
on higher-dimensional objects, like slices of a volume. To efficiently code these operations
we make use of vigra's multi-math facility and it's bindAt array method, which makes
these subarrays easily available.
TODO: while this is convenient, it's not too fast, as it's neither multithreaded nor
vectorized. Still in most 'normal' scenarios the execution time is negligible...
TODO: there are 'pathological' cases where one brace is larger than the other brace
and the width of the core together. These cases can't be handled for all bracing modes
and will result in an exception.
*/
#ifndef VSPLINE_BRACE_H
#define VSPLINE_BRACE_H
#include <vigra/multi_array.hxx>
#include <vigra/multi_iterator.hxx>
#include "common.h"
namespace vspline {
/// class bracer encodes the entire bracing process. Note that contrary
/// to my initial implementation, class bracer is now used exclusively
/// for populating the frame around a core area of data. It has no code
/// to determine which size a brace/frame should have. This is now
/// determined in class bspline, see especially class bspline's methods
/// get_left_brace_size(), get_right_brace_size() and setup_metrics().
template < class view_type >
struct bracer
{
typedef typename view_type::value_type value_type ;
enum { dimension = view_type::actual_dimension } ;
/// for spherical images, we require special treatment for two-dimensional
/// input data, because we need to shift the values by 180 degrees, or half
/// the margin's width. But to compile, we also have to give a procedure
/// for the other cases (not 2D), so this is first:
template < typename value_type >
void shift_assign ( value_type target , value_type source )
{
// should not ever get used, really...
}
/// specialized routine for the 2D case (the slice itself is 1D)
template < typename value_type >
void shift_assign ( MultiArrayView < 1 , value_type > target ,
MultiArrayView < 1 , value_type > source )
{
// bit sloppy here, with pathological data (very small source.size()) this will
// be imprecise for odd sizes, for even sizes it's always fine. But then full
// sphericals always have size 2N * N, so odd sizes should not occur at all for dim 0
// TODO what if it's not a full spherical?
target = source ;
return ;
auto si = source.begin() + source.size() / 2 ;
auto se = source.end() ;
for ( auto& ti : target )
{
ti = *si ;
++si ;
if ( si >= se )
si = source.begin() ;
}
}
/// apply the bracing to the array, performing the required copy/arithmetic operations
/// to the 'frame' around the core. This routine performs the operation along axis dim.
/// This is also the routine to be used for explicitly extrapolating a signal:
/// you place the data into the center of a larger array, and pass in the sizes of the
/// 'empty' space which is to be filled with the extrapolated data.
///
/// the bracing is done one-left-one-right, to avoid corner cases as best as posible.
void apply ( view_type & a , // containing array
bc_code bc , // boundary condition code
int lsz , // space to the left which needs to be filled
int rsz , // ditto, to the right
int axis ) // axis along which to apply bracing
{
int w = a.shape ( axis ) ; // width of containing array along axis 'axis'
int m = w - ( lsz + rsz ) ; // width of 'core' array
if ( m < 1 ) // has to be at least 1
throw shape_mismatch ( "combined brace sizes must be at least one less than container size" ) ;
if ( ( lsz > m + rsz )
|| ( rsz > m + lsz ) )
{
// not enough data to fill brace
if ( bc == PERIODIC || bc == NATURAL || bc == MIRROR || bc == REFLECT )
throw std::out_of_range ( "each brace must be smaller than the sum of it's opposite brace and the core's width" ) ;
}
int l0 = lsz - 1 ; // index of innermost empty slice on the left; like begin()
int r0 = lsz + m ; // ditto, on the right
int lp = l0 + 1 ; // index of leftmost occupied slice (p for pivot)
int rp = r0 - 1 ; // index of rightmost occupied slice
int l1 = -1 ; // index one before outermost empty slice to the left
int r1 = w ; // index one after outermost empty slice on the right; like end()
int lt = l0 ; // index to left target slice
int rt = r0 ; // index to right target slice ;
int ls , rs ; // indices to left and right source slice, will be set below
int ds = 1 ; // step for source index, +1 == forẃard, used for all mirroring modes
// for periodic bracing, it's set to -1.
switch ( bc )
{
case PERIODIC :
{
ls = l0 + m ;
rs = r0 - m ;
ds = -1 ; // step through source in reverse direction
break ;
}
case NATURAL :
case MIRROR :
{
ls = l0 + 2 ;
rs = r0 - 2 ;
break ;
}
case CONSTANT :
case SPHERICAL :
case REFLECT :
{
ls = l0 + 1 ;
rs = r0 - 1 ;
break ;
}
case ZEROPAD :
{
break ;
}
case IDENTITY :
{
// these modes perform no bracing, return prematurely
return ;
}
default:
{
cerr << "bracing for BC code " << bc_name[bc] << " is not supported" << endl ;
break ;
}
}
for ( int i = max ( lsz , rsz ) ; i > 0 ; --i )
{
if ( lt > l1 )
{
switch ( bc )
{
case PERIODIC :
case MIRROR :
case REFLECT :
{
// with these three bracing modes, we simply copy from source to target
a.bindAt ( axis , lt ) = a.bindAt ( axis , ls ) ;
break ;
}
case NATURAL :
{
// here, we subtract the source slice from twice the 'pivot'
// easiest would be:
// a.bindAt ( axis , lt ) = a.bindAt ( axis , lp ) * value_type(2) - a.bindAt ( axis , ls ) ;
// but this fails in 1D TODO: why?
auto target = a.bindAt ( axis , lt ) ; // get a view to left target slice
target = a.bindAt ( axis , lp ) ; // assign value of left pivot slice
target *= value_type(2) ; // double that
target -= a.bindAt ( axis , ls ) ; // subtract left source slice
break ;
}
case CONSTANT :
{
// here, we repeat the 'pivot' slice
a.bindAt ( axis , lt ) = a.bindAt ( axis , lp ) ;
break ;
}
case ZEROPAD :
{
// fill with 0
a.bindAt ( axis , lt ) = value_type() ;
break ;
}
case SPHERICAL : // needs special treatment
{
shift_assign ( a.bindAt ( axis , lt ) , a.bindAt ( axis , ls ) ) ;
break ;
}
default :
// default: leave untouched
break ;
}
--lt ;
ls += ds ;
}
if ( rt < r1 )
{
// essentially the same, but with rs instead of ls, etc.
switch ( bc )
{
case PERIODIC :
case MIRROR :
case REFLECT :
{
// with these three bracing modes, we simply copy from source to target
a.bindAt ( axis , rt ) = a.bindAt ( axis , rs ) ;
break ;
}
case NATURAL :
{
// here, we subtract the source slice from twice the 'pivot'
// the easiest would be:
// a.bindAt ( axis , rt ) = a.bindAt ( axis , rp ) * value_type(2) - a.bindAt ( axis , rs ) ;
// but this fails in 1D TODO: why?
auto target = a.bindAt ( axis , rt ) ; // get a view to right targte slice
target = a.bindAt ( axis , rp ) ; // assign value of pivot slice
target *= value_type(2) ; // double that
target -= a.bindAt ( axis , rs ) ; // subtract source slice
break ;
}
case CONSTANT :
{
// here, we repeat the 'pivot' slice
a.bindAt ( axis , rt ) = a.bindAt ( axis , rp ) ;
break ;
}
case ZEROPAD :
{
// fill with 0
a.bindAt ( axis , rt ) = value_type() ;
break ;
}
case SPHERICAL : // needs special treatment
{
shift_assign ( a.bindAt ( axis , rt ) , a.bindAt ( axis , rs ) ) ;
break ;
}
default :
// default: leave untouched
break ;
}
++rt ;
rs -= ds ;
}
}
}
/// This variant of apply braces along all axes in one go.
static void apply
( view_type& a , ///< target array, containing the core and (empty) frame
vigra::TinyVector < bc_code , dimension > bcv , ///< boundary condition codes
vigra::TinyVector < int , dimension > left_corner , ///< sizes of left braces
vigra::TinyVector < int , dimension > right_corner ) ///< sizes of right braces
{
for ( int dim = 0 ; dim < dimension ; dim++ )
apply ( a , bcv[dim] , left_corner[dim] , right_corner[dim] , dim ) ;
}
} ;
} ; // end of namespace vspline
#endif // VSPLINE_BRACE_H
|