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
44namespace 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
54template <typename T>
55class 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
81typedef Filter2<float> Filter2f;
82typedef Filter2<double> Filter2d;
83
84
85//
86// 2D box filter.
87//
88
89template <typename T>
90class 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
104template <typename T>
105class 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
119template <typename T>
120class 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
138template <typename T>
139class 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
166template <typename T>
167class 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
189template <typename T>
190class 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
217template <typename T>
218class 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
230template <typename T>
231class 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.
249template <typename Filter>
250typename 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
259template <typename T>
260inline 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
268template <typename T>
269inline T Filter2<T>::get_xradius() const
270{
271 return m_xradius;
272}
273
274template <typename T>
275inline T Filter2<T>::get_yradius() const
276{
277 return m_yradius;
278}
279
280
281//
282// BoxFilter2 class implementation.
283//
284
285template <typename T>
286inline BoxFilter2<T>::BoxFilter2(const T xradius, const T yradius)
287 : Filter2<T>(xradius, yradius)
288{
289}
290
291template <typename T>
292inline 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
302template <typename T>
303inline TriangleFilter2<T>::TriangleFilter2(const T xradius, const T yradius)
304 : Filter2<T>(xradius, yradius)
305{
306}
307
308template <typename T>
309inline 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
321template <typename T>
322inline 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
332template <typename T>
333inline 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
344template <typename T>
345APPLESEED_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
358template <typename T>
359inline 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
369template <typename T>
370inline 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
381template <typename T>
382APPLESEED_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
396template <typename T>
397inline 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
414template <typename T>
415inline 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
445template <typename T>
446inline 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
455template <typename T>
456inline 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
463template <typename T>
464APPLESEED_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
470template <typename T>
471APPLESEED_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
481template <typename T>
482inline BlackmanHarrisFilter2<T>::BlackmanHarrisFilter2(const T xradius, const T yradius)
483 : Filter2<T>(xradius, yradius)
484{
485}
486
487template <typename T>
488inline 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
495template <typename T>
496APPLESEED_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
510template <typename T>
511inline FastBlackmanHarrisFilter2<T>::FastBlackmanHarrisFilter2(const T xradius, const T yradius)
512 : Filter2<T>(xradius, yradius)
513{
514}
515
516template <typename T>
517inline 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
524template <typename T>
525APPLESEED_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
539template <typename Filter>
540typename 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