Bezier lengths w/o lots of resources

What chip?

This is an article about fast computation of curve sections using fixed point. The discussion is quite interesting.

George

Yes. Perhaps there is a better word than naive. It's poor on rounding error control and converges slowly. As a result there are a lot of calculations that are thown away.

OK, I've had chance to try this out in the case of a cubic. Start with just the end points then subdivide in smaller steps of the control parameter. For simplicity one line segment then two, four etc. Then printed out the summed length and the romberg estimate. In the example below the path happens to be exactly 2.0, making it easy to see the errors.

Control points (0.000000,0.000000) (0.000000,-1.000000) (1.000000,1.000000) (1.000000,0.000000)

Control points (0.000000,0.000000) (0.000000,1.000000) (1.000000,1.000000) (1.000000,0.000000)

sqrts length Romberg 1 1.000000 2 1.802776 2.070368 4 1.950719 2.000034 8 1.987716 2.000048 16 1.996931 2.000003 32 1.999233 2.000000 64 1.999808 2.000000 128 1.999952 2.000000 256 1.999988 2.000000 512 1.999997 2.000000

1024 1.999999 2.000000

so to get to 1 part per million the straightforward ('naive') method uses

1024 square roots for the result, and (assume no better way to start than from a single line segment) 1023 sqrt used getting there. The romberg method gets to a slightly better result using 48 square roots (32+16) and wastes

1. That seems like a worthwhile acceleration to me. The convergence power of h is 2 for the naive method and 4 for the romberg.

The control points are hardly much of an upper bound in the case of a small number of control points. In my example, above, the upper board is 3.

It might help to know the context. What is the purpose. What degree of polynomials are you talking about. I think you said 2-D in another post but... What hardware have you in mind?

peter

I'll assume the above is NOT what you used but, rather, the one which follows (?)

I'm not sure we're on the same page, here. :-( (see below) I think you are playing solely with "t" (or whatever you are using as the control variable) and subdividing *that* (?)

[did you run this in Maple/Mathematica/etc.? If so, can you post your code?]

No, you update the upper -- and lower! -- bounds as you subdivide the curve.

E.g., the first iteration of "my" algorithm has the length of the segments in the control polygon as 3 (= 1+1+1) and the length of the *chord* as 1 (obvious). So, my initial estimate for length is 2 ((3+1)/2).

Next, I subdivide the curve creating two curves:

(0,0) (0,0.5) (0.25,0.75) (0.5,0.75)

and

(0.5,0.75) (0.75,0.75) (1,0.5) (1,0)

(in this example, they are mirror images of each other)

[apologies if I screw this up... I'm doing it in my head as I type :-/ ]

The sum of the control polygon segments is ~1.103 (1/2 + sqrt(2)/4 + 1/4) and the length of the chord is ~0.901 (sqrt(0.5^2 + 0.75^2)). The estimate for the length of *this* curve (which is just the "left half" of the original curve) would the be ~1.002 ((1.103 + 0.901)/2). Since the left and right halves of the original curve are mirror images, the length of the *right* curve will be the same leaving the estimate for the total curve to be ~2.004

[sorry, going any further than that without a calculator and pen/paper is too hard for my little brain! :> ]

Note that many of the calculations performed in the first curve's assessment can be salvaged in the assessment of the second curve. (I haven't sat down and counted how many sqrts are *required* at each subdivision of the curve and how many can be "inferred").

The problem with all of these is coming up with criteria that you can automagically use to stop the recursion. E.g., I would simply use something like |estimate[i+1]| < |estimate[i]| - epsilon But, this assumes the process converges monotonically...

[attributions elided]

Agreed. Hence my suggestion to leave results "demoralized" ^^^ ;-) at times.

Of course, "integer" is really "fixed point" as X(i) and Y(i) are real numbers.

I don't think a naive "lets use this operator on all values" approach will work effectively. E.g., in regions of low curvature you can get away with bigger segment lengths (eats up a bigger portion of the workload quicker!). And, in high curvature, you quickly fall into the mud trying to get nice tight approximations to sharp curves.

So, maybe use floats for the big numbers and a fixed point scheme for the small increments. (or, two different fixed point representations conveniently scaled :> )

