Constrained Cubic Polynomial Fitting: A Bezier Approach?

[hotvette]

Anybody know the math/theory behind linear least squares where the curve is forced to go through the first and last data points?I'm specifically dealing with cubic polynomials.

In standard linear least squares formulation (i.e. A^TAc = A^Ty) the curve doesn't, in general, go through any of the data points. I've actually managed to come up with a formulation based on Bezier math that allows standard least squares treatment to be used, but the development is problem dependent and rather long and tedious.

Just wondering if there is an easiser way. Ideas?

 

[CRGreathouse]

A cubic polynomial is of the form y = (x - a)(x - b)(x - c)d. If you want it to go through two particular points, that constrains two of the variables -- just adjust the others as needed.

 

[hotvette]

Hmmm, not sure I follow. Picking a and b will set 2 of the zeros of the cubic, meaning it goes though the points (a,0) and (b,0), but I'm trying to force it to go through (x₀, y₀) and (x₁, y₁).

After a-l-o-t of work I was able to come up with:

y = (c(3x³-27x²+72x-48) + d(-3x³+18x²-27x+12) + x³+6x²-42x+62)/27

for (x₀, y₀) = (1,1) and (x₁, y₁) = (4,2).

Let me re-word the question. Is there a relatively straightforward way to construct an equation of a cubic polynomial that goes though 2 specific points and has 2 parameters to adjust? I need it in this form to evaluate ∂y/∂c and ∂y/∂d for least squares treatment. Unless there is another way...

 

[AlephZero]

A cubic polynomial is of the form y = (x - a)(x - b)(x - c)d. If you want it to go through two particular points, that constrains two of the variables -- just adjust the others as needed.

Yes... but if two of the constants a, b, c are complex, then it's not a very good place to start from. For example, x^3 + x gives a = 0, b = i, c = -i.

Bezier or Bernstein polymonials (apologies, can't remember which is the right name without looking it up) seem a better idea.

If your cubic goes through $$(x_0,y_0)$$ and $$(x_1,y_1)$$

$$y = a(x-x_0)^3 + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + d(x-x_1)^3$$

$$y_0$$ and $$y_1$$ fix the values of a and d. Find b and c by least squares.

 

[hotvette]

If your cubic goes through $$(x_0,y_0)$$ and $$(x_1,y_1)$$

$$y = a(x-x_0)^3 + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + d(x-x_1)^3$$

Perfect, thanks! I somehow frequently manage to come up with the most complicated way to do things...

$$y = a(x-1)^2(x-4) + b(x-1)(x-4)^2 + \frac{(x+2)}{3}$$

is much easier.
-----------------

Question - where'd you find that form of a cubic poly?

 

[AlephZero]

Question - where'd you find that form of a cubic poly?

Try Googling for the B-splines or Bernstein polynomials, (the B in B-splines probably stands for "Bernstein"). Bezier splines are very closely related. The "slope" control point at one end of a Bezier spline depends on a and b, and at the other on c and d.

I first saw this in some tutorial notes on splines, long before the Internet existed - I can't remember exactly where.

 

[nezard]

$$y = a(x-x_0)^3 + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + d(x-x_1)^3$$

$$y_0$$ and $$y_1$$ fix the values of a and d. Find b and c by least squares.

I am sitting on the same problem for several day already and follow this discussion with the great interest :)

Now, if you fix values of a and d (with previous statement), how do you apply/formulate LSS for remaining points?

Thx

 

[hotvette]

Welcome to the forums.

There may be multiple ways to approach it, but here is what I did.

1. Pick arbitrary values for b and c
2. Based on the values of a,b,c,d calculate f(x) and (f-y) for each point
3. Form J by calculating ∂y/∂b and ∂y/∂c for each point
4. Solve J^TJs = J^T(f-y) for s
5. Update values for b & c: [b' c'] = [b c] - s

It only takes 1 iteration because the problem is linear. Net is to apply non-linear LS technique. Check out the attachment in the last post of this thread:

https://www.physicsforums.com/showthread.php?t=181532
.

 

[nezard]

Thanks for your answer hotvette, but this seems to be little to advanced for me

To be honest, I also do not see how did you evolve here

$$y = a(x-1)^2(x-4) + b(x-1)(x-4)^2 + \frac{(x+2)}{3}$$

