Bezier lengths w/o lots of resources [Embedded]

Prev: Simple hack to get $600 to your home.
Next: Applied Computing 2010: 1st call extension until 25 June 2010

From: Boudewijn Dijkstra on 1 Jul 2010 10:17

Op Tue, 15 Jun 2010 03:23:44 +0200 schreef D Yuniskis
<not.going.to.be(a)seen.com>:
> Not the *ideal* group for this post -- but the "right"
> group would undoubtedly throw *way* more resources at
> the solution... resources that I don't *have*! :<
>
> I'm looking for a reasonably efficient (time + space)
> technique at approximating Bezier curve lengths from
> folks who may have already "been there"...
>
> I was thinking of subdividing the curve per De Casteljau.
> Then, computing the lengths of the individual segments
> obtained therefrom. Summing these as a lower limit
> on the curve length with the length of the control
> polygon acting as an upper limit.
>
> But, that's still a lot of number crunching! If I work
> in a quasi fixed point representation of the space, it
> still leaves me with lots of sqrt() operations (which
> are hard to reduce to table lookups in all but a trivial
> space! :< )
>
> So, obviously, I need a "trick". Something to throw away
> that gives me a "good enough" answer -- the sort of thing
> the "ideal" group wouldn't consider settling for ;-)
>
> Suggestions? Pointers?

Last week it dawned to me that any function, whether real or complex,
travels a path of which the distance can be calculated. (Except at
vertical asymptotes.) After some dabbling I came up with the following
definition.

The distance travelled by a function f(t) from t=t0 to t=t1 is given by:
length(f(t), t0, t1) =
int(abs(diff(f(t),t)),t=t0..t1)
where int=integrate, abs=absolute, diff=derivative

Note 1: the use of abs() yields sqrt() operations if f(t) is complex.

Note 2: for a complex function z(t)=x(t)+j*y(t):
length(z(t), t0, t1) != length(x(t), t0, t1) + length(y(t), t0, t1)
but it is usually a good guess for the upper bound.

Note 3: for the same complex function z(t):
length(z(t), t0, t1) != sqrt(length(x(t), t0, t1)^2 + length(y(t), t0,
t1)^2)
but it is usually a good ballpark estimate.

With the cubic Bezier curve function:

B[3](t) =
(1-t)^3*P[0]+3*(1-t)^2*t*P[1]+3*(1-t)*t^2*P[2]+t^3*P[3]

We get the following results when combining:

length(B[3](t), 0, 1) =
Int(abs(3*(1-t)^2*P[0]+6*(1-t)*t*P[1]-3*(1-t)^2*P[1]+3*t^2*P[2]-6*(1-t)*t*P[2]-3*t^2*P[3]),
t=0..1)

Unfortunately my old copy of Maple wasn't able to simplify this any
further. But the result can be calculated using numerical integration.

Note 4: if you know the roots of the x- and y-function (where the
derivative is zero), you can set up interpolation points.

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

From: Albert van der Horst on 1 Jul 2010 12:19

From: Peter Dickerson on 2 Jul 2010 06:15

"D Yuniskis" <not.going.to.be(a)seen.com> wrote in message
news:i0hl35$amk$1(a)speranza.aioe.org...
> Hi Peter,
>
> Peter Dickerson wrote:
>> "D Yuniskis" <not.going.to.be(a)seen.com> wrote in message
>> news:i0derm$en9$1(a)speranza.aioe.org...
>>
>>> OK, I made a stab at trying to "adjust" this for
>>> a more efficient algorithm (2.5/iteration).
>>
>> Sorry, I don't understand where the 2.5 sqrts/iteration comes from. When
>> you split the line and have two sets of four control points you need
>> eight distances (sqrt) but none of those have been computed at the
>> previous level. The only place I see five appear is that the are five new
>> points to calculate, one of which is used twice (giving six) and the
>> start and end points from the original (now we have eight). None of that
>> uses sqrts and the time to calculate it is potentially insignificant
>> relative to the sqrts (e.g. it can be done with integers of sufficient
>> size by rescaling by eight at each nesting).
>
> Given a Start point (S), departing control point (D), arriving
> control point (A), and End point (E), [i.e., these are P0 - P3]:
>
> The first iteration computes the lengths of |SD|, |DA|, |AE| and
> compares them against |SE| (i.e., 4 sqrts).
>
> Then the curve is subdivided into two half-curves ("Left" and
> "Right"). 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
>
> Looking at the Left half curve L0, L1, L2, L3:
>
> Compute |L0L1|, |L1L2| and |L2L3| and compare against |L0L3|.
>
> But, |L0L1| is actually just half of |AD| -- which we already
> have from the first iteration!