Or, do it all fixed point and just waste bits. E.g., using serial adders and multipliers means you just need to spend money on flip-flops.

Correct. The same is true when doing things like computing the length of "nearly vertical/horizontal" segments (the shorter axis gets lost when squared and added to the squared *longer* axis)

You can't even (naively) do that! I.e., you have to look at the magnitude of the segment length in order to decide where it will least disturb your error (i.e., if one of those "samples" ends up being large, then accumulating it with the 99 other "small ones" can dwarf their contributions)

It depends. If you can get rid of the sqrt, then it might be worth it to for a grid approach (pretending to be drawing the curve), and then add up all the single-pixel straight and diagonal (times 0.5*sqrt(2)) lengths. But, if you're not careful, you might end up zigzagging, adding significant amounts of error.

--
Gemaakt met Opera's revolutionaire e-mailprogramma:  
http://www.opera.com/mail/
(remove the obvious prefix to reply by mail)

I was thinking that for very small (dx,dy) you could probably do a table lookup or other approximation. But, even that suffers from not being able to put a lid on the error E.g., what if I have my control/end -points: (0.0001,0.0001) (0.1025,0.3077) (0.7442,0.0330) (0.0040,0.9908)

The "pixel" analogy doesn't work for me. There's nothing magical about a displacement of "1" (e.g., pixel). It could just as easily be "10" or "0.1" (don't think in terms of displays)

Addition in floating point representation must be done with equal exponents, so the smaller value needs to be denormalized to have the same exponent as the larger.

In decimal notation 9.97E03+5.0E01 = 9.97E03+0.05E03 = 10.02E3 =

1.002E04. Of course, you could postpone the overflow normalization, if other additions/subtractions will immediately follow and if the arithmetic unit has sufficient headroom bits.

Fixed point add/sub is the same as integer, in multiplication, you just have to shift right the result by the number of fractional bits etc.

As long as you avoid too frequent float/fixed or fixed/float conversions in a single expression since (ne)normalization is expensive using traditional instruction sets.

Addition, subtraction, multiplication and denormalization can be done with right shifting shift registers, however, normalization (and division) would also require left shifting shift registers.

_Bit_ serial computing was used a lot in the tube era (and in the HP-35 pocket scientific calculator), but I haven't seen it used more recently.

These days, a huge number of bit serial processors might make sense if each serial processor is controlling a single pixel or block of pixels and some of these processors could be completely shut down to reduce power consumption when no update is required. But otherwise, for equal total throughput, both serial and parallel systems would require about the same amount of arithmetic elements, but the more complex control logic on serial processors would have to be duplicated on more processors.

Sorry, yes. I pasted in one example and then decided that it would be clearer if I used an example with simple result.

Probably not. Yes I'm using the control parameter, which is commonly called t.

It's etc. I though the premise was a resource challenged system. So I wrote it in C - well its actually C++ but in C style. Pasted at the end, but only one line for the Romberg bit.

This is the bit I'm missing. How do you get the control points for the two halves from the control points of the whole? I haven't thought about this so perhaps its trivial, in which case it's clearly a good way to go. Even so, you can still use the Romberg approach on the results of two subdivisions. Its all about getting a better estimate based on two other estimates and knowledge of how the algorithm converges (and assumption that it does converge in a regular manner, but we're integrating a curve that is continuous and differentable)

Well the line segment summation is monotonic and bounded above, so that process must converge. The straight line segment between two points is always the shortest, so breaking it into two segments to better follow the curve must be longer (or equal).

the code... sorry this is going to have the formatting destroyed... note that I am making no attempt to optimize the Bezier point calculations since that wasn't the point I was making. Beziers can be calculated using just additions once the paraeters are convertd to the relevent form. Nor do I use integers. I simply fill out an array with the points befor hand to separate that from the integration. The Rombergy bit is just one line

-------------------------------------- // Bezier.cpp : Defines the entry point for the console application. //

#include "stdio.h" #include "stdlib.h" #include "math.h"

#define POINT_SHIFT 10 #define POINT_MAX (1

[snip]

OK.

De Casteljau's algorithm.

