Revert "Gauss filter calculation"

This reverts commit 53ec7dc7cb.

Reason for revert: Segv on very specific machines.

Original change's description:
> Gauss filter calculation
> 
> Change-Id: I921ef815d4f788c312aa729f353b6ea154140555
> Reviewed-on: https://skia-review.googlesource.com/67723
> Commit-Queue: Herb Derby <herb@google.com>
> Reviewed-by: Robert Phillips <robertphillips@google.com>

TBR=herb@google.com,robertphillips@google.com

Change-Id: I15164809d081dee0076e815b40fbfdbc6374cfba
No-Presubmit: true
No-Tree-Checks: true
No-Try: true
Reviewed-on: https://skia-review.googlesource.com/69641
Reviewed-by: Herb Derby <herb@google.com>
Commit-Queue: Herb Derby <herb@google.com>
This commit is contained in:
Herb Derby 2017-11-09 22:39:51 +00:00 committed by Skia Commit-Bot
parent 77e487dfc0
commit 66918078bb
5 changed files with 0 additions and 284 deletions

View File

@ -134,8 +134,6 @@ 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",

View File

@ -212,7 +212,6 @@ tests_sources = [
"$_tests/SkColor4fTest.cpp",
"$_tests/SkDOMTest.cpp",
"$_tests/SkFixed15Test.cpp",
"$_tests/SkGaussFilterTest.cpp",
"$_tests/SkImageTest.cpp",
"$_tests/SkLiteDLTest.cpp",
"$_tests/SkNxTest.cpp",

View File

@ -1,152 +0,0 @@
/*
* 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;
}

View File

@ -1,41 +0,0 @@
/*
* 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

View File

@ -1,88 +0,0 @@
/*
* 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);
}
}