PS Bezier Curve through Fuzzy Data!
By Don Lancaster                                                                   Version 3.2 February 14, 1998
Copyright c. 1998 by Don Lancaster and Synergetics, Box 809, Thatcher AZ, 85552
(520) 428-4073. synergetics@tinaja.com All commercial rights and all electronic media
rights are *fully* reserved. Linking welcome. Reposting is expressly forbidden.

Consulting services available.
Further support on www.tinaja.com/info01.html  

(The following is believed correct. Please report any errors or differing experiences.)

%

% SUMMARY: Guru Don Lancaster gives us a tutorial and heavy duty set of
% utilities that allow you to convert noisy data points into
% smooth cubic splines. This is a preliminary release.
%
% This is especially useful in capturing fonts from bitmaps,
% plotting engineering curves, or smoothing scanned artwork.
%
% Yes, the key constraint problems of preserving the input
% or angle and staying accurate on the end points are specifically
% handled.
%

%%%%%%%%%%%% (A) INSTALL AND RUN GONZO UTILITIES %%%%%%%%%%%%%%%%%%

% The demo uses my gonzo utilities from
% http://www.tinaja.com/psutils/gonzo20.ps

(C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities

gonzo begin % activate the gonzo utilities
ps.util.1 begin
% printerror
nuisance begin

%%%%%%%%%%%% (B) INTRO AND THEORY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% How do you take a bunch of noisy dots and run a smooth curve through
% them? This problem crops up often when designing typogrophy
% or doing engineering charts and graphs.

% Even trickier, you can't just tie cubic splines together. More
% often than not, you want to add CONSTRAINTS. Usually, adjacent splines
% should have matching endpoints. For best appearance, the entry and exit
% angles should match. Ideally, the rate of change of curvature should
% also match at the joints.

% One VERY handy tool: given four or more noisy data points within
% the range of what a single cubic spline can handle, find a very good
% Bezier curve that precisely goes through the end points, while smoothly
% coming reasonably near all of the others. AND starts off at a precisely
% selected angle, letting you smoothly continue your curve.

% This key tool would solve the general problem of taking a noisy
% character shape or an engineering graph, and drawing an efficient
% group-of-splines PostScript curve through everything.

% One simple older tool set involved the CURVETRACE routines found in
% GONZO20.PS and in the older utilities. See our classic MEOWWRRR.PS
% for a use example. The curvetracing takes a series of data values
% of x0-y0-ang0, x1-y1-ang1, up through ... xn-yn-angn. They draws some
% "fairly weak" splines through the selected points, using one spline per
% data triad. Slopes and endpoints all match.

% While very useful for cartoons and digitized signatures, my curvetrace
% routines are intolerant of fuzzy or noisy data. Although they form
% sparse data set, curvetracing uses lots of weak splines, instead of
% a few good ones. Angles also have to be specified for every point.

% Now, the CORRECT way to handle the problem appears in CURVEFITTING WITH
% PIECEWISE PARAMETIC CUBICS by Michael Plass and Maureen Stone (with a
% little help by some guy named Warnock). This appeared in the July 1983
% issue of COMPUTER GRAPHICS on pages 229-239, and otherwise known as the
% 1983 SIGGRAPH PROCEEDINGS.

% The only little problem I had with this paper was doing the vector
% Newton-Raphson iterations of the Hermite basis functions. Seriously, it
% was painfully obvious that bunches of time and effort would be required
% to master the postgraduate math involved. Especially while translating
% everything to convenient PostScript end user procs.

% So, what I thought we would do here is some dumb slow tools that may
% let you do the same stuff in a more intuitive fashion. What follows
% is open ended, intended as both tools you can use now to smooth fuzzy
% data, and as intermediate steps to fancier code that uses more
% elegant math.

% Naturally, once you smooth your data once, you steal the plans and just
% use the final splines for all repeat copies of the same job.

% We will attack three different constraint problems: Type A holds
% the entry angle and the end points. Type B holds the exit angle and
% the end points. Type C holds the end points only.

% Now let's pick a bunch of slightly noisy data points, arranged
% as [x0 y0 x1 y1 ... xn yn]. We want at least four points, but can have
% as many as you need. The assumption is made that one spline can handle
% this portion of the curve we are trying to fit. One spline can produce
% a smoothly flowing curve that has at most one cusp, one loop, or two
% inflection points. For brevity, we'll use whole numbers here ...

/rawdata [2 3 5 8 8 14 11 17 14 17 16 15 18 11] def

% Next, let's draw a fancy gray grid...

100 100 10 setgrid 20 20 showgrid

% do a simple Gonzo dot-plotter...

/dot { currentpoint newpath 0.150 0 360 arc fill } def

/plotdots {/mat exch def 0 2 mat length 1 sub {/posn exch def mat
posn get mat posn 1 add get mt dot} for} def

rawdata plotdots

% and optionally look at them...

% copypage

% Let's assume our previous spline left us with a required entrance angle % of sixty degrees...

/entang 60 def

% The key to any successive approximation scheme is to start out with
% something reasonable as our first guess. As a wild guess, and
% remembering our constraints, let's try a rule of...

% We'll keep things simple and obvious here for a while, rather than
% going for speed or compact code...

/mat {rawdata} def

/delta1 {mat 2 get mat 0 get sub dup mul mat 3 get mat 1 get sub
dup mul add sqrt} def

/initsplineguess{

/x0a {mat 0 get} def /y0a {mat 1 get} def

/x1a {delta1 entang cos mul 2 mul x0a add } def /y1a {delta1 entang
sin mul 2 mul y0a add} def

/x2a {mat dup length 4 sub get dup mat dup length 2 sub get sub
2 mul add} def /y2a {mat dup length 3 sub get dup mat dup length
1 sub get sub 2 mul add} def

/x3a {mat dup length 2 sub get} def /y3a {mat dup length 1 sub
get} def } def

initsplineguess

% Let's see what our first cut looks like...

line1

/showspline {x0a y0a moveto x1a y1a x2a y2a x3a y3a curveto stroke} def

showspline

grestore % get off grid for now

showpage % show initial guess on first page

% Which ain't half bad for amateurs.

%%%%%%%%%%%%%%% (C) POINT OPTIMIZATION %%%%%%%%%%%%%%%%%%%

% Now, instead of ultra fancy math,
% let's use some slow and intuitive successive approximations to get
% us started understanding what's coming down.

% The game plan from here: We find how close the curve comes to each one
% of the intermediate points. Then we find the total sum of the errors
% squared for the points. Then we adjust x1, y1, x2, and y2 incrementally
% to minimize this error. In this case, x1 and Y1 are locked together by
% entang. With some luck, things will rapidly converge for resonable
% types of curves that are likely to occur.

% But first a review: There are two forms for a Bezier curve. These are
% the graph space and the equation space.

% In graph space, a Bezier curve starts at x0,y0 and ends at x3 y3.
% It leaves x0,y0 with a direction and an enthuasiasm set by the first
% control point x1,y1 and enters x3,y3 with a direction and an
% enthuasiasm set by the second control point x2,y2.

% Here are the same Bezier formulas in the algebraic equation space:
%
%          x(t) = A(t^3) + B(t^2) + C(t) + x0
%          y(t) = D(t^3) + E(t^2) + F(t) + y0
%
% Note that (t) varies from 0 to 1 from the % beginning to the end of your spline curve.

% You can think of a Bezier spline as a three dimensional spiral in
% x, y, and z space. As you go from the initial values of x and y
% to the final values of x and y, t goes exactly from 0 to 1. It is
% sometimes useful to think of t as TIME, starting at t=0 and ending
% with t=1.

% As the red book tells us, here is how to get from equation space
% to graph space:
%
%     x1 = x0 + C/3                 y1 = y0 + F/3
%     x2 = x1 + C/3 + B/3           y2 = y1 + F/3 + E/3
%     x3 = x0 + C + B + A           y3 = y0 + F + E + D

% You can use simple algebra to get the far more useful inverse form.
% Here is how to get from graph space back to equation space ..
%
%     A = x3 - 3x2 + 3x1 + x0       D = y3 - 3y2 + 3 y1 + y0
%     B = 3x2 - 6x1 + 3x0           E = 3y2 - 6y1 + 3y0
%     C = 3x1 - 3x0                 F = 3y1 - 3y0

% This grabs six graph values off the stack...

/grabgraf {/y3 exch def /x3 exch def /y2 exch def /x2 exch def
/y1 exch def /x1 exch def /y0 exch def /x0 exch def} def

% Having predefined all four x and y values this finds x(t) given tt in
% the range of 0 to 1 ...

/xtt {x3 x2 3 mul sub x1 3 mul add x0 sub tt 3 exp mul x2 3 mul x1 6 mul
neg add x0 3 mul add tt dup mul mul add x1 3 mul x0 3 mul neg add tt
mul add x0 add} def

% this finds y(t) given tt in the range of 0 to 1 ...

/ytt {y3 y2 3 mul sub y1 3 mul add y0 sub tt 3 exp mul y2 3 mul y1
6 mul neg add y0 3 mul add tt dup mul mul add y1 3 mul y0 3 mul neg
add tt mul add y0 add} def

% this finds the magnitude of the error distance between a given tt and
% points x, y on the stack

/finderrdist {ytt sub exch xtt sub dup mul exch dup mul add sqrt} def

% In general, tt does NOT travel uniformly from 0 to 1. It tends to
% advance more violently on the "more bent" portions of the curve.

% Let's guess a tt value for each data point, initially assuming an
% equal tt spacing...

/initttguess {/ttguess [ 0 mat length 2 div dup 1 exch div /inc exch
def 2 sub cvi {dup inc add} repeat 1 ] def} def

initttguess

% Now, let's improve each tt value by incrementing it in distinc steps to
% find the minimum distance. What we are doing here is sliding along the
% t space until we are as close to each noisy point as we can get.

% We will try half percent tt steps initially. This seems good enough for
% now...

/distinc .005 def

% We first try a direction to see if the error gets larger or smaller.
% Then we continue making the error smaller till we overshoot. After
% an overshoot, we use the "parabola fit interpolation" of Appendix A
% to put a better guess into the ttguess array.

/improvettguess { /tt curtt def /curinc distinc def mat curpos 2 mul
get /curx exch def mat curpos 2 mul 1 add get /cury exch def mark
curx cury finderrdist /tt tt curinc add def curx cury finderrdist
exch dup 2 index sub 3 1 roll 2 index 0 gt {pop}{exch pop exch neg
exch /curinc curinc neg def /tt tt curinc add def} ifelse /tt tt
curinc add def {curx cury finderrdist dup 2 index exch sub 3 -1 roll
pop exch 1 index 0 gt {/tt tt curinc add def}{pop exit} ifelse}
loop curinc 0 lt {exch /curinc curinc neg def /tt tt curinc 2 mul
add def} if interpolate tt add distinc sub dup /tt exch def
ttguess exch curpos exch put cleartomark} def

% Here's the interpolation used (see Appendix A) ...

/interpolate {neg exch div dup 1 exch sub exch 1 add dup 0 eq
{0.00000001} if div distinc mul 2 div} def

% this takes care of optomizing all the inside values, given the eight
% graphspace Bezier points...

/setnearest { grabgraf 1 1 ttguess length 2 sub {/curpos exch def
ttguess curpos get /curtt exch def improvettguess} for} def

x0a y0a x1a y1a x2a y2a x3a y3a setnearest

% At this point, the ttguess array has values that accurate approximate
% the nearest point on the x(t), y(t) curve to the points to be fitted

% finderrors measures the error distances for all the intermediate points
% and returns an array of those actual distances

/finderrors {mark 1 1 ttguess length 2 sub {/posn exch def /tt ttguess
posn get def mat posn 2 mul get mat posn 2 mul 1 add get finderrdist }
for] /currenterrors exch def} def