Given a Start point (S), departing control point (D), arriving control point (A), and End point (E), [i.e., these are P0 - P3], the control points for the two half curve are: M = (D+A)/2

L0 = S R3 = E

L1 = (S+D)/2 R2 = (A+E)/2

L2 = (L1+M)/2 R1 = (M+R2)/2

L3 = R0 = (L2+R1)/2

(it's easier to visualize if you draw it out)

Note that L3 (=R0) is a point on the original curve that you have just "created" at the expense of just a few adds and shifts. L0 and R3 are the original endpoints (S & E).

Yes.

OK, I'll try to take a look through it tomorrow. Thanks!

--don

[snip]

OK, so it is trivial for cubics. So, you can start with the first estimate as the average of the control point path and the start to end path. Then subdivide (binary, say) and get a new estimate. Then combine these estimates using the Romberg method. Assuming you know the convergence power you can also estimate the error. If the error is OK keep the result and pass it up, else recurse again. Somebody somewhere has to decide what good enough is, perhaps in fractional terms. The benefit of your scheme is that it doesn't need to recurse so deepy for some parts of the curve as others.

Peter

b0 = a0 b1 = (a0 + a1)/2 b2 = (a0 + 2*a1 + a2)/4

c0 = b3 = (a0 + 3*a1 + 3*a2 + a3)/8

c1 = (a3 + 2*a2 + a1)/4 c2 = (a3 + a2)/2 c3 = a3

Or, using midpoints only:

b0 = a0 c3 = a3

b1 = (a1 + a0)/2 c2 = (a3 + a2)/2

m = (a1 + a2)/2

b2 = (m + b1)/2 c1 = (m + c2)/2

c0 = b3 = (c1 + b2)/2

Yes, I'm familiar with the issues (I've had to write many floating point "support" libraries over the years). I think that this would be a bad fit for a conventional "one-size-fits-all" floating point implementation. Too many operations are simple arithmetic (add/sub/div2) and using floats makes those disproportionately expensive.