Not sure that's so. You mean |SD|, I guess.

> However, we do have to compute |L1L2|. And |L0L3|.
> There's two sqrts.
>
> OTOH, when we compute |L2L3|, this is really half of |L2R1|.
> So, whatever result we get there we can use for the |R0R1|
> computation! So, that's a 0.5 sqrt.
>
> (Did I do that right?)

OK. L2 is the middle of L2 and R1. An aweful lot of degneracy here.

So, you're say that you need five additional sqrts to compute one line, not
five sqrts. In my tables the sqrts column is the number of sqrts needed to
calculate the second column. But in you case you need to include *some of*
the counts higher up the table. Then when it comes to using Romberg on your
method you need all the sqrts done for the previous row whereas you only
need some to compute the current. So not directly comparable. In particular
if we have some method to estimate the depth needed ahead of time then my
method would not need to start from the top (although its a convenient way
to calculate the sub-curve endpoints). So I would only do sqrts that are
needed where as you'd need to do some earlier ones and pass them down.
Still, the naive method+Romberg is very similar work to your's without, and
botgh are O(h^4).

>>> Trusting your "results", I assume the table would be:
>>>
>>> sqrts DY length DY Romberg
>>> 4 1.000000000
>>> 5 1.002470605 1.002635312
>>> 10 1.000128079 0.999971911
>>> 20 1.000007988 0.999999982
>>> 40 1.000000499 1.000000000
>>>

see above.

>>> I.e., at this point, I have invested *79* sqrts (4+5+10+...).
>>> Note, however, that some of those might never occur if the
>>> estimate for a particular "subcurve" (ick!) looks "good
>>> enough" (it need not be further subdivided).
>>
>> You've still lost e.
>
> Assume I am correct with 2.5 sqrts per "half curve".
> And 4 on the first iteration.
>
> So, second iteration is 5 (2.5 + 2.5) sqrts (assuming
> I am not looking at the closeness of fit).
>
> Third iteration cuts each of the two half curves from
> the second iteration in half using 5 sqrts for each
> of the second iteration half curves -- 10 sqrts total.
>
> 20 for the next. etc.
>
> Since my algorithm starts with the *first* ITERATION
> and conditionally subdivides, it takes on the costs of
> each iteration as it subdivides the curve(s). So,
> my total cost is the sum of the costs of *every*
> iteration.
>
> I.e., my code computes "sum of segments" and "endpoint
> length"; checks to see how close they are to each other
> (or, does your magic Romberg estimate and compares
> *that* to these lengths to see how closely *it*
> fits); and, if the fit isn;t good enough, *then*
> it subdivides the curve and does another iteration of
> chord length computations. (but, it has already
> incurred the cost of the computations leading up to
> that decision to recurse!)

I'm just doing the calculation to a fixed depth to see how fast it
converges. There is no reason why I couldn't compare row against the
previous and stop when they are close enough. Close enough means meeting
some nearness criterion for the non-romberg results. That nearness will be
much larger than the actual error after Romberg.

>>> If I blindly force N subdivisions, it changes the mix.
>>>
>>> I'll try to make some time to code up a "real" algorithm
>>> and instrument SQRT() to do the counting for me with
>>> various "estimate tolerances".
>>
>> Surely you can identify any 'sharp' regions of the curve from the radius
>> of curvature (or at least its square, otherwise you'll need a sqrt for
>> that). A proxy might be dx/dt*d2y/dt2 - dy/dt*d2x/dt2 being small.

Peter

From: Albert van der Horst on 2 Jul 2010 14:15