finderrors

% Executing finderrors give us a value of .146415 for the first point
% first pass, which "looks about right" and converges to the .146321
% actual distance found at higher resolution. Other values of .569515,
% .846671, .789305, and .442002 also "look about right" on the graph.

% To find the total sum-of-squares error, each error value is squared
% and added...

/findsoserror { 0 currenterrors {dup mul add} forall /currenterror
exch def } def

% findsoserror

% Our current sum-of-squares error is 1.86033.

% The key tools are now in place. We can find the error for any given
% Bezier curve fit to our noisy data. Now the fun begins...

% In general, there are eight variables in the Bezier graph space for
% any given spline. Fixing the end points fixes four of these variables.
% fixing the entang entry angle locks x1 and y1 to each other. We
% thus have three variables left that we can adjust to try and get the
% best fit we can.

% The trial increment is the normalized percentage change we are going
% to make for each test improvement...

/trialinc .005 def % use half percent increment steps

% This finds the length of the end-to-end vector. It will be handy to
% normalize the incremental test steps...

/chord { mat 0 get mat dup length 2 sub get sub dup mul mat 1 get mat
dup length 1 sub get sub dup mul add sqrt} def

% This combines the previous utilities into a single error finder

/findcurrenterror {x0a y0a x1a y1a x2a y2a x3a y3a setnearest
finderrors findsoserror currenterror} def