if you've had this

$$y = a(x-x_0)^3 + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + d(x-x_1)^3$$

BTW, it would be more then wonderful if you could give us parts of your example that you did on Bezier CPs.

many thanks

 

[hotvette]

Thanks for your answer hotvette, but this seems to be little to advanced for me

To be honest, I also do not see how did you evolve here...

Yeah, it's a bit convoluted. What I wanted was a form where y = y₁ at x = x₁ and y = y₀ at x = x₀. If you evaluate the expression at x = x₀ and x = x₁ you get:

$$y_1 = a(x_1-x_0)^3 \ \ and\ \ y_0 = d(x_0-x_1)^3$$

which means

$$a = \frac{y_1}{(x_1-x_0)^3} \ \ and\ \ d = \frac{y_0}{(x_0-x_1)^3}$$

Substituting those expressions:

$$y = y_1\frac{(x-x_0)^3}{(x_1-x_0)^3} + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + y_0\frac{(x-x_1)^3}{(x_0-x_1)^3}$$

I then discovered (thanks to HallsofIvy) that the cubes aren't really needed on the 1st and last terms:

$$y = y_1\frac{(x-x_0)}{(x_1-x_0)} + b(x-x_0)^2(x-x_1) + c(x-x_0)(x-x_1)^2 + y_0\frac{(x-x_1)}{(x_0-x_1)}$$

Substituting:

$$x_0 = 1\ \ \ y_0 = 1\ \ \ x_1 = 4\ \ \ y_1 = 2$$

and simplifying will yield the final form.

Re your second question, what are Bezier CPs? Cubic Polynomials? What specifically would you like to know?

 

[nezard]

I really appreciate your reply hotvette and I see your point here.

However, in the meantime I managed to discover another way to approximate “my” curve, yielding less computation overhead if compared to the least squares. To put it straight, my input is not discrete set of data (points), but the full mathematical expression y=y(x), or even parametric y=y(t) & x=x(t). The goal of an approximation here was both performance and visual-quality driven, e.g. to avoid splitting the curve into arbitrary number of very short segments interconnected with the straight line, giving overall appearance of one smooth curve, where each point is calculated by a dedicated call to above stated functions. Instead, I wanted limited number of fine cubic splines, or more precisely Bezier cubic splines, to approximate everything. Bezier simply therefore, it is an integral part of java native graphics library (there is no other curve built in java) and it comes with the neat built in antialiasing mechanism for smoothness.

The problem unexpectedly arose by defining control points for Bezier, since they are not on the curve itself. As I was trying to find solution I stumbled upon this forum and several of your posts related to this particular problem :)

BTW, the CP actually meant Control Point for a cubic Bezier :)

 

[hotvette]

OK, sounds like your task is to find Bezier control points to represent certain functions, like y = e^x or y = sin(x), as accurately as possible.

Have you figured out how to do it, or do you need help? If you need help, tell is some specific function examples you are working with.

 

[nezard]

Well, I do have a solution here, but as someone said “there is always room for improvement”

I'll describe briefly both, problem and the actual solution, and you are welcome to comment on it. Although I am not going into details, they can always be provided if needed. Also in case you would choose quite different approach, please feel invited to present it

(1) Assume function R = a + b cos(t) in polar coordinates, or in parametric form
x = a cos(t) + b cos(t)^2, y = a sin(t) + b cos(t)sin(t), which you want to approximate as close as possible with Bezier, subdivided in N sections, where N is usually in range 8-16 (to get one smooth and continuous curve).

(2) Obviously, every point on the curve is known, so start-point and end-point for Bezier is not an issue, regardless of how the overall range t=0-2Pi is subdivided into sections. The only problem here is how to find control points

(3) As the curve function is known, its derivative is known as well, which suits well the basic idea of Hermite approximation. A typical Hermite is defined with four points, having control points defined as the first derivative in its start, and in its end point, or more precisely B = {P1, P2, P1', P2'} for any arbitrary segment.

(4) Once having the Hermite approximation Points it is not a problem to transform them into Bezier of four points.

