Fix Bezier Curve Errors with Colinear Points

[danieldouglas92]

Add the functionality within bezier_curve.cc to determine the total arc-length of the bezier curve across all user-specified coordinates, and the length of the bezier curve between each incremental pair of user-specified coordinates.

[github-actions]

Benchmark	Main	Feature	Difference (99.9% CI)
Slab interpolation simple none	1.114 ± 0.009 (s=386)	1.117 ± 0.009 (s=424)	+0.1% .. +0.4%
Slab interpolation curved simple none	1.116 ± 0.011 (s=393)	1.116 ± 0.008 (s=416)	-0.1% .. +0.3%
Spherical slab interpolation simple none	1.094 ± 0.008 (s=407)	1.095 ± 0.010 (s=418)	-0.1% .. +0.3%
Slab interpolation simple curved CMS	1.165 ± 0.013 (s=404)	1.168 ± 0.015 (s=370)	-0.0% .. +0.5%
Spherical slab interpolation simple CMS	1.442 ± 0.012 (s=331)	1.421 ± 0.013 (s=300)	-1.6% .. -1.2%
Spherical fault interpolation simple none	1.100 ± 0.009 (s=385)	1.103 ± 0.009 (s=435)	+0.1% .. +0.4%
Cartesian min max surface	2.535 ± 0.017 (s=181)	2.545 ± 0.020 (s=176)	+0.1% .. +0.6%
Spherical min max surface	7.155 ± 0.066 (s=69)	7.179 ± 0.080 (s=59)	-0.3% .. +1.0%

[danieldouglas92]

At this point I think the code that I have is correct, there was an issue with the existing implementation of the Bezier curve when dealing with points that were co-linear.

If you have points P1 and P2, to compute the Bezier curve two additional control points CP1 and CP2 must be determined which control the shape of the curve. Sometimes, the control points would be outside of the bounds of Points P1 and P2, as illustrated below:

x----o--------x----o
P1 CP1 P2 CP2

However, the control points should lie within on the vector P1P2, not beyond it, like this:
x----o--------o----x
P1 CP1 CP2 P2

This is important because when the control points lie beyond the vector P1P2, the arclength of the Bezier curve is artificially made longer.

Currently, some of the tests are still failing but they are mostly related to tests which utilize coordinates that are colinear (with the exception of the test slab_interpolation_simple_curved_CMS.dat
@MFraters this is ready for a look when you have the time!

[MFraters]

This looks good to me. thanks for finding and fixing the bug! I have a few small comments, but otherwise I think it is good to go.

[danieldouglas92]

Thanks for the review! I've addressed your comments, it turns out it wasn't slab_interpolation_simple_curved_CMS that was failing, but actually slab_interpolation_simple_CMS that was failing. This means that all the tests that are actually failing currently involve coordinates that are colinear. If you agree that it is ok to update the test results then I will do that in a follow up commit and hopefully that will wrap this up.

[danieldouglas92]

After talking to @MFraters I've tentatively updated the test results. The test results should be more correct now, because the tests that were failing involved colinear points. However, we want to do a little more checking to make sure that the new code still produces correct results, so do not merge this quite yet.

[danieldouglas92]

Finally returning to this to check the output from the tests that were failing to confirm that things look good with the changes in this branch. Good news is that the changes in this branch actually improve the output! Below are screenshots comparing the output from 4 tests that were failing on the main branch, and on the changes in this branch.

Here is the test dat_cartesian_slab_x_and_y_direction

Here is the test dat_cartesian_fault_x_and_y_direction

Here is the test dat_slab_interpolation_simple_CMS. This one is trickier to determine if the updates in this branch are "correct", the first pair of screenshots shows the top view where it seems like the main branch performs better.

However, from a side view this is how the main branch compares with the changes in this branch.

Lastly, here is the comparison for dat_smooth_composition_fault.

For 3 of the 4 tests, the tests prior to this change produce features that are disjointed/discontinuous, and the changes within this branch certainly fix that. For the dat_slab_interpolation_simple_CMS test, it's less clear, though I suspect this might just be due to the trench coordinates used in the test. The coordinates jump from values of single meters to values of 1000's of meters:

[0,0],[1,2],[2,3],[10,10],[11,11],[50e3,50e3],[125e3,75e3]

Thoughts on this and these comparisons @MFraters?

[danieldouglas92]

@MFraters Given the results of Issue #861, and also the way that the tests look with/without this feature, I think it's looking like the changes added here should be implemented soon. Thoughts?

[danieldouglas92]

I've fixed the merge conflicts and I'll fix the broken tests sometime tomorrow.

[MFraters]

Thanks! It does seem to fix #861, but I need to look at bit more at it to remember what this pull request is doing. I will have another look at it this Friday.

[MFraters]

For the dat_slab_interpolation_simple_CMS test, it's less clear, though I suspect this might just be due to the trench coordinates used in the test. The coordinates jump from values of single meters to values of 1000's of meters:

So, does the issue disappear when you use values on the same order of magnitude and/or change it to:

[0,0],[50e3,50e3],[125e3,75e3]

[MFraters]

Generally looks good, although these are kind of two different pull requests in one. I am fine with that, but please split it up in the changelog.

Can you remind me why you need the arclength?

Can you maybe add a test of a bezier curve of which you know the arc length, to make sure the function actually computes a reasonable length?

[danieldouglas92]

@MFraters Thanks for taking a look at this again. The bezier curve length was something that was needed for the contours implementation that we were working on at the Hackathon (PR #829), and since (as you said) it actually doesn't do anything in this PR I think I'll remove the the total arclength from this PR and save that for a follow up PR.

[danieldouglas92]

@MFraters The slab gap introduced in the slab_interpolation_simple_CMS seems to only happen when there are trench points within a certain configuration. I'll try to layout what I found clearly below:

For the dat_slab_interpolation_simple_CMS test, it's less clear, though I suspect this might just be due to the trench coordinates used in the test. The coordinates jump from values of single meters to values of 1000's of meters:

So, does the issue disappear when you use values on the same order of magnitude and/or change it to:

[0,0],[50e3,50e3],[125e3,75e3]

this set of coordinates works perfectly fine for both the current main and this branch, since these points are not colinear. Just as a reference, this is what the slab looks like when the trench is [0,0],[50e3,50e3],[125e3,75e3]:

If I make the multiply some of the coordinates by 1000 so that they are all on the order of km's instead of m's:

[0,0],[1e3,2e3],[2e3,3e3],[10e3,10e3],[11e3,11e3],[50e3,50e3],[125e3,75e3]

then the slab gap is still present:

The slab gap only seems to be removed depending on the placement of trench coordinates around the trench-kink point at [50e3,50e3]. For x<50e3, the trench is (mostly) a straight line with a slope of 1, but for x>50e3 the trench is a straight line with a slope of 3. If there are several colinear coordinates between [0,0] and [50e3,50e3], then there is a slab gap. However, there can be several colinear points between [50e3,50e3] and [125e3,75e3] and the slab does not have slab gap.

I need to double check the code, but I suspect that this is due to the changes here overconstraining the bezier curve on either side of this kink point.