In article <i0a96e$6dq$1(a)news.eternal-september.org>,
Peter Dickerson <first.last(a)tiscali.invalid> wrote:
>"Albert van der Horst" <albert(a)spenarnc.xs4all.nl> wrote in message
>news:l4qams.fv1(a)spenarnc.xs4all.nl...
>> In article <i01bmh$lgv$1(a)news.eternal-september.org>,
>> Peter Dickerson <first.last(a)tiscali.invalid> wrote:
>>>"D Yuniskis" <not.going.to.be(a)seen.com> wrote in message
>> <SNIP>
>>>
>>>> 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)
>>>
>>>There aren't any hard curves. It's a cubic.
>>
>> 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.
>
>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.

What did you do? I inserted a ds = sqrt( x'*x' +y'*y')dt
and got the following result for the 5 point (10th order) Gaussian:

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

Derive control points
(0.000000,1.000000)
(1.000000,0.000000)
(0.000000,-1.000000)

#ival length Romberg Gauss
1 1.00000000 ------------ 2.00000000
2 1.80277564 2.07036752 2.00000000
4 1.95071911 2.00003360 2.00000000
8 1.98771584 2.00004808 2.00000000
16 1.99693127 2.00000308 2.00000000
32 1.99923296 2.00000019 2.00000000
64 1.99980825 2.00000001 2.00000000
128 1.99995206 2.00000000 2.00000000
256 1.99998802 2.00000000 2.00000000
512 1.99999700 2.00000000 2.00000000
1024 1.99999925 2.00000000 2.00000000

Sorry, if the above looks like a fake ;-)
Run the program below if you don't believe it.
So doing integration of the whole interval at the
expense of 10 sqrt gives 8 digits precision,
good enough for government work. Splitting up the interval
doesn't improve or deteriorate precision.

Tabulating the function integrated shows a cusp,
the sqrt causes the derivative to change sign.
This should be bad for Gauss, but in this case the
discontinuity is at an interval boundary, but not
for 1 interval. Did we get lucky?

Well here is the code. No warranty of any kind.

---------------------8<----------------------------
//bezier.c : Defines the entry point for the console application.
// Trying to integrate the length of a cubic spline by different
// methods. Author : Albert van der Horst e.a. july 2010

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

/* From numerical recipes with improved accuracy constants
from LGPL-ed internet */
double qgauss( double (*func)(double), double a, double b )
{
int j;
double xr,xm,dx,s ;
static double x[] = {
0.0,
0.1488743389816312108848260,
0.4333953941292471907992659,
0.6794095682990244062343274,
0.8650633666889845107320967,
0.9739065285171717200779640
};
static double w[] = {
0.0,
0.2955242247147528701738930,
0.2692667193099963550912269,
0.2190863625159820439955349,
0.1494513491505805931457763,
0.0666713443086881375935688
};
xm=0.5*(b+a);
xr=0.5*(b-a);
s=0;
for (j=1; j<=5; j++)
{
dx = xr*x[j];
s += w[j] * ( (*func)(xm+dx) + (*func)(xm-dx) );
}
return s *= xr;
}

double qgauss2( double (*func)(double), double a, double b, int div )
{
double s = 0;
for (int i=0; i<div; i++)
s += qgauss( func, a+i*(b-a)/div, a+(i+1)*(b-a)/div );
return s;
}

#define POINT_SHIFT 10
#define POINT_MAX (1<<POINT_SHIFT)

struct POINT { double x; double y; };

POINT control[4];
POINT control_der[3];
POINT curve[POINT_MAX+1];
POINT deriv[POINT_MAX+1];
double lengths[POINT_SHIFT+1];

void fill_control()
{
// get control points
control[0].x = 0.0;
control[0].y = 0.0;

control[1].x = 0.0;
control[1].y = 1.0;

control[2].x = 1.0;
control[2].y = 1.0;

control[3].x = 1.0;
control[3].y = 0.0;
}

/* Fill control points for the derivative. */
void fill_control_der()
{
for (int i=0; i<3; i++)
{

// get control points
control_der[i].x = control[i+1].x - control[i].x ;
control_der[i].y = control[i+1].y - control[i].y ;
}
}

