raster/gaussian_blur.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/raster/gaussian_blur.h" | ||
| 2 | |||
| 3 | #include "terraingen/util/checks.h" | ||
| 4 | |||
| 5 | #include <cmath> | ||
| 6 | #include <cstddef> | ||
| 7 | #include <stdexcept> | ||
| 8 | #include <vector> | ||
| 9 | |||
| 10 | namespace terraingen { | ||
| 11 | |||
| 12 | namespace { | ||
| 13 | |||
| 14 | /// Unnormalized Gaussian weights of standard deviation `sigma` cells out to three sigma. | ||
| 15 | 24 | std::vector<double> GaussianKernel(double sigma) { | |
| 16 | 24 | if (sigma == 0.0) { | |
| 17 | 3 | return {1.0}; | |
| 18 | } | ||
| 19 | 23 | const auto radius = static_cast<std::ptrdiff_t>(std::ceil(3.0 * sigma)); | |
| 20 | 23 | std::vector<double> kernel(static_cast<std::size_t>(2 * radius + 1)); | |
| 21 | 336 | for (std::ptrdiff_t i = -radius; i <= radius; ++i) { | |
| 22 | 313 | const auto x = static_cast<double>(i); | |
| 23 | 313 | kernel[static_cast<std::size_t>(i + radius)] = std::exp(-0.5 * x * x / (sigma * sigma)); | |
| 24 | } | ||
| 25 | 23 | return kernel; | |
| 26 | 23 | } | |
| 27 | |||
| 28 | /// Convolves `values` along one axis of a `width` x `height` grid with `kernel`, truncated at the | ||
| 29 | /// grid's edges. | ||
| 30 | 48 | std::vector<double> Convolve(const std::vector<double>& values, std::size_t width, | |
| 31 | std::size_t height, const std::vector<double>& kernel, | ||
| 32 | bool alongRows) { | ||
| 33 | 48 | const auto radius = static_cast<std::ptrdiff_t>(kernel.size() / 2); | |
| 34 | 48 | const auto length = static_cast<std::ptrdiff_t>(alongRows ? width : height); | |
| 35 | 48 | std::vector<double> result(values.size(), 0.0); | |
| 36 | 1008 | for (std::size_t row = 0; row < height; ++row) { | |
| 37 | 32252 | for (std::size_t col = 0; col < width; ++col) { | |
| 38 | 31292 | const auto position = static_cast<std::ptrdiff_t>(alongRows ? col : row); | |
| 39 | 31292 | double sum = 0.0; | |
| 40 | 538184 | for (std::ptrdiff_t k = -radius; k <= radius; ++k) { | |
| 41 | 506892 | const std::ptrdiff_t at = position + k; | |
| 42 | 506892 | if (at < 0 || at >= length) { | |
| 43 | 68844 | continue; | |
| 44 | } | ||
| 45 | 438048 | const std::size_t index = alongRows ? row * width + static_cast<std::size_t>(at) | |
| 46 | 214086 | : static_cast<std::size_t>(at) * width + col; | |
| 47 | 438048 | sum += kernel[static_cast<std::size_t>(k + radius)] * values[index]; | |
| 48 | } | ||
| 49 | 31292 | result[row * width + col] = sum; | |
| 50 | } | ||
| 51 | } | ||
| 52 | 48 | return result; | |
| 53 | } | ||
| 54 | |||
| 55 | } // namespace | ||
| 56 | |||
| 57 | 10 | Raster<float> GaussianBlur(const Raster<float>& raster, double sigmaPixels) { | |
| 58 | 10 | return GaussianBlur(raster, sigmaPixels, sigmaPixels); | |
| 59 | } | ||
| 60 | |||
| 61 | 16 | Raster<float> GaussianBlur(const Raster<float>& raster, double sigmaColumns, double sigmaRows) { | |
| 62 | 16 | if (!(IsNonNegative(sigmaColumns) && IsNonNegative(sigmaRows))) { | |
| 63 | 2 | throw std::invalid_argument("Gaussian sigma must be finite and not negative"); | |
| 64 | } | ||
| 65 | 14 | Raster<float> result = raster; | |
| 66 | 14 | if ((sigmaColumns == 0.0 && sigmaRows == 0.0) || raster.Empty()) { | |
| 67 | 2 | return result; | |
| 68 | } | ||
| 69 | 12 | const std::size_t width = raster.Width(); | |
| 70 | 12 | const std::size_t height = raster.Height(); | |
| 71 | 24 | std::vector<double> weighted(raster.Size()); | |
| 72 | 12 | std::vector<double> mask(raster.Size()); | |
| 73 | 7835 | for (std::size_t i = 0; i < raster.Size(); ++i) { | |
| 74 | 7823 | const float value = raster.Values()[i]; | |
| 75 | 7823 | const bool valid = !raster.IsNoData(value); | |
| 76 | 7823 | mask[i] = valid ? 1.0 : 0.0; | |
| 77 | 7823 | weighted[i] = valid ? static_cast<double>(value) : 0.0; | |
| 78 | } | ||
| 79 | 12 | const std::vector<double> across = GaussianKernel(sigmaColumns); | |
| 80 | 12 | const std::vector<double> down = GaussianKernel(sigmaRows); | |
| 81 | const std::vector<double> numerator = | ||
| 82 | 12 | Convolve(Convolve(weighted, width, height, across, true), width, height, down, false); | |
| 83 | const std::vector<double> denominator = | ||
| 84 | 12 | Convolve(Convolve(mask, width, height, across, true), width, height, down, false); | |
| 85 | 7835 | for (std::size_t i = 0; i < raster.Size(); ++i) { | |
| 86 | 7823 | if (mask[i] > 0.0) { | |
| 87 | 7821 | result.Values()[i] = static_cast<float>(numerator[i] / denominator[i]); | |
| 88 | } | ||
| 89 | } | ||
| 90 | 12 | return result; | |
| 91 | 12 | } | |
| 92 | |||
| 93 | } // namespace terraingen | ||
| 94 |