% findcurrenterror % uncomment to activate

% these three routines are used to optimize x2a. x2a is incremented once
% to determine a direction. It is then incremented in the correct
% direction until a minimum is found...

/incx1a {/x1a x1a trialinc chord mul add def} def /decx1a {/x1a x1a
trialinc chord mul sub def} def

/optimizex1 {mark findcurrenterror incx1a findcurrenterror 2 copy gt
{{incx1a findcurrenterror 2 copy lt{exit} if} loop decx1a } { exch
decx1a {decx1a findcurrenterror 2 copy lt{exit} if} loop incx1a}
ifelse cleartomark} def

% optimizex1 % uncomment to activate

% these three routines are used to optimize x2a. x2a is incremented once
% to determine a direction. It is then incremented in the correct
% direction until a minimum is found...

/incx2a {/x2a x2a trialinc chord mul add def} def /decx2a {/x2a
x2a trialinc chord mul sub def} def

/optimizex2 {mark findcurrenterror incx2a findcurrenterror 2 copy gt
{{incx2a findcurrenterror 2 copy lt{exit} if} loop decx2a }
{ exch decx2a {decx2a findcurrenterror 2 copy lt{exit} if} loop
incx2a} ifelse cleartomark} def

% optimizex2 % uncomment to activate

% these three routines are used to optimize y1a. Y1a is incremented once
% to determine a direction. It is then incremented in the correct
% direction until a minimum is found...

