%!
%
% ************************************************************************
% ************************************************************************
%
% BEZIER CURVE THROUGH FUZZY DATA
%
% ************************************************************************
% 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
%
% 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 over and over again when designing typogrophy
% or doing engineering charts and graphs.
% Even trickier, you can't just tie a bunch of cubic splines together. More
% often than not, you want to add CONSTRAINTS. For instance, adjacent splines
% should have matching endpoints. For best appearance, the entry and exit
% angles should match. And ideally, the rate of change of curvature should
% also match at the joints.
% One VERY handy tool: given four or more noisy data points that are 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 from your previous spline.
% This key tool would solve the general problem of taking a noisy
% character shape or an engineering graph, and drawing a smooth and efficient
% group-of-splines PostScript curve through everything.
% One simple older tool set involved the CURVETRACE routines found in my
% GONZO20.PS and in the older gonzo utilities. See our classic MEOWWRRR.PS
% for a use example. The curvetracing routines take a series of data values
% of x0-y0-ang0, x1-y1-ang1, up through ... xn-yn-angn. They then draw 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 curvetracing
% routines are intolerant of fuzzy or noisy data. Although they form a very
% sparse data set, curvetracing also 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 space
% 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 and slow tools that may
% let you do much of 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 and faster 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.
% We'll use my gonzo utilities to simplify our study here. Lets start off
% by doing a persistent download of GONZO15.PTL and then activating the
% guru tools by saying, of all things...
% Now let's pick a bunch of slightly noisy data points, carefully 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 at first...
/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 % put two errors on stack; what we
incx1ay1a % have and what we could get
findcurrenterror
2 copy gt % save these errors If inc is high
{ {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 required 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