diff options
author | Herb Derby <herb@google.com> | 2017-11-03 13:36:55 -0400 |
---|---|---|
committer | Skia Commit-Bot <skia-commit-bot@chromium.org> | 2017-11-10 19:58:57 +0000 |
commit | 66498bc41604fd1c9c43580c74a542813b97b549 (patch) | |
tree | 01f38bfb3c151739ecd5ba1daa67414dc56b064a | |
parent | c0ae2c8a6272bec976a6b617408f98856c79862d (diff) |
Try 2 for Gauss filter calculation
Originally reviewed at:
https://skia-review.googlesource.com/c/skia/+/67723
Change-Id: Ie62d81f818899f3a79df888c1594d3fbccf6d414
Reviewed-on: https://skia-review.googlesource.com/69681
Commit-Queue: Herb Derby <herb@google.com>
Reviewed-by: Greg Daniel <egdaniel@google.com>
-rw-r--r-- | gn/core.gni | 2 | ||||
-rw-r--r-- | gn/tests.gni | 1 | ||||
-rw-r--r-- | src/core/SkGaussFilter.cpp | 152 | ||||
-rw-r--r-- | src/core/SkGaussFilter.h | 41 | ||||
-rw-r--r-- | tests/SkGaussFilterTest.cpp | 97 |
5 files changed, 293 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..2fdf4ac573 --- /dev/null +++ b/tests/SkGaussFilterTest.cpp @@ -0,0 +1,97 @@ +/* + * 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); + auto check = [&](double sigma, SkGaussFilter::Type type) { + 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, SkGaussFilter::Type::Gaussian); + } + + check(maxSigma, SkGaussFilter::Type::Gaussian); + } + + { + + for (double sigma = 0.0; sigma < 2.0; sigma += 0.1) { + check(sigma, SkGaussFilter::Type::Bessel); + } + + check(maxSigma, SkGaussFilter::Type::Bessel); + } +} |