Polybezier Curve Fitting

[GWS]

OK .. now for a bit of technical stuff ..    courtesy of Abel Nov 2004 ..

"The API contains a graphics function, PolyBezier, that can be most useful to those who would scheme and plot. A cursory search of these forums and sample programs didn't reveal any examples and I thought a brief expo might be of interest.

The essential problem is that we have a display of a set of data points that we wish to connect. In zeroth order, we might assume a constant value over each interval, a bar chart. In first order, straight lines. PolyBezier lets us connect these points with cubic splines or 3rd order polynomials. These suffice for an aesthetic rendering to the human eye of a 'smooth' curve.

For simplicity, we'll assume that the points are equally spaced along the x-axis, e.g. x(i) = 0, 1, 2, ...
Within each interval with x normalized, y(x) is described by a cubic polynomial,
y(x) = a + b*x + c*x^2 + d*x^3 (0<x<1)

Requiring that the curve pass through the data endpoints of each interval sets two conditions on a solution for a,b,c,d. A third condition we'll adopt is that the slope be continuous at each endpoint so that the initial slope for the (i+1)th interval matches the final slope of the (i)th interval. We still need a 4th condition for a solution, so lets posit that extension of the curve to x(i+2) also passes through y(i+2). For a complete solution we require additional assumptions about the first and last intervals, and we'll assume zero slopes at the plot's extrema. Note that how we choose to depict the curve between any two points is a matter of aesthetics and personal taste - there isn't any 'right' way (we could as well specified that c=d=0 and gotten straight lines).

The solution for a-d is thus given by:
y(i) = a
y(i+1) = a + b + c + d
y(i+2) = a + 2*b + 4*c + 8*d
b = b(i-1) + 2*c(i-1) + 3*d(i-1)
For the last equation, the right-hand values are taken from the solution for the preceding interval.

Implementation of the PolyBezier functions requires a pointer to an array of integer values sequenced x(0),y(0),x(1),y(1),... Each 3rd point will be a true data point, i.e. i=0,3,6,... and the two intervening points are called control points whose coordinates are obtained from solutions for a,b,c,d. As all data points are equally spaced on the x-axis, the normalized Bezier solutions for the two control points are:
x = x(i) + 1/3, y = y(i) + b/3
x = x(i) + 2/3, y = y(i) + (2*b + c)/3

Finally, a little IBasic trick to construct the necessary array:
DEF graph[2,(3*n-2)]:int

where n is the number of data points and
graph[0,i]=x(i)
graph[1,i]=y(i)

Abel "

Code Select Expand


AUTODEFINE "OFF"
DECLARE "gdi32",PolyBezier(hdc:int,points[]:int,count:int),int
DEF n,max,run,hdc:int
n=20:'Number of points
max=30:'Maximum for plot
'Assign n arbitrary values less than max
DEF data[n]:double
LET data[0], 00.0,02.5,05.0,07.5,10.0
LET data[5], 06.6,15.5,27.7,11.4,10.0
LET data[10],08.0,06.0,04.0,02.0,00.0
LET data[15],07.0,06.0,05.0,04.0,00.0
DEF graph[2,(3*n-2)]:int
DEF i,j,L,T,W,H:int
DEF dx,dy,a,b,c,d:double
b=0:' initial slope
DEF wnd:WINDOW
WINDOW wnd, 100, 100, 400, 400, @CAPTION, 0, "Bezier Plotting", procedure
SETWINDOWCOLOR wnd, 0x00
FRONTPEN wnd,0xFF
GETCLIENTSIZE wnd,L,T,W,H
dx=W/(3*n-3):'Horizontal pixel scale
dy=H/max:'    Vertical pixel scale
'All x coordinates, control points included, equally spaced
FOR i=0 to 3*n-3
  graph[0,i]=dx*i
NEXT i
'Solve for a,b,c,d in the interval x(i)-> x(i+1)
'  where y(x) = a + b*x +c*x^2 +d*x^3
'Cubic spline constraints are:
'  1. The solution passes through points y(i), y(i+1), and y(i+2)
'  2. The slope is continuous across intermediate points
'  3. The slope is zero at the end points
j=0
FOR i=0 TO n-1
  graph[1,j]=H-dy*data[i]
  circle wnd,graph[0,j],graph[1,j],3
  a=data[i]
  SELECT 1
    CASE (i<(n-2))
      c=0.25*(-7*a+8*data[i+1]-data[i+2]-6*b)
      d=data[i+1]-a-b-c
      graph[1,j+1]=H-dy*(a+b/3):'       Control point 1
      graph[1,j+2]=H-dy*(a+2*b/3+c/3):' Control point 2
    CASE (i=n-2):'Special case for last interval forcing zero slope
      c=-3*a-2*b+3*data[i+1]
      graph[1,j+1]=H-dy*(a+b/3)
      graph[1,j+2]=H-dy*(a+2*b/3+c/3)
  ENDSELECT
  b=b+2*c+3*d
  j=j+3
NEXT i
hdc=GETHDC(wnd)
PolyBezier(hdc,graph[0,0],3*n-2)
ReleaseHDC wnd,hdc
run=1:WAITUNTIL run=0:CLOSEWINDOW wnd:END
SUB procedure
SELECT @CLASS
  CASE @IDCLOSEWINDOW
    run=0
ENDSELECT
RETURN



"Those who've had the time to play with the code above may have discovered that, when plotting discontinuous data such as a step function or delta function (a single spike above the baseline), a ringing trails the discontinuity. As discussed, four conditions are needed to calculate the shape of each plotted segment. Three are spent specifying that the curve pass through the segment endpoints and maintain a continuous slope across segment boundaries.