(5) The problem with Hermite approach however, is approximation accuracy, for the curve (segment) has good shape, but it still does not fit the original curve exactly, except in start and end points. Here I do the trick, by calculating the needed scale factor, required to multiply derivatives, for the curve to pass through one arbitrary point on the segment. As a result, we are getting acceptable approximation through three points per segment. There is of course still some minor remaining deviation between the original shape and the approximated one, but its the one I can well live with.

 

[hotvette]

Your method sounds perfectly valid. Whether or not there's a 'better' approach depends on a number of things:

- whether you are enforcing C2 continuity at the connection points
- whether the calcs are done on the fly and therefore need to be computationally efficient
- I'm sure there are more

Just for kicks, I used a=1 and b=3. Also, because of symmetry, I only looked at t = 0 to pi. I broke the function into 4 segments, which is plenty. I used points where the function has vertical or horzontal slope as the Bezier connect points. After a fair amount of playing around, I came up with the following control points.

segment 1:
P0 = (4, 0)
P1 = (4, 1.234536617)
P2 = (3.013679806, 2.24436971)
P3 = (1.814333489, 2.24436971)

segment 2:
P0 = (1.814333489, 2.24436971)
P1 = (0.614987172, 2.24436971)
P2 = (-0.083333333, 1.314059632)
P3 = (-0.083333333, 0.493006649)

segment 3:
P0 =(-0.083333333, 0.493006649)
P1=(-0.083333333, -0.328046334)
P2 = (0.534044209, -0.840154844)
P3 = (1.102333177, -0.840154844)

segment 4:
P0 = (1.102333177, -0.840154844)
P1 = (1.670622146, -0.840154844)
P2 = (2, -0.406642476)
P3 = (2, 0)

Visually, I cannot tell the different between the connected Bezier segments and the function. Plot it and see for yourself. I'll tell you the details of what I did tomorrow. It's late and I'm tired.
:zzz:

 

[nezard]

For illustration only, I enclosed a comparison of your approximation and the real function of Limacon of Pascal, showing a fairly good match. The approximation, in magenta, also shows all Bezier points. This (a=1 and b=3) however, is not very difficult shape. The real “fun” begins when a=4, 5, 6 or 7 having b=3. Then you might “wish” to have more then four segments in range 0-Pi.

As far as conditions are considered, yes, it should be calculated on the fly (many times in a loop per second) and therefore the performance is an issue. And for the other, I haven't tried forcing the continuity at connection points, for I consider each segment separately, but somehow I imagine having difficulties with it on the span of four segments anyway.

I am looking forward seeing the details of your method

 

[hotvette]

Since computation performance is important, I don't think there is anything I can suggest that is better than what you are currently doing. What I did was computationally intensive and doesn't appear to have resulted in anything better. I'm surprised at the accuracy you are trying to achieve if the example I gave is considered a 'fairly good match'.

Anyway, here is what I did for a=1 and b=3:

1. Split the curve for t=0 to pi into 4 segments.

2. The boundary of each segment is a point at which the slope is either horizontal or vertical.

3. Enforce C1 and C2 continuity at segment intersections. This results in 2 control points to adjust in the first segment and 1 control point to adjust in each of the remaining segments

4. Sum perpendicular distances at 10 equally distributed points in each segment

5. Using a simplex optimization program, adjust the control points in each successive segment by minimizing the qty calculated in 4.

It sounds like what you are doing is enforcing C1 (and C0 obviously) continuity at each segment boundary, which results in 2 control points to adjust for each segment. The control points are chosen to force an exact match at some point within the segment (probably around t = 0.5). I've seen articles on fitting a Bezier segment to a quarter circle that use the exact same approach. One article insists that the method is the only correct one...

One thing I'm curious about is what method you are using to judge quality of fit. I can think of a least 4 ways. In the case of the 1st segment (a=1, b=3), I came up with 4 different sets of control points depending on the quality of fit criteria:

Px0 = Px1 = 4
Px3 = 1.814333489

Py0 = 0
Py2 = Py3 = 2.24436971

Case 1:
Px2 = 3.13839852
Py1 = 1.110593575

Case 2:
Px2 = 3.00408392
Py1 = 1.244587413

Case 3:
Px2 = 3.013679808
Py1 = 1.234536616

Case 4:
Px2 = 2.931618476
Py1 = 1.319853858

I'm curious as to which you consider a better fit and why.

 

[nezard]

