Fix Bezier Curve Errors with Colinear Points

CHANGELOG.md

@@ -23,6 +23,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 - Added a Euler pole velocity field (in spherical coordinates) for the continental plate, oceanic plate and mantel layer. \[Menno Fraters; 2025-08-25; [#850](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/850)\]
 
 ### Changed
+- Modify the implementation of the Bezier curve to account for the special case of co-linear points. \[Daniel Douglas; 2025-01-23; [#799](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/799)\]
 - The tian2019 composition model now returns a mass fraction instead of a mass percentage. \[Daniel Douglas; 2024-11-12; [#767](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/767)\]
 - Only link to MPI libraries if the cmake variable USE_MPI has been set. No longer automatically link to MPI if MPI is found. \[Rene Gassmoeller; 2025-01-20; [#792](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/792)\]
 - Change the Doxygen documentation design using the Doxygen Awesome theme. Also fix the main README logo so it appears in the doxygen start page. \[Rene Gassmoeller; 2025-01-21; [#807](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/807)\]

include/world_builder/objects/bezier_curve.h

@@ -66,10 +66,12 @@ namespace WorldBuilder
         ClosestPointOnCurve closest_point_on_curve_segment(const Point<2> &p, const bool verbose = false) const;
 
         /**
-         * @brief
+         * @brief Returns a point that lies on the bezier curve at some interval x between coordinate i and coordinate
+         *        i + 1. If x = 0, returns point i, if x = 1, returns point i + 1.
          *
-         * @param i
-         * @param x
+         * @param i The index of the coordinate that defines the bezier curve.
+         * @param x The value used to determine additional points that lie on the bezier curve between coordinates i
+         *          and i + 1
          * @return Point<2>
          */
         Point<2> operator()(const size_t i, const double x) const;
@@ -79,7 +81,6 @@ namespace WorldBuilder
         std::vector<std::array<Point<2>,2 > > control_points;
         std::vector<double> lengths;
         std::vector<double> angles;
-
     };
   }
 

source/world_builder/objects/bezier_curve.cc

@@ -46,6 +46,13 @@ namespace WorldBuilder
       std::vector<double> angle_constraints = angle_constraints_input;
       angle_constraints.resize(n_points,NaN::DQNAN);
 
+      // Whether points on the line are co-linear. If they are, this is a special case for the bezier curve.
+      // We determine if the points are co-linear by calculating the area of the triangle formed by the points.
+      // If the area is less than epsilon, the points are co-linear.
+      bool points_are_colinear = false;
+      const unsigned int max_arclength_discretization = 10;
+      const double epsilon = 1e-9;
+
       // if no angle is provided, compute the angle as the average angle between the previous and next point.
       // The first angle points at the second point and the last angle points at the second to last point.
       // The check points are set at a distance of 1/10th the line length from the point in the direction of the angle.
@@ -71,6 +78,13 @@ namespace WorldBuilder
               // Calculate the line between the current point and the following point
               const Point<2> P3P2 = points[p_i+1]-points[p_i];
 
+              // Check if the points are co-linear by determining the area of the triangle
+              // formed by the 3 points. This is a special case of the bezier curve.
+              if ( std::abs(points[p_i-1][0] * (points[p_i][1] - points[p_i+1][1]) +
+                            points[p_i][0] * (points[p_i+1][1] - points[p_i-1][1]) +
+                            points[p_i+1][0] * (points[p_i-1][1] - points[p_i][1])) < epsilon)
+                points_are_colinear = true;
+
               // Calculate the angles of the two lines determined above
               const double angle_p1p2 = atan2(P1P2[1],P1P2[0]);
               const double angle_p3p1 = atan2(P3P2[1],P3P2[0]);
@@ -115,6 +129,7 @@ namespace WorldBuilder
             control_points[0][0][1] = sin(angles[0])*length*fraction_of_length+p1[1];
             control_points[0][1][0] = cos(angles[1])*length*fraction_of_length+p2[0];
             control_points[0][1][1] = sin(angles[1])*length*fraction_of_length+p2[1];
+
             {
               // Determine which side of the line the control points lie on
               const bool side_of_line_1 =  (p1[0] - p2[0]) * (control_points[0][1][1] - p1[1])
@@ -123,12 +138,42 @@ namespace WorldBuilder
               const bool side_of_line_2 =  (p1[0] - p2[0]) * (p3[1] - p1[1])
                                            - (p1[1] - p2[1]) * (p3[0] - p1[0])
                                            < 0;
-              if (side_of_line_1 == side_of_line_2)
+
+              // The points are co-linear, so we need to check if the control points are within the line p1p2
+              if (points_are_colinear)
+                {
+                  const bool cp_1_within_p1p2 = (std::min(p1[0], p2[0]) - epsilon <= control_points[0][0][0] && control_points[0][0][0] <= std::max(p1[0], p2[0]) + epsilon) &&
+                                                (std::min(p1[1], p2[1]) - epsilon <= control_points[0][0][1] && control_points[0][0][1] <= std::max(p1[1], p2[1]) + epsilon);
+                  const bool cp_2_within_p1p2 = (std::min(p1[0], p2[0]) - epsilon <= control_points[0][1][0] && control_points[0][1][0] <= std::max(p1[0], p2[0]) + epsilon) &&
+                                                (std::min(p1[1], p2[1]) - epsilon <= control_points[0][1][1] && control_points[0][1][1] <= std::max(p1[1], p2[1]) + epsilon);
+                  // If the control points are not within the line p1p2, we need to move them within the line p1p2
+                  if (!cp_1_within_p1p2)
+                    {
+                      control_points[0][0][0] = cos(angles[0]+Consts::PI)*length*fraction_of_length+p1[0];
+                      control_points[0][0][1] = sin(angles[0]+Consts::PI)*length*fraction_of_length+p1[1];
+                    }
+                  if (!cp_2_within_p1p2)
+                    {
+                      control_points[0][1][0] = cos(angles[0]+Consts::PI)*length*fraction_of_length+p1[0];
+                      control_points[0][1][1] = sin(angles[0]+Consts::PI)*length*fraction_of_length+p1[1];
+                    }
+                }
+
+              else if (side_of_line_1 == side_of_line_2)
                 {
                   // use a 180 degree rotated angle to create this control_point
                   control_points[0][1][0] = cos(angles[1]+Consts::PI)*length*fraction_of_length+p2[0];
                   control_points[0][1][1] = sin(angles[1]+Consts::PI)*length*fraction_of_length+p2[1];
                 }
+
+              // There is no closed-form analytic way to express the arc-length of a cubic bezier curve. We approximate
+              // the arc-length by dividing the curve into 10 points and piecewise linearly connect them. We also store the
+              // length of the bezier curve within each of these intervals. We calculate the points that lie on the bezier
+              // curve using the operator function below.
+              for (unsigned int t_value = 1; t_value <= max_arclength_discretization; ++t_value)
+                {
+                  lengths[0] = (operator()(0, static_cast <double> (t_value)/max_arclength_discretization) - operator()(0, static_cast <double> (t_value - 1)/max_arclength_discretization)).norm();
+                }
             }
           }
 
@@ -139,7 +184,6 @@ namespace WorldBuilder
               const double length = (points[p_i]-points[p_i+1]).norm(); // can be squared
               control_points[p_i][0][0] = cos(angles[p_i])*length*fraction_of_length+p1[0];
               control_points[p_i][0][1] = sin(angles[p_i])*length*fraction_of_length+p1[1];
-
               {
                 // Determine which side of the line the control points lie on
                 const bool side_of_line_1 =  (p1[0] - p2[0]) * (control_points[p_i-1][1][1] - p1[1])
@@ -148,16 +192,43 @@ namespace WorldBuilder
                 const bool side_of_line_2 =  (p1[0] - p2[0]) * (control_points[p_i][0][1] - p1[1])
                                              - (p1[1] - p2[1]) * (control_points[p_i][0][0] - p1[0])
                                              < 0;
-                if (side_of_line_1 == side_of_line_2)
+
+                // The points are co-linear, so we need to check if the control points are within the line p1p2
+                if (points_are_colinear)
+                  {
+                    const bool cp_1_within_p1p2 = (std::min(p1[0], p2[0]) <= control_points[p_i][0][0] && control_points[p_i][0][0] <= std::max(p1[0], p2[0])) &&
+                                                  (std::min(p1[1], p2[1]) <= control_points[p_i][0][1] && control_points[p_i][0][1] <= std::max(p1[1], p2[1]));
+                    const bool cp_2_within_p1p2 = (std::min(p1[0], p2[0]) <= control_points[p_i][1][0] && control_points[p_i][1][0] <= std::max(p1[0], p2[0])) &&
+                                                  (std::min(p1[1], p2[1]) <= control_points[p_i][1][1] && control_points[p_i][1][1] <= std::max(p1[1], p2[1]));
+                    // If the control points are not within the line p1p2, we need to move them within the line p1p2
+                    if (!cp_1_within_p1p2)
+                      {
+                        control_points[p_i][0][0] = cos(angles[p_i]+Consts::PI)*length*fraction_of_length+p1[0];
+                        control_points[p_i][0][1] = sin(angles[p_i]+Consts::PI)*length*fraction_of_length+p1[1];
+                      }
+                    if (!cp_2_within_p1p2)
+                      {
+                        control_points[p_i][1][0] = cos(angles[p_i]+Consts::PI)*length*fraction_of_length+p1[0];
+                        control_points[p_i][1][1] = sin(angles[p_i]+Consts::PI)*length*fraction_of_length+p1[1];
+                      }
+                  }
+
+                // Check to see if the angles are different. If the angles are the same, points p1, p2, and p3
+                // are co-linear, and therefore the control points will also be co-linear with p1, p2 and p3. This
+                // makes determining which 'side' the control points lie meaningless.
+                else if (side_of_line_1 == side_of_line_2)
                   {
                     // use a 180 degree rotated angle to create this control_point
                     control_points[p_i][0][0] = cos(angles[p_i]+Consts::PI)*length*fraction_of_length+p1[0];
                     control_points[p_i][0][1] = sin(angles[p_i]+Consts::PI)*length*fraction_of_length+p1[1];
                   }
               }
 
-              control_points[p_i][1][0] = cos(angles[p_i+1])*length*fraction_of_length+points[p_i+1][0];
-              control_points[p_i][1][1] = sin(angles[p_i+1])*length*fraction_of_length+points[p_i+1][1];
+              if (!points_are_colinear)
+                {
+                  control_points[p_i][1][0] = cos(angles[p_i+1])*length*fraction_of_length+p2[0];
+                  control_points[p_i][1][1] = sin(angles[p_i+1])*length*fraction_of_length+p2[1];
+                }
 
               if (p_i+1 < n_points-1)
                 {
@@ -168,15 +239,53 @@ namespace WorldBuilder
                   const bool side_of_line_2 =  (p1[0] - p2[0]) * (p3[1] - p1[1])
                                                - (p1[1] - p2[1]) * (p3[0] - p1[0])
                                                < 0;
-                  if (side_of_line_1 == side_of_line_2)
+                  // The points are co-linear, so we need to check if the control points are within the line p1p2
+                  if (points_are_colinear)
+                    {
+                      const bool cp_1_within_p1p2 = (std::min(p1[0], p2[0]) <= control_points[p_i][0][0] && control_points[p_i][0][0] <= std::max(p1[0], p2[0])) &&
+                                                    (std::min(p1[1], p2[1]) <= control_points[p_i][0][1] && control_points[p_i][0][1] <= std::max(p1[1], p2[1]));
+                      const bool cp_2_within_p1p2 = (std::min(p1[0], p2[0]) <= control_points[p_i][1][0] && control_points[p_i][1][0] <= std::max(p1[0], p2[0])) &&
+                                                    (std::min(p1[1], p2[1]) <= control_points[p_i][1][1] && control_points[p_i][1][1] <= std::max(p1[1], p2[1]));
+                      // If the control points are not within the line p1p2, we need to move them within the line p1p2
+                      if (!cp_1_within_p1p2)
+                        {
+                          control_points[p_i][0][0] = cos(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p1[0];
+                          control_points[p_i][0][1] = sin(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p1[1];
+                        }
+                      if (!cp_2_within_p1p2)
+                        {
+                          control_points[p_i][1][0] = cos(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p1[0];
+                          control_points[p_i][1][1] = sin(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p1[1];
+                        }
+                    }
+
+                  // Check to see if the angles are different. If the angles are the same, points p1, p2, and p3
+                  // are co-linear, and therefore the control points will also be co-linear with p1, p2 and p3. This
+                  // makes determining which 'side' the control points lie meaningless.
+                  else if (side_of_line_1 == side_of_line_2)
                     {
                       // use a 180 degree rotated angle to create this control_point
                       control_points[p_i][1][0] = cos(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p2[0];
                       control_points[p_i][1][1] = sin(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p2[1];
                     }
                 }
+
+              // There is no closed-form analytic way to express the arc-length of a cubic bezier curve. We approximate
+              // the arc-length by dividing the curve into 10 points and piecewise linearly connect them. We also store the
+              // length of the bezier curve within each of these intervals. We calculate the points that lie on the bezier
+              // curve using the operator function below.
+              for (unsigned int t_value = 1; t_value <= max_arclength_discretization; ++t_value)
+                {
+                  lengths[p_i] = (operator()(p_i, static_cast <double> (t_value)/max_arclength_discretization) -
+                                  operator()(p_i, static_cast <double> ((t_value - 1))/max_arclength_discretization)).norm();
+                }
             }
         }
+
+      else
+        {
+          lengths[0] = (points[0]-points[1]).norm();
+        }
     }
 
 

tests/gwb-dat/cartesian_fault_x_and_y_direction/screen-output.log

@@ -1,5 +1,5 @@
 # x y z d g T vx vy vz c0 c1 tag 
--20000 90000 990000 10000 20 6.44612 7.44612 8.44612 0 1 0
+-20000 90000 990000 10000 20 6.34463 7.34463 8.34463 0 1 0
 -20000 580000 990000 10000 20 4 5 6 0 1 0
 -20000 910000 990000 10000 20 4 5 6 0 1 0
 70000 -30000 990000 10000 10 1 2 3 1 0 0

tests/gwb-dat/cartesian_slab_x_and_y_direction/screen-output.log

@@ -1,5 +1,5 @@
 # x y z d g T vx vy vz c0 c1 tag 
--20000 90000 990000 10000 20 6.44612 7.44612 8.44612 0 1 0
+-20000 90000 990000 10000 20 6.34463 7.34463 8.34463 0 1 0
 -20000 580000 990000 10000 20 4 5 6 0 1 0
 -20000 910000 990000 10000 20 4 5 6 0 1 0
 70000 -30000 990000 10000 10 1 2 3 1 0 0

tests/gwb-dat/slab_interpolation_simple_CMS/screen-output.log

@@ -56,6 +56,6 @@
 -104000 120000 150000 50e3 1622.56 0 0 0 -1
 -14000 70000 150000 50e3 1622.56 0 0 0 -1
 -163600 124000 -10 200010 1692.16 0 0 0 -1
-90000 200000 0 200000 1692.16 0 0 0 -1
+90000 200000 0 200000 600 0 0 0 0
 90000 180000 0 200000 600 0 0 0 0
 86000 0 0 200000 1692.16 0 0 0 -1

tests/gwb-dat/smooth_composition_fault/screen-output.log

@@ -1,9 +1,9 @@
 # x y z d g T vx vy vz c0 c1 tag 
-6371000 -24.5 -1. 0 1600 0 0 0 0.145198 0 0
+6371000 -24.5 -1. 0 1600 0 0 0 0.999675 0 0
 6365000 -23. -1.3 6e3 1602.64 0 0 0 0.728204 0 0
 6361000 -23.3 -1.3 10e3 1604.4 0 0 0 0.879318 0 0
 6371000 -23. -1. 0 1600 0 0 0 0.999489 0 0
-6370000 -23.4 -1.66 1e3 1600.44 0 0 0 0.000138979 0 0
+6370000 -23.4 -1.66 1e3 1600.44 0 0 0 0 0 -1
 6370000 -23.4 -1.45 1e3 1600.44 0 0 0 0.271888 0 0
 6370000 -23.4 -1.25 1e3 1600.44 0 0 0 0.968088 0 0
 6370000 -23.4 -1.0 1e3 1600.44 0 0 0 0.999865 0 0

tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Features__Distance_to_Feature_Plane.txt

@@ -3,10 +3,10 @@ Test
 
 [0] = 7141.78
 [1] = 70.7107
-[2] = 482.865
-[3] = 153.282
-[4] = 468.712
-[5] = 139.151
+[2] = 481.055
+[3] = 155.091
+[4] = 466.902
+[5] = 140.96
 [6] = 1070.96
 [7] = 6141.53
 [8] = inf