Improve math for Bezier approximation

[johnzeringue]

@noahmclean, you had given me this to use to approximate parametric curves with Beziers, and Dr. Bowring and I had quite the time with it today.

We were working under the assumption that t₀ and t₁ were t-values in our application (more specifically the min and max t-values for each Bezier). So the kinda stuff we did for the first pass looked something like this:

Let f(t) = (x(t), y(t)).

x_P1 = x(t₀) + ((t₁ - t₀) / 3) * f'(t₀) y_P1 = y(t₀) + ((t₁ - t₀) / 3) * f'(t₀)

x_P2 = x(t₁) - ((t₁ - t₀) / 3) * f'(t₁) y_P2 = y(t₁) - ((t₁ - t₀) / 3) * f'(t₁)

After closer inspection, this is clearly nonsense, because in some cases (particularly ours) x, y, and t will be of drastically different magnitudes. So, after quite a bit of frustration, I just went ahead and wrote something to this effect myself:

Let f(t) = (x(t), y(t)).

x_P1 = x(t₀) + (|f(t₁) - f(t₀)| / 3) * cos(tan^-1(f'(t₀))) y_P1 = y(t₀) + (|f(t₁) - f(t₀)| / 3) * sin(tan^-1(f'(t₀)))

x_P2 = x(t₁) - (|f(t₁) - f(t₀)| / 3) * cos(tan^-1(f'(t₁))) y_P2 = y(t₁) - (|f(t₁) - f(t₀)| / 3) * sin(tan^-1(f'(t₁)))

Really the only issue with it is that the degree of concavity is constant, but I don't know any better way to do it. So I just wanted to pass this along to you to give you an opportunity to either point out what I missed in the text that you provided or improve upon what I came up with.

No one will ever be able to tell the difference for our purposes right now. I just wonder if we'll ever need to improve this to handle curvier curves in the future.

[johnzeringue]

Actually, now that I think about this, it's probably pretty straightforward to solve for better values using a point or two on the curve. I'm not going to bother with it myself for a while, though. These sorts of things'll drive you crazy if you let them.

[noahmclean]

@johnzeringue, it looks like you were on the right track with the first set of equations above, then got a bit off track with the trigonometry. That first set of equations looks like the ones I pointed out here: http://www.math.ubc.ca/~cass/graphics/manual/pdf/ch6.pdf . In that article, the variables represented by B, B', and P are vectors, with an x- and a y-component, and t is a scalar. You can do the math by writing out separate equations for x and y, as above, or by leaving them together and adding vectors and multiplying them by scalars--your choice. I'll divide them up, as you have done.

To make the math a bit more clear, let's let x = f(t) and y = g(t). For the Wetherill concordia, that means that x = r207_235r = exp(lambda235_t) - 1 = f(t) and y = r206_238r = exp(lambda238_t) - 1 = g(t)

The next thing we need to calculate is the derivatives of each function with respect to t: f'(t) and g'(t). This is pretty easy for the exponential functions involved, and will in general be possible even for more complicated cases like the r207_206r axis you'll meet in the Tera-Wasserburg concordia plot. For Wetherill, f'(t) = lambda235_exp(lambda235_t) and g'(t) = lambda238_exp(lambda238_t)

Armed with this information, we can calculate the four control points needed for a cubic Bezier curve segment. You've got the correct equations above, but let me substitute g'(t) into the y-coordinate equations:

xP0 = f(t0) yP0 = g(t0)

xP1 = f(t0) + ((t1 - t0) / 3) * f'(t0) yP1 = g(t0) + ((t1 - t0) / 3) * g'(t0)

xP2 = f(t1) - ((t1 - t0) / 3) * f'(t1) yP2 = g(t1) - ((t1 - t0) / 3) * g'(t1)

xP3 = f(t1) yP3 = g(t1)

The t0 and t1 in the equations above are indeed the same 't' that you use to label the concordia plots--in units of years. Thus, drawing a segment from zero to ten million years means that t0 = 0, t1 = 10e6, and the derivatives f'(t) and g'(t) above are evaluated with those values of t.

