1 #include <mbgl/util/tiny_sdf.hpp>
2 
3 #include <mbgl/util/math.hpp>
4 
5 #include <algorithm>
6 
7 namespace mbgl {
8 namespace util {
9 
10 namespace tinysdf {
11 
12 static const double INF = 1e20;
13 
14 // 1D squared distance transform
edt1d(std::vector<double> & f,std::vector<double> & d,std::vector<int16_t> & v,std::vector<double> & z,uint32_t n)15 void edt1d(std::vector<double>& f,
16            std::vector<double>& d,
17            std::vector<int16_t>& v,
18            std::vector<double>& z,
19            uint32_t n) {
20     v[0] = 0;
21     z[0] = -INF;
22     z[1] = +INF;
23 
24     for (uint32_t q = 1, k = 0; q < n; q++) {
25         double s = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * q - 2 * v[k]);
26         while (s <= z[k]) {
27             k--;
28             s = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * q - 2 * v[k]);
29         }
30         k++;
31         v[k] = q;
32         z[k] = s;
33         z[k + 1] = +INF;
34     }
35 
36     for (uint32_t q = 0, k = 0; q < n; q++) {
37         while (z[k + 1] < q) k++;
38         d[q] = (q - v[k]) * (q - v[k]) + f[v[k]];
39     }
40 }
41 
42 
43 // 2D Euclidean distance transform by Felzenszwalb & Huttenlocher https://cs.brown.edu/~pff/dt/
edt(std::vector<double> & data,uint32_t width,uint32_t height,std::vector<double> & f,std::vector<double> & d,std::vector<int16_t> & v,std::vector<double> & z)44 void edt(std::vector<double>& data,
45          uint32_t width,
46          uint32_t height,
47          std::vector<double>& f,
48          std::vector<double>& d,
49          std::vector<int16_t>& v,
50          std::vector<double>& z) {
51     for (uint32_t x = 0; x < width; x++) {
52         for (uint32_t y = 0; y < height; y++) {
53             f[y] = data[y * width + x];
54         }
55         edt1d(f, d, v, z, height);
56         for (uint32_t y = 0; y < height; y++) {
57             data[y * width + x] = d[y];
58         }
59     }
60     for (uint32_t y = 0; y < height; y++) {
61         for (uint32_t x = 0; x < width; x++) {
62             f[x] = data[y * width + x];
63         }
64         edt1d(f, d, v, z, width);
65         for (uint32_t x = 0; x < width; x++) {
66             data[y * width + x] = std::sqrt(d[x]);
67         }
68     }
69 }
70 
71 } // namespace tinysdf
72 
transformRasterToSDF(const AlphaImage & rasterInput,double radius,double cutoff)73 AlphaImage transformRasterToSDF(const AlphaImage& rasterInput, double radius, double cutoff) {
74     uint32_t size = rasterInput.size.width * rasterInput.size.height;
75     uint32_t maxDimension = std::max(rasterInput.size.width, rasterInput.size.height);
76 
77     AlphaImage sdf(rasterInput.size);
78 
79     // temporary arrays for the distance transform
80     std::vector<double> gridOuter(size);
81     std::vector<double> gridInner(size);
82     std::vector<double> f(maxDimension);
83     std::vector<double> d(maxDimension);
84     std::vector<double> z(maxDimension + 1);
85     std::vector<int16_t> v(maxDimension);
86 
87     for (uint32_t i = 0; i < size; i++) {
88         double a = double(rasterInput.data[i]) / 255; // alpha value
89         gridOuter[i] = a == 1.0 ? 0.0 : a == 0.0 ? tinysdf::INF : std::pow(std::max(0.0, 0.5 - a), 2.0);
90         gridInner[i] = a == 1.0 ? tinysdf::INF : a == 0.0 ? 0.0 : std::pow(std::max(0.0, a - 0.5), 2.0);
91     }
92 
93     tinysdf::edt(gridOuter, rasterInput.size.width, rasterInput.size.height, f, d, v, z);
94     tinysdf::edt(gridInner, rasterInput.size.width, rasterInput.size.height, f, d, v, z);
95 
96     for (uint32_t i = 0; i < size; i++) {
97         double distance = gridOuter[i] - gridInner[i];
98         sdf.data[i] = std::max(0l, std::min(255l, ::lround(255.0 - 255.0 * (distance / radius + cutoff))));
99     }
100 
101     return sdf;
102 }
103 
104 } // namespace util
105 } // namespace mbgl
106