BLI: Refactor vector types & functions to use templates

This patch implements the vector types (i.e:float2) by making heavy
usage of templating. All vector functions are now outside of the vector
classes (inside the blender::math namespace) and are not vector size
dependent for the most part.

In the ongoing effort to make shaders less GL centric, we are aiming
to share more code between GLSL and C++ to avoid code duplication.

Motivations:
- We are aiming to share UBO and SSBO structures between GLSL and C++.
  This means we will use many of the existing vector types and others we
  currently don't have (uintX, intX). All these variations were asking
  for many more code duplication.
- Deduplicate existing code which is duplicated for each vector size.
- We also want to share small functions. Which means that vector functions
  should be static and not in the class namespace.
- Reduce friction to use these types in new projects due to their
  incompleteness.
- The current state of the BLI_(float|double|mpq)(2|3|4).hh is a bit of a
  let down. Most clases are incomplete, out of sync with each others with
  different codestyles, and some functions that should be static are not
  (i.e: float3::reflect()).

Upsides:
- Still support .x, .y, .z, .w for readability.
- Compact, readable and easilly extendable.
- All of the vector functions are available for all the vectors types and
  can be restricted to certain types. Also template specialization let us
  define exception for special class (like mpq).
- With optimization ON, the compiler unroll the loops and performance is
  the same.

Downsides:
- Might impact debugability. Though I would arge that the bugs are rarelly
  caused by the vector class itself (since the operations are quite trivial)
  but by the type conversions.
- Might impact compile time. I did not saw a significant impact since the
  usage is not really widespread.
- Functions needs to be rewritten to support arbitrary vector length. For
  instance, one can't call len_squared_v3v3 in math::length_squared() and
  call it a day.
- Type cast does not work with the template version of the math:: vector
  functions. Meaning you need to manually cast float * and (float *)[3] to
  float3 for the function calls.
  i.e: math::distance_squared(float3(nearest.co), positions[i]);
- Some parts might loose in readability:
  float3::dot(v1.normalized(), v2.normalized())
  becoming
  math::dot(math::normalize(v1), math::normalize(v2))
  But I propose, when appropriate, to use
  using namespace blender::math; on function local or file scope to
  increase readability. dot(normalize(v1), normalize(v2))

Consideration:
- Include back .length() method. It is quite handy and is more C++
  oriented.
- I considered the GLM library as a candidate for replacement.
  It felt like too much for what we need and would be difficult to
  extend / modify to our needs.
- I used Macros to reduce code in operators declaration and potential
  copy paste bugs. This could reduce debugability and could be reverted.
- This touches delaunay_2d.cc and the intersection code. I would like to
  know @Howard Trickey (howardt) opinion on the matter.
- The noexcept on the copy constructor of mpq(2|3) is being removed.
  But according to @Jacques Lucke (JacquesLucke) it is not a real problem
  for now.

I would like to give a huge thanks to @Jacques Lucke (JacquesLucke) who
helped during this and pushed me to reduce the duplication further.

Reviewed By: brecht, sergey, JacquesLucke

Differential Revision: http://developer.blender.org/D13791

1 files changed, 43 insertions, 45 deletions

diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index 1f150137ba3..982759ffcff 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -30,15 +30,13 @@
 #  include "BLI_array.hh"
 #  include "BLI_assert.h"
 #  include "BLI_delaunay_2d.h"
-#  include "BLI_double3.hh"
-#  include "BLI_float3.hh"
 #  include "BLI_hash.hh"
 #  include "BLI_kdopbvh.h"
 #  include "BLI_map.hh"
 #  include "BLI_math_boolean.hh"
 #  include "BLI_math_mpq.hh"
-#  include "BLI_mpq2.hh"
-#  include "BLI_mpq3.hh"
+#  include "BLI_math_vec_mpq_types.hh"
+#  include "BLI_math_vec_types.hh"
 #  include "BLI_polyfill_2d.h"
 #  include "BLI_set.hh"
 #  include "BLI_span.hh"
@@ -198,14 +196,14 @@ void Face::populate_plane(bool need_exact)
       for (int i : index_range()) {
         co[i] = vert[i]->co_exact;
       }
-      normal_exact = mpq3::cross_poly(co);
+      normal_exact = math::cross_poly(co.as_span());
     }
     else {
       mpq3 tr02 = vert[0]->co_exact - vert[2]->co_exact;
       mpq3 tr12 = vert[1]->co_exact - vert[2]->co_exact;
-      normal_exact = mpq3::cross(tr02, tr12);
+      normal_exact = math::cross(tr02, tr12);
     }
-    mpq_class d_exact = -mpq3::dot(normal_exact, vert[0]->co_exact);
+    mpq_class d_exact = -math::dot(normal_exact, vert[0]->co_exact);
     plane = new Plane(normal_exact, d_exact);
   }
   else {
@@ -215,14 +213,14 @@ void Face::populate_plane(bool need_exact)
       for (int i : index_range()) {
         co[i] = vert[i]->co;
       }
-      normal = double3::cross_poly(co);
+      normal = math::cross_poly(co.as_span());
     }
     else {
       double3 tr02 = vert[0]->co - vert[2]->co;
       double3 tr12 = vert[1]->co - vert[2]->co;
-      normal = double3::cross_high_precision(tr02, tr12);
+      normal = math::cross(tr02, tr12);
     }
