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 review of the fundamentals
of
% cubic spline math, along with several new utilities that let
% you precisely find individual x and y values at any point
% along a cubic spline, curveto, or Bezier curve.
%
% Note: This is primary a tutorial from which you edit and extract
% PostScript routines for your own uses. There is no printed output.
% 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 modified code AS A STANDARD ASCII TEXTFILE with
% unique newnewfilename.
%
% (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.
%
% I had to COMPLETELY solve lots of points along a cubic spline
for an
% upcoming RETOUCH.PS universal gamma corrector and photo improver. So
% I thought I'd try to summarize the key spline math here.
% A cubic spline is two parametric equations in t (or time) space...
% x = At^3 + Bt^2 + Ct + D
% y = Dt^3 + Et^2 + Ft + G
%
% The t variable goes precisely from 0 to 1 as you move along the curve.
% This can generate a straight line, a curved line, a curve with a cusp
% in it, a curve with an inflection point (such as a sinewave) in it, or
% certain curves with simple loops in them.
% Cubing can be a pain, so the above equations can be rewritten
in a
% "cubeless" form that calculates quickly...
% x = ((At + B)t + C)t + D
% y = ((Dt + E)t + F)t + G
% In graph space, the cubic spline is defined by eight points.
A pair of
% initial points x0 and y0. A pair of end points x3 and y3. A pair of
% first influence points x1 and y1. A second pair of influence points
% x2 and y2.
%
% The first influence point sets the SLOPE and the ENTHUASIASM with
% which the initial point is exited. The second influence point sets
% the SLOPE and the ENTHUASIASM with which the final point is entered.
% Other names for the enthuasiasm are the TENSION or then VELOCITY.
%
% Simple algebra gets from the graph space to the equation space...
%
% A = x3 - 3x2 + 3x1 -x0 E = y3 - 3y2 + 3y1 -y0
% B = 3x2 - 6x1 + 3x0 F = 3y2 - 6y1 + 3y0
% C = 3x1 - 3x0 G = 3y1 - 3y0
% D = x0 H = y0
%
% You solve backwards to get from equation space to graph space...
%
% x0 = D y0 = H
% x1 = D + C/3 y1 = H + G/3
% x2 = D + 2C/3 + B/3 y2 = Y + 2G/3 + F/3
% x3 = D + C + B + A y3 = H + G + F + E
%
% The length of a cubic spline is not obvious, and I have yet to prove
a
% closed form expression exists. For now, you chop the spline up into
% small straight line approximations and add them up.
%
% Intermediate points along a spline can be found several ways.
% One route is to use FLATTENPATH to break up
% the curve into straight line segments. A second is CONSTANT T to
% break the curve up into equal changes in t. And a third crude method
% breaks the curve into CONSTANT S (or constant length) segments.
%
% I wasn't happy with the latter crude approximation. So I
% decided to look further. The key problem is to find which y goes with
% what x. In RETOUCH.PS, we need to relate x values to y values for one
% or more gamma correction and "photo improvement" tables.
%
% This is not trivial. An x can have one y value associated with it,
% three y values, or infinite y values. You get one solution
% if x is monotonic. You get three solutions if x reverses on itself for
% a sinewave or a loop. And you get infinite solutions if the curve is
% really a vertical line of infinite x slope.
%
% To further complicate matters, two of your three x values can end up
% "outside" of the t = 0 to t = 1 range that defines the curve.
% This is easy enough to test for, but still another detail to watch.
%
% We can usually ignore the extremely rare vertical straight line case.
% If you have got to test for it, it only occurs when x0 = x1 = x2 = x3.
% You end up with a slope of y = mX + b where m is infinite.
%
% I don't know how to find a closed solution to finding y given
% x. You first have to use x to solve for t and then you solve t for y.
%
% x = At^3 + Bt^2 + Ct + D
%
% ...is an example of a CUBIC equation.
It will
% have either ONE real root or else THREE real roots. There is no "cubic
% formula" comparable to the "quadratic formula". Let's
worry about
% finding ONE real root that always exists. This may be all you need.
% Otherwise, we can test and proceed from there.
%
% One useful way to do this is to take a guess for a t value. See what
x
% you get. Note the error. Reduce the error and try again. Keep this up
% till you have a root to acceptable accuracy.
%
% A good first guess is to normalize x so it ranges from 0 to 1 and then
% simply guess that x = t. This will be close for curves that aren't
% bent very much. And a useful guess for ALL spline curves.
%
% Let's use a better ploy to get our approximation to close quickly.
% It is called the NEWTON-RAPHELSON method, but is much simpler than it
% sounds. Say we get an error of x - x1. x1 is the current x for our
% current guess. At x1, our spline curve has a slope. Find the slope.
% The slope is expressed as rise/run.
%
% Now, on any triangle...
%
% rise = run x (rise/run)
%
% This gives us a very good improvement for our next approximation. It
% turns out that the "adjust for slope" method converges very
rapidly.
% Three passes are usually good enough.
%
% If our curve has an equation of...
%
% x
= At^3 + Bt^2 + Ct + D
%
% ...its
slope will be...
%
% x'
= 3At^2 +2Bt + C
%
% And the dt/dx slope will be its inverse or 1/(3At^2 + 2Bt +C)
%
% This is easily calculated. We'll have a code example in just a bit.
% % The next guess will be...
%
% nextguess = currentt + (curentx - x)(currentslope)
%
% In the case of our upcoming RETOUCH.PS gamma corrector, the spline
% is almost always monotonic and a single root is all we need.
%
% But what if we need to know for sure whether there are multiple y
% values for any given x? First note that the equation...
%
% At^3
+ Bt^2 + Ct + D = 0
%
% ...has at least one solution of (t - K) where K is the root
% we just found by successive approximation.
%
% Divide (t - K) into At^3 + Bt^2 + Ct + D to see what is left. It
% will be this quadratic in t...
%
%
% At^2
+ (B - KA)t + (KKA - BK +C)
% ___________________________________________________
% |
% (t+K) | At^3 + Bt^2
+ Ct +
D
% At^3 + KAt^2
% ==================================================
% (B
- KA)t^2 Ct D
% (B
- KA)t^2 (BK -KKA)t
% ==================================================
% (C
- BK + KAA)t D
% (C
- BK + KAA)t KKKA + BKK + CK
% ========================================
% 0
0
%
% Note that AK^3 + BK^2 + CK has to equal D if K is a valid root.
%
%
% Which gives us...
%
% (A)t^2 + (B - KA)t + (AK^2
- BK + C) = 0
%
% This will be of form...
%
%
ax^2 + bx + c = 0
%
% As with any quadratic, if b^2 - 4ac is positive or zero, there are
% TWO real roots, equal to the second and third root of the cubic.
% If b^2 - 4ac is negative, there will be a complex conjugate pair of
% imaginary roots. Meaning that you have one real cubic spline root.
%
% For RETOUCH.PS, monotonic x values are all I needed. The implicit
% assumption is that there is always only y value for a given x value.
% Let's find this root first. We can use this to develop fancier code
% that handles multiple roots in the t = 0 to t = 1 space..
% Let's develop some PostScript code utilities...
% grabcoeffs takes the eight spline xy values and changes them
into
% some intermediate variables handy to calculate x and y as a function
% of t...
/grabcoeffs {
/y3a exch def %
final y
/x3a exch def %
final x
/y2a exch def %
2nd influence y
/x2a exch def %
2nd influence x
/y1a exch def %
1st influence y
/x1a exch def %
1st influence x
/y0a exch def %
initial y
/x0a exch def %
initial x
/A x3a x1a x2a sub 3 mul add x0a sub store %
as in A(x)t^3
/B x2a x1a dup add sub x0a add 3 mul store %
as in B(x)t^2
/C x1a x0a sub 3 mul store %
as in C(x)t
/D y3a y1a y2a sub 3 mul add y0a sub store %
as in D(y)t^3
/E y2a y1a dup add sub y0a add 3 mul store %
as in E(y)t^2
/F y1a y0a sub 3 mul store %
as in F(y)t
} def
% And this is the "magic" code to find x and y given t from 0 to 1..
/ytt { dup dup D mul E add mul F add mul y0a add} def /xtt {
dup dup A
mul B add mul C add mul x0a add} def
% This code evaluates the dt/dx slope at a value of t...
/slopet { dup A mul 3 mul B 2 mul add mul C add dup 0 eq {0.000001}
if 1
exch div} def
% Here is some code to find the first real root of the cubic spline...
% /tforx solves for one possible value of t for the selected x...
/tforx {/curx exch def x3a x0a sub curx x0a sub exch div /tguess
exch def
4 {tguess dup dup xtt curx sub exch slopet mul sub /tguess exch store}
repeat tguess} def
% Let's try it. We'll first use the spline plotted in HACK62.PS
2 3 3 7.5 7 7 8 4 grabcoeffs
% Find the t value at x = 3, ferinstance...
3 tforx
% Plotting the intermediate guesses, our first t guess gives
an x value
% of 2.72222. The first adjustment gives 3.01453. The second adjusment
% gives 3.00003. The third and higher adjustments are integer 3 beyond
% PostScript precision. The series very rapidly converges. At least for
% this example.
% Find the corresponding y for the t left on the stack by tforx...
ytt
% Which in this case turns out to be 5.23226.
% If you are doing this on your own, you can print the stack
values back
% to the host, or else install a printing error trapper and purposely
% place a "zowie" in your file anytime you want a stack dump.
% At this point in our demo, we have a y value on the stack.
Pop it
% to continue...
pop
% Single roots are all I needed for TOUCHUP.PS. To get fancy,
/find3roots
% returns all real roots as an array.
% NOTE THAT UNNEEDED REAL ROOTS OFTEN OCCUR OUTSIDE OF THE
% t = 0 TO t = 1 RANGE.
% Dividing the found root t = K out of the cubic leaves this
quadratic...
%
% (A)t^2 + (B - KA)t
+ (AK^2 - BK + C) = 0
%
% Where K in this case equals our best tguess.
/find3roots { mark /a A def /b B tguess neg A mul sub def /c
tguess dup
mul A mul tguess neg B mul sub C add def tguess b dup mul a c mul 4 mul
sub dup 0 ge {sqrt a 2 mul div b neg a 2 mul div exch 2 copy add 3 1 roll
sub } if ] } def
find3roots
% An array of length 1 or length 3 now appears on the stack top.
It shows
% t values, some of which may be outside of the allowable 0 to 1 range.
% ////////////// THE FINAL RESULTS //////////////////////
% /find1y returns the y value for a given x AFTER the spline
has been set
% Use this fast code for MONOTONIC functions without loops or cusps...
/find1y {tforx ytt} def
% demo...
2 3 3 7.5 7 7 8 4 grabcoeffs % input the spline
3 find1y % find y for x = 3
% This returns the correct value of 5.23226. Any other y value
is
% easily found.
% pop the answer off the stack
pop
% /find3y returns all possible y values for a given x AFTER the
% spline has been set. Use this slower code if you suspect loops and/or
% cusps.
/find3y {tforx pop find3roots mark exch {dup dup 1 gt exch 0
lt or not
{ytt}{pop} ifelse} forall ] } def
% In this loop demo, we ask for y values for 17 different x values..
2 1 11 5 -1 5 8 1 grabcoeffs
% spline having a loop 2 0.5 8.01 {dup find3y} for
% grab x,y data values
% the following x point and y value arrays are returned to the stack...
% 2.0 [1.0]
% 2.5 [1.2281]
% 3.0 [1.46983]
% 3.5 [1.67319]
% 4.0 [1.58829]
% 4.5 [3.78729 3.31007 2.3312]
% 5.0 [4.0 2.71429 2.71429]
% 5.5 [3.78729 2.3312 3.31007]
% 6.0 [1.58829]
% 6.5 [1.67319]
% 7.0 [1.46983]
% 7.5 [1.2281]
% 8.0 [1.0]
% In the "loopless" region, you have only a single
y value. And, for
% that matter only a single real root, with a complex pair for the other
% two. In the "loop" region, there are 3 possible y values for
every x.
% Finding which of the 3 y values to actually use is left as
an exercise
% for the student. Generally, you will switch from one y value to three
% as you cross the edge of a loop or cusp. Your IMMEDIATE PREVIOUS value
% will be a single value and will be much closer to one of the three.
% Ferinstance, in the above example the most likely y value in
the
% [3.78729 3.31007 2.3312] array to follow a y value of [1.58829]
% would be 2.3312. Closer data points would make this more obvious.
% Naturally, if you were just interested in drawing the curve,
you would
% use -curveto- instead. The methods shown here are needed when you want
% to extract numeric y values from the spline for one reason or another.
% See TOUCHUP.PS for a stunning example of what all this is good for...
% Have fun...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Click here for the Downloadable
PostScript Sourcecode 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. | Find out what a tinaja is. | |
Find research solutions. | View recommended books. | |
Visit the acrobat library. | Return to home page. | |