aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--gn/core.gni2
-rw-r--r--gn/tests.gni1
-rw-r--r--src/core/SkGaussFilter.cpp152
-rw-r--r--src/core/SkGaussFilter.h41
-rw-r--r--tests/SkGaussFilterTest.cpp88
5 files changed, 284 insertions, 0 deletions
diff --git a/gn/core.gni b/gn/core.gni
index 21dcff4f44..358c3875c8 100644
--- a/gn/core.gni
+++ b/gn/core.gni
@@ -134,6 +134,8 @@ skia_core_sources = [
"$_src/core/SkFindAndPlaceGlyph.h",
"$_src/core/SkArenaAlloc.cpp",
"$_src/core/SkArenaAlloc.h",
+ "$_src/core/SkGaussFilter.cpp",
+ "$_src/core/SkGaussFilter.h",
"$_src/core/SkFlattenable.cpp",
"$_src/core/SkFlattenableSerialization.cpp",
"$_src/core/SkFont.cpp",
diff --git a/gn/tests.gni b/gn/tests.gni
index b051326d8d..3f8283ddad 100644
--- a/gn/tests.gni
+++ b/gn/tests.gni
@@ -212,6 +212,7 @@ tests_sources = [
"$_tests/SkColor4fTest.cpp",
"$_tests/SkDOMTest.cpp",
"$_tests/SkFixed15Test.cpp",
+ "$_tests/SkGaussFilterTest.cpp",
"$_tests/SkImageTest.cpp",
"$_tests/SkLiteDLTest.cpp",
"$_tests/SkNxTest.cpp",
diff --git a/src/core/SkGaussFilter.cpp b/src/core/SkGaussFilter.cpp
new file mode 100644
index 0000000000..548ff4398d
--- /dev/null
+++ b/src/core/SkGaussFilter.cpp
@@ -0,0 +1,152 @@
+/*
+ * Copyright 2017 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+
+#include "SkGaussFilter.h"
+
+#include <cmath>
+#include "SkTypes.h"
+
+static constexpr double kPi = 3.14159265358979323846264338327950288;
+
+// The value when we can stop expanding the filter. The spec implies that 3% is acceptable, but
+// we just use 1%.
+static constexpr double kGoodEnough = 1.0 / 100.0;
+
+// Normalize the values of gauss to 1.0, and make sure they add to one.
+// NB if n == 1, then this will force gauss[0] == 1.
+static void normalize(int n, double* gauss) {
+ // Carefully add from smallest to largest to calculate the normalizing sum.
+ double sum = 0;
+ for (int i = n-1; i >= 1; i--) {
+ sum += 2 * gauss[i];
+ }
+ sum += gauss[0];
+
+ // Normalize gauss.
+ for (int i = 0; i < n; i++) {
+ gauss[i] /= sum;
+ }
+
+ // The factors should sum to 1. Take any remaining slop, and add it to gauss[0]. Add the
+ // values in such a way to maintain the most accuracy.
+ sum = 0;
+ for (int i = n - 1; i >= 1; i--) {
+ sum += 2 * gauss[i];
+ }
+
+ gauss[0] = 1 - sum;
+}
+
+static int calculate_bessel_factors(double sigma, double *gauss) {
+ auto var = sigma * sigma;
+
+ // The two functions below come from the equations in "Handbook of Mathematical Functions"
+ // by Abramowitz and Stegun. Specifically, equation 9.6.10 on page 375. Bessel0 is given
+ // explicitly as 9.6.12
+ // BesselI_0 for 0 <= sigma < 2.
+ // NB the k = 0 factor is just sum = 1.0.
+ auto besselI_0 = [](double t) -> double {
+ auto tSquaredOver4 = t * t / 4.0;
+ auto sum = 1.0;
+ auto factor = 1.0;
+ auto k = 1;
+ // Use a variable number of loops. When sigma is small, this only requires 3-4 loops, but
+ // when sigma is near 2, it could require 10 loops. The same holds for BesselI_1.
+ while(factor > 1.0/1000000.0) {
+ factor *= tSquaredOver4 / (k * k);
+ sum += factor;
+ k += 1;
+ }
+ return sum;
+ };
+ // BesselI_1 for 0 <= sigma < 2.
+ auto besselI_1 = [](double t) -> double {
+ auto tSquaredOver4 = t * t / 4.0;
+ auto sum = t / 2.0;
+ auto factor = sum;
+ auto k = 1;
+ while (factor > 1.0/1000000.0) {
+ factor *= tSquaredOver4 / (k * (k + 1));
+ sum += factor;
+ k += 1;
+ }
+ return sum;
+ };
+
+ // The following formula for calculating the Gaussian kernel is from
+ // "Scale-Space for Discrete Signals" by Tony Lindeberg.
+ // gauss(n; var) = besselI_n(var) / (e^var)
+ auto d = std::exp(var);
+ double b[6] = {besselI_0(var), besselI_1(var)};
+ gauss[0] = b[0]/d;
+ gauss[1] = b[1]/d;
+
+ int n = 1;
+ // The recurrence relation below is from "Numerical Recipes" 3rd Edition.
+ // Equation 6.5.16 p.282
+ while (gauss[n] > kGoodEnough) {
+ b[n+1] = -(2*n/var) * b[n] + b[n-1];
+ gauss[n+1] = b[n+1] / d;
+ n += 1;
+ }
+
+ normalize(n, gauss);
+
+ return n;
+}
+
+static int calculate_gauss_factors(double sigma, double* gauss) {
+ SkASSERT(0 <= sigma && sigma < 2);
+
+ // From the SVG blur spec: 8.13 Filter primitive <feGaussianBlur>.
+ // H(x) = exp(-x^2/ (2s^2)) / sqrt(2π * s^2)
+ auto var = sigma * sigma;
+ auto expGaussDenom = -2 * var;
+ auto normalizeDenom = std::sqrt(2 * kPi) * sigma;
+
+ // Use the recursion relation from "Incremental Computation of the Gaussian" by Ken
+ // Turkowski in GPUGems 3. Page 877.
+ double g0 = 1.0 / normalizeDenom;
+ double g1 = std::exp(1.0 / expGaussDenom);
+ double g2 = g1 * g1;
+
+ gauss[0] = g0;
+ g0 *= g1;
+ g1 *= g2;
+ gauss[1] = g0;
+
+ int n = 1;
+ while (gauss[n] > kGoodEnough) {
+ g0 *= g1;
+ g1 *= g2;
+ gauss[n+1] = g0;
+ n += 1;
+ }
+
+ normalize(n, gauss);
+
+ return n;
+}
+
+SkGaussFilter::SkGaussFilter(double sigma, Type type) {
+ SkASSERT(0 <= sigma && sigma < 2);
+
+ if (type == Type::Bessel) {
+ fN = calculate_bessel_factors(sigma, fBasis);
+ } else {
+ fN = calculate_gauss_factors(sigma, fBasis);
+ }
+}
+
+int SkGaussFilter::filterDouble(double* values) const {
+ for (int i = 0; i < fN; i++) {
+ values[i] = fBasis[i];
+ }
+ return fN;
+}
+
diff --git a/src/core/SkGaussFilter.h b/src/core/SkGaussFilter.h
new file mode 100644
index 0000000000..9af45c875b
--- /dev/null
+++ b/src/core/SkGaussFilter.h
@@ -0,0 +1,41 @@
+/*
+ * Copyright 2017 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#ifndef SkGaussFilter_DEFINED
+#define SkGaussFilter_DEFINED
+
+#include <cstdint>
+
+// Define gaussian filters for values of sigma < 2. Produce values good to 1 part in 1,000,000.
+// Gaussian produces values as defined in the SVG 1.1 spec:
+// https://www.w3.org/TR/SVG/filters.html#feGaussianBlurElement
+// Bessel produces values as defined in "Scale-Space for Discrete Signals" by Tony Lindeberg
+class SkGaussFilter {
+public:
+ enum class Type : bool {
+ Gaussian,
+ Bessel
+ };
+
+ // Type selects which method is used to calculate the gaussian factors.
+ SkGaussFilter(double sigma, Type type);
+
+ int radius() const { return fN - 1; }
+ int width() const { return 2 * this->radius() + 1; }
+
+ // Take an array of values where the gaussian factors will be placed. Return the number of
+ // values filled.
+ int filterDouble(double values[5]) const;
+
+private:
+ double fBasis[5];
+ int fN;
+};
+
+#endif // SkGaussFilter_DEFINED
+
+
diff --git a/tests/SkGaussFilterTest.cpp b/tests/SkGaussFilterTest.cpp
new file mode 100644
index 0000000000..958e811ce5
--- /dev/null
+++ b/tests/SkGaussFilterTest.cpp
@@ -0,0 +1,88 @@
+/*
+ * Copyright 2017 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "SkGaussFilter.h"
+
+#include <cmath>
+#include <tuple>
+#include <vector>
+#include "Test.h"
+
+// one part in a million
+static constexpr double kEpsilon = 0.000001;
+
+static double careful_add(int n, double* gauss) {
+ // Sum smallest to largest to retain precision.
+ double sum = 0;
+ for (int i = n - 1; i >= 1; i--) {
+ sum += 2.0 * gauss[i];
+ }
+ sum += gauss[0];
+ return sum;
+}
+
+DEF_TEST(SkGaussFilterCommon, r) {
+ using Test = std::tuple<double, SkGaussFilter::Type, std::vector<double>>;
+
+ auto golden_check = [&](const Test& test) {
+ double sigma; SkGaussFilter::Type type; std::vector<double> golden;
+ std::tie(sigma, type, golden) = test;
+ SkGaussFilter filter{sigma, type};
+ double result[5];
+ size_t n = filter.filterDouble(result);
+ REPORTER_ASSERT(r, n == golden.size());
+ double sum = careful_add(n, result);
+ REPORTER_ASSERT(r, sum == 1.0);
+ for (size_t i = 0; i < golden.size(); i++) {
+ REPORTER_ASSERT(r, std::abs(golden[i] - result[i]) < kEpsilon);
+ }
+ };
+
+ // The following two sigmas account for about 85% of all sigmas used for masks.
+ // Golden values generated using Mathematica.
+ auto tests = {
+ // 0.788675 - most common mask sigma.
+ // GaussianMatrix[{{Automatic}, {.788675}}, Method -> "Gaussian"]
+ Test{0.788675, SkGaussFilter::Type::Gaussian, {0.506205, 0.226579, 0.0203189}},
+
+ // GaussianMatrix[{{Automatic}, {.788675}}]
+ Test{0.788675, SkGaussFilter::Type::Bessel, {0.593605, 0.176225, 0.0269721}},
+
+ // 1.07735 - second most common mask sigma.
+ // GaussianMatrix[{{Automatic}, {1.07735}}, Method -> "Gaussian"]
+ Test{1.07735, SkGaussFilter::Type::Gaussian, {0.376362, 0.244636, 0.0671835}},
+
+ // GaussianMatrix[{{4}, {1.07735}}, Method -> "Bessel"]
+ Test{1.07735, SkGaussFilter::Type::Bessel, {0.429537, 0.214955, 0.059143, 0.0111337}},
+ };
+
+ for (auto& test : tests) {
+ golden_check(test);
+ }
+}
+
+DEF_TEST(SkGaussFilterSweep, r) {
+ // The double just before 2.0.
+ const double maxSigma = nextafter(2.0, 0.0);
+ for (auto type : {SkGaussFilter::Type::Gaussian, SkGaussFilter::Type::Bessel}) {
+
+ auto check = [&](double sigma) {
+ SkGaussFilter filter{sigma, type};
+ double result[5];
+ int n = filter.filterDouble(result);
+ REPORTER_ASSERT(r, n <= 5);
+ double sum = careful_add(n, result);
+ REPORTER_ASSERT(r, sum == 1.0);
+ };
+
+ for (double sigma = 0.0; sigma < 2.0; sigma += 0.1) {
+ check(sigma);
+ }
+
+ check(maxSigma);
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-07 08:27:52 +0000