It appears that part of the confusion stems from your conclusion that "After closer inspection, this is clearly nonsense, because in some cases (particularly ours) x, y, and t will be of drastically different magnitudes." It helps if, for the purposes of factor analysis (getting the units right), you can think of the derivatives as fractions. In the right hand side of the equation for xP1, there are two terms: f(t0) is in the same units as the x-axis, then the second is a term in units of t multiplied by dx/dt, If you think of dt as 'a small change in t' and (t1-t0)/3 as 'a small distance in t', then the two units cancel out, leaving you with the units of 'small change in x'. Combined with the 'small change in y' from the second term in yP0, the first control point P1 is just P0 shifted over and up in the direction that the curve was travelling at P0, which is quantified by its derivative ( f'(t0), g'(t0) ). Hopefully a little experimentation on your part should demonstrate that the x- and y-variables having very different magnitudes doesn't matter--they don't appear in the same equations above.

The Java curveTo and the svg Path commands do a nice job of eliminating redundant information when making the piecewise continuous curves. Since the P3 of one segment will be the P0 of the next, the P0 can get left out of subsequent curveTo statements or 'C' (cubic Bezier) path segments.

In general, Topsoil will be making lots of these paths for different plots--I see Tera-Wasserburg is up next. For the general case of curves in plots with two measured variables (like concordia curves, regression lines, isochrons or other iso-curves, etc, and their uncertainty envelopes), we will be able to analytically calculate the derivatives as above, though it might take a bit more math. For one-variable plots (e.g. KDE plots, weighted mean plots, MSWD distributions), the algorithm outlined above can be generalized so that x = f(t) = t, and f'(t) = 1.

Once this works, we can later add in the ability to test for multiple intersections with the plot box, and to choose appropriate endpoints for the curve segments so that they all have roughly the same length.

[johnzeringue]

Thanks, this was a silly oversight on our end. I'll reimplement it sometime soon.

I was confused over the units, because I was thinking that B' was dy/dx instead of (dx/dt, dy/dt).

Would you mind helping me with the derivatives for the Tera-Wasserburg, then?

[noahmclean]

For the Tera-Wasserburg derivatives: x = f(t) = r238_206r = 1/( exp(lambda238_t) - 1) y = g(t) = r207_206r = ( exp(lambda235_t) - 1 )/( r238_235s * (exp(lambda238*t) - 1) )

f'(t) = - lambda238 * exp(lambda238_t) / ( exp(lambda238_t) - 1 )^2 g'(t) = - (lambda235 * exp(lambda235_t) - exp(lambda238t)(lambda238 + lambda235 * exp(lambda235_t) - lambda238 * exp(lambda235_t)))/(r238_235s * (exp(lambda238_t) - 1)^2)

Don't forget the order of operations in the expression for f'(t) -- make sure to square the denominator rather than the entire expression.

[johnzeringue]

Does the math on page eleven still apply to 3D points?

[noahmclean]

Yup, everything still works. In more than two dimensions, it helps to think of (and calculate using) the points as vectors and the derivatives as gradients. Does javafx support 3d 'space curves'?

[johnzeringue]

JavaFX supports triangle meshes, so yes and no. They won't be true curves, but approximations of the intended surfaces using triangles.

[noahmclean]

Using the Bezier curve equations above will give you a one-dimensional object (a curve or line) embedded in a 2D or 3D space, not a surface (two-dimensional object).

[johnzeringue]

Is that the fancy way to say that we'l have a one-dimensional codomain from either a one- or two-dimensional domain?

[noahmclean]

Not quite, related but different. The equations discussed in section 6.5 of http://www.math.ubc.ca/~cass/graphics/manual/pdf/ch6.pdf generate what this Wikipedia article calls curves: http://en.wikipedia.org/wiki/Curve. Can we use the words 'curve' and 'surface' and not get confused?

[johnzeringue]

Ah, I see. You just meant that a curve only has one-dimension of "size", because it's infinitely thin.

I'm on track now. The term "space curve" confused me earlier.

As far as JavaFX support for 3D curves, it's not especially great out of the box. Options are (1) create cylindrical-ish meshes to create curves with visible thickness or (2) project the curves ourselves into 2D Beziers (which is probably the better solution). Note that anything that lies in some 2D plane is a special case that's easier to handle, since we should be able to recreate it using simple transformations.

[johnzeringue]

This was implemented in the latest release.