/incy1a {/y1a y1a trialinc chord mul add def} def /decy1a {/y1a y1a
trialinc chord mul sub def} def

/optimizey1 {mark findcurrenterror incy1a findcurrenterror 2 copy
gt {{incy1a findcurrenterror 2 copy lt{exit} if} loop decy1a }
{ exch decy1a {decy1a findcurrenterror 2 copy lt{exit} if} loop
incy1a} ifelse cleartomark} def

% optimizey1 % uncomment to activate

% these three routines are used to optimize y2a. Y2a is incremented once
% to determine a direction. It is then incremented in the correct
% direction until a minimum is found...

/incy2a {/y2a y2a trialinc chord mul add def} def /decy2a {/y2a y2a
trialinc chord mul sub def} def

/optimizey2 {mark findcurrenterror incy2a findcurrenterror 2 copy gt
{{incy2a findcurrenterror 2 copy lt{exit} if} loop decy2a } { exch
decy2a {decy2a findcurrenterror 2 copy lt{exit} if} loop incy2a}
ifelse cleartomark} def

% optimizey2 % uncomment to activate

% these three routines are used to optimize x1a and y1a, WHEN THEY ARE
% CONSTRAINED TO LIE ALONG THE ENTANG AXIS. They get incremented once
% to determine a direction. They are then incremented in the correct
% direction until a minimum is found...