POINT curve_cal( double t )
{
double t1 = 1.0-t;
POINT p;

p.x = control[0].x*t1*t1*t1 + control[1].x*3.0*t1*t1*t +
control[2].x*3.0*t1*t*t + control[3].x*t*t*t;
p.y = control[0].y*t1*t1*t1 + control[1].y*3.0*t1*t1*t +
control[2].y*3.0*t1*t*t + control[3].y*t*t*t;
return p;
}

POINT deriv_cal( double t )
{
double t1 = 1.0-t;
POINT p;

p.x = control_der[0].x*3.0*t1*t1 + control_der[1].x*6.0*t1*t + control_der[2].x*3.0*t*t ;
p.y = control_der[0].y*3.0*t1*t1 + control_der[1].y*6.0*t1*t + control_der[2].y*3.0*t*t ;
return p;
}

double deriv_len( double t )
{
POINT p = deriv_cal( t );
/* return sqrt(fabs(p.x*p.y)); */
return sqrt(p.x*p.x + p.y*p.y);
}

int main(int argc, char* argv[])
{
fill_control();
fill_control_der();

if ( argc <= 4 )
{
for (int i=1; i<argc; i++)
{
if ( sscanf( argv[i], "(%f,%f)", &control[i].x, &control[i].y )
!= 2 )
{
printf( "Invalid argument: %s\n", argv[i] );
}
}

}

printf( "Control points\n" );
for (int i=0; i<4; i++)
printf( " (%f,%f)\n", control[i].x, control[i].y );
printf( "\n" );

printf( "Derive control points\n" );
for (int i=0; i<3; i++)
printf( " (%f,%f)\n", control_der[i].x, control_der[i].y );
printf( "\n" );

// fill in point array
for (int i=0; i<=POINT_MAX; i++)
{
double t = (double)i/POINT_MAX;
curve[i] = curve_cal( t );
}

// fill in deriv array
for (int i=0; i<=POINT_MAX; i++)
{
double t = (double)i/POINT_MAX;
deriv[i] = deriv_cal( t );
}

printf( "sqrts length Romberg Gauss \n" );
// compute lengths
for (int j=0; j<=POINT_SHIFT; j++)
{
int step = 1<<(POINT_SHIFT-j);
double length = 0.0;
for (int i=0; i<POINT_MAX; i += step)
{
double dx = curve[i+step].x - curve[i].x;
double dy = curve[i+step].y - curve[i].y;
length += sqrt( dx*dx+dy*dy );
}
lengths[j] = length;
printf( "%4d %12.8f", 1<<j, length );

// romberg improvement uses previous result, so not appricable to
// first time round
if ( j > 0 )
printf( " %12.8f", (lengths[j]*4 - lengths[j-1])/3.0 );
else
printf( " ------------" );
printf( "%12.8f\n", qgauss2( deriv_len, 0.0, 1.0, 1<<j ));
printf( "\n" );
}
/* Comment in if you want to tabulate the function integrated by Gauss.
for (int i=0; i<POINT_MAX; i += 32)
{
double t = i; t /= POINT_MAX;
printf( "%8.4f %12.8f", t, deriv[i].x*deriv[i].y );
printf( "\n" );
}
*/
}
---------------------8<----------------------------

For experimenting it would be better to allow for easier filling
in of the control points. Oh well. There must be something left
to do for others.

>
>Peter
>
>

Groetjes Albert

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

From: Peter Dickerson on 2 Jul 2010 15:03

