Cleanup: Subdiv: Remove common_ prefix
[blender.git] / source / blender / blenlib / BLI_math_vector.hh
blobc7d9d8f2dc4aacefde28559779de83ecabac14be
1 /* SPDX-FileCopyrightText: 2022 Blender Authors
3 * SPDX-License-Identifier: GPL-2.0-or-later */
5 #pragma once
7 /** \file
8 * \ingroup bli
9 */
11 #include <type_traits>
13 #include "BLI_math_base.hh"
14 #include "BLI_math_vector_types.hh"
15 #include "BLI_span.hh"
16 #include "BLI_utildefines.h"
18 namespace blender::math {
20 /**
21 * Returns true if the given vectors are equal within the given epsilon.
22 * The epsilon is scaled for each component by magnitude of the matching component of `a`.
24 template<typename T, int Size>
25 [[nodiscard]] inline bool almost_equal_relative(const VecBase<T, Size> &a,
26 const VecBase<T, Size> &b,
27 const T &epsilon_factor)
29 for (int i = 0; i < Size; i++) {
30 const float epsilon = epsilon_factor * math::abs(a[i]);
31 if (math::distance(a[i], b[i]) > epsilon) {
32 return false;
35 return true;
38 template<typename T, int Size> [[nodiscard]] inline VecBase<T, Size> abs(const VecBase<T, Size> &a)
40 VecBase<T, Size> result;
41 for (int i = 0; i < Size; i++) {
42 result[i] = a[i] >= 0 ? a[i] : -a[i];
44 return result;
47 /**
48 * Returns -1 if \a a is less than 0, 0 if \a a is equal to 0, and +1 if \a a is greater than 0.
50 template<typename T, int Size>
51 [[nodiscard]] inline VecBase<T, Size> sign(const VecBase<T, Size> &a)
53 BLI_UNROLL_MATH_VEC_OP_VEC(math::sign, a);
56 template<typename T, int Size>
57 [[nodiscard]] inline VecBase<T, Size> min(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
59 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::min, a, b);
62 template<typename T, int Size>
63 [[nodiscard]] inline VecBase<T, Size> max(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
65 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::max, a, b);
68 template<typename T, int Size>
69 [[nodiscard]] inline VecBase<T, Size> clamp(const VecBase<T, Size> &a,
70 const VecBase<T, Size> &min,
71 const VecBase<T, Size> &max)
73 VecBase<T, Size> result = a;
74 for (int i = 0; i < Size; i++) {
75 result[i] = math::clamp(result[i], min[i], max[i]);
77 return result;
80 template<typename T, int Size>
81 [[nodiscard]] inline VecBase<T, Size> clamp(const VecBase<T, Size> &a, const T &min, const T &max)
83 VecBase<T, Size> result = a;
84 for (int i = 0; i < Size; i++) {
85 result[i] = math::clamp(result[i], min, max);
87 return result;
90 template<typename T, int Size>
91 [[nodiscard]] inline VecBase<T, Size> step(const VecBase<T, Size> &edge,
92 const VecBase<T, Size> &value)
94 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::step, edge, value);
97 template<typename T, int Size>
98 [[nodiscard]] inline VecBase<T, Size> step(const T &edge, const VecBase<T, Size> &value)
100 VecBase<T, Size> result = value;
101 for (int i = 0; i < Size; i++) {
102 result[i] = math::step(edge, result[i]);
104 return result;
107 template<typename T, int Size>
108 [[nodiscard]] inline VecBase<T, Size> mod(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
110 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::mod, a, b);
113 template<typename T, int Size>
114 [[nodiscard]] inline VecBase<T, Size> mod(const VecBase<T, Size> &a, const T &b)
116 BLI_assert(b != 0);
117 VecBase<T, Size> result;
118 for (int i = 0; i < Size; i++) {
119 result[i] = math::mod(a[i], b);
121 return result;
125 * Safe version of mod(a, b) that returns 0 if b is equal to 0.
127 template<typename T, int Size>
128 [[nodiscard]] inline VecBase<T, Size> safe_mod(const VecBase<T, Size> &a,
129 const VecBase<T, Size> &b)
131 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::safe_mod, a, b);
135 * Safe version of mod(a, b) that returns 0 if b is equal to 0.
137 template<typename T, int Size>
138 [[nodiscard]] inline VecBase<T, Size> safe_mod(const VecBase<T, Size> &a, const T &b)
140 if (b == 0) {
141 return VecBase<T, Size>(0);
143 VecBase<T, Size> result;
144 for (int i = 0; i < Size; i++) {
145 result[i] = math::mod(a[i], b);
147 return result;
151 * Return the value of x raised to the y power.
152 * The result is undefined if x < 0 or if x = 0 and y <= 0.
154 template<typename T, int Size>
155 [[nodiscard]] inline VecBase<T, Size> pow(const VecBase<T, Size> &x, const T &y)
157 VecBase<T, Size> result;
158 for (int i = 0; i < Size; i++) {
159 result[i] = math::pow(x[i], y);
161 return result;
165 * Return the value of x raised to the y power.
166 * The result is x if x < 0 or if x = 0 and y <= 0.
168 template<typename T, int Size>
169 [[nodiscard]] inline VecBase<T, Size> safe_pow(const VecBase<T, Size> &x, const T &y)
171 VecBase<T, Size> result;
172 for (int i = 0; i < Size; i++) {
173 result[i] = math::safe_pow(x[i], y);
175 return result;
179 * Return the value of x raised to the y power.
180 * The result is the given fallback if x < 0 or if x = 0 and y <= 0.
182 template<typename T, int Size>
183 [[nodiscard]] inline VecBase<T, Size> fallback_pow(const VecBase<T, Size> &x,
184 const T &y,
185 const VecBase<T, Size> &fallback)
187 VecBase<T, Size> result;
188 for (int i = 0; i < Size; i++) {
189 result[i] = math::fallback_pow(x[i], y, fallback[i]);
191 return result;
195 * Return the value of x raised to the y power.
196 * The result is undefined if x < 0 or if x = 0 and y <= 0.
198 template<typename T, int Size>
199 [[nodiscard]] inline VecBase<T, Size> pow(const VecBase<T, Size> &x, const VecBase<T, Size> &y)
201 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::pow, x, y);
204 /** Per-element square. */
205 template<typename T, int Size>
206 [[nodiscard]] inline VecBase<T, Size> square(const VecBase<T, Size> &a)
208 BLI_UNROLL_MATH_VEC_OP_VEC(math::square, a);
211 /* Per-element exponent. */
212 template<typename T, int Size> [[nodiscard]] inline VecBase<T, Size> exp(const VecBase<T, Size> &x)
214 BLI_UNROLL_MATH_VEC_OP_VEC(math::exp, x);
218 * Returns \a a if it is a multiple of \a b or the next multiple or \a b after \b a .
219 * In other words, it is equivalent to `divide_ceil(a, b) * b`.
220 * It is undefined if \a a is negative or \b b is not strictly positive.
222 template<typename T, int Size>
223 [[nodiscard]] inline VecBase<T, Size> ceil_to_multiple(const VecBase<T, Size> &a,
224 const VecBase<T, Size> &b)
226 VecBase<T, Size> result;
227 for (int i = 0; i < Size; i++) {
228 BLI_assert(a[i] >= 0);
229 BLI_assert(b[i] > 0);
230 result[i] = ((a[i] + b[i] - 1) / b[i]) * b[i];
232 return result;
236 * Integer division that returns the ceiling, instead of flooring like normal C division.
237 * It is undefined if \a a is negative or \b b is not strictly positive.
239 template<typename T, int Size>
240 [[nodiscard]] inline VecBase<T, Size> divide_ceil(const VecBase<T, Size> &a,
241 const VecBase<T, Size> &b)
243 VecBase<T, Size> result;
244 for (int i = 0; i < Size; i++) {
245 BLI_assert(a[i] >= 0);
246 BLI_assert(b[i] > 0);
247 result[i] = (a[i] + b[i] - 1) / b[i];
249 return result;
252 template<typename T, int Size>
253 inline void min_max(const VecBase<T, Size> &vector, VecBase<T, Size> &min, VecBase<T, Size> &max)
255 min = math::min(vector, min);
256 max = math::max(vector, max);
260 * Returns 0 if denominator is exactly equal to 0.
262 template<typename T, int Size>
263 [[nodiscard]] inline VecBase<T, Size> safe_divide(const VecBase<T, Size> &a,
264 const VecBase<T, Size> &b)
266 BLI_UNROLL_MATH_VEC_FUNC_VEC_VEC(math::safe_divide, a, b);
270 * Returns 0 if denominator is exactly equal to 0.
272 template<typename T, int Size>
273 [[nodiscard]] inline VecBase<T, Size> safe_divide(const VecBase<T, Size> &a, const T &b)
275 return (b != 0) ? a / b : VecBase<T, Size>(0.0f);
278 template<typename T, int Size>
279 [[nodiscard]] inline VecBase<T, Size> floor(const VecBase<T, Size> &a)
281 BLI_UNROLL_MATH_VEC_OP_VEC(math::floor, a);
284 template<typename T, int Size>
285 [[nodiscard]] inline VecBase<T, Size> round(const VecBase<T, Size> &a)
287 BLI_UNROLL_MATH_VEC_OP_VEC(math::round, a);
290 template<typename T, int Size>
291 [[nodiscard]] inline VecBase<T, Size> ceil(const VecBase<T, Size> &a)
293 BLI_UNROLL_MATH_VEC_OP_VEC(math::ceil, a);
297 * Per-element square root.
298 * Negative elements are evaluated to NaN.
300 template<typename T, int Size>
301 [[nodiscard]] inline VecBase<T, Size> sqrt(const VecBase<T, Size> &a)
303 BLI_UNROLL_MATH_VEC_OP_VEC(math::sqrt, a);
307 * Per-element square root.
308 * Negative elements are evaluated to zero.
310 template<typename T, int Size>
311 [[nodiscard]] inline VecBase<T, Size> safe_sqrt(const VecBase<T, Size> &a)
313 VecBase<T, Size> result;
314 for (int i = 0; i < Size; i++) {
315 result[i] = a[i] >= T(0) ? math ::sqrt(a[i]) : T(0);
317 return result;
321 * Per-element inverse.
322 * Zero elements are evaluated to NaN.
324 template<typename T, int Size> [[nodiscard]] inline VecBase<T, Size> rcp(const VecBase<T, Size> &a)
326 BLI_UNROLL_MATH_VEC_OP_VEC(math::rcp, a);
330 * Per-element inverse.
331 * Zero elements are evaluated to zero.
333 template<typename T, int Size>
334 [[nodiscard]] inline VecBase<T, Size> safe_rcp(const VecBase<T, Size> &a)
336 BLI_UNROLL_MATH_VEC_OP_VEC(math::safe_rcp, a);
339 template<typename T, int Size>
340 [[nodiscard]] inline VecBase<T, Size> fract(const VecBase<T, Size> &a)
342 BLI_UNROLL_MATH_VEC_OP_VEC(math::fract, a);
346 * Dot product between two vectors.
347 * Equivalent to component wise multiplication followed by summation of the result.
348 * Equivalent to the cosine of the angle between the two vectors if the vectors are normalized.
349 * \note prefer using `length_manhattan(a)` than `dot(a, vec(1))` to get the sum of all components.
351 template<typename T, int Size>
352 [[nodiscard]] inline T dot(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
354 T result = a[0] * b[0];
355 for (int i = 1; i < Size; i++) {
356 result += a[i] * b[i];
358 return result;
362 * Returns summation of all components.
364 template<typename T, int Size> [[nodiscard]] inline T length_manhattan(const VecBase<T, Size> &a)
366 T result = math::abs(a[0]);
367 for (int i = 1; i < Size; i++) {
368 result += math::abs(a[i]);
370 return result;
373 template<typename T, int Size> [[nodiscard]] inline T length_squared(const VecBase<T, Size> &a)
375 return dot(a, a);
378 template<typename T, int Size> [[nodiscard]] inline T length(const VecBase<T, Size> &a)
380 return math::sqrt(length_squared(a));
383 /** Return true if each individual column is unit scaled. Mainly for assert usage. */
384 template<typename T, int Size> [[nodiscard]] inline bool is_unit_scale(const VecBase<T, Size> &v)
386 /* Checks are flipped so NAN doesn't assert because we're making sure the value was
387 * normalized and in the case we don't want NAN to be raising asserts since there
388 * is nothing to be done in that case. */
389 const T test_unit = math::length_squared(v);
390 return (!(math::abs(test_unit - T(1)) >= AssertUnitEpsilon<T>::value) ||
391 !(math::abs(test_unit) >= AssertUnitEpsilon<T>::value));
394 template<typename T, int Size>
395 [[nodiscard]] inline T distance_manhattan(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
397 return length_manhattan(a - b);
400 template<typename T, int Size>
401 [[nodiscard]] inline T distance_squared(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
403 return length_squared(a - b);
406 template<typename T, int Size>
407 [[nodiscard]] inline T distance(const VecBase<T, Size> &a, const VecBase<T, Size> &b)
409 return length(a - b);
412 template<typename T, int Size>
413 [[nodiscard]] inline VecBase<T, Size> reflect(const VecBase<T, Size> &incident,
414 const VecBase<T, Size> &normal)
416 BLI_assert(is_unit_scale(normal));
417 return incident - 2.0 * dot(normal, incident) * normal;
420 template<typename T, int Size>
421 [[nodiscard]] inline VecBase<T, Size> refract(const VecBase<T, Size> &incident,
422 const VecBase<T, Size> &normal,
423 const T &eta)
425 float dot_ni = dot(normal, incident);
426 float k = 1.0f - eta * eta * (1.0f - dot_ni * dot_ni);
427 if (k < 0.0f) {
428 return VecBase<T, Size>(0.0f);
430 return eta * incident - (eta * dot_ni + sqrt(k)) * normal;
434 * Project \a p onto \a v_proj .
435 * Returns zero vector if \a v_proj is degenerate (zero length).
437 template<typename T, int Size>
438 [[nodiscard]] inline VecBase<T, Size> project(const VecBase<T, Size> &p,
439 const VecBase<T, Size> &v_proj)
441 if (UNLIKELY(is_zero(v_proj))) {
442 return VecBase<T, Size>(0.0f);
444 return v_proj * (dot(p, v_proj) / dot(v_proj, v_proj));
447 template<typename T, int Size>
448 [[nodiscard]] inline VecBase<T, Size> normalize_and_get_length(const VecBase<T, Size> &v,
449 T &out_length)
451 out_length = length_squared(v);
452 /* A larger value causes normalize errors in a scaled down models with camera extreme close. */
453 constexpr T threshold = std::is_same_v<T, double> ? 1.0e-70 : 1.0e-35f;
454 if (out_length > threshold) {
455 out_length = sqrt(out_length);
456 return v / out_length;
458 /* Either the vector is small or one of it's values contained `nan`. */
459 out_length = 0.0;
460 return VecBase<T, Size>(0.0);
463 template<typename T, int Size>
464 [[nodiscard]] inline VecBase<T, Size> normalize(const VecBase<T, Size> &v)
466 T len;
467 return normalize_and_get_length(v, len);
471 * \return cross perpendicular vector to \a a and \a b.
472 * \note Return zero vector if \a a and \a b are collinear.
473 * \note The length of the resulting vector is equal to twice the area of the triangle between \a a
474 * and \a b ; and it is equal to the sine of the angle between \a a and \a b if they are
475 * normalized.
476 * \note Blender 3D space uses right handedness: \a a = thumb, \a b = index, return = middle.
478 template<typename T>
479 [[nodiscard]] inline VecBase<T, 3> cross(const VecBase<T, 3> &a, const VecBase<T, 3> &b)
481 return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
485 * Same as `cross(a, b)` but use double float precision for the computation.
487 [[nodiscard]] inline VecBase<float, 3> cross_high_precision(const VecBase<float, 3> &a,
488 const VecBase<float, 3> &b)
490 return {float(double(a.y) * double(b.z) - double(a.z) * double(b.y)),
491 float(double(a.z) * double(b.x) - double(a.x) * double(b.z)),
492 float(double(a.x) * double(b.y) - double(a.y) * double(b.x))};
496 * \param poly: Array of points around a polygon. They don't have to be co-planar.
497 * \return Best fit plane normal for the given polygon loop or zero vector if point
498 * array is too short. Not normalized.
500 template<typename T> [[nodiscard]] inline VecBase<T, 3> cross_poly(Span<VecBase<T, 3>> poly)
502 /* Newell's Method. */
503 int nv = int(poly.size());
504 if (nv < 3) {
505 return VecBase<T, 3>(0, 0, 0);
507 const VecBase<T, 3> *v_prev = &poly[nv - 1];
508 const VecBase<T, 3> *v_curr = &poly[0];
509 VecBase<T, 3> n(0, 0, 0);
510 for (int i = 0; i < nv;) {
511 n[0] = n[0] + ((*v_prev)[1] - (*v_curr)[1]) * ((*v_prev)[2] + (*v_curr)[2]);
512 n[1] = n[1] + ((*v_prev)[2] - (*v_curr)[2]) * ((*v_prev)[0] + (*v_curr)[0]);
513 n[2] = n[2] + ((*v_prev)[0] - (*v_curr)[0]) * ((*v_prev)[1] + (*v_curr)[1]);
514 v_prev = v_curr;
515 ++i;
516 if (i < nv) {
517 v_curr = &poly[i];
520 return n;
524 * Return normal vector to a triangle.
525 * The result is not normalized and can be degenerate.
527 template<typename T>
528 [[nodiscard]] inline VecBase<T, 3> cross_tri(const VecBase<T, 3> &v1,
529 const VecBase<T, 3> &v2,
530 const VecBase<T, 3> &v3)
532 return cross(v1 - v2, v2 - v3);
536 * Return normal vector to a triangle.
537 * The result is normalized but can still be degenerate.
539 template<typename T>
540 [[nodiscard]] inline VecBase<T, 3> normal_tri(const VecBase<T, 3> &v1,
541 const VecBase<T, 3> &v2,
542 const VecBase<T, 3> &v3)
544 return normalize(cross_tri(v1, v2, v3));
548 * Per component linear interpolation.
549 * \param t: interpolation factor. Return \a a if equal 0. Return \a b if equal 1.
550 * Outside of [0..1] range, use linear extrapolation.
552 template<typename T, typename FactorT, int Size>
553 [[nodiscard]] inline VecBase<T, Size> interpolate(const VecBase<T, Size> &a,
554 const VecBase<T, Size> &b,
555 const FactorT &t)
557 return a * (1 - t) + b * t;
561 * \return Point halfway between \a a and \a b.
563 template<typename T, int Size>
564 [[nodiscard]] inline VecBase<T, Size> midpoint(const VecBase<T, Size> &a,
565 const VecBase<T, Size> &b)
567 return (a + b) * 0.5;
571 * Return `vector` if `incident` and `reference` are pointing in the same direction.
573 template<typename T, int Size>
574 [[nodiscard]] inline VecBase<T, Size> faceforward(const VecBase<T, Size> &vector,
575 const VecBase<T, Size> &incident,
576 const VecBase<T, Size> &reference)
578 return (dot(reference, incident) < 0) ? vector : -vector;
582 * \return Index of the component with the greatest magnitude.
584 template<typename T> [[nodiscard]] inline int dominant_axis(const VecBase<T, 3> &a)
586 VecBase<T, 3> b = abs(a);
587 return ((b.x > b.y) ? ((b.x > b.z) ? 0 : 2) : ((b.y > b.z) ? 1 : 2));
591 * \return the maximum component of a vector.
593 template<typename T, int Size> [[nodiscard]] inline T reduce_max(const VecBase<T, Size> &a)
595 T result = a[0];
596 for (int i = 1; i < Size; i++) {
597 if (a[i] > result) {
598 result = a[i];
601 return result;
605 * \return the minimum component of a vector.
607 template<typename T, int Size> [[nodiscard]] inline T reduce_min(const VecBase<T, Size> &a)
609 T result = a[0];
610 for (int i = 1; i < Size; i++) {
611 if (a[i] < result) {
612 result = a[i];
615 return result;
619 * \return the sum of the components of a vector.
621 template<typename T, int Size> [[nodiscard]] inline T reduce_add(const VecBase<T, Size> &a)
623 T result = a[0];
624 for (int i = 1; i < Size; i++) {
625 result += a[i];
627 return result;
631 * \return the product of the components of a vector.
633 template<typename T, int Size> [[nodiscard]] inline T reduce_mul(const VecBase<T, Size> &a)
635 T result = a[0];
636 for (int i = 1; i < Size; i++) {
637 result *= a[i];
639 return result;
643 * \return the average of the components of a vector.
645 template<typename T, int Size> [[nodiscard]] inline T average(const VecBase<T, Size> &a)
647 return reduce_add(a) * (T(1) / T(Size));
651 * Calculates a perpendicular vector to \a v.
652 * \note Returned vector can be in any perpendicular direction.
653 * \note Returned vector might not the same length as \a v.
655 template<typename T> [[nodiscard]] inline VecBase<T, 3> orthogonal(const VecBase<T, 3> &v)
657 const int axis = dominant_axis(v);
658 switch (axis) {
659 case 0:
660 return {-v.y - v.z, v.x, v.x};
661 case 1:
662 return {v.y, -v.x - v.z, v.y};
663 case 2:
664 return {v.z, v.z, -v.x - v.y};
666 return v;
670 * Calculates a perpendicular vector to \a v.
671 * \note Returned vector can be in any perpendicular direction.
673 template<typename T> [[nodiscard]] inline VecBase<T, 2> orthogonal(const VecBase<T, 2> &v)
675 return {-v.y, v.x};
679 * Returns true if vectors are equal within the given epsilon.
681 template<typename T, int Size>
682 [[nodiscard]] inline bool is_equal(const VecBase<T, Size> &a,
683 const VecBase<T, Size> &b,
684 const T epsilon = T(0))
686 for (int i = 0; i < Size; i++) {
687 if (math::abs(a[i] - b[i]) > epsilon) {
688 return false;
691 return true;
695 * Return true if the absolute values of all components are smaller than given epsilon (0 by
696 * default).
698 * \note Does not compute the actual length of the vector, for performance.
700 template<typename T, int Size>
701 [[nodiscard]] inline bool is_zero(const VecBase<T, Size> &a, const T epsilon = T(0))
703 for (int i = 0; i < Size; i++) {
704 if (math::abs(a[i]) > epsilon) {
705 return false;
708 return true;
712 * Returns true if at least one component is exactly equal to 0.
714 template<typename T, int Size> [[nodiscard]] inline bool is_any_zero(const VecBase<T, Size> &a)
716 for (int i = 0; i < Size; i++) {
717 if (a[i] == T(0)) {
718 return true;
721 return false;
725 * Return true if the squared length of the vector is (almost) equal to 1 (with a
726 * `10 * std::numeric_limits<T>::epsilon()` epsilon error by default).
728 template<typename T, int Size>
729 [[nodiscard]] inline bool is_unit(const VecBase<T, Size> &a,
730 const T epsilon = T(10) * std::numeric_limits<T>::epsilon())
732 const T length = length_squared(a);
733 return math::abs(length - T(1)) <= epsilon;
736 /** Intersections. */
738 template<typename T> struct isect_result {
739 enum {
740 LINE_LINE_COLINEAR = -1,
741 LINE_LINE_NONE = 0,
742 LINE_LINE_EXACT = 1,
743 LINE_LINE_CROSS = 2,
744 } kind;
745 typename T::base_type lambda;
748 template<typename T, int Size>
749 [[nodiscard]] isect_result<VecBase<T, Size>> isect_seg_seg(const VecBase<T, Size> &v1,
750 const VecBase<T, Size> &v2,
751 const VecBase<T, Size> &v3,
752 const VecBase<T, Size> &v4);
754 } // namespace blender::math