/incx1ay1a {/x1a x1a trialinc chord mul entang cos mul add def
/y1a y1a trialinc chord mul entang sin mul add def} def

/decx1ay1a {/x1a x1a trialinc chord mul entang cos mul sub def /y1a
y1a trialinc chord mul entang sin mul sub def} def

/optimizex1y1 {mark findcurrenterror incx1ay1a findcurrenterror 2 copy gt

{ {incx1ay1a findcurrenterror 2 copy lt{exit} if} loop decx1ay1a }

{ exch decx1ay1a {decx1ay1a findcurrenterror 2 copy lt{exit} if}
loop incx1ay1a

} ifelse cleartomark } def

% these three routines are used to optimize x2a and y2a, WHEN THEY ARE
% CONSTRAINED TO LIE ALONG THE EXTANG AXIS. They get incremented once
% to determine a direction. They are then incremented in the correct
% direction until a minimum is found...

/incx2ay2a {/x2a x2a trialinc chord mul entang cos mul add def /y2a
y2a trialinc chord mul entang sin mul add def} def

/decx2ay2a {/x2a x2a trialinc chord mul entang cos mul sub def /y2a
y2a trialinc chord mul entang sin mul sub def} def

/optimizex2y2 {mark findcurrenterror incx2ay2a findcurrenterror 2 copy
gt {{incx2ay2a findcurrenterror 2 copy lt{exit} if} loop decx2ay2a}
{ exch decx2ay2a {decx2ay2a findcurrenterror 2 copy lt{exit} if}
loop incx2ay2a } ifelse cleartomark} def

/extang 0 def % default

% this reports the currentspline to host for recording...

/makestring {10 string cvs print ( ) print} def

/reportsplinetohost {(\r\r) print flush x0a makestring y0a makestring
(moveto\r) print flush x1a makestring y1a makestring x2a makestring
y2a makestring x3a makestring y3a makestring (curveto\r\r) print flush
(The total least squares error is ) print flush findcurrenterror
makestring (.\r\r) print flush (The entry angle is ) print y1a y0a
sub x1a x0a sub dup 0 eq {pop 0.000001} if atan 100 mul round 100 div
dup 180 gt {360 sub} if makestring (degrees.\r\r\r) print flush
(The exit angle is ) print y3a y2a sub x3a x2a sub dup 0 eq
{pop 0.000001} if atan 100 mul round 100 div dup 180 gt {360 sub}
if makestring (degrees.\r\r\r) print flush} def

% reportsplinetohost % uncomment to use

% this combines the optimization process...

%%%%%%%%%(D) HOLD ENTRY ANGLE, ENTRY AND EXIT POINTS %%%%%%%%%%%%%%

% "A" spline constraints are as follows: Entry angle specified; entry
% and exit points are exact and fixed

/findbestAspline {/mat exch def initsplineguess initttguess optimizey2
optimizex2 optimizex1y1} def

/tryAagain {optimizey2 optimizex2 optimizex1y1} def

%%%%%%%%%%%% (D) HOLD EXIT ANGLE, ENTRY AND EXIT POINTS %%%%%%%%%%%%

% "B" spline constraints are as follows: Exit angle specified; entry
% and exit points are exact and fixed

/findbestBspline {/mat exch def initsplineguess initttguess optimizey1
optimizex1 optimizex2y2} def

/tryBagain {optimizey1 optimizex1 optimizex2y2} def

