1 #include <algorithm>
2 #include <cfloat>
3 #include <climits>
4 #include <cmath>
5 #include <cstddef>
6 #include <cstdlib>
7 #include <stdexcept>
8 #include <vector>
9 #include "common/except.h"
10 #include "common/libm_wrapper.h"
11 #include "common/matrix.h"
12 #include "common/zassert.h"
13 #include "filter.h"
14
15 namespace zimg {
16 namespace resize {
17
18 namespace {
19
20 constexpr double PI = 3.14159265358979323846;
21
sinc(double x)22 double sinc(double x) noexcept
23 {
24 // Guaranteed to not yield division by zero on IEEE machine with accurate sin(x).
25 return x == 0.0 ? 1.0 : zimg_x_sin(x * PI) / (x * PI);
26 }
27
poly3(double x,double c0,double c1,double c2,double c3)28 double poly3(double x, double c0, double c1, double c2, double c3) noexcept
29 {
30 return c0 + x * (c1 + x * (c2 + x * c3));
31 }
32
round_halfup(double x)33 double round_halfup(double x) noexcept
34 {
35 /* When rounding on the pixel grid, the invariant
36 * round(x - 1) == round(x) - 1
37 * must be preserved. This precludes the use of modes such as
38 * half-to-even and half-away-from-zero.
39 */
40 return x < 0 ? std::floor(x + 0.5) : std::floor(x + 0.49999999999999994);
41 }
42
43
matrix_to_filter(const RowMatrix<double> & m)44 FilterContext matrix_to_filter(const RowMatrix<double> &m)
45 {
46 size_t width = 0;
47
48 for (size_t i = 0; i < m.rows(); ++i) {
49 width = std::max(width, m.row_right(i) - m.row_left(i));
50 }
51 zassert_d(width, "empty matrix");
52
53 if (width > floor_n(UINT_MAX, AlignmentOf<uint16_t>::value))
54 error::throw_<error::OutOfMemory>();
55 if (width > floor_n(UINT_MAX, AlignmentOf<float>::value))
56 error::throw_<error::OutOfMemory>();
57
58 FilterContext e{};
59
60 try {
61 e.filter_width = static_cast<unsigned>(width);
62 e.filter_rows = static_cast<unsigned>(m.rows());
63 e.input_width = static_cast<unsigned>(m.cols());
64 e.stride = static_cast<unsigned>(ceil_n(width, AlignmentOf<float>::value));
65 e.stride_i16 = static_cast<unsigned>(ceil_n(width, AlignmentOf<uint16_t>::value));
66
67 if (e.filter_rows > UINT_MAX / e.stride || e.filter_rows > UINT_MAX / e.stride_i16)
68 error::throw_<error::OutOfMemory>();
69
70 e.data.resize(static_cast<size_t>(e.stride) * e.filter_rows);
71 e.data_i16.resize(static_cast<size_t>(e.stride_i16) * e.filter_rows);
72 e.left.resize(e.filter_rows);
73 } catch (const std::length_error &) {
74 error::throw_<error::OutOfMemory>();
75 }
76
77 for (size_t i = 0; i < m.rows(); ++i) {
78 unsigned left = static_cast<unsigned>(std::min(m.row_left(i), m.cols() - width));
79 double f32_err = 0.0f;
80 double i16_err = 0;
81
82 double f32_sum = 0.0;
83 int16_t i16_sum = 0;
84 int16_t i16_greatest = 0;
85 size_t i16_greatest_idx = 0;
86
87 /* Dither filter coefficients when rounding them to their storage format.
88 * This minimizes accumulation of error and ensures that the filter
89 * continues to sum as close to 1.0 as possible after rounding.
90 */
91 for (size_t j = 0; j < width; ++j) {
92 double coeff = m[i][left + j];
93
94 double coeff_expected_f32 = coeff - f32_err;
95 double coeff_expected_i16 = coeff * (1 << 14) - i16_err;
96
97 float coeff_f32 = static_cast<float>(coeff_expected_f32);
98 int16_t coeff_i16 = static_cast<int16_t>(std::lrint(coeff_expected_i16));
99
100 f32_err = static_cast<double>(coeff_f32) - coeff_expected_f32;
101 i16_err = static_cast<double>(coeff_i16) - coeff_expected_i16;
102
103 if (std::abs(coeff_i16) > i16_greatest) {
104 i16_greatest = coeff_i16;
105 i16_greatest_idx = j;
106 }
107
108 f32_sum += coeff_f32;
109 i16_sum += coeff_i16;
110
111 e.data[i * e.stride + j] = coeff_f32;
112 e.data_i16[i * e.stride_i16 + j] = coeff_i16;
113 }
114
115 /* The final sum may still be off by a few ULP. This can not be fixed for
116 * floating point data, since the error is dependent on summation order,
117 * but for integer data, the error can be added to the greatest coefficient.
118 */
119 zassert_d(1.0 - f32_sum <= FLT_EPSILON, "error too great");
120 zassert_d(std::abs((1 << 14) - i16_sum) <= 1, "error too great");
121
122 e.data_i16[i * e.stride_i16 + i16_greatest_idx] += (1 << 14) - i16_sum;
123
124 e.left[i] = left;
125 }
126
127 return e;
128 }
129
130 } // namespace
131
132
133 Filter::~Filter() = default;
134
support() const135 unsigned PointFilter::support() const { return 0; }
136
operator ()(double x) const137 double PointFilter::operator()(double x) const { return 1.0; }
138
139
support() const140 unsigned BilinearFilter::support() const { return 1; }
141
operator ()(double x) const142 double BilinearFilter::operator()(double x) const
143 {
144 return std::max(1.0 - std::abs(x), 0.0);
145 }
146
147
BicubicFilter(double b,double c)148 BicubicFilter::BicubicFilter(double b, double c) :
149 p0{ ( 6.0 - 2.0 * b ) / 6.0 },
150 p2{ (-18.0 + 12.0 * b + 6.0 * c) / 6.0 },
151 p3{ ( 12.0 - 9.0 * b - 6.0 * c) / 6.0 },
152 q0{ ( 8.0 * b + 24.0 * c) / 6.0 },
153 q1{ ( -12.0 * b - 48.0 * c) / 6.0 },
154 q2{ ( 6.0 * b + 30.0 * c) / 6.0 },
155 q3{ ( -b - 6.0 * c) / 6.0 }
156 {}
157
support() const158 unsigned BicubicFilter::support() const { return 2; }
159
operator ()(double x) const160 double BicubicFilter::operator()(double x) const
161 {
162 x = std::abs(x);
163
164 if (x < 1.0)
165 return poly3(x, p0, 0.0, p2, p3);
166 else if (x < 2.0)
167 return poly3(x, q0, q1, q2, q3);
168 else
169 return 0.0;
170 }
171
172
support() const173 unsigned Spline16Filter::support() const { return 2; }
174
operator ()(double x) const175 double Spline16Filter::operator()(double x) const
176 {
177 x = std::abs(x);
178
179 if (x < 1.0) {
180 return poly3(x, 1.0, -1.0 / 5.0, -9.0 / 5.0, 1.0);
181 } else if (x < 2.0) {
182 x -= 1.0;
183 return poly3(x, 0.0, -7.0 / 15.0, 4.0 / 5.0, -1.0 / 3.0);
184 } else {
185 return 0.0;
186 }
187 }
188
189
support() const190 unsigned Spline36Filter::support() const { return 3; }
191
operator ()(double x) const192 double Spline36Filter::operator()(double x) const
193 {
194 x = std::abs(x);
195
196 if (x < 1.0) {
197 return poly3(x, 1.0, -3.0 / 209.0, -453.0 / 209.0, 13.0 / 11.0);
198 } else if (x < 2.0) {
199 x -= 1.0;
200 return poly3(x, 0.0, -156.0 / 209.0, 270.0 / 209.0, -6.0 / 11.0);
201 } else if (x < 3.0) {
202 x -= 2.0;
203 return poly3(x, 0.0, 26.0 / 209.0, -45.0 / 209.0, 1.0 / 11.0);
204 } else {
205 return 0.0;
206 }
207 }
208
209
support() const210 unsigned Spline64Filter::support() const { return 4; }
211
operator ()(double x) const212 double Spline64Filter::operator()(double x) const
213 {
214 x = std::abs(x);
215
216 if (x < 1.0) {
217 return poly3(x, 1.0, -3.0 / 2911.0, -6387.0 / 2911.0, 49.0 / 41.0);
218 } else if (x < 2.0) {
219 x -= 1.0;
220 return poly3(x, 0.0, -2328.0 / 2911.0, 4032.0 / 2911.0, -24.0 / 41.0);
221 } else if (x < 3.0) {
222 x -= 2.0;
223 return poly3(x, 0.0, 582.0 / 2911.0, -1008.0 / 2911.0, 6.0 / 41.0);
224 } else if (x < 4.0) {
225 x -= 3.0;
226 return poly3(x, 0.0, -97.0 / 2911.0, 168.0 / 2911.0, -1.0 / 41.0);
227 } else {
228 return 0.0;
229 }
230 }
231
232
LanczosFilter(unsigned taps)233 LanczosFilter::LanczosFilter(unsigned taps) : taps{ taps }
234 {
235 if (taps <= 0)
236 error::throw_<error::IllegalArgument>("lanczos tap count must be positive");
237 }
238
support() const239 unsigned LanczosFilter::support() const { return taps; }
240
operator ()(double x) const241 double LanczosFilter::operator()(double x) const
242 {
243 x = std::abs(x);
244 return x < taps ? sinc(x) * sinc(x / taps) : 0.0;
245 }
246
247
compute_filter(const Filter & f,unsigned src_dim,unsigned dst_dim,double shift,double width)248 FilterContext compute_filter(const Filter &f, unsigned src_dim, unsigned dst_dim, double shift, double width)
249 {
250 double scale = static_cast<double>(dst_dim) / width;
251 double step = std::min(scale, 1.0);
252 double support = static_cast<double>(f.support()) / step;
253 unsigned filter_size = std::max(static_cast<unsigned>(std::ceil(support)) * 2U, 1U);
254
255 if (support > static_cast<unsigned>(UINT_MAX / 2))
256 error::throw_<error::ResamplingNotAvailable>("filter width too great");
257
258 try {
259 RowMatrix<double> m{ dst_dim, src_dim };
260
261 for (unsigned i = 0; i < dst_dim; ++i) {
262 // Position of output sample on input grid.
263 double pos = (i + 0.5) / scale + shift;
264 double begin_pos = round_halfup(pos - filter_size / 2.0) + 0.5;
265
266 double total = 0.0;
267 for (unsigned j = 0; j < filter_size; ++j) {
268 double xpos = begin_pos + j;
269 total += f((xpos - pos) * step);
270 }
271
272 size_t left = SIZE_MAX;
273
274 for (unsigned j = 0; j < filter_size; ++j) {
275 double xpos = begin_pos + j;
276 double real_pos;
277
278 // Mirror the position if it goes beyond image bounds.
279 if (xpos < 0.0)
280 real_pos = -xpos;
281 else if (xpos >= src_dim)
282 real_pos = 2.0 * src_dim - xpos;
283 else
284 real_pos = xpos;
285
286 // Clamp the position if it is still out of bounds.
287 real_pos = std::min(std::max(real_pos, 0.0), std::nextafter(src_dim, -INFINITY));
288
289 size_t idx = static_cast<size_t>(std::floor(real_pos));
290 m[i][idx] += f((xpos - pos) * step) / total;
291 left = std::min(left, idx);
292 }
293
294 // Force allocating an entry to keep the left offset table sorted.
295 if (m[i][left] == 0.0) {
296 m[i][left] = DBL_EPSILON;
297 m[i][left] = 0.0;
298 }
299 }
300
301 return matrix_to_filter(m);
302 } catch (const std::length_error &) {
303 error::throw_<error::OutOfMemory>();
304 }
305 }
306
307 } // namespace resize
308 } // namespace zimg
309