-    double d = -double3::dot(normal, vert[0]->co);
+    double d = -math::dot(normal, vert[0]->co);
     plane = new Plane(normal, d);
   }
 }
@@ -1098,15 +1096,15 @@ static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis)
  */
 static double supremum_dot_cross(const double3 &a, const double3 &b)
 {
-  double3 abs_a = double3::abs(a);
-  double3 abs_b = double3::abs(b);
+  double3 abs_a = math::abs(a);
+  double3 abs_b = math::abs(b);
   double3 c;
   /* This is dot(cross(a, b), cross(a,b)) but using absolute values for a and b
    * and always using + when operation is + or -. */
   c[0] = abs_a[1] * abs_b[2] + abs_a[2] * abs_b[1];
   c[1] = abs_a[2] * abs_b[0] + abs_a[0] * abs_b[2];
   c[2] = abs_a[0] * abs_b[1] + abs_a[1] * abs_b[0];
-  return double3::dot(c, c);
+  return math::dot(c, c);
 }
 
 /* The index of dot when inputs are plane_coords with index 1 is much higher.
@@ -1143,11 +1141,11 @@ static int filter_plane_side(const double3 &p,
                              const double3 &abs_plane_p,
                              const double3 &abs_plane_no)
 {
-  double d = double3::dot(p - plane_p, plane_no);
+  double d = math::dot(p - plane_p, plane_no);
   if (d == 0.0) {
     return 0;
   }
-  double supremum = double3::dot(abs_p + abs_plane_p, abs_plane_no);
+  double supremum = math::dot(abs_p + abs_plane_p, abs_plane_no);
   double err_bound = supremum * index_plane_side * DBL_EPSILON;
   if (fabs(d) > err_bound) {
     return d > 0 ? 1 : -1;
@@ -1178,9 +1176,9 @@ static inline mpq3 tti_interp(
   ab -= b;
   ac = a;
   ac -= c;
-  mpq_class den = mpq3::dot_with_buffer(ab, n, dotbuf);
+  mpq_class den = math::dot_with_buffer(ab, n, dotbuf);
   BLI_assert(den != 0);
-  mpq_class alpha = mpq3::dot_with_buffer(ac, n, dotbuf) / den;
+  mpq_class alpha = math::dot_with_buffer(ac, n, dotbuf) / den;
   return a - alpha * ab;
 }
 
@@ -1209,7 +1207,7 @@ static inline int tti_above(const mpq3 &a,
   n.y = ba.z * ca.x - ba.x * ca.z;
   n.z = ba.x * ca.y - ba.y * ca.x;
 
-  return sgn(mpq3::dot_with_buffer(ad, n, dotbuf));
+  return sgn(math::dot_with_buffer(ad, n, dotbuf));
 }
 
 /**
@@ -1428,11 +1426,11 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   const double3 &d_r2 = vr2->co;
   const double3 &d_n2 = tri2.plane->norm;
 
-  const double3 &abs_d_p1 = double3::abs(d_p1);
-  const double3 &abs_d_q1 = double3::abs(d_q1);
-  const double3 &abs_d_r1 = double3::abs(d_r1);
-  const double3 &abs_d_r2 = double3::abs(d_r2);
-  const double3 &abs_d_n2 = double3::abs(d_n2);
+  const double3 &abs_d_p1 = math::abs(d_p1);
+  const double3 &abs_d_q1 = math::abs(d_q1);
+  const double3 &abs_d_r1 = math::abs(d_r1);
+  const double3 &abs_d_r2 = math::abs(d_r2);
+  const double3 &abs_d_n2 = math::abs(d_n2);
 
   int sp1 = filter_plane_side(d_p1, d_r2, d_n2, abs_d_p1, abs_d_r2, abs_d_n2);
   int sq1 = filter_plane_side(d_q1, d_r2, d_n2, abs_d_q1, abs_d_r2, abs_d_n2);
@@ -1448,9 +1446,9 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   }
 
   const double3 &d_n1 = tri1.plane->norm;
-  const double3 &abs_d_p2 = double3::abs(d_p2);
-  const double3 &abs_d_q2 = double3::abs(d_q2);
-  const double3 &abs_d_n1 = double3::abs(d_n1);
+  const double3 &abs_d_p2 = math::abs(d_p2);
+  const double3 &abs_d_q2 = math::abs(d_q2);
+  const double3 &abs_d_n1 = math::abs(d_n1);
 
   int sp2 = filter_plane_side(d_p2, d_r1, d_n1, abs_d_p2, abs_d_r1, abs_d_n1);
   int sq2 = filter_plane_side(d_q2, d_r1, d_n1, abs_d_q2, abs_d_r1, abs_d_n1);
@@ -1477,17 +1475,17 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   if (sp1 == 0) {
     buf[0] = p1;
     buf[0] -= r2;
-    sp1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
+    sp1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
   }
   if (sq1 == 0) {
     buf[0] = q1;
     buf[0] -= r2;
-    sq1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
+    sq1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
   }
   if (sr1 == 0) {
     buf[0] = r1;
     buf[0] -= r2;
-    sr1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
+    sr1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
   }
 
   if (dbg_level > 1) {
@@ -1509,17 +1507,17 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   if (sp2 == 0) {
     buf[0] = p2;
     buf[0] -= r1;
-    sp2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
+    sp2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
   }
   if (sq2 == 0) {
     buf[0] = q2;
     buf[0] -= r1;
-    sq2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
+    sq2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
   }
   if (sr2 == 0) {
     buf[0] = r2;
     buf[0] -= r1;
-    sr2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
+    sr2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
   }
 
   if (dbg_level > 1) {
@@ -1721,7 +1719,7 @@ static CDT_data prepare_cdt_input(const IMesh &tm, int t, const Vector<ITT_value
   BLI_assert(tm.face(t)->plane_populated());
   ans.t_plane = tm.face(t)->plane;
   BLI_assert(ans.t_plane->exact_populated());
-  ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact);
+  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
   prepare_need_tri(ans, tm, t);
   for (const ITT_value &itt : itts) {
     switch (itt.kind) {
@@ -1757,7 +1755,7 @@ static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm,
   BLI_assert(tm.face(t0)->plane_populated());
   ans.t_plane = tm.face(t0)->plane;
   BLI_assert(ans.t_plane->exact_populated());
-  ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact);
+  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
   for (const int t : cl) {
     prepare_need_tri(ans, tm, t);
   }
@@ -2004,9 +2002,9 @@ static bool is_quad_flip_first_third(const double3 &v1,
                                      const double3 &normal)
 {
   double3 dir_v3v1 = v3 - v1;
-  double3 tangent = double3::cross_high_precision(dir_v3v1, normal);
-  double dot = double3::dot(v1, tangent);
-  return (double3::dot(v4, tangent) >= dot) || (double3::dot(v2, tangent) <= dot);
+  double3 tangent = math::cross(dir_v3v1, normal);
+  double dot = math::dot(v1, tangent);
+  return (math::dot(v4, tangent) >= dot) || (math::dot(v2, tangent) <= dot);
 }
 
 /**
@@ -2124,7 +2122,7 @@ static Array<Face *> exact_triangulate_poly(Face *f, IMeshArena *arena)
     f->populate_plane(false);
   }
   const double3 &poly_normal = f->plane->norm;
-  int axis = double3::dominant_axis(poly_normal);
+  int axis = math::dominant_axis(poly_normal);
   /* If project down y axis as opposed to x or z, the orientation
    * of the polygon will be reversed.
    * Yet another reversal happens if the poly normal in the dominant
@@ -2203,15 +2201,15 @@ static bool face_is_degenerate(const Face *f)
   }
   double3 da = v2->co - v0->co;
   double3 db = v2->co - v1->co;
-  double3 dab = double3::cross_high_precision(da, db);
-  double dab_length_squared = dab.length_squared();
+  double3 dab = math::cross(da, db);
+  double dab_length_squared = math::length_squared(dab);
   double err_bound = supremum_dot_cross(dab, dab) * index_dot_cross * DBL_EPSILON;
   if (dab_length_squared > err_bound) {
     return false;
   }
   mpq3 a = v2->co_exact - v0->co_exact;
   mpq3 b = v2->co_exact - v1->co_exact;
-  mpq3 ab = mpq3::cross(a, b);
+  mpq3 ab = math::cross(a, b);
   if (ab.x == 0 && ab.y == 0 && ab.z == 0) {
     return true;
   }
@@ -2231,8 +2229,8 @@ static bool any_degenerate_tris_fast(const Array<Face *> triangulation)
     }
     double3 da = v2->co - v0->co;
     double3 db = v2->co - v1->co;
-    double da_length_squared = da.length_squared();
-    double db_length_squared = db.length_squared();
+    double da_length_squared = math::length_squared(da);
+    double db_length_squared = math::length_squared(db);
     if (da_length_squared == 0.0 || db_length_squared == 0.0) {
       return true;
     }
@@ -2240,8 +2238,8 @@ static bool any_degenerate_tris_fast(const Array<Face *> triangulation)
      * The triangle is almost degenerate if sin t is almost 0.
      * sin^2 t = |da x db|^2 / (|da|^2 |db|^2)
      */
-    double3 dab = double3::cross_high_precision(da, db);
-    double dab_length_squared = dab.length_squared();
+    double3 dab = math::cross(da, db);
+    double dab_length_squared = math::length_squared(dab);
     double sin_squared_t = dab_length_squared / (da_length_squared * db_length_squared);
     if (sin_squared_t < 1e-8) {
       return true;