1 /** @file homkermap.c
2 ** @brief Homogeneous kernel map - Definition
3 ** @author Andrea Vedaldi
4 **/
5
6 /*
7 Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
8 Copyright (C) 2013 Andrea Vedaldi.
9 All rights reserved.
10
11 This file is part of the VLFeat library and is made available under
12 the terms of the BSD license (see the COPYING file).
13 */
14
15 /** @file homkermap.h
16
17 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
18 @page homkermap Homogeneous kernel map
19 @author Andrea Vedaldi
20 @tableofcontents
21 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
22
23 @ref homkermap.h implements the homogeneous kernel maps introduced in
24 @cite{vedaldi10efficient},@cite{vedaldi12efficient}. Such maps are
25 efficient linear representations of popular kernels such as the
26 intersection, $\chi^2$, and Jensen-Shannon ones.
27
28 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
29 @section homkermap-starting Getting started
30 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
31
32 The homogeneous kernel map is implemented as an object of type
33 ::VlHomogeneousKernelMap. To use thois object, first create an
34 instance by using ::vl_homogeneouskernelmap_new, then use
35 ::vl_homogeneouskernelmap_evaluate_d or
36 ::vl_homogeneouskernelmap_evaluate_f (depdening on whether the data is
37 @c double or @c float) to compute the feature map $ \Psi(x)
38 $. When done, dispose of the object by calling
39 ::vl_homogeneouskernelmap_delete.
40
41 @code
42 double gamma = 1.0 ;
43 int order = 1 ;
44 double period = -1 ; // use default
45 double psi [3] ;
46 vl_size psiStride = 1 ;
47 double x = 0.5 ;
48 VlHomogeneousKernelMap * hom = vl_homogeneouskernelmap_new(
49 VlHomogeneousKernelChi2, gamma, order, period,
50 VlHomogeneousKernelMapWindowRectangular) ;
51 vl_homogeneouskernelmap_evaluate_d(hom, psi, psiStride, x) ;
52 vl_homogeneouskernelmap_delete(x) ;
53 @endcode
54
55 The constructor ::vl_homogeneouskernelmap_new takes the kernel type @c
56 kernel (see ::VlHomogeneousKernelType), the homogeneity order @c gamma
57 (use one for the standard $1$-homogeneous kernels), the approximation
58 order @c order (usually order one is enough), the period @a period
59 (use a negative value to use the default period), and a window type @c
60 window (use ::VlHomogeneousKernelMapWindowRectangular if unsure). The
61 approximation order trades off the quality and dimensionality of the
62 approximation. The resulting feature map $ \Psi(x) $, computed by
63 ::vl_homogeneouskernelmap_evaluate_d or
64 ::vl_homogeneouskernelmap_evaluate_f , is <code>2*order+1</code>
65 dimensional.
66
67 The code pre-computes the map $ \Psi(x) $ for efficient
68 evaluation. The table spans values of $ x $ in the range
69 $[2^{-20}, 2^{8}) $. In particular, values smaller than $
70 2^{-20} $ are treated as zeroes (which results in a null feature).
71
72 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
73 @section homkermap-fundamentals Fundamentals
74 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
75
76 The <em>homogeneous kernel map</em> is a finite dimensional linear
77 approximation of homogeneous kernels, including the intersection,
78 $\chi^2$, and Jensen-Shannon kernels. These kernels are frequently
79 used in computer vision applications because they are particular
80 suited to data in the format of histograms, which includes many common
81 visual descriptors.
82
83 Let $x,y \in \mathbb{R}_+$ be non-negative scalars and let $k(x,y) \in
84 \mathbb{R}$ be an homogeneous kernel such as the $\chi^2$ and or the
85 intersection ones:
86
87 @f[
88 k_{\mathrm{inters}}(x,y) = \min\{x, y\},
89 \quad
90 k_{\chi^2}(x,y) = 2 \frac{(x - y)^2}{x+y}.
91 @f]
92
93 For vectorial data $ \mathbf{x},\mathbf{y} \in \mathbb{R}_+^d $, the
94 homogeneous kernels is defined as an <em>additive combination</em> of
95 scalar kernels $K(\mathbf{x},\mathbf{y}) = \sum_{i=1}^d k(x_i,y_i)$.
96
97 The <em>homogeneous kernel map</em> of order $n$ is a vector function
98 $\Psi(x) \in \mathbb{R}^{2n+1}$ such that, for any choice of $x, y \in
99 \mathbb{R}_+$, the following approximation holds:
100
101 @f[
102 k(x,y) \approx \langle \Psi(x), \Psi(y) \rangle.
103 @f]
104
105 Given the feature map for the scalar case, the corresponding feature
106 map $\Psi(\mathbf{x})$ for the vectorial case is obtained by stacking
107 $[\Psi(x_1), \dots, \Psi(x_n)]$. Note that the stacked feature
108 $\Psi(\mathbf{x})$ has dimension $d(2n+1)$.
109
110 Using linear analysis tools (e.g. a linear support vector machine)
111 on top of dataset that has been encoded by the homogeneous kernel map
112 is therefore approximately equivalent to using a method based
113 on the corresponding non-linear kernel.
114
115 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
116 @subsection homkermap-overview-negative Extension to the negative reals
117 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
118
119 Any positive (semi-)definite kernel $k(x,y)$ defined on the
120 non-negative reals $x,y \in \mathbb{R}_+$ can be extended to the
121 entire real line by using the definition:
122
123 @f[
124 k_\pm(x,y) = \operatorname{sign}(x) \operatorname{sign}(y) k(|x|,|y|).
125 @f]
126
127 The homogeneous kernel map implements this extension by defining
128 $\Psi_\pm(x) = \operatorname{sign}(x) \Psi(|x|)$. Note that other
129 extensions are possible, such as
130
131 @f[
132 k_\pm(x,y) = H(xy) \operatorname{sign}(y) k(|x|,|y|)
133 @f]
134
135 where $H$ is the Heaviside function, but may result in higher
136 dimensional feature maps.
137
138 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
139 @subsection homkermap-overview-homogeneity Homogeneity degree
140 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
141
142 Any (1-)homogeneous kernel $k_1(x,y)$ can be extended to a so called
143 $\gamma$-homgeneous kernel $k_\gamma(x,y)$ by the definition
144
145 @f[
146 k_\gamma(x,y) = (xy)^{\frac{\gamma}{2}} \frac{k_1(x,y)}{\sqrt{xy}}
147 @f]
148
149 Smaller values of $\gamma$ enhance the kernel non-linearity and are
150 sometimes beneficial in applications (see
151 @cite{vedaldi10efficient},@cite{vedaldi12efficient} for details).
152
153 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
154 @subsection homkermap-overview-window Windowing and period
155 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
156
157 This section discusses aspects of the homogeneous kernel map which are
158 more technical and may be skipped. The homogeneous kernel map
159 approximation is based on periodizing the kernel; given the kernel
160 signature
161
162 @f[
163 \mathcal{K}(\lambda) = k(e^{\frac{\lambda}{2}}, e^{-\frac{\lambda}{2}})
164 @f]
165
166 the homogeneous kernel map is a feature map for the windowed and
167 periodized kernel whose signature is given by
168
169 @f[
170 \hat{\mathcal{K}}(\lambda)
171 =
172 \sum_{i=-\infty}^{+\infty} \mathcal{K}(\lambda + k \Lambda) W(\lambda + k \Lambda)
173 @f]
174
175 where $W(\lambda)$ is a windowing function and $\Lambda$ is the
176 period. This implementation of the homogeneous kernel map supports the
177 use of a <em>uniform window</em> ($ W(\lambda) = 1 $) or of a
178 <em>rectangular window</em> ($ W(\lambda) =
179 \operatorname{rect}(\lambda/\Lambda) $). Note that $ \lambda =
180 \log(y/x) $ is equal to the logarithmic ratio of the arguments of the
181 kernel. Empirically, the rectangular window seems to have a slight
182 edge in applications.
183
184 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
185 @section homkermap-details Implementation details
186 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
187
188 This implementation uses the expressions given in
189 @cite{vedaldi10efficient},@cite{vedaldi11efficient} to compute in
190 closed form the maps $\Psi(x)$ for the supported kernel types. For
191 efficiency reasons, it precomputes $\Psi(x)$ for a large range of
192 values of the argument when the homogeneous kernel map object is
193 created.
194
195 The internal table stores $\Psi(x) \in \mathbb{R}^{2n+1}$ by sampling
196 $x\geq 0$. This uses the internal decomposition of IEEE floating point
197 representations (@c float and @c double) in mantissa and exponent:
198 <pre>
199 x = mantissa * (2**exponent),
200 minExponent <= exponent <= maxExponent,
201 1 <= matnissa < 2.
202 </pre>
203 Each octave is further sampled in @c numSubdivisions sublevels.
204
205 When the map $\Psi(x)$ is evaluated, @c x is decomposed again into
206 exponent and mantissa to index the table. The output is obtained by
207 bilinear interpolation from the appropriate table entries.
208
209 **/
210
211 /* ---------------------------------------------------------------- */
212 #ifndef VL_HOMKERMAP_INSTANTIATING
213 /* ---------------------------------------------------------------- */
214
215 #include "homkermap.h"
216 #include "mathop.h"
217 #include <math.h>
218
219 struct _VlHomogeneousKernelMap
220 {
221 VlHomogeneousKernelType kernelType ;
222 double gamma ;
223 VlHomogeneousKernelMapWindowType windowType ;
224 vl_size order ;
225 double period ;
226 vl_size numSubdivisions ;
227 double subdivision ;
228 vl_index minExponent ;
229 vl_index maxExponent ;
230 double * table ;
231 } ;
232
233 /** @internal @brief Sample the kernel specturm
234 ** @param self homogeneous kernel map.
235 ** @param omega sampling frequency.
236 ** @return the spectrum sampled at @a omega.
237 **/
238
239 VL_INLINE double
vl_homogeneouskernelmap_get_spectrum(VlHomogeneousKernelMap const * self,double omega)240 vl_homogeneouskernelmap_get_spectrum (VlHomogeneousKernelMap const * self, double omega)
241 {
242 assert (self) ;
243 switch (self->kernelType) {
244 case VlHomogeneousKernelIntersection:
245 return (2.0 / VL_PI) / (1 + 4 * omega*omega) ;
246 case VlHomogeneousKernelChi2:
247 return 2.0 / (exp(VL_PI * omega) + exp(-VL_PI * omega)) ;
248 case VlHomogeneousKernelJS:
249 return (2.0 / log(4.0)) *
250 2.0 / (exp(VL_PI * omega) + exp(-VL_PI * omega)) /
251 (1 + 4 * omega*omega) ;
252 default:
253 abort() ;
254 }
255 }
256
257 /* helper */
sinc(double x)258 VL_INLINE double sinc(double x)
259 {
260 if (x == 0.0) return 1.0 ;
261 return sin(x) / x ;
262 }
263
264 /** @internal @brief Sample the smoothed kernel spectrum
265 ** @param self homogeneous kernel map.
266 ** @param omega sampling frequency.
267 ** @return the spectrum sampled at @a omega after smoothing.
268 **/
269
270 VL_INLINE double
vl_homogeneouskernelmap_get_smooth_spectrum(VlHomogeneousKernelMap const * self,double omega)271 vl_homogeneouskernelmap_get_smooth_spectrum (VlHomogeneousKernelMap const * self, double omega)
272 {
273 double kappa_hat = 0 ;
274 double omegap ;
275 double epsilon = 1e-2 ;
276 double const omegaRange = 2.0 / (self->period * epsilon) ;
277 double const domega = 2 * omegaRange / (2 * 1024.0 + 1) ;
278 assert (self) ;
279 switch (self->windowType) {
280 case VlHomogeneousKernelMapWindowUniform:
281 kappa_hat = vl_homogeneouskernelmap_get_spectrum(self, omega) ;
282 break ;
283 case VlHomogeneousKernelMapWindowRectangular:
284 for (omegap = - omegaRange ; omegap <= omegaRange ; omegap += domega) {
285 double win = sinc((self->period/2.0) * omegap) ;
286 win *= (self->period/(2.0*VL_PI)) ;
287 kappa_hat += win * vl_homogeneouskernelmap_get_spectrum(self, omegap + omega) ;
288 }
289 kappa_hat *= domega ;
290 /* project on the postivie orthant (see PAMI) */
291 kappa_hat = VL_MAX(kappa_hat, 0.0) ;
292 break ;
293 default:
294 abort() ;
295 }
296 return kappa_hat ;
297 }
298
299 /* ---------------------------------------------------------------- */
300 /* Constructors and destructors */
301 /* ---------------------------------------------------------------- */
302
303 /** @brief Create a new homgeneous kernel map
304 ** @param kernelType type of homogeneous kernel.
305 ** @param gamma kernel homogeneity degree.
306 ** @param order approximation order.
307 ** @param period kernel period.
308 ** @param windowType type of window used to truncate the kernel.
309 ** @return the new homogeneous kernel map.
310 **
311 ** The function intializes a new homogeneous kernel map for the
312 ** specified kernel type, homogeneity degree, approximation order,
313 ** period, and truncation window. See @ref homkermap-fundamentals for
314 ** details.
315 **
316 ** The homogeneity degree @c gamma must be positive (the standard
317 ** kernels are obtained by setting @c gamma to 1). When unsure, set
318 ** @c windowType to ::VlHomogeneousKernelMapWindowRectangular. The @c
319 ** period should be non-negative; specifying a negative or null value
320 ** causes the function to switch to a default value.
321 **
322 ** The function returns @c NULL if there is not enough free memory.
323 **/
324
325 VlHomogeneousKernelMap *
vl_homogeneouskernelmap_new(VlHomogeneousKernelType kernelType,double gamma,vl_size order,double period,VlHomogeneousKernelMapWindowType windowType)326 vl_homogeneouskernelmap_new (VlHomogeneousKernelType kernelType,
327 double gamma,
328 vl_size order,
329 double period,
330 VlHomogeneousKernelMapWindowType windowType)
331 {
332 int tableWidth, tableHeight ;
333 VlHomogeneousKernelMap * self = vl_malloc(sizeof(VlHomogeneousKernelMap)) ;
334 if (! self) return NULL ;
335
336 assert(gamma > 0) ;
337
338 assert(kernelType == VlHomogeneousKernelIntersection ||
339 kernelType == VlHomogeneousKernelChi2 ||
340 kernelType == VlHomogeneousKernelJS) ;
341
342 assert(windowType == VlHomogeneousKernelMapWindowUniform ||
343 windowType == VlHomogeneousKernelMapWindowRectangular) ;
344
345 if (period < 0) {
346 switch (windowType) {
347 case VlHomogeneousKernelMapWindowUniform:
348 switch (kernelType) {
349 case VlHomogeneousKernelChi2: period = 5.86 * sqrt(order + 0) + 3.65 ; break ;
350 case VlHomogeneousKernelJS: period = 6.64 * sqrt(order + 0) + 7.24 ; break ;
351 case VlHomogeneousKernelIntersection: period = 2.38 * log(order + 0.8) + 5.6 ; break ;
352 }
353 break ;
354 case VlHomogeneousKernelMapWindowRectangular:
355 switch (kernelType) {
356 case VlHomogeneousKernelChi2: period = 8.80 * sqrt(order + 4.44) - 12.6 ; break ;
357 case VlHomogeneousKernelJS: period = 9.63 * sqrt(order + 1.00) - 2.93; break ;
358 case VlHomogeneousKernelIntersection: period = 2.00 * log(order + 0.99) + 3.52 ; break ;
359 }
360 break ;
361 }
362 period = VL_MAX(period, 1.0) ;
363 }
364
365 self->kernelType = kernelType ;
366 self->windowType = windowType ;
367 self->gamma = gamma ;
368 self->order = order ;
369 self->period = period ;
370 self->numSubdivisions = 8 + 8*order ;
371 self->subdivision = 1.0 / self->numSubdivisions ;
372 self->minExponent = -20 ;
373 self->maxExponent = 8 ;
374
375 tableHeight = (int) (2*self->order + 1) ;
376 tableWidth = (int) (self->numSubdivisions * (self->maxExponent - self->minExponent + 1)) ;
377 self->table = vl_malloc (sizeof(double) *
378 (tableHeight * tableWidth + 2*(1+self->order))) ;
379 if (! self->table) {
380 vl_free(self) ;
381 return NULL ;
382 }
383
384 {
385 vl_index exponent ;
386 vl_uindex i, j ;
387 double * tablep = self->table ;
388 double * kappa = self->table + tableHeight * tableWidth ;
389 double * freq = kappa + (1+self->order) ;
390 double L = 2.0 * VL_PI / self->period ;
391
392 /* precompute the sampled periodicized spectrum */
393 j = 0 ;
394 i = 0 ;
395 while (i <= self->order) {
396 freq[i] = j ;
397 kappa[i] = vl_homogeneouskernelmap_get_smooth_spectrum(self, j * L) ;
398 ++ j ;
399 if (kappa[i] > 0 || j >= 3*i) ++ i ;
400 }
401
402 /* fill table */
403 for (exponent = self->minExponent ;
404 exponent <= self->maxExponent ; ++ exponent) {
405
406 double x, Lxgamma, Llogx, xgamma ;
407 double sqrt2kappaLxgamma ;
408 double mantissa = 1.0 ;
409
410 for (i = 0 ; i < self->numSubdivisions ;
411 ++i, mantissa += self->subdivision) {
412 x = ldexp(mantissa, (int)exponent) ;
413 xgamma = pow(x, self->gamma) ;
414 Lxgamma = L * xgamma ;
415 Llogx = L * log(x) ;
416
417 *tablep++ = sqrt(Lxgamma * kappa[0]) ;
418 for (j = 1 ; j <= self->order ; ++j) {
419 sqrt2kappaLxgamma = sqrt(2.0 * Lxgamma * kappa[j]) ;
420 *tablep++ = sqrt2kappaLxgamma * cos(freq[j] * Llogx) ;
421 *tablep++ = sqrt2kappaLxgamma * sin(freq[j] * Llogx) ;
422 }
423 } /* next mantissa */
424 } /* next exponent */
425 }
426 return self ;
427 }
428
429 /** @brief Delete an object instance.
430 ** @param self object.
431 ** The function deletes the specified map object.
432 **/
433
434 void
vl_homogeneouskernelmap_delete(VlHomogeneousKernelMap * self)435 vl_homogeneouskernelmap_delete (VlHomogeneousKernelMap * self)
436 {
437 vl_free(self->table) ;
438 self->table = NULL ;
439 vl_free(self) ;
440 }
441
442 /* ---------------------------------------------------------------- */
443 /* Retrieve data and parameters */
444 /* ---------------------------------------------------------------- */
445
446 /** @brief Get the map order.
447 ** @param self object.
448 ** @return the map order.
449 **/
450
451 vl_size
vl_homogeneouskernelmap_get_order(VlHomogeneousKernelMap const * self)452 vl_homogeneouskernelmap_get_order (VlHomogeneousKernelMap const * self)
453 {
454 assert(self) ;
455 return self->order ;
456 }
457
458 /** @brief Get the map dimension.
459 ** @param self object.
460 ** @return the map dimension (2 @c order +1).
461 **/
462
463 vl_size
vl_homogeneouskernelmap_get_dimension(VlHomogeneousKernelMap const * self)464 vl_homogeneouskernelmap_get_dimension (VlHomogeneousKernelMap const * self)
465 {
466 assert(self) ;
467 return 2 * self->order + 1 ;
468 }
469
470 /** @brief Get the kernel type.
471 ** @param self object.
472 ** @return kernel type.
473 **/
474
475 VlHomogeneousKernelType
vl_homogeneouskernelmap_get_kernel_type(VlHomogeneousKernelMap const * self)476 vl_homogeneouskernelmap_get_kernel_type (VlHomogeneousKernelMap const * self)
477 {
478 assert(self) ;
479 return self->kernelType ;
480 }
481
482 /** @brief Get the window type.
483 ** @param self object.
484 ** @return window type.
485 **/
486
487 VlHomogeneousKernelMapWindowType
vl_homogeneouskernelmap_get_window_type(VlHomogeneousKernelMap const * self)488 vl_homogeneouskernelmap_get_window_type (VlHomogeneousKernelMap const * self)
489 {
490 assert(self) ;
491 return self->windowType ;
492 }
493
494 /* ---------------------------------------------------------------- */
495 /* Process data */
496 /* ---------------------------------------------------------------- */
497
498 /** @fn ::vl_homogeneouskernelmap_evaluate_d(VlHomogeneousKernelMap const*,double*,vl_size,double)
499 ** @brief Evaluate map
500 ** @param self map object.
501 ** @param destination output buffer.
502 ** @param stride stride of the output buffer.
503 ** @param x value to expand.
504 **
505 ** The function evaluates the feature map on @a x and stores the
506 ** resulting <code>2*order+1</code> dimensional vector to
507 ** @a destination[0], @a destination[stride], @a destination[2*stride], ....
508 **/
509
510 /** @fn ::vl_homogeneouskernelmap_evaluate_f(VlHomogeneousKernelMap const*,float*,vl_size,double)
511 ** @copydetails ::vl_homogeneouskernelmap_evaluate_d(VlHomogeneousKernelMap const*,double*,vl_size,double)
512 **/
513
514 #define FLT VL_TYPE_FLOAT
515 #define VL_HOMKERMAP_INSTANTIATING
516 #include "homkermap.c"
517
518 #define FLT VL_TYPE_DOUBLE
519 #define VL_HOMKERMAP_INSTANTIATING
520 #include "homkermap.c"
521
522 /* VL_HOMKERMAP_INSTANTIATING */
523 #endif
524
525 /* ---------------------------------------------------------------- */
526 #ifdef VL_HOMKERMAP_INSTANTIATING
527 /* ---------------------------------------------------------------- */
528
529 #include "float.h"
530
531 void
VL_XCAT(vl_homogeneouskernelmap_evaluate_,SFX)532 VL_XCAT(vl_homogeneouskernelmap_evaluate_,SFX)
533 (VlHomogeneousKernelMap const * self,
534 T * destination,
535 vl_size stride,
536 double x)
537 {
538 /* break value into exponent and mantissa */
539 int exponent ;
540 int unsigned j ;
541 double mantissa = frexp(x, &exponent) ;
542 double sign = (mantissa >= 0.0) ? +1.0 : -1.0 ;
543 mantissa *= 2*sign ;
544 exponent -- ;
545
546 if (mantissa == 0 ||
547 exponent <= self->minExponent ||
548 exponent >= self->maxExponent) {
549 for (j = 0 ; j < 2*self->order+1 ; ++j) {
550 *destination = (T) 0.0 ;
551 destination += stride ;
552 }
553 return ;
554 }
555 {
556 vl_size featureDimension = 2*self->order + 1 ;
557 double const * v1 = self->table +
558 (exponent - self->minExponent) * self->numSubdivisions * featureDimension ;
559 double const * v2 ;
560 double f1, f2 ;
561
562 mantissa -= 1.0 ;
563 while (mantissa >= self->subdivision) {
564 mantissa -= self->subdivision ;
565 v1 += featureDimension ;
566 }
567 v2 = v1 + featureDimension ;
568 for (j = 0 ; j < featureDimension ; ++j) {
569 f1 = *v1++ ;
570 f2 = *v2++ ;
571 *destination = (T) sign * ((f2 - f1) * (self->numSubdivisions * mantissa) + f1) ;
572 destination += stride ;
573 }
574 }
575 }
576
577 #undef FLT
578 #undef VL_HOMKERMAP_INSTANTIATING
579 /* VL_HOMKERMAP_INSTANTIATING */
580 #endif
581