The 4th condition leaves room for creativity. For simplicity, it was originally assumed that the curve segment calculated to pass through y(i) and y(i+1) should also pass through y(i+2). This leads to a forward-looking asymmetry that may account for some ringing. From a more balanced perspective, we should have a constraint that the curve also pass through y(i-1), but this requires quartic splines and these aren't built into the API.

Instead, we might assume something about the distances by which the curve segment misses these adjacent points, i.e.
{y(i-1) - a + b - c + d} and {y(i+2) - a -2b -4c -8d}.

When one assumes the algebraic sum of these two differences is zero (codex infra), the ringing in plots of simple discontinuous functions is suppressed and symmetric and asymmetric data are more cleanly represented. The rationale for such an assumption is that, at large distances from a given segment, a cubic spline diverges towards opposing infinities and setting this sum to zero attempts a balance of the divergencies left and right of the given segment. A statistician might see opportunity for better weighting including points yet further removed. This approximation may suffice for 99% of plotters but, for the remaining 1%, perhaps there's room for more thorough analysis.

Abel
"The purpose of computing is insight, not numbers." - R. Hamming "

Revised code:

Code Select Expand


AUTODEFINE "OFF"
DECLARE "gdi32",PolyBezier(hdc:int,points[]:int,count:int),int
DEF n,max,run,hdc:int
n=20:'Number of points
max=30:'Maximum for plot
'Assign n arbitrary values less than max
DEF data[n]:double
LET data[0], 00.0,02.5,05.0,07.5,10.0
LET data[5], 06.6,15.5,27.7,11.4,10.0
LET data[10],08.0,06.0,04.0,02.0,00.0
LET data[15],07.0,06.0,05.0,04.0,00.0
DEF graph[2,(3*n-2)]:int
DEF i,j,L,T,W,H:int
DEF dx,dy,a,b,c,d:double
b=0:' initial slope
DEF wnd:WINDOW
WINDOW wnd, 100, 100, 400, 400, @CAPTION, 0, "Bezier Plotting", procedure
SETWINDOWCOLOR wnd, 0x00
FRONTPEN wnd,0xFF
GETCLIENTSIZE wnd,L,T,W,H
dx=W/(3*n-3):'Horizontal pixel scale
dy=H/max:'    Vertical pixel scale
'All x coordinates, control points included, equally spaced
FOR i=0 to 3*n-3
  graph[0,i]=dx*i
NEXT i
'Solve for a,b,c,d in the interval x(i)-> x(i+1)
'  where y(x) = a + b*x +c*x^2 +d*x^3
'Cubic spline constraints are:
'  1. The solution passes through points y(i), y(i+1), and y(i+2)
'  2. The slope is continuous across intermediate points
'  3. The slope is zero at the end points
'Delta and Step
LET data[0], 03.0,03.0,03.0,03.0,03.0
LET data[5], 25.0,03.0,03.0,03.0,03.0
LET data[10],03.0,03.0,25.0,25.0,25.0
LET data[15],25.0,25.0,25.0,25.0,25.0
'Replacement for the original loop
'The original algorithm is included but commented out to facilitate comparison.
j=0
FOR i=0 TO n-1
  graph[1,j]=H-dy*data[i]
  circle wnd,graph[0,j],graph[1,j],3
  a=data[i]
  SELECT 1
    CASE (i=0):'Special case for first segment
      c=0.25*(-7*a+8*data[i+1]-data[i+2]-6*b)
    CASE (i<(n-2))
'     c=0.25*(-7*a+8*data[i+1]-data[i+2]-6*b):' former algorithm
      c=a-3*b+0.5*(7*(data[i+1]-a)-data[i-1]-data[i+2])
    CASE (i=n-2):'Special case for last segment
      c=-3*a-2*b+3*data[i+1]
  ENDSELECT
  IF i<(n-1)
    d=data[i+1]-a-b-c
    graph[1,j+1]=H-dy*(a+b/3):'       Control point 1
    graph[1,j+2]=H-dy*(a+2*b/3+c/3):' Control point 2
  ENDIF
  b=b+2*c+3*d
  j=j+3
NEXT i

hdc=GETHDC(wnd)
PolyBezier(hdc,graph[0,0],3*n-2)
ReleaseHDC wnd,hdc
run=1:WAITUNTIL run=0:CLOSEWINDOW wnd:END
SUB procedure
SELECT @CLASS
  CASE @IDCLOSEWINDOW
    run=0
ENDSELECT
RETURN



Phew ..  .. that should keep you busy ..

best wishes,

Graham

[aurelCB]

Very good Graham,one question why you use LET in:

Code Select Expand

DEF data[n]:double
LET data[0], 00.0,02.5,05.0,07.5,10.0
LET data[5], 06.6,15.5,27.7,11.4,10.0
LET data[10],08.0,06.0,04.0,02.0,00.0
LET data[15],07.0,06.0,05.0,04.0,00.0
regards
zlatko

[GWS]

It's an alternative assignment method - no longer mentioned in the user's guide.

Originally, Paul suggested it was just a little faster than the a = b method.

From Paul's Asteroids program :

'load our embedded explosion sound
'before anyone asks....
'LET a,b is the same as a = b
'its just a little faster.
SUB initsound
DEF snd[4780],value:char
let snd,82,73,70,70,164,18,0,0,87,65
let snd[10],86,69,102,109,116,32,16,0,0
let snd[19],0,1,0,1,0,17,43,0,0

best wishes,

Graham

[aurelCB]

Yeah,i now that is little faster but I use this method for clear strings in parser in AB becose I use
Autodefine On and i think that work better this hard coded way or maby i'm wrong?