%%%%%%%%%%%% (E) HOLD ENTRY AND EXIT POINTS ONLY %%%%%%%%%%%%%%%%%%

% "C" spline constraints are as follows: Entry and exit points
% are exact and fixed; angles may vary

/findbestCspline {/mat exch def initsplineguess initttguess optimizey1
optimizex2 optimizex1 optimizey2} def

/tryCagain {optimizey1 optimizex2 optimizex1 optimizey2} def

% and finally a control point showing utility...

/showcontrolpoints{gsave x0a y0a moveto 1.5 dup scale dot grestore gsave
x1a y1a moveto 2 dup scale dot grestore gsave x2a y2a moveto 2.5 dup
scale dot grestore gsave x3a y3a moveto 3 dup scale dot grestore} def

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% And that's where we are for now. The code currently takes around
% five seconds per typical spline optimization on Distiller. The code
% can obviously be dramatically simplified and improved. But hey, it
% seems to work.

% And I have run roughshod over traditional math here. There's no
% guarantee that this will converge or will converge independently
% of how the points are adjusted.

% Obvious extensions are to design constraints for end points held
% only and for exit angle fixed.

% Let me know what you want to see further here.

%%%%%%%%%%%%% (F) SOME EXAMPLES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% (A) Our original problem...

% Run the Gonzo utilities (already on in this example, so commented)

