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_VECTOR_H
31#define APPLESEED_FOUNDATION_MATH_VECTOR_H
32
33// appleseed.foundation headers.
34#include "foundation/math/scalar.h"
35#include "foundation/utility/poison.h"
36
37// Imath headers.
38#ifdef APPLESEED_ENABLE_IMATH_INTEROP
39#include "foundation/platform/exrheaderguards.h"
40BEGIN_EXR_INCLUDES
41#include "OpenEXR/ImathVec.h"
42END_EXR_INCLUDES
43#endif
44
45// Standard headers.
46#include <algorithm>
47#include <cassert>
48#include <cmath>
49#include <cstddef>
50
51namespace foundation
52{
53
54//
55// N-dimensional vector class of arbitrary type.
56//
57
58template <typename T, size_t N>
59class Vector
60{
61 public:
62 // Types.
63 typedef T ValueType;
64 typedef Vector<T, N> VectorType;
65
66 // Dimension.
67 static const size_t Dimension = N;
68
69 // Constructors.
70 Vector(); // leave all components uninitialized
71 explicit Vector(const ValueType* rhs); // initialize with array of N scalars
72 explicit Vector(const ValueType val); // set all components to 'val'
73
74 // Construct a vector from another vector of a different type.
75 template <typename U>
76 explicit Vector(const Vector<U, N>& rhs);
77
78 // Unchecked array subscripting.
79 ValueType& operator[](const size_t i);
80 const ValueType& operator[](const size_t i) const;
81
82 private:
83 // Vector components.
84 ValueType m_comp[N];
85};
86
87// Poisoning.
88template <typename T, size_t N>
89class PoisonImpl<Vector<T, N> >
90{
91 public:
92 static void do_poison(Vector<T, N>& v);
93};
94
95// Exact inequality and equality tests.
96template <typename T, size_t N> bool operator!=(const Vector<T, N>& lhs, const Vector<T, N>& rhs);
97template <typename T, size_t N> bool operator==(const Vector<T, N>& lhs, const Vector<T, N>& rhs);
98
99// Approximate equality tests.
100template <typename T, size_t N> bool feq(const Vector<T, N>& lhs, const Vector<T, N>& rhs);
101template <typename T, size_t N> bool feq(const Vector<T, N>& lhs, const Vector<T, N>& rhs, const T eps);
102
103// Approximate zero tests.
104template <typename T, size_t N> bool fz(const Vector<T, N>& v);
105template <typename T, size_t N> bool fz(const Vector<T, N>& v, const T eps);
106
107// Vector arithmetic.
108template <typename T, size_t N> Vector<T, N> operator+ (const Vector<T, N>& lhs, const Vector<T, N>& rhs);
109template <typename T, size_t N> Vector<T, N> operator- (const Vector<T, N>& lhs, const Vector<T, N>& rhs);
110template <typename T, size_t N> Vector<T, N> operator- (const Vector<T, N>& lhs);
111template <typename T, size_t N> Vector<T, N> operator* (const Vector<T, N>& lhs, const T rhs);
112template <typename T, size_t N> Vector<T, N> operator* (const T lhs, const Vector<T, N>& rhs);
113template <typename T, size_t N> Vector<T, N> operator/ (const Vector<T, N>& lhs, const T rhs);
114template <typename T, size_t N> Vector<T, N> operator* (const Vector<T, N>& lhs, const Vector<T, N>& rhs);
115template <typename T, size_t N> Vector<T, N> operator/ (const Vector<T, N>& lhs, const Vector<T, N>& rhs);
116template <typename T, size_t N> Vector<T, N>& operator+=(Vector<T, N>& lhs, const Vector<T, N>& rhs);
117template <typename T, size_t N> Vector<T, N>& operator-=(Vector<T, N>& lhs, const Vector<T, N>& rhs);
118template <typename T, size_t N> Vector<T, N>& operator*=(Vector<T, N>& lhs, const T rhs);
119template <typename T, size_t N> Vector<T, N>& operator/=(Vector<T, N>& lhs, const T rhs);
120template <typename T, size_t N> Vector<T, N>& operator*=(Vector<T, N>& lhs, const Vector<T, N>& rhs);
121template <typename T, size_t N> Vector<T, N>& operator/=(Vector<T, N>& lhs, const Vector<T, N>& rhs);
122
123// Dot product.
124template <typename T, size_t N> T dot(const Vector<T, N>& lhs, const Vector<T, N>& rhs);
125
126// Vector square norm, norm and normalization.
127template <typename T, size_t N> T square_norm(const Vector<T, N>& v);
128template <typename T, size_t N> T norm(const Vector<T, N>& v);
129template <typename T, size_t N> Vector<T, N> normalize(const Vector<T, N>& v);
130template <typename T, size_t N> Vector<T, N> safe_normalize(const Vector<T, N>& v, const Vector<T, N>& fallback);
131template <typename T, size_t N> Vector<T, N> safe_normalize(const Vector<T, N>& v);
132
133// Bring the norm of a nearly-unit vector closer to 1 by performing a single Newton-Raphson step.
134template <typename T, size_t N>
135Vector<T, N> improve_normalization(const Vector<T, N>& v);
136
137// Bring the norm of a nearly-unit vector closer to 1 by performing a set number of Newton-Raphson steps.
138template <size_t Steps, typename T, size_t N>
139Vector<T, N> improve_normalization(const Vector<T, N>& v);
140
141// Return true if a vector is normalized (unit-length), false otherwise.
142template <typename T, size_t N> bool is_normalized(const Vector<T, N>& v);
143template <typename T, size_t N> bool is_normalized(const Vector<T, N>& v, const T eps);
144
145// Return n if dot(n, i) < 0, -n otherwise.
146template <typename T, size_t N>
147Vector<T, N> faceforward(
148 const Vector<T, N>& n,
149 const Vector<T, N>& i);
150
151// Return v if it is in the same hemisphere as ref, -v otherwise.
152template <typename T, size_t N>
153Vector<T, N> flip_to_same_hemisphere(
154 const Vector<T, N>& v,
155 const Vector<T, N>& ref);
156
157// Return v with the component along n zeroed.
158template <typename T, size_t N>
159Vector<T, N> project(
160 const Vector<T, N>& v,
161 const Vector<T, N>& n);
162
163// Return the reflection vector given an incoming vector i and a unit-length normal vector n.
164template <typename T, size_t N>
165Vector<T, N> reflect(
166 const Vector<T, N>& i,
167 const Vector<T, N>& n);
168
169//
170// Sets t to the refracted vector given a unit-length incoming vector i, a unit-length normal
171// vector n and a relative index of refraction eta. eta is the ratio of the index of
172// refraction in the volume containing the incoming vector to that of the volume being entered.
173// Returns false in the case of total internal reflection.
174//
175// Geometry:
176//
177// i n eta_i
178// eta = -----
179// ^ ^ eta_t
180// \ |
181// \ |
182// \ |
183// \ | IOR = eta_i
184// \ |
185// \|
186// ---------------+--------------- interface
187// |
188// |
189// | IOR = eta_t
190// |
191// |
192// |
193// v
194//
195// t = refract(i, n, eta_i / eta_t)
196//
197
198template <typename T, size_t N>
199bool refract(
200 const Vector<T, N>& i,
201 const Vector<T, N>& n,
202 const T eta,
203 Vector<T, N>& t);
204
205// Clamp the argument to [min, max].
206template <typename T, size_t N>
207Vector<T, N> clamp(
208 const Vector<T, N>& v,
209 const T min,
210 const T max);
211
212// Clamp the argument to [0,1].
213template <typename T, size_t N> Vector<T, N> saturate(const Vector<T, N>& v);
214
215// Return the smallest or largest signed component of a vector.
216template <typename T, size_t N> T min_value(const Vector<T, N>& v);
217template <typename T, size_t N> T max_value(const Vector<T, N>& v);
218
219// Return the index of the smallest or largest signed component of a vector.
220template <typename T, size_t N> size_t min_index(const Vector<T, N>& v);
221template <typename T, size_t N> size_t max_index(const Vector<T, N>& v);
222
223// Return the index of the smallest or largest component of a vector, in absolute value.
224template <typename T, size_t N> size_t min_abs_index(const Vector<T, N>& v);
225template <typename T, size_t N> size_t max_abs_index(const Vector<T, N>& v);
226
227// Component-wise min/max of two vectors.
228template <typename T, size_t N> Vector<T, N> component_wise_min(const Vector<T, N>& lhs, const Vector<T, N>& rhs);
229template <typename T, size_t N> Vector<T, N> component_wise_max(const Vector<T, N>& lhs, const Vector<T, N>& rhs);
230
231
232//
233// 2-dimensional vector class of arbitrary type.
234//
235
236template <typename T>
237class Vector<T, 2>
238{
239 public:
240 // Types.
241 typedef T ValueType;
242 typedef Vector<T, 2> VectorType;
243
244 // Dimension.
245 static const size_t Dimension = 2;
246
247 // Vector components.
248 ValueType x, y;
249
250 // Constructors.
251 Vector(); // leave all components uninitialized
252 explicit Vector(const ValueType* rhs); // initialize with array of 2 scalars
253 explicit Vector(const ValueType val); // set all components to 'val'
254 Vector( // set individual components
255 const ValueType x,
256 const ValueType y);
257
258 // Construct a vector from another vector of a different type.
259 template <typename U>
260 explicit Vector(const Vector<U, 2>& rhs);
261
262#ifdef APPLESEED_ENABLE_IMATH_INTEROP
263
264 // Implicit construction from an Imath::Vec2.
265 Vector(const Imath::Vec2<T>& rhs);
266
267 // Reinterpret this vector as an Imath::Vec2.
268 operator Imath::Vec2<T>&();
269 operator const Imath::Vec2<T>&() const;
270
271#endif
272
273 // Unchecked array subscripting.
274 ValueType& operator[](const size_t i);
275 const ValueType& operator[](const size_t i) const;
276};
277
278// Determinant of the 2D matrix whose first column is lhs and second column is rhs.
279template <typename T> T det(const Vector<T, 2>& lhs, const Vector<T, 2>& rhs);
280
281
282//
283// 3-dimensional vector class of arbitrary type.
284//
285
286template <typename T>
287class Vector<T, 3>
288{
289 public:
290 // Types.
291 typedef T ValueType;
292 typedef Vector<T, 3> VectorType;
293
294 // Dimension.
295 static const size_t Dimension = 3;
296
297 // Vector components.
298 ValueType x, y, z;
299
300 // Constructors.
301 Vector(); // leave all components uninitialized
302 explicit Vector(const ValueType* rhs); // initialize with array of 3 scalars
303 explicit Vector(const ValueType val); // set all components to 'val'
304 Vector( // set individual components
305 const ValueType x,
306 const ValueType y,
307 const ValueType z);
308
309 // Construct a vector from another vector of a different type.
310 template <typename U>
311 explicit Vector(const Vector<U, 3>& rhs);
312
313#ifdef APPLESEED_ENABLE_IMATH_INTEROP
314
315 // Implicit construction from an Imath::Vec3.
316 Vector(const Imath::Vec3<T>& rhs);
317
318 // Reinterpret this vector as an Imath::Vec3.
319 operator Imath::Vec3<T>&();
320 operator const Imath::Vec3<T>&() const;
321
322#endif
323
324 // Build a unit vector from two angles.
325 static VectorType make_unit_vector(
326 const ValueType theta, // angle with Y basis vector, in radians
327 const ValueType phi); // angle with X basis vector, in radians
328 static VectorType make_unit_vector(
329 const ValueType cos_theta, // cosine of angle with Y basis vector
330 const ValueType sin_theta, // sine of angle with Y basis vector
331 const ValueType cos_phi, // cosine of angle with X basis vector
332 const ValueType sin_phi); // sine of angle with X basis vector
333
334 // Unchecked array subscripting.
335 ValueType& operator[](const size_t i);
336 const ValueType& operator[](const size_t i) const;
337};
338
339// Cross product.
340template <typename T> Vector<T, 3> cross(const Vector<T, 3>& lhs, const Vector<T, 3>& rhs);
341
342
343//
344// 4-dimensional vector class of arbitrary type.
345//
346
347template <typename T>
348class Vector<T, 4>
349{
350 public:
351 // Types.
352 typedef T ValueType;
353 typedef Vector<T, 4> VectorType;
354
355 // Dimension.
356 static const size_t Dimension = 4;
357
358 // Vector components.
359 ValueType x, y, z, w;
360
361 // Constructors.
362 Vector(); // leave all components uninitialized
363 explicit Vector(const ValueType* rhs); // initialize with array of 4 scalars
364 explicit Vector(const ValueType val); // set all components to 'val'
365 Vector( // set individual components
366 const ValueType x,
367 const ValueType y,
368 const ValueType z,
369 const ValueType w);
370
371 // Construct a vector from another vector of a different type.
372 template <typename U>
373 explicit Vector(const Vector<U, 4>& rhs);
374
375 // Unchecked array subscripting.
376 ValueType& operator[](const size_t i);
377 const ValueType& operator[](const size_t i) const;
378};
379
380
381//
382// Full specializations for 1D, 2D, 3D and 4D vectors of type int, float and double.
383//
384
385typedef Vector<int, 1> Vector1i;
386typedef Vector<size_t, 1> Vector1u;
387typedef Vector<float, 1> Vector1f;
388typedef Vector<double, 1> Vector1d;
389
390typedef Vector<int, 2> Vector2i;
391typedef Vector<size_t, 2> Vector2u;
392typedef Vector<float, 2> Vector2f;
393typedef Vector<double, 2> Vector2d;
394
395typedef Vector<int, 3> Vector3i;
396typedef Vector<size_t, 3> Vector3u;
397typedef Vector<float, 3> Vector3f;
398typedef Vector<double, 3> Vector3d;
399
400typedef Vector<int, 4> Vector4i;
401typedef Vector<size_t, 4> Vector4u;
402typedef Vector<float, 4> Vector4f;
403typedef Vector<double, 4> Vector4d;
404
405
406//
407// Utility class to find the Imath equivalent of a given foundation::Vector type.
408//
409
410#ifdef APPLESEED_ENABLE_IMATH_INTEROP
411
412template <typename T, size_t N>
413struct ImathVecEquivalent; // no equivalent in the general case
414
415template <typename T>
416struct ImathVecEquivalent<T, 2>
417{
418 typedef Imath::Vec2<T> Type;
419};
420
421template <typename T>
422struct ImathVecEquivalent<T, 3>
423{
424 typedef Imath::Vec3<T> Type;
425};
426
427#endif
428
429
430//
431// N-dimensional vector implementation.
432//
433
434template <typename T, size_t N>
435inline Vector<T, N>::Vector()
436{
437}
438
439template <typename T, size_t N>
440inline Vector<T, N>::Vector(const ValueType* rhs)
441{
442 assert(rhs);
443
444 for (size_t i = 0; i < N; ++i)
445 m_comp[i] = rhs[i];
446}
447
448template <typename T, size_t N>
449inline Vector<T, N>::Vector(const ValueType val)
450{
451 for (size_t i = 0; i < N; ++i)
452 m_comp[i] = val;
453}
454
455template <typename T, size_t N>
456template <typename U>
457inline Vector<T, N>::Vector(const Vector<U, N>& rhs)
458{
459 for (size_t i = 0; i < N; ++i)
460 m_comp[i] = static_cast<ValueType>(rhs[i]);
461}
462
463template <typename T, size_t N>
464inline T& Vector<T, N>::operator[](const size_t i)
465{
466 assert(i < Dimension);
467 return m_comp[i];
468}
469
470template <typename T, size_t N>
471inline const T& Vector<T, N>::operator[](const size_t i) const
472{
473 assert(i < Dimension);
474 return m_comp[i];
475}
476
477template <typename T, size_t N>
478void PoisonImpl<Vector<T, N> >::do_poison(Vector<T, N>& v)
479{
480 for (size_t i = 0; i < N; ++i)
481 poison(v[i]);
482}
483
484template <typename T, size_t N>
485inline bool operator!=(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
486{
487 for (size_t i = 0; i < N; ++i)
488 {
489 if (lhs[i] != rhs[i])
490 return true;
491 }
492
493 return false;
494}
495
496template <typename T, size_t N>
497inline bool operator==(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
498{
499 return !(lhs != rhs);
500}
501
502template <typename T, size_t N>
503inline bool feq(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
504{
505 for (size_t i = 0; i < N; ++i)
506 {
507 if (!feq(lhs[i], rhs[i]))
508 return false;
509 }
510
511 return true;
512}
513
514template <typename T, size_t N>
515inline bool feq(const Vector<T, N>& lhs, const Vector<T, N>& rhs, const T eps)
516{
517 for (size_t i = 0; i < N; ++i)
518 {
519 if (!feq(lhs[i], rhs[i], eps))
520 return false;
521 }
522
523 return true;
524}
525
526template <typename T, size_t N>
527inline bool fz(const Vector<T, N>& v)
528{
529 for (size_t i = 0; i < N; ++i)
530 {
531 if (!fz(v[i]))
532 return false;
533 }
534
535 return true;
536}
537
538template <typename T, size_t N>
539inline bool fz(const Vector<T, N>& v, const T eps)
540{
541 for (size_t i = 0; i < N; ++i)
542 {
543 if (!fz(v[i], eps))
544 return false;
545 }
546
547 return true;
548}
549
550template <typename T, size_t N>
551inline Vector<T, N> operator+(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
552{
553 Vector<T, N> result;
554
555 for (size_t i = 0; i < N; ++i)
556 result[i] = lhs[i] + rhs[i];
557
558 return result;
559}
560
561template <typename T, size_t N>
562inline Vector<T, N> operator-(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
563{
564 Vector<T, N> result;
565
566 for (size_t i = 0; i < N; ++i)
567 result[i] = lhs[i] - rhs[i];
568
569 return result;
570}
571
572template <typename T, size_t N>
573inline Vector<T, N> operator-(const Vector<T, N>& lhs)
574{
575 Vector<T, N> result;
576
577 for (size_t i = 0; i < N; ++i)
578 result[i] = -lhs[i];
579
580 return result;
581}
582
583template <typename T, size_t N>
584inline Vector<T, N> operator*(const Vector<T, N>& lhs, const T rhs)
585{
586 Vector<T, N> result;
587
588 for (size_t i = 0; i < N; ++i)
589 result[i] = lhs[i] * rhs;
590
591 return result;
592}
593
594template <typename T, size_t N>
595inline Vector<T, N> operator*(const T lhs, const Vector<T, N>& rhs)
596{
597 return rhs * lhs;
598}
599
600template <typename T, size_t N>
601inline Vector<T, N> operator/(const Vector<T, N>& lhs, const T rhs)
602{
603 Vector<T, N> result;
604
605 for (size_t i = 0; i < N; ++i)
606 result[i] = lhs[i] / rhs;
607
608 return result;
609}
610
611template <size_t N>
612inline Vector<float, N> operator/(const Vector<float, N>& lhs, const float rhs)
613{
614 return lhs * (1.0f / rhs);
615}
616
617template <size_t N>
618inline Vector<double, N> operator/(const Vector<double, N>& lhs, const double rhs)
619{
620 return lhs * (1.0 / rhs);
621}
622
623template <size_t N>
624inline Vector<long double, N> operator/(const Vector<long double, N>& lhs, const long double rhs)
625{
626 return lhs * (1.0L / rhs);
627}
628
629template <typename T, size_t N>
630inline Vector<T, N> operator*(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
631{
632 Vector<T, N> result;
633
634 for (size_t i = 0; i < N; ++i)
635 result[i] = lhs[i] * rhs[i];
636
637 return result;
638}
639
640template <typename T, size_t N>
641inline Vector<T, N> operator/(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
642{
643 Vector<T, N> result;
644
645 for (size_t i = 0; i < N; ++i)
646 result[i] = lhs[i] / rhs[i];
647
648 return result;
649}
650
651template <typename T, size_t N>
652inline Vector<T, N>& operator+=(Vector<T, N>& lhs, const Vector<T, N>& rhs)
653{
654 for (size_t i = 0; i < N; ++i)
655 lhs[i] += rhs[i];
656
657 return lhs;
658}
659
660template <typename T, size_t N>
661inline Vector<T, N>& operator-=(Vector<T, N>& lhs, const Vector<T, N>& rhs)
662{
663 for (size_t i = 0; i < N; ++i)
664 lhs[i] -= rhs[i];
665
666 return lhs;
667}
668
669template <typename T, size_t N>
670inline Vector<T, N>& operator*=(Vector<T, N>& lhs, const T rhs)
671{
672 for (size_t i = 0; i < N; ++i)
673 lhs[i] *= rhs;
674
675 return lhs;
676}
677
678template <typename T, size_t N>
679inline Vector<T, N>& operator/=(Vector<T, N>& lhs, const T rhs)
680{
681 for (size_t i = 0; i < N; ++i)
682 lhs[i] /= rhs;
683
684 return lhs;
685}
686
687template <size_t N>
688inline Vector<float, N>& operator/=(Vector<float, N>& lhs, const float rhs)
689{
690 return lhs *= 1.0f / rhs;
691}
692
693template <size_t N>
694inline Vector<double, N>& operator/=(Vector<double, N>& lhs, const double rhs)
695{
696 return lhs *= 1.0 / rhs;
697}
698
699template <size_t N>
700inline Vector<long double, N>& operator/=(Vector<long double, N>& lhs, const long double rhs)
701{
702 return lhs *= 1.0L / rhs;
703}
704
705template <typename T, size_t N>
706inline Vector<T, N>& operator*=(Vector<T, N>& lhs, const Vector<T, N>& rhs)
707{
708 for (size_t i = 0; i < N; ++i)
709 lhs[i] *= rhs[i];
710
711 return lhs;
712}
713
714template <typename T, size_t N>
715inline Vector<T, N>& operator/=(Vector<T, N>& lhs, const Vector<T, N>& rhs)
716{
717 for (size_t i = 0; i < N; ++i)
718 lhs[i] /= rhs[i];
719
720 return lhs;
721}
722
723template <typename T, size_t N>
724inline T dot(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
725{
726 T result = T(0.0);
727
728 for (size_t i = 0; i < N; ++i)
729 result += lhs[i] * rhs[i];
730
731 return result;
732}
733
734template <typename T, size_t N>
735inline T square_norm(const Vector<T, N>& v)
736{
737 return dot(v, v);
738}
739
740template <typename T, size_t N>
741inline T norm(const Vector<T, N>& v)
742{
743 return std::sqrt(square_norm(v));
744}
745
746template <typename T, size_t N>
747inline Vector<T, N> normalize(const Vector<T, N>& v)
748{
749 const T n = norm(v);
750 assert(n > T(0.0));
751 return v / n;
752}
753
754template <typename T, size_t N>
755inline Vector<T, N> safe_normalize(const Vector<T, N>& v, const Vector<T, N>& fallback)
756{
757 assert(is_normalized(fallback));
758 const T n = norm(v);
759 return n > 0.0 ? v / n : fallback;
760}
761
762template <typename T, size_t N>
763inline Vector<T, N> safe_normalize(const Vector<T, N>& v)
764{
765 Vector<T, N> result = v;
766
767 const T n = norm(result);
768
769 if (n > 0.0)
770 result /= n;
771 else result[0] = T(1.0);
772
773 assert(is_normalized(result));
774
775 return result;
776}
777
778template <typename T, size_t N>
779inline Vector<T, N> improve_normalization(const Vector<T, N>& v)
780{
781 return improve_normalization<1, T, N>(v);
782}
783
784template <size_t Steps, typename T, size_t N>
785inline Vector<T, N> improve_normalization(const Vector<T, N>& v)
786{
787 Vector<T, N> result = v;
788
789 for (size_t i = 0; i < Steps; ++i)
790 result *= (T(3.0) - square_norm(result)) * T(0.5);
791
792 return result;
793}
794
795template <typename T, size_t N>
796inline bool is_normalized(const Vector<T, N>& v)
797{
798 return feq(square_norm(v), T(1.0), make_eps<T>(1.0e-4f, 1.0e-5));
799}
800
801template <typename T, size_t N>
802inline bool is_normalized(const Vector<T, N>& v, const T eps)
803{
804 return feq(square_norm(v), T(1.0), eps);
805}
806
807template <typename T, size_t N>
808inline Vector<T, N> faceforward(
809 const Vector<T, N>& n,
810 const Vector<T, N>& i)
811{
812 return dot(n, i) < T(0.0) ? n : -n;
813}
814
815template <typename T, size_t N>
816inline Vector<T, N> flip_to_same_hemisphere(
817 const Vector<T, N>& v,
818 const Vector<T, N>& ref)
819{
820 return dot(v, ref) < T(0.0) ? -v : v;
821}
822
823template <typename T, size_t N>
824inline Vector<T, N> project(
825 const Vector<T, N>& v,
826 const Vector<T, N>& n)
827{
828 return v - dot(v, n) * n;
829}
830
831template <typename T, size_t N>
832inline Vector<T, N> reflect(
833 const Vector<T, N>& i,
834 const Vector<T, N>& n)
835{
836 assert(is_normalized(n));
837
838 const T dot_in = dot(i, n);
839 return (dot_in + dot_in) * n - i;
840}
841
842template <typename T, size_t N>
843inline bool refract(
844 const Vector<T, N>& i,
845 const Vector<T, N>& n,
846 const T eta,
847 Vector<T, N>& t)
848{
849 //
850 // Reference: http://en.wikipedia.org/wiki/Snell's_law
851 //
852
853 assert(is_normalized(i));
854 assert(is_normalized(n));
855
856 const T cos_theta_i = dot(i, n);
857 assert(cos_theta_i >= T(0.0));
858
859 const T sin_theta_i2 = T(1.0) - square(cos_theta_i);
860 const T sin_theta_t2 = sin_theta_i2 * square(eta);
861 const T cos_theta_t2 = T(1.0) - sin_theta_t2;
862
863 if (cos_theta_t2 < T(0.0))
864 {
865 // Total internal reflection.
866 return false;
867 }
868
869 const T cos_theta_t = std::sqrt(cos_theta_t2);
870
871 t = (eta * cos_theta_i - cos_theta_t) * n - eta * i;
872 assert(is_normalized(t));
873
874 return true;
875}
876
877template <typename T, size_t N>
878inline Vector<T, N> clamp(
879 const Vector<T, N>& v,
880 const T min,
881 const T max)
882{
883 Vector<T, N> result;
884
885 for (size_t i = 0; i < N; ++i)
886 result[i] = clamp(v[i], min, max);
887
888 return result;
889}
890
891template <typename T, size_t N>
892inline Vector<T, N> saturate(const Vector<T, N>& v)
893{
894 Vector<T, N> result;
895
896 for (size_t i = 0; i < N; ++i)
897 result[i] = saturate(v[i]);
898
899 return result;
900}
901
902template <typename T, size_t N>
903inline T min_value(const Vector<T, N>& v)
904{
905 T value = v[0];
906
907 for (size_t i = 1; i < N; ++i)
908 {
909 if (value > v[i])
910 value = v[i];
911 }
912
913 return value;
914}
915
916template <typename T, size_t N>
917inline T max_value(const Vector<T, N>& v)
918{
919 T value = v[0];
920
921 for (size_t i = 1; i < N; ++i)
922 {
923 if (value < v[i])
924 value = v[i];
925 }
926
927 return value;
928}
929
930template <typename T, size_t N>
931inline size_t min_index(const Vector<T, N>& v)
932{
933 size_t index = 0;
934 T value = v[0];
935
936 for (size_t i = 1; i < N; ++i)
937 {
938 const T x = v[i];
939
940 if (value > x)
941 {
942 value = x;
943 index = i;
944 }
945 }
946
947 return index;
948}
949
950template <typename T, size_t N>
951inline size_t max_index(const Vector<T, N>& v)
952{
953 size_t index = 0;
954 T value = v[0];
955
956 for (size_t i = 1; i < N; ++i)
957 {
958 const T x = v[i];
959
960 if (value < x)
961 {
962 value = x;
963 index = i;
964 }
965 }
966
967 return index;
968}
969
970template <typename T, size_t N>
971inline size_t min_abs_index(const Vector<T, N>& v)
972{
973 size_t index = 0;
974 T value = std::abs(v[0]);
975
976 for (size_t i = 1; i < N; ++i)
977 {
978 const T x = std::abs(v[i]);
979
980 if (value > x)
981 {
982 value = x;
983 index = i;
984 }
985 }
986
987 return index;
988}
989
990template <typename T, size_t N>
991inline size_t max_abs_index(const Vector<T, N>& v)
992{
993 size_t index = 0;
994 T value = std::abs(v[0]);
995
996 for (size_t i = 1; i < N; ++i)
997 {
998 const T x = std::abs(v[i]);
999
1000 if (value < x)
1001 {
1002 value = x;
1003 index = i;
1004 }
1005 }
1006
1007 return index;
1008}
1009
1010template <typename T, size_t N>
1011inline Vector<T, N> component_wise_min(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
1012{
1013 Vector<T, N> result;
1014
1015 for (size_t i = 0; i < N; ++i)
1016 result[i] = std::min(lhs[i], rhs[i]);
1017
1018 return result;
1019}
1020
1021template <typename T, size_t N>
1022inline Vector<T, N> component_wise_max(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
1023{
1024 Vector<T, N> result;
1025
1026 for (size_t i = 0; i < N; ++i)
1027 result[i] = std::max(lhs[i], rhs[i]);
1028
1029 return result;
1030}
1031
1032
1033//
1034// 2D vector implementation.
1035//
1036
1037template <typename T>
1038inline Vector<T, 2>::Vector()
1039{
1040}
1041
1042template <typename T>
1043inline Vector<T, 2>::Vector(const ValueType* rhs)
1044{
1045 assert(rhs);
1046 x = rhs[0];
1047 y = rhs[1];
1048}
1049
1050template <typename T>
1051inline Vector<T, 2>::Vector(const ValueType val)
1052 : x(val)
1053 , y(val)
1054{
1055}
1056
1057template <typename T>
1058inline Vector<T, 2>::Vector(
1059 const ValueType x_,
1060 const ValueType y_)
1061 : x(x_)
1062 , y(y_)
1063{
1064}
1065
1066template <typename T>
1067template <typename U>
1068inline Vector<T, 2>::Vector(const Vector<U, 2>& rhs)
1069 : x(static_cast<ValueType>(rhs.x))
1070 , y(static_cast<ValueType>(rhs.y))
1071{
1072}
1073
1074#ifdef APPLESEED_ENABLE_IMATH_INTEROP
1075
1076template <typename T>
1077inline Vector<T, 2>::Vector(const Imath::Vec2<T>& rhs)
1078 : x(rhs.x)
1079 , y(rhs.y)
1080{
1081}
1082
1083template <typename T>
1084inline Vector<T, 2>::operator Imath::Vec2<T>&()
1085{
1086 return reinterpret_cast<Imath::Vec2<T>&>(*this);
1087}
1088
1089template <typename T>
1090inline Vector<T, 2>::operator const Imath::Vec2<T>&() const
1091{
1092 return reinterpret_cast<const Imath::Vec2<T>&>(*this);
1093}
1094
1095#endif
1096
1097template <typename T>
1098inline T& Vector<T, 2>::operator[](const size_t i)
1099{
1100 assert(i < Dimension);
1101 return (&x)[i];
1102}
1103
1104template <typename T>
1105inline const T& Vector<T, 2>::operator[](const size_t i) const
1106{
1107 assert(i < Dimension);
1108 return (&x)[i];
1109}
1110
1111template <typename T>
1112inline T det(const Vector<T, 2>& lhs, const Vector<T, 2>& rhs)
1113{
1114 return lhs.x * rhs.y - lhs.y * rhs.x;
1115}
1116
1117
1118//
1119// 3D vector implementation.
1120//
1121
1122template <typename T>
1123inline Vector<T, 3>::Vector()
1124{
1125}
1126
1127template <typename T>
1128inline Vector<T, 3>::Vector(const ValueType* rhs)
1129{
1130 assert(rhs);
1131 x = rhs[0];
1132 y = rhs[1];
1133 z = rhs[2];
1134}
1135
1136template <typename T>
1137inline Vector<T, 3>::Vector(const ValueType val)
1138 : x(val)
1139 , y(val)
1140 , z(val)
1141{
1142}
1143
1144template <typename T>
1145inline Vector<T, 3>::Vector(
1146 const ValueType x_,
1147 const ValueType y_,
1148 const ValueType z_)
1149 : x(x_)
1150 , y(y_)
1151 , z(z_)
1152{
1153}
1154
1155template <typename T>
1156template <typename U>
1157inline Vector<T, 3>::Vector(const Vector<U, 3>& rhs)
1158 : x(static_cast<ValueType>(rhs.x))
1159 , y(static_cast<ValueType>(rhs.y))
1160 , z(static_cast<ValueType>(rhs.z))
1161{
1162}
1163
1164#ifdef APPLESEED_ENABLE_IMATH_INTEROP
1165
1166template <typename T>
1167inline Vector<T, 3>::Vector(const Imath::Vec3<T>& rhs)
1168 : x(rhs.x)
1169 , y(rhs.y)
1170 , z(rhs.z)
1171{
1172}
1173
1174template <typename T>
1175inline Vector<T, 3>::operator Imath::Vec3<T>&()
1176{
1177 return reinterpret_cast<Imath::Vec3<T>&>(*this);
1178}
1179
1180template <typename T>
1181inline Vector<T, 3>::operator const Imath::Vec3<T>&() const
1182{
1183 return reinterpret_cast<const Imath::Vec3<T>&>(*this);
1184}
1185
1186#endif
1187
1188template <typename T>
1189inline Vector<T, 3> Vector<T, 3>::make_unit_vector(
1190 const ValueType theta,
1191 const ValueType phi)
1192{
1193 const ValueType sin_theta = std::sin(theta);
1194
1195 return
1196 Vector<T, 3>(
1197 std::cos(phi) * sin_theta,
1198 std::cos(theta),
1199 std::sin(phi) * sin_theta);
1200}
1201
1202template <typename T>
1203inline Vector<T, 3> Vector<T, 3>::make_unit_vector(
1204 const ValueType cos_theta,
1205 const ValueType sin_theta,
1206 const ValueType cos_phi,
1207 const ValueType sin_phi)
1208{
1209 return
1210 Vector<T, 3>(
1211 cos_phi * sin_theta,
1212 cos_theta,
1213 sin_phi * sin_theta);
1214}
1215
1216template <typename T>
1217inline T& Vector<T, 3>::operator[](const size_t i)
1218{
1219 assert(i < Dimension);
1220 return (&x)[i];
1221}
1222
1223template <typename T>
1224inline const T& Vector<T, 3>::operator[](const size_t i) const
1225{
1226 assert(i < Dimension);
1227 return (&x)[i];
1228}
1229
1230template <typename T>
1231inline Vector<T, 3> cross(const Vector<T, 3>& lhs, const Vector<T, 3>& rhs)
1232{
1233 Vector<T, 3> result;
1234 result.x = lhs.y * rhs.z - rhs.y * lhs.z;
1235 result.y = lhs.z * rhs.x - rhs.z * lhs.x;
1236 result.z = lhs.x * rhs.y - rhs.x * lhs.y;
1237 return result;
1238}
1239
1240
1241//
1242// 4D vector implementation.
1243//
1244
1245template <typename T>
1246inline Vector<T, 4>::Vector()
1247{
1248}
1249
1250template <typename T>
1251inline Vector<T, 4>::Vector(const ValueType* rhs)
1252{
1253 assert(rhs);
1254 x = rhs[0];
1255 y = rhs[1];
1256 z = rhs[2];
1257 w = rhs[3];
1258}
1259
1260template <typename T>
1261inline Vector<T, 4>::Vector(const ValueType val)
1262 : x(val)
1263 , y(val)
1264 , z(val)
1265 , w(val)
1266{
1267}
1268
1269template <typename T>
1270inline Vector<T, 4>::Vector(
1271 const ValueType x_,
1272 const ValueType y_,
1273 const ValueType z_,
1274 const ValueType w_)
1275 : x(x_)
1276 , y(y_)
1277 , z(z_)
1278 , w(w_)
1279{
1280}
1281
1282template <typename T>
1283template <typename U>
1284inline Vector<T, 4>::Vector(const Vector<U, 4>& rhs)
1285 : x(static_cast<ValueType>(rhs.x))
1286 , y(static_cast<ValueType>(rhs.y))
1287 , z(static_cast<ValueType>(rhs.z))
1288 , w(static_cast<ValueType>(rhs.w))
1289{
1290}
1291
1292template <typename T>
1293inline T& Vector<T, 4>::operator[](const size_t i)
1294{
1295 assert(i < Dimension);
1296 return (&x)[i];
1297}
1298
1299template <typename T>
1300inline const T& Vector<T, 4>::operator[](const size_t i) const
1301{
1302 assert(i < Dimension);
1303 return (&x)[i];
1304}
1305
1306} // namespace foundation
1307
1308#endif // !APPLESEED_FOUNDATION_MATH_VECTOR_H
1309