1 /*
2 * Copyright 2017 Google Inc.
3 *
4 * Use of this source code is governed by a BSD-style license that can be
5 * found in the LICENSE file.
6 */
7
8
9 #include "include/core/SkTypes.h"
10 #include "include/private/SkFloatingPoint.h"
11 #include "src/core/SkGaussFilter.h"
12 #include <cmath>
13
14 // The value when we can stop expanding the filter. The spec implies that 3% is acceptable, but
15 // we just use 1%.
16 static constexpr double kGoodEnough = 1.0 / 100.0;
17
18 // Normalize the values of gauss to 1.0, and make sure they add to one.
19 // NB if n == 1, then this will force gauss[0] == 1.
normalize(int n,double * gauss)20 static void normalize(int n, double* gauss) {
21 // Carefully add from smallest to largest to calculate the normalizing sum.
22 double sum = 0;
23 for (int i = n-1; i >= 1; i--) {
24 sum += 2 * gauss[i];
25 }
26 sum += gauss[0];
27
28 // Normalize gauss.
29 for (int i = 0; i < n; i++) {
30 gauss[i] /= sum;
31 }
32
33 // The factors should sum to 1. Take any remaining slop, and add it to gauss[0]. Add the
34 // values in such a way to maintain the most accuracy.
35 sum = 0;
36 for (int i = n - 1; i >= 1; i--) {
37 sum += 2 * gauss[i];
38 }
39
40 gauss[0] = 1 - sum;
41 }
42
calculate_bessel_factors(double sigma,double * gauss)43 static int calculate_bessel_factors(double sigma, double *gauss) {
44 auto var = sigma * sigma;
45
46 // The two functions below come from the equations in "Handbook of Mathematical Functions"
47 // by Abramowitz and Stegun. Specifically, equation 9.6.10 on page 375. Bessel0 is given
48 // explicitly as 9.6.12
49 // BesselI_0 for 0 <= sigma < 2.
50 // NB the k = 0 factor is just sum = 1.0.
51 auto besselI_0 = [](double t) -> double {
52 auto tSquaredOver4 = t * t / 4.0;
53 auto sum = 1.0;
54 auto factor = 1.0;
55 auto k = 1;
56 // Use a variable number of loops. When sigma is small, this only requires 3-4 loops, but
57 // when sigma is near 2, it could require 10 loops. The same holds for BesselI_1.
58 while(factor > 1.0/1000000.0) {
59 factor *= tSquaredOver4 / (k * k);
60 sum += factor;
61 k += 1;
62 }
63 return sum;
64 };
65 // BesselI_1 for 0 <= sigma < 2.
66 auto besselI_1 = [](double t) -> double {
67 auto tSquaredOver4 = t * t / 4.0;
68 auto sum = t / 2.0;
69 auto factor = sum;
70 auto k = 1;
71 while (factor > 1.0/1000000.0) {
72 factor *= tSquaredOver4 / (k * (k + 1));
73 sum += factor;
74 k += 1;
75 }
76 return sum;
77 };
78
79 // The following formula for calculating the Gaussian kernel is from
80 // "Scale-Space for Discrete Signals" by Tony Lindeberg.
81 // gauss(n; var) = besselI_n(var) / (e^var)
82 auto d = std::exp(var);
83 double b[SkGaussFilter::kGaussArrayMax] = {besselI_0(var), besselI_1(var)};
84 gauss[0] = b[0]/d;
85 gauss[1] = b[1]/d;
86
87 // The code below is tricky, and written to mirror the recursive equations from the book.
88 // The maximum spread for sigma == 2 is guass[4], but in order to know to stop guass[5]
89 // is calculated. At this point n == 5 meaning that gauss[0..4] are the factors, but a 6th
90 // element was used to calculate them.
91 int n = 1;
92 // The recurrence relation below is from "Numerical Recipes" 3rd Edition.
93 // Equation 6.5.16 p.282
94 while (gauss[n] > kGoodEnough) {
95 b[n+1] = -(2*n/var) * b[n] + b[n-1];
96 gauss[n+1] = b[n+1] / d;
97 n += 1;
98 }
99
100 normalize(n, gauss);
101
102 return n;
103 }
104
SkGaussFilter(double sigma)105 SkGaussFilter::SkGaussFilter(double sigma) {
106 SkASSERT(0 <= sigma && sigma < 2);
107
108 fN = calculate_bessel_factors(sigma, fBasis);
109 }
110