Home | History | Annotate | Download | only in core 
      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 "SkGaussFilter.h"
     10 
     11 #include <cmath>
     12 #include "SkTypes.h"
     13 
     14 static constexpr double kPi = 3.14159265358979323846264338327950288;
     15 
     16 // The value when we can stop expanding the filter. The spec implies that 3% is acceptable, but
     17 // we just use 1%.
     18 static constexpr double kGoodEnough = 1.0 / 100.0;
     19 
     20 // Normalize the values of gauss to 1.0, and make sure they add to one.
     21 // NB if n == 1, then this will force gauss[0] == 1.
     22 static void normalize(int n, double* gauss) {
     23     // Carefully add from smallest to largest to calculate the normalizing sum.
     24     double sum = 0;
     25     for (int i = n-1; i >= 1; i--) {
     26         sum += 2 * gauss[i];
     27     }
     28     sum += gauss[0];
     29 
     30     // Normalize gauss.
     31     for (int i = 0; i < n; i++) {
     32         gauss[i] /= sum;
     33     }
     34 
     35     // The factors should sum to 1. Take any remaining slop, and add it to gauss[0]. Add the
     36     // values in such a way to maintain the most accuracy.
     37     sum = 0;
     38     for (int i = n - 1; i >= 1; i--) {
     39         sum += 2 * gauss[i];
     40     }
     41 
     42     gauss[0] = 1 - sum;
     43 }
     44 
     45 static int calculate_bessel_factors(double sigma, double *gauss) {
     46     auto var = sigma * sigma;
     47 
     48     // The two functions below come from the equations in "Handbook of Mathematical Functions"
     49     // by Abramowitz and Stegun. Specifically, equation 9.6.10 on page 375. Bessel0 is given
     50     // explicitly as 9.6.12
     51     // BesselI_0 for 0 <= sigma < 2.
     52     // NB the k = 0 factor is just sum = 1.0.
     53     auto besselI_0 = [](double t) -> double {
     54         auto tSquaredOver4 = t * t / 4.0;
     55         auto sum = 1.0;
     56         auto factor = 1.0;
     57         auto k = 1;
     58         // Use a variable number of loops. When sigma is small, this only requires 3-4 loops, but
     59         // when sigma is near 2, it could require 10 loops. The same holds for BesselI_1.
     60         while(factor > 1.0/1000000.0) {
     61             factor *= tSquaredOver4 / (k * k);
     62             sum += factor;
     63             k += 1;
     64         }
     65         return sum;
     66     };
     67     // BesselI_1 for 0 <= sigma < 2.
     68     auto besselI_1 = [](double t) -> double {
     69         auto tSquaredOver4 = t * t / 4.0;
     70         auto sum = t / 2.0;
     71         auto factor = sum;
     72         auto k = 1;
     73         while (factor > 1.0/1000000.0) {
     74             factor *= tSquaredOver4 / (k * (k + 1));
     75             sum += factor;
     76             k += 1;
     77         }
     78         return sum;
     79     };
     80 
     81     // The following formula for calculating the Gaussian kernel is from
     82     // "Scale-Space for Discrete Signals" by Tony Lindeberg.
     83     // gauss(n; var) = besselI_n(var) / (e^var)
     84     auto d = std::exp(var);
     85     double b[SkGaussFilter::kGaussArrayMax] = {besselI_0(var), besselI_1(var)};
     86     gauss[0] = b[0]/d;
     87     gauss[1] = b[1]/d;
     88 
     89     // The code below is tricky, and written to mirror the recursive equations from the book.
     90     // The maximum spread for sigma == 2 is guass[4], but in order to know to stop guass[5]
     91     // is calculated. At this point n == 5 meaning that gauss[0..4] are the factors, but a 6th
     92     // element was used to calculate them.
     93     int n = 1;
     94     // The recurrence relation below is from "Numerical Recipes" 3rd Edition.
     95     // Equation 6.5.16 p.282
     96     while (gauss[n] > kGoodEnough) {
     97         b[n+1] = -(2*n/var) * b[n] + b[n-1];
     98         gauss[n+1] = b[n+1] / d;
     99         n += 1;
    100     }
    101 
    102     normalize(n, gauss);
    103 
    104     return n;
    105 }
    106 
    107 static int calculate_gauss_factors(double sigma, double* gauss) {
    108     SkASSERT(0 <= sigma && sigma < 2);
    109 
    110     // From the SVG blur spec: 8.13 Filter primitive <feGaussianBlur>.
    111     // H(x) = exp(-x^2/ (2s^2)) / sqrt(2 * s^2)
    112     auto var = sigma * sigma;
    113     auto expGaussDenom = -2 * var;
    114     auto normalizeDenom = std::sqrt(2 * kPi) * sigma;
    115 
    116     // Use the recursion relation from "Incremental Computation of the Gaussian"  by Ken
    117     // Turkowski in GPUGems 3. Page 877.
    118     double g0 = 1.0 / normalizeDenom;
    119     double g1 = std::exp(1.0 / expGaussDenom);
    120     double g2 = g1 * g1;
    121 
    122     gauss[0] = g0;
    123     g0 *= g1;
    124     g1 *= g2;
    125     gauss[1] = g0;
    126 
    127     int n = 1;
    128     while (gauss[n] > kGoodEnough) {
    129         g0 *= g1;
    130         g1 *= g2;
    131         gauss[n+1] = g0;
    132         n += 1;
    133     }
    134 
    135     normalize(n, gauss);
    136 
    137     return n;
    138 }
    139 
    140 SkGaussFilter::SkGaussFilter(double sigma, Type type) {
    141     SkASSERT(0 <= sigma && sigma < 2);
    142 
    143     if (type == Type::Bessel) {
    144         fN = calculate_bessel_factors(sigma, fBasis);
    145     } else {
    146         fN = calculate_gauss_factors(sigma, fBasis);
    147     }
    148 }
    149