% (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities

% gonzo begin % activate the gonzo utilities
% ps.util.1 begin
% % printerror
% nuisance begin

100 100 10 setgrid             % start a gonzo grid
20 20 showgrid                 % and draw part of it
line1                          % a nice line width

/entang 60 def                 % entry angle in degrees
/distinc .005 def              % distance increment as size fraction
/trialinc .005 def             % approximation increment as size fraction

[2 3 5 8 8 14 11 17 14 17 16 15 18 11] % the data set, 4 points minimum

dup plotdots % show the data points on the graph

findbestAspline               % find the best A constrained spline
showspline                    % draw it
reportsplinetohost            % report key values to host

grestore % get off grid

showpage % show results

% (B) Spline with inflection point and exit angle of 90 degrees...

% Run the Gonzo utilities (already on in this example, so commented)

% (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities

% gonzo begin     % activate the gonzo utilities
% ps.util.1 begin
% % printerror
% nuisance begin

100 100 10 setgrid        % start a gonzo grid
20 20 showgrid            % and draw part of it
line1                     % a nice line width

/entang 70 def            % entry angle in degrees
/distinc .05 def          % distance increment as size fraction
/trialinc .05 def         % approximation increment as size fraction

[2 10 5 12 8 11 11 8 14 6 17 5 19 10] % the data set, 4 points minimum

dup plotdots % show the data points on the graph

/distinc 0.1 def
/trialinc 0.1 def

findbestAspline % find the best A constrained spline
showspline

/distinc 0.05 def
/trialinc 0.05 def

tryAagain % make a second pass showspline

/distinc 0.005 def
/trialinc 0.005

tryAagain % and a third showspline
tryAagain % and a fourth

showspline % draw it

reportsplinetohost % report key values to host

grestore % get off grid

showpage % show results

% Note: If the optimization run seems to take forever, drop the
% distinc and trialinc values as shown here. In this case, three
% passes were needed to converge on a good solution. If the
% results are not close enough, this means it is unrealistic to
% try and fit the data points with a single spline. Try splitting
% the data in half and try again.

% Here is the error versus the number of passes for this curve...

% 1st trial ...... 2.80568
% 2nd trial ...... 1.92304
% 3rd trial ...... 0.88609
% 4th trial ...... 0.08672
% 5th trial ...... 0.08671

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% (C) Spline with a loop ...

% A cubic spline can have a cusp discontinuity or even a loop in it.
% The algorithm used here only can work towards a nearby minimum,
% which is not necessarily the best possible minimum.

% The code only works properly IF YOUR FIRST GUESS IS A GOOD ONE.
% If you want to do a loop, you have to use a good initial guess
% that is loop specific. Like so...

% Run the Gonzo utilities (already on in this example, so commented)

% (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities

% gonzo begin       % activate the gonzo utilities
% ps.util.1 begin
% % printerror
% nuisance begin

save /loopsnap exch def % save because of a custom guess

100 100 10 setgrid         % start a gonzo grid
20 20 showgrid             % and draw part of it
line1                      % a nice line width

/entang 30 def             % entry angle in degrees
/distinc .005 def          % distance increment as size fraction
/trialinc .005 def         % approximation increment as size fraction

% This special guess is data specific....

/initsplineguess{

/x0a 2 def /y0a 5 def /x1a 25 def /y1a 15 def /x2a -5 def /y2a 15
def /x3a 18 def /y3a 5 def

} def

initsplineguess

showcontrolpoints

showspline

[2 5 5 7 8 8 12 12 13 14 12 17 10 18 8 17
7 14 8 12 12 8 15 7 18 5] % the data set, 4 points minimum

dup plotdots % show the data points on the graph

/distinc 0.05 def /trialinc 0.05 def

findbestCspline % find the best C constrained spline

(still at it...) print flush

4 {tryCagain} repeat showspline

reportsplinetohost % report key values to host

grestore % get off grid

showpage % show results

loopsnap restore % get rid of special guess

% As you can see, it is far more accurate to split the loop up into at
% least a pair of splines. It is unrealistic to try and do all but the
% simplest of loops, especially if the input or exit angle must be held.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SUPPLEMENT A - INTERPOLATING ALONG A PARABOLA;
% =============================================

% Interpolating along a parabola to find a minimum distance to an
% off-curve point.

% This scheme finds the minimum of a parabola given three points near
% that minimum. It is exact if the parabola is exact with no higher order
% terms. Most very short pieces of most cubic splines can be accurately
% approximated by a perfect parabola.

% Let y1, y2, and y3 be three points of equal x spacing d along a
% parabola near its minimum. Let y2 be an x spacing k away from the
% origin.

% Let the parabola be of form y = ax^2 + b. Then

% take some differences...

% y2 - y1 = ak^2 + b - ak^2 +2adk -ad^2 - b = a (2dk - d^2)

% y3 - y2 = ak^2 + 2akd +2ad^2 +b -ak^2 - b = a (2dk + d^2)

% and divide them, calling D the ratio of the y differences...

% D = (y2 - y1)/(y3 - y2) = ad (2k - d) / ad (2k + d) = (2k - d)/(2k + d)

% solve for k, the offset between the minimum and the y1 value...

% 2kD + dD = 2k - d % 2kD - 2k = -dD - d % 2k(1 - D) = dD + d = d (D + 1)

% k = d (1 - D) / 2(1 + D)

% k is the distance from the x distance of y2 to the minimum value.

% Curiously, the math is independent of the a and b scale and offset
% values of the parabola itself. "If you've seen one parabola, you've
% seen them all."

% Assume that the stack holds (y3 - y2) on top of (y2 - y1) and that
% inca is the x spacing. Adjustment from the x position of y2 will be...

% /interpolate {neg exch div dup 1 exch sub exch 1 add
% dup 0 eq {0.00000001} if div inca mul 2 div} def

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright c. 1998 by Don Lancaster and Synergetics, Box 809, % Thatcher AZ, 85552 (520) 428-4073. synergetics@tinaja.com % All commercial rights and all electronic media rights are % *fully* reserved. Linking welcome. Reposting is expressly forbidden.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%EOF

Copyright c. 1998 by Don Lancaster and Synergetics, Box 809, Thatcher AZ, 85552
(520) 428-4073. synergetics@tinaja.com All commercial rights and all electronic media
rights are *fully* reserved. Linking welcome. Reposting is expressly forbidden.

Consulting services available.
Further support on www.tinaja.com/info01.html  

Please click here to...