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.
%
% The GENERAL PROCEDURE for using PostScript-as-language is as
% follows:
% (1) Carefully READ the entire sourcecode file inside your
% favorite WP or Editor as an ASCII textfile.
%
% (2) Carefully MODIFY the filenames and variables and such
% so they apply to you.
%
% (3) SAVE the modified code AS A STANDARD ASCII TEXTFILE
% using a unique filename.
%
% (4) RUN the modified code by dragging and dropping it into
% Acrobat Distiller. GhostScript may alternately
be used.
%
% (5) TEST, INTERPRET, and MODIFY your results.
%
%%%%%%%%%%%% (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...
% FIRST influence point at TWICE the distance to the second
% data point ALONG THE ENTANG AXIS.
%
% SECOND influence point at TWICE the distance from the next
% to the last data point IN THE DIRECTION OF THAT DATA POINT.
% 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
% y1 = a(k-d)^2 + b
% y2 = ak^2 + b
% y3 = a(k+d)^2 + b
% 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
Click here for the Downloadable
PostScript Sourcecode for this file.
Click here for the Acrobat
PDF Graphic Demo Figure for this file.
Remember to first READ the file in an editor, MODIFY it
to match your
needs, SAVE it as an ORDINARY ASCII TEXTFILE, and
then RUN
it by dragging and dropping into Acrobat Distiller.
Consulting services available.
Further support on www.tinaja.com/info01.html
Please click here to...
Get a Synergetics catalog. | Send Don Lancaster email. | |
Start your own tech venture. | Pick up surplus bargains. | |
Sponsor a display banner. |