SkGaussFilter.cpp 3.5 KB

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