Exactly. You *know* you are going to compute so come up with a numerical representation that fits the complexity of the calculations regardless of how "clean" it is (gee, why bother supporting subnormals, NaN's, etc.?)

For "commodity" processors and/or semicustom logic, I think a fixed point (though not "integer") scheme is probably going to be the winner. I need to look at how to cut costs in the basic algorithm(s). Then, see if there are any tricks I can employ to further constrain the domain. E.g., perhaps scaling the coordinates of the points before running the algorithm and then "de"-scaling the results...

I like using serial adders/multipliers when designing dedicated processors -- i.e., where the algorithm is (literally?) hard-wired. It's small, easy to debug and scale, etc. E.g., handling 32 bit values is just as easy as 64 (just add more bits of "data store"). By contrast, going to wider parallel implementations means having to add *lots* of logic as you add each new bit.

Sorry, you've lost me (easy to do, this early in the morning :< )

Well, I know that the error can't be more than half of the difference between the "polygon length" and chord length. I.e., once they are within 2*error of each other, I *know* the length fits to within "one error" (?)

Well, it *could* not-need to recurse as deeply. It depends on how tight the curve is. E.g., a (pseudo) circle is probably the worst "overall" curve to follow as there are *no* "straight runs".

The problem witht he algorithm is that the cost of recursion can quickly escalate (i.e., the cost of the computations for each level of recursion is high). And, there is no way of bounding this a priori. :<

Oh! I thought I was reiterating what you were proposing to see if I had it right.

Yes, but that bound is quite loose.

Well, the curve can't be tight unless the poly is high order. The order of the poly, or number of control points hasn't been specified except that examples have been cubic. Anyway, if you specify that you want the final error be some small fraction of the initial upper bound then you can use that fraction on each division of the curve. Or you can assign half the limit to each half of the division. Not perfect because either way you'll probably do more work than strictly necessary (but you wouldn't know that until you'd done all the work...).

Sorry, my confusion revolves around your Romberg statement. I.e., how would that apply since I am no longer operating with "t"?

See above.

Yes, cubic. (sorry, I guess I had assumed as much)

Initial upper bound is essentially meaningless. I will probably express error (tolerance?) as some finite value. Or, proportional to the actual length (since very small curves may need tighter absolute errors)

Exactly.

Sorry, by "the cost of recursion can quickly escalate", I meant the cost of each "iteration" of the algorithm grows quickly. (i.e., "the cost of the computations for each level of recursion is high").

In other words, if the "current" estimate isn't good enough, you (I) have to subdivide to get the *next* better estimate. I can't get "just a little bit better"... I have to take the full step of a subdivision and computing all of the "lengths" that this entails. (note that some of the computations are trivial based on the lengths computed for the "current" subdivision/estimate)

I need to look at the (bezier) math and try to identify criteria that make "easy" curves vs. "hard" curves. Perhaps there is something that can be gleaned from that to condition the algorithm employed at each step. (e.g., look at the tangents at the various control/end -points wrt the chords themselves)

--don

It's got little to do with t pre se. The Romberg approach in a nutshell goes like this: Given a method of estimating a value (the pathlength) that is parameterised by a small value (the 'stepsize' h) such that error in the estimate decreases with decreasing h, and the estimate converges to the exact value as h -> 0. Then one can use estimates from two different values of h (say h and h/2) to produce a new estimate that converges more quickly. This assumes that the dominant error in the estimate can be expressed as a power of h. Implicitly we presume that the amount of work computing th estimate goes like 1/h or less, but that's not strictly necessary. The method doesn't say anything about how the estimates are computed.

In the your case of splitting the curve in half the first estimate can be thought of as having a stepsize 1, then splitting the curve and estimating the pathlength of each half and adding together has stepsixe 1/2. Repeating with another round of splitting the stepsize is 1/4. Now the Romberg estimate would be (4*estimate(1/2) - estimate(1))/3 assuming the error term is O(h^2).

Something concrete to assume...

OK, like so many ppm of the true answer or so many ppm of the original upper bound. I'm failrly sure that

1 >= length of curve/distance along the control points >= 0.5

I see this as the cost of the calculation not the recursion since the main difference between this method and stepping along the curve is the order that things are done, not how many.

Equivalent to restricting yourself to binary splits.

There aren't any hard curves. It's a cubic.

Peter

I did try out Gaussian quad but was disappointed. The problem is that the length element ds/dt that you need to integrat has a sqrt in it while Gauss does a good job for poly fits. I only tried a five point for lack of coefficent values.

Peter

Exactly. They are very smooth. It gives me the impression that somehow Gaussian quadrature would do a terrific job. They forego the need for small steps, with all the round off errors that result from it.

It requires to express the length as a function of one variable, something I don't know how to do for bezier curves. See 4.5 Gaussian Quadrature in Numerical Recipes.

Groetjes Albert

--

--
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

Ah, OK. So, it just tries to look-ahead "better" based on history (?). E.g., my estimate[i] are averages of polygon/chord length and endpoint length. Track successive estimate[i]'s and, at each step, compute an ESTIMATE[i] from Romberg? This would have no "resale value" in the next iteration of my algorithm but would (hopefully) be a better representation of the *real* "curve length" -- better than my "average of chords and endpoints" (?)

In some pathological cases, that's not true. (e.g., look at your original example, I think -- sorry, early in the morning to be doing this without pencil and paper...)

If you walk up the curve (t 0->1), you have to take a fixed dt for each pass. Even in regions of low curvature.

I can spend those calculations elsewhere -- in places where the curvature demands more work for a particular degree of fitness.

The cost of the actual recursion is tiny (another stack frame and function invocation) compared to all the number crunching that it (hopefully) avoids (or, spends in better places).

Sorry, I was trying to describe the (total) work required for particular "types" of curves.

E.g., (making up numbers entirely in my head :> )

(0,0) (0,0) (1000,1000) (1000,1000)

is an *easy* curve. The estimate converges instantly (regardless of how you -- rationally -- create that estimate).

OTOH,

(0,0) (0,400) (10,0) (0,0)

is a "hard" curve. The estimate requires multiple iterations to get a given "goodness".

The point of my statement (above) was to try to come up with some simple/algorithmic criteria (short of actually *solving* the curve) that would let you determine how hard you have to work on a particular curve (i.e., set of control points) to get a particular degree of fit.