1 | |
2 | // |
3 | // This source file is part of appleseed. |
4 | // Visit http://appleseedhq.net/ for additional information and resources. |
5 | // |
6 | // This software is released under the MIT license. |
7 | // |
8 | // Copyright (c) 2010-2013 Francois Beaune, Jupiter Jazz Limited |
9 | // Copyright (c) 2014-2017 Francois Beaune, The appleseedhq Organization |
10 | // |
11 | // Permission is hereby granted, free of charge, to any person obtaining a copy |
12 | // of this software and associated documentation files (the "Software"), to deal |
13 | // in the Software without restriction, including without limitation the rights |
14 | // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
15 | // copies of the Software, and to permit persons to whom the Software is |
16 | // furnished to do so, subject to the following conditions: |
17 | // |
18 | // The above copyright notice and this permission notice shall be included in |
19 | // all copies or substantial portions of the Software. |
20 | // |
21 | // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
22 | // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
23 | // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
24 | // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
25 | // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
26 | // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
27 | // THE SOFTWARE. |
28 | // |
29 | |
30 | #ifndef APPLESEED_FOUNDATION_MATH_FILTER_H |
31 | #define APPLESEED_FOUNDATION_MATH_FILTER_H |
32 | |
33 | // appleseed.foundation headers. |
34 | #include "foundation/math/fastmath.h" |
35 | #include "foundation/math/qmc.h" |
36 | #include "foundation/math/scalar.h" |
37 | #include "foundation/math/vector.h" |
38 | #include "foundation/platform/compiler.h" |
39 | |
40 | // Standard headers. |
41 | #include <cmath> |
42 | #include <cstddef> |
43 | |
44 | namespace foundation |
45 | { |
46 | |
47 | // |
48 | // Base class for 2D reconstruction filters. |
49 | // |
50 | // The filters are not normalized (they don't integrate to 1 over their domain). |
51 | // The return value of evaluate() is undefined if (x, y) is outside the filter's domain. |
52 | // |
53 | |
54 | template <typename T> |
55 | class Filter2 |
56 | { |
57 | public: |
58 | typedef T ValueType; |
59 | |
60 | Filter2(const T xradius, const T yradius); |
61 | |
62 | virtual ~Filter2() {} |
63 | |
64 | T get_xradius() const; |
65 | T get_yradius() const; |
66 | |
67 | virtual T evaluate(const T x, const T y) const = 0; |
68 | |
69 | protected: |
70 | const T m_xradius; |
71 | const T m_yradius; |
72 | const T m_rcp_xradius; |
73 | const T m_rcp_yradius; |
74 | }; |
75 | |
76 | |
77 | // |
78 | // Specializations. |
79 | // |
80 | |
81 | typedef Filter2<float> Filter2f; |
82 | typedef Filter2<double> Filter2d; |
83 | |
84 | |
85 | // |
86 | // 2D box filter. |
87 | // |
88 | |
89 | template <typename T> |
90 | class BoxFilter2 |
91 | : public Filter2<T> |
92 | { |
93 | public: |
94 | BoxFilter2(const T xradius, const T yradius); |
95 | |
96 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
97 | }; |
98 | |
99 | |
100 | // |
101 | // 2D triangle filter. |
102 | // |
103 | |
104 | template <typename T> |
105 | class TriangleFilter2 |
106 | : public Filter2<T> |
107 | { |
108 | public: |
109 | TriangleFilter2(const T xradius, const T yradius); |
110 | |
111 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
112 | }; |
113 | |
114 | |
115 | // |
116 | // 2D Gaussian filter. |
117 | // |
118 | |
119 | template <typename T> |
120 | class GaussianFilter2 |
121 | : public Filter2<T> |
122 | { |
123 | public: |
124 | GaussianFilter2( |
125 | const T xradius, |
126 | const T yradius, |
127 | const T alpha); |
128 | |
129 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
130 | |
131 | private: |
132 | const T m_alpha; |
133 | const T m_shift; |
134 | |
135 | static T gaussian(const T x, const T alpha); |
136 | }; |
137 | |
138 | template <typename T> |
139 | class FastGaussianFilter2 |
140 | : public Filter2<T> |
141 | { |
142 | public: |
143 | FastGaussianFilter2( |
144 | const T xradius, |
145 | const T yradius, |
146 | const T alpha); |
147 | |
148 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
149 | |
150 | private: |
151 | const T m_alpha; |
152 | const T m_shift; |
153 | |
154 | static T gaussian(const T x, const T alpha); |
155 | }; |
156 | |
157 | |
158 | // |
159 | // 2D Mitchell-Netravali filter. |
160 | // |
161 | // Reference: |
162 | // |
163 | // http://www.cs.utexas.edu/~fussell/courses/cs384g/lectures/mitchell/Mitchell.pdf |
164 | // |
165 | |
166 | template <typename T> |
167 | class MitchellFilter2 |
168 | : public Filter2<T> |
169 | { |
170 | public: |
171 | MitchellFilter2( |
172 | const T xradius, |
173 | const T yradius, |
174 | const T b, |
175 | const T c); |
176 | |
177 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
178 | |
179 | private: |
180 | T m_a3, m_a2, m_a0; |
181 | T m_b3, m_b2, m_b1, m_b0; |
182 | }; |
183 | |
184 | |
185 | // |
186 | // 2D Lanczos filter. |
187 | // |
188 | |
189 | template <typename T> |
190 | class LanczosFilter2 |
191 | : public Filter2<T> |
192 | { |
193 | public: |
194 | LanczosFilter2( |
195 | const T xradius, |
196 | const T yradius, |
197 | const T tau); |
198 | |
199 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
200 | |
201 | private: |
202 | const T m_rcp_tau; |
203 | |
204 | static T lanczos(const T x, const T rcp_tau); |
205 | static T sinc(const T x); |
206 | }; |
207 | |
208 | |
209 | // |
210 | // 2D Blackman-Harris filter. |
211 | // |
212 | // Reference: |
213 | // |
214 | // http://en.wikipedia.org/wiki/Window_function#Higher-order_generalized_cosine_windows |
215 | // |
216 | |
217 | template <typename T> |
218 | class BlackmanHarrisFilter2 |
219 | : public Filter2<T> |
220 | { |
221 | public: |
222 | BlackmanHarrisFilter2(const T xradius, const T yradius); |
223 | |
224 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
225 | |
226 | private: |
227 | static T blackman(const T x); |
228 | }; |
229 | |
230 | template <typename T> |
231 | class FastBlackmanHarrisFilter2 |
232 | : public Filter2<T> |
233 | { |
234 | public: |
235 | FastBlackmanHarrisFilter2(const T xradius, const T yradius); |
236 | |
237 | virtual T evaluate(const T x, const T y) const APPLESEED_OVERRIDE; |
238 | |
239 | private: |
240 | static T blackman(const T x); |
241 | }; |
242 | |
243 | |
244 | // |
245 | // Utilities. |
246 | // |
247 | |
248 | // Compute the normalization factor for a given filter. |
249 | template <typename Filter> |
250 | typename Filter::ValueType compute_normalization_factor( |
251 | const Filter& filter, |
252 | const size_t sample_count = 1024); |
253 | |
254 | |
255 | // |
256 | // Filter2 class implementation. |
257 | // |
258 | |
259 | template <typename T> |
260 | inline Filter2<T>::Filter2(const T xradius, const T yradius) |
261 | : m_xradius(xradius) |
262 | , m_yradius(yradius) |
263 | , m_rcp_xradius(T(1.0) / xradius) |
264 | , m_rcp_yradius(T(1.0) / yradius) |
265 | { |
266 | } |
267 | |
268 | template <typename T> |
269 | inline T Filter2<T>::get_xradius() const |
270 | { |
271 | return m_xradius; |
272 | } |
273 | |
274 | template <typename T> |
275 | inline T Filter2<T>::get_yradius() const |
276 | { |
277 | return m_yradius; |
278 | } |
279 | |
280 | |
281 | // |
282 | // BoxFilter2 class implementation. |
283 | // |
284 | |
285 | template <typename T> |
286 | inline BoxFilter2<T>::BoxFilter2(const T xradius, const T yradius) |
287 | : Filter2<T>(xradius, yradius) |
288 | { |
289 | } |
290 | |
291 | template <typename T> |
292 | inline T BoxFilter2<T>::evaluate(const T x, const T y) const |
293 | { |
294 | return T(1.0); |
295 | } |
296 | |
297 | |
298 | // |
299 | // TriangleFilter2 class implementation. |
300 | // |
301 | |
302 | template <typename T> |
303 | inline TriangleFilter2<T>::TriangleFilter2(const T xradius, const T yradius) |
304 | : Filter2<T>(xradius, yradius) |
305 | { |
306 | } |
307 | |
308 | template <typename T> |
309 | inline T TriangleFilter2<T>::evaluate(const T x, const T y) const |
310 | { |
311 | const T nx = x * Filter2<T>::m_rcp_xradius; |
312 | const T ny = y * Filter2<T>::m_rcp_yradius; |
313 | return (T(1.0) - std::abs(nx)) * (T(1.0) - std::abs(ny)); |
314 | } |
315 | |
316 | |
317 | // |
318 | // GaussianFilter2 class implementation. |
319 | // |
320 | |
321 | template <typename T> |
322 | inline GaussianFilter2<T>::GaussianFilter2( |
323 | const T xradius, |
324 | const T yradius, |
325 | const T alpha) |
326 | : Filter2<T>(xradius, yradius) |
327 | , m_alpha(alpha) |
328 | , m_shift(gaussian(T(1.0), alpha)) |
329 | { |
330 | } |
331 | |
332 | template <typename T> |
333 | inline T GaussianFilter2<T>::evaluate(const T x, const T y) const |
334 | { |
335 | const T nx = x * Filter2<T>::m_rcp_xradius; |
336 | const T ny = y * Filter2<T>::m_rcp_yradius; |
337 | |
338 | const T fx = gaussian(nx, m_alpha) - m_shift; |
339 | const T fy = gaussian(ny, m_alpha) - m_shift; |
340 | |
341 | return fx * fy; |
342 | } |
343 | |
344 | template <typename T> |
345 | APPLESEED_FORCE_INLINE T GaussianFilter2<T>::gaussian(const T x, const T alpha) |
346 | { |
347 | return |
348 | static_cast<T>( |
349 | std::exp( |
350 | static_cast<float>(-alpha * x * x))); |
351 | } |
352 | |
353 | |
354 | // |
355 | // FastGaussianFilter2 class implementation. |
356 | // |
357 | |
358 | template <typename T> |
359 | inline FastGaussianFilter2<T>::FastGaussianFilter2( |
360 | const T xradius, |
361 | const T yradius, |
362 | const T alpha) |
363 | : Filter2<T>(xradius, yradius) |
364 | , m_alpha(alpha) |
365 | , m_shift(gaussian(T(1.0), alpha)) |
366 | { |
367 | } |
368 | |
369 | template <typename T> |
370 | inline T FastGaussianFilter2<T>::evaluate(const T x, const T y) const |
371 | { |
372 | const T nx = x * Filter2<T>::m_rcp_xradius; |
373 | const T ny = y * Filter2<T>::m_rcp_yradius; |
374 | |
375 | const T fx = gaussian(nx, m_alpha) - m_shift; |
376 | const T fy = gaussian(ny, m_alpha) - m_shift; |
377 | |
378 | return fx * fy; |
379 | } |
380 | |
381 | template <typename T> |
382 | APPLESEED_FORCE_INLINE T FastGaussianFilter2<T>::gaussian(const T x, const T alpha) |
383 | { |
384 | // Use foundation::fast_exp() because foundation::faster_exp() is way too inaccurate. |
385 | return |
386 | static_cast<T>( |
387 | fast_exp( |
388 | static_cast<float>(-alpha * x * x))); |
389 | } |
390 | |
391 | |
392 | // |
393 | // MitchellFilter2 class implementation. |
394 | // |
395 | |
396 | template <typename T> |
397 | inline MitchellFilter2<T>::MitchellFilter2( |
398 | const T xradius, |
399 | const T yradius, |
400 | const T b, |
401 | const T c) |
402 | : Filter2<T>(xradius, yradius) |
403 | { |
404 | m_a3 = T(1.0 / 6.0) * (T(12.0) - T(9.0) * b - T(6.0) * c); |
405 | m_a2 = T(1.0 / 6.0) * (T(-18.0) + T(12.0) * b + T(6.0) * c); |
406 | m_a0 = T(1.0 / 6.0) * (T(6.0) - T(2.0) * b); |
407 | |
408 | m_b3 = T(1.0 / 6.0) * (-b - T(6.0) * c); |
409 | m_b2 = T(1.0 / 6.0) * (T(6.0) * b + T(30.0) * c); |
410 | m_b1 = T(1.0 / 6.0) * (T(-12.0) * b - T(48.0) * c); |
411 | m_b0 = T(1.0 / 6.0) * (T(8.0) * b + T(24.0) * c); |
412 | } |
413 | |
414 | template <typename T> |
415 | inline T MitchellFilter2<T>::evaluate(const T x, const T y) const |
416 | { |
417 | const T nx = x * Filter2<T>::m_rcp_xradius; |
418 | const T x1 = std::abs(nx + nx); |
419 | const T x2 = x1 * x1; |
420 | const T x3 = x2 * x1; |
421 | |
422 | const T fx = |
423 | x1 < T(1.0) |
424 | ? m_a3 * x3 + m_a2 * x2 + m_a0 |
425 | : m_b3 * x3 + m_b2 * x2 + m_b1 * x1 + m_b0; |
426 | |
427 | const T ny = y * Filter2<T>::m_rcp_yradius; |
428 | const T y1 = std::abs(ny + ny); |
429 | const T y2 = y1 * y1; |
430 | const T y3 = y2 * y1; |
431 | |
432 | const T fy = |
433 | y1 < T(1.0) |
434 | ? m_a3 * y3 + m_a2 * y2 + m_a0 |
435 | : m_b3 * y3 + m_b2 * y2 + m_b1 * y1 + m_b0; |
436 | |
437 | return fx * fy; |
438 | } |
439 | |
440 | |
441 | // |
442 | // LanczosFilter2 class implementation. |
443 | // |
444 | |
445 | template <typename T> |
446 | inline LanczosFilter2<T>::LanczosFilter2( |
447 | const T xradius, |
448 | const T yradius, |
449 | const T tau) |
450 | : Filter2<T>(xradius, yradius) |
451 | , m_rcp_tau(tau) |
452 | { |
453 | } |
454 | |
455 | template <typename T> |
456 | inline T LanczosFilter2<T>::evaluate(const T x, const T y) const |
457 | { |
458 | const T nx = x * Filter2<T>::m_rcp_xradius; |
459 | const T ny = y * Filter2<T>::m_rcp_yradius; |
460 | return lanczos(nx, m_rcp_tau) * lanczos(ny, m_rcp_tau); |
461 | } |
462 | |
463 | template <typename T> |
464 | APPLESEED_FORCE_INLINE T LanczosFilter2<T>::lanczos(const T x, const T rcp_tau) |
465 | { |
466 | const T theta = Pi<T>() * x; |
467 | return theta == T(0.0) ? T(1.0) : sinc(theta * rcp_tau) * sinc(theta); |
468 | } |
469 | |
470 | template <typename T> |
471 | APPLESEED_FORCE_INLINE T LanczosFilter2<T>::sinc(const T x) |
472 | { |
473 | return std::sin(x) / x; |
474 | } |
475 | |
476 | |
477 | // |
478 | // BlackmanHarrisFilter2 class implementation. |
479 | // |
480 | |
481 | template <typename T> |
482 | inline BlackmanHarrisFilter2<T>::BlackmanHarrisFilter2(const T xradius, const T yradius) |
483 | : Filter2<T>(xradius, yradius) |
484 | { |
485 | } |
486 | |
487 | template <typename T> |
488 | inline T BlackmanHarrisFilter2<T>::evaluate(const T x, const T y) const |
489 | { |
490 | const T nx = T(0.5) * (T(1.0) + x * Filter2<T>::m_rcp_xradius); |
491 | const T ny = T(0.5) * (T(1.0) + y * Filter2<T>::m_rcp_yradius); |
492 | return blackman(nx) * blackman(ny); |
493 | } |
494 | |
495 | template <typename T> |
496 | APPLESEED_FORCE_INLINE T BlackmanHarrisFilter2<T>::blackman(const T x) |
497 | { |
498 | return |
499 | T(0.35875) |
500 | - T(0.48829) * std::cos((T(2.0) * Pi<T>()) * x) |
501 | + T(0.14128) * std::cos((T(4.0) * Pi<T>()) * x) |
502 | - T(0.01174) * std::cos((T(6.0) * Pi<T>()) * x); // original coefficient is 0.01168, modified to ensure 0 at borders |
503 | } |
504 | |
505 | |
506 | // |
507 | // FastBlackmanHarrisFilter2 class implementation. |
508 | // |
509 | |
510 | template <typename T> |
511 | inline FastBlackmanHarrisFilter2<T>::FastBlackmanHarrisFilter2(const T xradius, const T yradius) |
512 | : Filter2<T>(xradius, yradius) |
513 | { |
514 | } |
515 | |
516 | template <typename T> |
517 | inline T FastBlackmanHarrisFilter2<T>::evaluate(const T x, const T y) const |
518 | { |
519 | const T nx = T(0.5) * (T(1.0) + x * Filter2<T>::m_rcp_xradius); |
520 | const T ny = T(0.5) * (T(1.0) + y * Filter2<T>::m_rcp_yradius); |
521 | return blackman(nx) * blackman(ny); |
522 | } |
523 | |
524 | template <typename T> |
525 | APPLESEED_FORCE_INLINE T FastBlackmanHarrisFilter2<T>::blackman(const T x) |
526 | { |
527 | return |
528 | T(0.35875) |
529 | - T(0.48829) * fast_cos_full_positive((T(2.0) * Pi<T>()) * x) |
530 | + T(0.14128) * fast_cos_full_positive((T(4.0) * Pi<T>()) * x) |
531 | - T(0.01174) * fast_cos_full_positive((T(6.0) * Pi<T>()) * x); // original coefficient is 0.01168, modified to ensure 0 at borders |
532 | } |
533 | |
534 | |
535 | // |
536 | // Utilities implementation. |
537 | // |
538 | |
539 | template <typename Filter> |
540 | typename Filter::ValueType compute_normalization_factor( |
541 | const Filter& filter, |
542 | const size_t sample_count) |
543 | { |
544 | typedef typename Filter::ValueType ValueType; |
545 | |
546 | const ValueType xradius = filter.get_xradius(); |
547 | const ValueType yradius = filter.get_yradius(); |
548 | |
549 | ValueType result(0.0); |
550 | |
551 | for (size_t i = 0; i < sample_count; ++i) |
552 | { |
553 | static const size_t Bases[1] = { 2 }; |
554 | |
555 | const Vector<ValueType, 2> s = |
556 | hammersley_sequence<ValueType, 2>(Bases, sample_count, i); |
557 | |
558 | const Vector<ValueType, 2> p( |
559 | xradius * (ValueType(2.0) * s.x - ValueType(1.0)), |
560 | yradius * (ValueType(2.0) * s.y - ValueType(1.0))); |
561 | |
562 | result += filter.evaluate(p.x, p.y); |
563 | } |
564 | |
565 | result *= ValueType(4.0) * xradius * yradius; |
566 | result /= static_cast<ValueType>(sample_count); |
567 | |
568 | return result; |
569 | } |
570 | |
571 | } // namespace foundation |
572 | |
573 | #endif // !APPLESEED_FOUNDATION_MATH_FILTER_H |
574 | |