Well, I did not mean to sound critical with 'fairly good match', but I'll rephrase my statement gladly into: For most practical use I can imagine it is indeed an excellent match. However, only if a < b, like a=1 and b=3. Since this is not always the case I enclose a snapshot of the curve when parameters are a=5 and b=3, which illustrate rather difficult shape to be approximated in the second segment, as well as dissatisfying approximation - with my method.

To be honest, seeing how they approximate 1/4 of the circle in java brought me to the idea to try doing the same with Limacon. Unfortunately, seeing it coded in java is not enough to be able to grasp an underlaying logic which makes it work the right way – in particular if you do not have strong or fresh math background like in my case. Luckily, I came across an excellent book of prof. David Salomon from CSU, which helped me gaining essential insight into the area.

Although not aware of articles you mention, yes I use point t=0.5 to adjust control points...

It was a pleasure of 'writing' to you - thanks

 

[nezard]

OK, from whatever reason I received only one half of you post at first and missed the rest, e.g. the question.

For this purpose there is no particular 'quality assurance' methodology I can think of. Visual observation and subjective judgment is the best I can offer at this point. Even though I would probably go for 3 or 4, all four cases would do for constant a=1 and b=3, as the misalignment does not seem to be that considerable (more then 1-2 pixel on the screen length) and all curves overall look smooth. Insofar I am curious about how do you evaluate quality of fit. BTW, have a look at the curve I am getting for four segments:

Segment 1:
P0 {4.,0.}
P1 {4.,1.26725}
P2 {2.98297,2.24437}
P3 {1.81433,2.24437}
Segment 2:
P0 {1.81433,2.24437}
P1 {0.695769,2.24437}
P2 {-0.0833333,1.38385}
P3 {-0.0833333,0.493007}
Segment 3:
P0 {-0.0833333,0.493007}
P1 {-0.0833333,-0.302675}
P2 {0.504082,-0.840155}
P3 {1.10233,-0.840155}
Segment 4:
P0 {1.10233,-0.840155}
P1 {1.62068,-0.840155}
P2 {2.,-0.451426}
P3 {2.,0.}


But the problem as I see it is not to choose from these four. It is rather to ensure the method to generate control points does at least as good results for other values of a and b. In other words, that you can relay on the method always to provide you with the 'fairly good match', in all ranges of a/b. If you have an algorithm for 1 and 3, give it a try with a=5 for instance.

 

[hotvette]

OK, I think I understand what you are trying to do. Seems like quite a challenge to develop an algorithm that will be adequate for all possible combinations of a and b. I need to ponder this a bit.

Question - what technique are you using to adjust the control point values to force the curve to go though a mid point? Some sort of iterative method?

 

[nezard]

The technique to adjust CPs is not iterative. An iteration would be to expensive I expect.

As I do not know how to edit formulas here, I enclosed a brief description of the method I use in pdf file.

 

[hotvette]

Thanks for the pdf. I haven't completely deciphered it yet, but did independently come up with what I believe to be the same equations (3 points and 2 derivatives). After playing around with 5 < a < 6, b = 3 and various schemes for distributing the Bezier segments, it seems to me there are two trouble spots, both near t = pi:

1. Sometimes the final segment had very strange shape.
2. If a segment is 'nearly' linear, the method seems to break down

In the end I came up with the following algorithm based on 8 segments for t = 0 to pi that seems to work pretty well:

1. If the curve has a vertical slope in the range pi < t < 0 (i.e 0 < b/a < 2), force a segment end point at that point. Use 3 segments after this point and 5 segments before it.

2. If no vertical slope between t = 0 and pi, then evenly distribute 8 segments

3. If a segment is 'nearly' linear (e.g. 0.95 < (y₀'*x₁')/(x₀'*y₁') < 1.05) then use four point interpolation for that segment.

I don't have the graphics tools (nor patience) to look at a bunch of plots, so I came up with an analytical method to judge quality of fit. I found the max perpendicular distance between the two curves using 100 points in each segment (probably an overkill but I wanted to make sure to get the max). After trying 10,000 combinations of a and b between 0 and 10, I found the max deviation between the curves divided by the max dimension of the entire curve to be 0.00005 which is less than 1 pixel on a 20" screen with 2000 pixels. Of course, it is possible I made a mistake somewhere, but it seems to work.

