More on Cubic Spline Math
By Don Lancaster                                                                   Version 1.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 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:

% 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...

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

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...