"Albert van der Horst" <albert(a)spenarnc.xs4all.nl> wrote in message
news:l4y01b.lcc(a)spenarnc.xs4all.nl...
> In article <i0a96e$6dq$1(a)news.eternal-september.org>,
> Peter Dickerson <first.last(a)tiscali.invalid> wrote:
>>"Albert van der Horst" <albert(a)spenarnc.xs4all.nl> wrote in message
>>news:l4qams.fv1(a)spenarnc.xs4all.nl...
>>> In article <i01bmh$lgv$1(a)news.eternal-september.org>,
>>> Peter Dickerson <first.last(a)tiscali.invalid> wrote:
>>>>"D Yuniskis" <not.going.to.be(a)seen.com> wrote in message
>>> <SNIP>
>>>>
>>>>> 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)
>>>>
>>>>There aren't any hard curves. It's a cubic.
>>>
>>> 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.
>>
>>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.
>
> What did you do? I inserted a ds = sqrt( x'*x' +y'*y')dt
> and got the following result for the 5 point (10th order) Gaussian:
>
> Control points
> (0.000000,0.000000)
> (0.000000,1.000000)
> (1.000000,1.000000)
> (1.000000,0.000000)
>
> Derive control points
> (0.000000,1.000000)
> (1.000000,0.000000)
> (0.000000,-1.000000)
>
> #ival length Romberg Gauss
> 1 1.00000000 ------------ 2.00000000
> 2 1.80277564 2.07036752 2.00000000
> 4 1.95071911 2.00003360 2.00000000
> 8 1.98771584 2.00004808 2.00000000
> 16 1.99693127 2.00000308 2.00000000
> 32 1.99923296 2.00000019 2.00000000
> 64 1.99980825 2.00000001 2.00000000
> 128 1.99995206 2.00000000 2.00000000
> 256 1.99998802 2.00000000 2.00000000
> 512 1.99999700 2.00000000 2.00000000
> 1024 1.99999925 2.00000000 2.00000000
>
> Sorry, if the above looks like a fake ;-)
> Run the program below if you don't believe it.
> So doing integration of the whole interval at the
> expense of 10 sqrt gives 8 digits precision,
> good enough for government work. Splitting up the interval
> doesn't improve or deteriorate precision.
>
> Tabulating the function integrated shows a cusp,
> the sqrt causes the derivative to change sign.
> This should be bad for Gauss, but in this case the
> discontinuity is at an interval boundary, but not
> for 1 interval. Did we get lucky?
>
> Well here is the code. No warranty of any kind.
>
> ---------------------8<----------------------------
> //bezier.c : Defines the entry point for the console application.
> // Trying to integrate the length of a cubic spline by different
> // methods. Author : Albert van der Horst e.a. july 2010
>
> #include "stdio.h"
> #include "stdlib.h"
> #include "math.h"
>
> /* From numerical recipes with improved accuracy constants
> from LGPL-ed internet */
> double qgauss( double (*func)(double), double a, double b )
> {
> int j;
> double xr,xm,dx,s ;
> static double x[] = {
> 0.0,
> 0.1488743389816312108848260,
> 0.4333953941292471907992659,
> 0.6794095682990244062343274,
> 0.8650633666889845107320967,
> 0.9739065285171717200779640
> };
> static double w[] = {
> 0.0,
> 0.2955242247147528701738930,
> 0.2692667193099963550912269,
> 0.2190863625159820439955349,
> 0.1494513491505805931457763,
> 0.0666713443086881375935688
> };
> xm=0.5*(b+a);
> xr=0.5*(b-a);
> s=0;
> for (j=1; j<=5; j++)
> {
> dx = xr*x[j];
> s += w[j] * ( (*func)(xm+dx) + (*func)(xm-dx) );
> }
> return s *= xr;
> }
>
> double qgauss2( double (*func)(double), double a, double b, int div )
> {
> double s = 0;
> for (int i=0; i<div; i++)
> s += qgauss( func, a+i*(b-a)/div, a+(i+1)*(b-a)/div );
> return s;
> }
>
> #define POINT_SHIFT 10
> #define POINT_MAX (1<<POINT_SHIFT)
>
> struct POINT { double x; double y; };
>
> POINT control[4];
> POINT control_der[3];
> POINT curve[POINT_MAX+1];
> POINT deriv[POINT_MAX+1];
> double lengths[POINT_SHIFT+1];
>
> void fill_control()
> {
> // get control points
> control[0].x = 0.0;
> control[0].y = 0.0;
>
> control[1].x = 0.0;
> control[1].y = 1.0;
>
> control[2].x = 1.0;
> control[2].y = 1.0;
>
> control[3].x = 1.0;
> control[3].y = 0.0;
> }
>
> /* Fill control points for the derivative. */
> void fill_control_der()
> {
> for (int i=0; i<3; i++)
> {
>
> // get control points
> control_der[i].x = control[i+1].x - control[i].x ;
> control_der[i].y = control[i+1].y - control[i].y ;
> }
> }
>
> POINT curve_cal( double t )
> {
> double t1 = 1.0-t;
> POINT p;
>
> p.x = control[0].x*t1*t1*t1 + control[1].x*3.0*t1*t1*t +
> control[2].x*3.0*t1*t*t + control[3].x*t*t*t;
> p.y = control[0].y*t1*t1*t1 + control[1].y*3.0*t1*t1*t +
> control[2].y*3.0*t1*t*t + control[3].y*t*t*t;
> return p;
> }
>
> POINT deriv_cal( double t )
> {
> double t1 = 1.0-t;
> POINT p;
>
> p.x = control_der[0].x*3.0*t1*t1 + control_der[1].x*6.0*t1*t +
> control_der[2].x*3.0*t*t ;
> p.y = control_der[0].y*3.0*t1*t1 + control_der[1].y*6.0*t1*t +
> control_der[2].y*3.0*t*t ;
> return p;
> }
>
> double deriv_len( double t )
> {
> POINT p = deriv_cal( t );
> /* return sqrt(fabs(p.x*p.y));
> */
> return sqrt(p.x*p.x + p.y*p.y);
> }
>
> int main(int argc, char* argv[])
> {
> fill_control();
> fill_control_der();
>
> if ( argc <= 4 )
> {
> for (int i=1; i<argc; i++)
> {
> if ( sscanf( argv[i], "(%f,%f)", &control[i].x, &control[i].y )
> != 2 )
> {
> printf( "Invalid argument: %s\n", argv[i] );
> }
> }
>
> }
>
> printf( "Control points\n" );
> for (int i=0; i<4; i++)
> printf( " (%f,%f)\n", control[i].x, control[i].y );
> printf( "\n" );
>
> printf( "Derive control points\n" );
> for (int i=0; i<3; i++)
> printf( " (%f,%f)\n", control_der[i].x, control_der[i].y );
> printf( "\n" );
>
>
> // fill in point array
> for (int i=0; i<=POINT_MAX; i++)
> {
> double t = (double)i/POINT_MAX;
> curve[i] = curve_cal( t );
> }
>
> // fill in deriv array
> for (int i=0; i<=POINT_MAX; i++)
> {
> double t = (double)i/POINT_MAX;
> deriv[i] = deriv_cal( t );
> }
>
> printf( "sqrts length Romberg Gauss \n" );
> // compute lengths
> for (int j=0; j<=POINT_SHIFT; j++)
> {
> int step = 1<<(POINT_SHIFT-j);
> double length = 0.0;
> for (int i=0; i<POINT_MAX; i += step)
> {
> double dx = curve[i+step].x - curve[i].x;
> double dy = curve[i+step].y - curve[i].y;
> length += sqrt( dx*dx+dy*dy );
> }
> lengths[j] = length;
> printf( "%4d %12.8f", 1<<j, length );
>
> // romberg improvement uses previous result, so not appricable to
> // first time round
> if ( j > 0 )
> printf( " %12.8f", (lengths[j]*4 - lengths[j-1])/3.0 );
> else
> printf( " ------------" );
> printf( "%12.8f\n", qgauss2( deriv_len, 0.0, 1.0, 1<<j ));
> printf( "\n" );
> }
> /* Comment in if you want to tabulate the function integrated by Gauss.
> for (int i=0; i<POINT_MAX; i += 32)
> {
> double t = i; t /= POINT_MAX;
> printf( "%8.4f %12.8f", t, deriv[i].x*deriv[i].y );
> printf( "\n" );
> }
> */
> }
> ---------------------8<----------------------------

obviously I did something wrong. I was expecting great things but got poor
results. Perhaps a typo in one of the constants. However, 5 points ought to
be 5 sqrts...

peter

First | Prev | Next | Last
Pages: 1 2 3 4 5 6 7 8 9 10 11 12
Prev: Simple hack to get $600 to your home.
Next: Applied Computing 2010: 1st call extension until 25 June 2010