Suggest you post what you consider to be the worst case values of a and b and I'll let you know the control points using my algorithm.

 

[nezard]

I am gladly accepting your invitation for parameters But, before I give you some a-s and b-s, I guess I shall explain how do I currently split the 0-Pi range and where I think the problem is hiding.

Heuristics:
1. The first segment is fixed to t=0,pi/3 and it is always approximated in excellent manner.
2. The range pi/3 < t < pi is variably approximated with n segments, where 5 < n < 10, equally subdivided over the whole range. Concerning approximation quality, it is as good as always producing a slip-up on one of the segments, for some combination of a/b and n.

This means, if i find n which satisfies one a/b, the slip-up will occur by slightly changed a/b and vice versa. One n is so far never valid for the whole range of a/b.

If I increase n to, let's say over 10, then I am getting more slip-ups, however each smaller in magnitude. This leads to impression of not so smooth surface.

As I am forced to move forward with other aspects of my development tasks I decided to revisit this problem later. However, I am still keen to try out your heuristics, if I figure out what point exactly did you mean with:
“1. If the curve has a vertical slope in the range pi < t < 0 (i.e 0 < b/a < 2), force a segment end point at that point.“
(a) is it any arbitrary point on the vertical or s-shaped slope, or
(b) is it some point of inflexion on one of the parametric functions for x, and y, or
(c) it is extreme point on one of the parametric functions?

The specific parameters b=5.8, a=3 and n=1+5 segments is for instance causing deviation from the proper curve shape. Interestingly, the deviation occurs in 4th segment, where x(t) is having inflexion, just around t=0.5 where I adjust the curve tension. So it is inflexion, either on parametric functions x(t) or on y(t), which I believe is responsible for slip-ups.

 

[hotvette]

Hmmm, I can't seem to re-produce your results for b=5.8, a=3 and n=1+5 segments. If the 1st segment is 0-pi/3 and the remainder (to pi) has 5 evenly spaced segments, then each subsequent segment is 2pi/15 long in t. This means the 4th segment has t=1.884955592 to 2.303834613, x = -0.373200267 to 0.589475637, and y = 1.148592317 to -0.65467902. So far so good? Trouble is, I don't see an inflection (where x' or y' changes sign within the segment) nor does the curve look strange. Are you sure you don't mean a=5.8, b=3? Even then, I don't see anything odd in the 4th segment. For b=5.8, a=3 and n=1+5 segments, the control points I get for the 4th segment using your method are:

Px, Py
-0.373200267, 1.148592317
-0.294311231, 0.351259409
0.088886107, -0.284957067
0.589475637, -0.65467902

Re the method I use, I find the point between 0 and pi where x' = 0 (if it exists), meaning dy/dx = infinity. In other words, where -0.5a/b = cos(t). This occurs at t=1.832390368 for b=5.8, a=3.

Still, I'd like to reproduce your result so we have a common starting point. Any ideas?

By the way, would you agree that b=5.8, a=3 gives the same result (shape of curve and same trouble spots as b=5.8/3, a=1? The only difference between the two curves is a scale factor of 3.

Also, please explain what you mean by 'slip-up'.

*************************************
Edit: OK, here are details of what I did for a=5,b=3, segments=8 for t=0,pi:

1. x'=0 at t=2.55590711
2. Five equal segments for t=0,2.55590711
3. Three equal segments for t=2.55590711,pi
4. 7th segment is nearly linear: f'(t0)/f'(t1) = 0.963364278, so use 4 point interpolation for that segment

The t-values for the segment boundaries are:

0
0.511181422
1.022362844
1.533544266
2.044725688
2.55590711
2.751135625
2.946364139
3.141592654

Worst deviation is in segment #4. Control points for this segment are:

Px, Py
0.19037845, 5.108183931
-0.711942681, 4.623221229
-1.31191621, 3.922223584
-1.657065624, 3.230657136

 

[hotvette]

This thread has come full circle. Turns out that the Bezier formulation using 3 points and dy/dx at both ends works nicely for least squares analysis because it couples the equations for x and y. In other words, y = f(x₀, x₁ , x₂ , y₀, y₁ , y₂ , dy₀/dx₀, dy₁/dx₁)

Let me know if you want to further collaborate on the function fitting. I do think I have an algorithm that works.