%!PS
% An improved Bezier Cubic Spline length finder and LINEAR(!) subdivider
% =======================================================================
% by Don Lancaster
% Copyright c 2006 by Don Lancaster & Synergetics, Box 809, Thatcher, AZ, 85552
% (928) 428-4073 Email: don@tinaja.com Website: http://www.tinaja.com
% Consulting services available http://www.tinaja.com/info01.html
% All commercial rights and all electronic media rights ~fully~ reserved.
% Linking usually welcome. Reposting expressly forbidden. Version 1.1
% IMPORTANT NOTE: Don Lancaster's file gonzo.ps is recommended but not needed
% for this program. This may be downloaded from http://www.tinaja.com/post01.asp#gonzo
% To activate after obvious location mods, uncomment ONE of the following two lines:
% (C:\\Documents and Settings\\don\\Desktop\\gonzo\\gonzo.ps) run % use internal gonzo
% (A:\\gonzo.ps) run % use external gonzo
% Plus the next line if Gonzo is to be activated...
% gonzo begin ps.util.1 begin nuisance begin
% NOTE THAT ALL PS FILENAME STRINGS !!!DEMAND!!! DOUBLE REVERSE SLASHES.
% A collection of utilities that let you find the length of a full or partial cubic
% spline using actual along-the-curve measurement for best possible accuracy.
% A "S-list" method provides accurate and constant spacing. Spacing is uniform,
% compared to "constant-T" methods which are closer together along the "more bent"
% portions of the cuve.
% Methods presented are accurate and conceptually simple but computationally intensive.
% Run modified text as PostScript ASCII textfile sent to Acrobat Distiller.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% TEMPORARY COPIES OF GONZO ROUTINES FOR DEMOS BELOW %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of "real" gonzo is strongly encouraged.
% Free download from http://www.tinaja.com/post01.asp#gonzo
/xrpt{gsave aload pop /trips exch def /dist exch def /rproc exch def
trips { gsave rproc grestore dist 0 translate } repeat grestore} def
/yrpt{gsave aload pop /trips exch def /dist exch def /rproc exch def
trips { gsave rproc grestore 0 dist translate } repeat grestore} def
/fat5 true def % 3 pixels wide every fifth line?
/fatter10 true def % 5 pixels wide every tenth line?
/fatterborder true def % 5 pixel gray outline border?
/setgrid { gsave /blocksize exch def translate
blocksize dup scale} def
/thingridlines 0 store
/ggray {60 0 {sub abs 0.21 lt {1}{0} ifelse} setscreen
gsave 0 60 0 {pop pop 1 add 1} setscreen grestore dup 25 eq
{pop 0.8}{100 eq {0.72}{0.765} ifelse} ifelse setgray
0.3 1 0.5 setrgbcolor % modded - not sure why
} def
/showgrid {gsave ggray /vblocks exch def /hblocks exch def
thingridlines setlinewidth
[{0 0 moveto 0 vblocks rlineto stroke} 1
hblocks 1 add] xrpt
[{0 0 moveto hblocks 0 rlineto stroke} 1
vblocks 1 add] yrpt
fatterborder { gsave newpath 0 0.96 blocksize div
dtransform round idtransform
setlinewidth pop 2 setlinecap
0 0 moveto hblocks 0 rlineto 0 vblocks rlineto
hblocks neg 0 rlineto closepath stroke grestore} if
fat5 { gsave newpath 0 0.48 blocksize div dtransform round
idtransform setlinewidth pop mark {5 0 moveto 0 vblocks rlineto
stroke} 5 hblocks 5 div cvi] xrpt mark {0 5 moveto hblocks 0 rlineto
stroke} 5 vblocks 5 div cvi] yrpt grestore} if
fatter10 { gsave newpath 0 0.96 blocksize div
dtransform round idtransform setlinewidth pop
mark {10 0 moveto 0 vblocks rlineto stroke} 10
hblocks 10 div cvi] xrpt mark {0 10 moveto hblocks 0 rlineto stroke} 10
vblocks 10 div cvi] yrpt grestore} if
grestore} def
/line1 {.06 dup setlinewidth 5 mul /erase exch def} def
/line2 {.12 dup setlinewidth 5 mul /erase exch def} def
/line3 {.18 dup setlinewidth 5 mul /erase exch def} def
/m {moveto} store
/dot { currentpoint newpath 0.150 0 360 arc fill } def
/mdot { m dot} def
/mergestr {2 copy length exch length add string dup dup 4 3 roll
4 index length exch putinterval 3 1 roll exch 0 exch putinterval} def
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% CUBIC SPLINE SERVICE UTILITIES %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% code to convert x0y0...x3y3 influence points to A-H equation constants.
% See http://www.tinaja.com/glib/bezmath.pdf for a derivation...
/findAH {
/A x3 x2 3 mul sub x1 3 mul add x0 sub store
/E y3 y2 3 mul sub y1 3 mul add y0 sub store
/B x2 3 mul x1 6 mul sub x0 3 mul add store
/F y2 3 mul y1 6 mul sub y0 3 mul add store
/C x1 3 mul x0 3 mul sub store
/G y1 3 mul y0 3 mul sub store
/D x0 store
/H y0 store
} store
% code to find current x and y given A-H and t and using "cubless" method
% input tt, exit with curx or cury defined but not on stack
/findx {/ttt exch store A ttt mul B add ttt mul C add ttt mul D add
/curx exch store} store
/findy {/ttt exch store E ttt mul F add ttt mul G add ttt mul H add
/cury exch store} store
% code to find curent s incremental chord length given curx cury prevx prevy...
% the usual sum-of squares of the x and y increment. Returns to curs but not on stack.
/finds {curx prevx sub dup mul cury prevy sub dup mul add sqrt
/curs exch store} store
% key to the routines are three arrays of identical length...
% The tlist is an array of Bezier t values, initially constant spaced.
% The slist is the RUNNING TOTAL of the chord lengths, ending with length.
% The nslist is a NORMALIZED 0-1 version of the slist.
% code to make a t list of numchords+1 samples. In the case of 100 samples,
% would be [0 0.01 0.02 .... 1.00]
/numchords 100 store % recommended default
/makecurtlist {/curtlist mark 0 numchords {dup 1 numchords div add} repeat ] store
} store
% code to make a slist from a tlist. Each entry is SUMMED distance to the tlist value;
% final value is Bezier length.
/makecurslist {
/prevx 0 store % initialize previous values
/prevy 0 store
/curslist mark 0 % start slist array
0 1 numchords {curtlist exch get /tt exch store % grab current
tt findx
tt findy % find curx and cury
finds % current chord
curs add dup % integrate running result and dup for next
/prevx curx store % save new chord end
/prevy cury store
} for
pop % undo final integration!
] store % save slist
} store
% /findsplinelen reads the last entry in the slist, equal to the total spline length.
/findsplinelen {curslist dup length 1 sub get /splinelen exch store} store
% /findtofs finds an interpolated t value for a given s distance. The fractional
% distance through a s array entry is calculated and converted to dt/ds. Once t
% is known, the corresponding x and y positions are easily found.
% findbezlen must be previously run.
/findtofs { /ss exch store % must be in absolute range of 0 to length.
1 1 curslist length 1 sub % scan each s value
{ /ii exch store
ss curslist ii get % view each accumulated s in turn till too big
le {exit} if } for % then exit with high s location
/ii ii 1 sub store % back up one so s is in box
curslist ii 1 add get
curslist ii get sub % current s delta
ss curslist ii get sub % excess into interval
exch div % current interpolated s fraction
curtlist ii 1 add get % t delta
curtlist ii get sub
mul % fractional t interpolation distance
curtlist ii get add % final t for selected s
} store
% Yeah, the above "count till overflow" is somewhat crude, but it seems more
% than fast enough. You can speed it up by a binary pretree (0, 0.25, 0.5 0.75 etc)
% Or else by "inverting" fot a linear curslist and nonlinear curtlist.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% HIGH LEVEL BEZIER LENGTH & SUBDIVISION PROCS %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% /findbezlen calculates the total Bezier spline length
% and returns it as /splinelen. Requires predefined x0y0...x3y3
% and predefined numchords. (100 recommended).
/findbezlen {findAH % calculate cubic coefficients
makecurtlist % create linear t array
makecurslist % create nonlinear s array
findsplinelen % report length as splinelen
} store
% /findxyofs finds the x and y position for a given sublength position.
% The sublength position is presently entered NORMALIZED in the 0 to 1.0
% range. findbezlen must be run previously.
/findxyofs {splinelen mul % denormalize
findtofs % find the interpolated t for the LINEAR s
dup findx % get x(t)
findy % get y(t)
curx cury % and return to stack
} store
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% DEMOS - REMOVE OR ALTER BEFORE REUSE %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%
%%%% (A) Find the length of a Bezier cubic spline %%%%%%
%%%%%%%%%%%%%
% Spline coefficients or defaults need entered here...
/testspline {
/x0 0 store /y0 0 store
/x1 3 store /y1 20 store
/x2 16 store /y2 4 store
/x3 20 store /y3 20 store
} store
testspline % enter values
/numchords 100 store % new entry here if not default=100
findbezlen % find length of spline
% and report to log file...
(\n\nThe length of the test spline is ) splinelen 20 string cvs mergestr
(.\n) mergestr print flush
%%%%%%%%%%%%%%
% Accuracy check of a spline of known length of 28.28427...
%%%%%%%%%%%%%%
/x0 0 store /y0 0 store
/x1 0.001 store /y1 0.001 store
/x2 19.999 store /y2 19.999 store
/x3 20 store /y3 20 store
findbezlen % find length of spline
(\nThe length of the 28.28427 known length spline is ) splinelen 20 string cvs mergestr
(.\n) mergestr print flush
%%%%%%%%%%%%%%%
% Accuracy check for increasing numchords. Values >> 10,000 may trip PS math limits.
%%%%%%%%%%%%%%%
testspline % enter values
/numchords 10 store % new entry here if not default=100
findbezlen % find length of spline
% and report to log file...
(\nThe length of the 10 sample test spline is ) splinelen 20 string cvs mergestr
(.\n) mergestr print flush
testspline % enter values
/numchords 100 store % new entry here if not default=100
findbezlen % find length of spline
% and report to log file...
(The length of the 100 sample test spline is ) splinelen 20 string cvs mergestr
(.\n) mergestr print flush
testspline % enter values
/numchords 1000 store % new entry here if not default=100
findbezlen % find length of spline
% and report to log file...
(The length of the 1000 sample test spline is ) splinelen 20 string cvs mergestr
(.\n) mergestr print flush
testspline % enter values
/numchords 10000 store % new entry here if not default=100
findbezlen % find length of spline
% and report to log file...
(The length of the 10000 sample test spline is ) splinelen 20 string cvs mergestr
(.\n\n) mergestr print flush
%%%%%%%%%%%%%%%
% Find x and y values a given distance into a spline
%%%%%%%%%%%%%%%
testspline % enter values
/numchords 100 store % new entry here if not default=100
findbezlen % find length of spline
0.271 findxyofs
(\n\nGoing 0.271 of the distance into testspline results in an x value of )
curx 20 string cvs mergestr ( and a y value of ) mergestr cury 20 string cvs
mergestr (.\n) mergestr print
% Plotted in red in next example below.
%%%%%%%%%%%%%%%
% Subdivide a cubic spline into 20 EQUAL chord length parts
%%%%%%%%%%%%%%%
% requires full or partial gonzo utilities
100 100 10 setgrid % draw a pretty graph
20 20 showgrid
testspline % enter values
/numchords 100 store % new entry here if not default=100
findbezlen % find length of spline
x0 y0 moveto
x1 y1 x2 y2 x3 y3 curveto
line1 stroke
0 0.05 1.001 { findxyofs m dot} for
1 0 0 setrgbcolor
2.67118 7.9798 m dot % plot previous example result
showpage
% eof