% ************************************************************************
% ************************************************************************
%
% MORE ON CUBIC SPLINE MATH
%
% ************************************************************************
%
% 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.
%
% Also has a cross reference of other PSRT cubic spline files.
%
% Copyright c 1993 by Don Lancaster. All rights reserved.
% Free help line and additional info: (602) 428-4073.
%
% ************************************************************************
% Name of textfile: BEZMATH.PS
% Source: SYNERGETICS
% Author: Don Lancaster
% Desc: More on cubic spline math
% Date: August 25, 1993
% Release: 1.0
% Status: Copyright c 1993 by Don Lancaster and Synergetics.
% 3860 West First Street, Thatcher, AZ. (602) 428-4073.
% All commercial rights reserved. Personal use permitted
% so long as this header remains present and intact.
% PostScript Secrets Book + Disk costs $39.50. VISA/MC.
% Approx length: 20K
%
% Keywords: PostScript, guru, spline, cubic, Bezier, root,
% approximation, path, cubic, length, equation
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X E % activate XON/XOFF if needed
% Z % values shown for Apple IIe.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: This is primary a tutorial from which you edit and extract individual
% PostScript routines for your own uses. There is no printed output.
% MORE ON CUBIC SPLINE MATH
% =========================
% by Don Lancaster
% I needed to COMPLETELY solve lots of points along a cubic spline for our
% upcoming RETOUCH.PS universal gamma corrector and photo improver. So
% I thought I'd try to summarize the key spline math here.
% One mathematician did send me an apparantly elegant and simple closed
% form method to find the length of a spline. I don't quite understand it
% and it is going out to others for an independent check. On the surface
% it looks brilliant. Let me know if you want to review this ahead of time.
% A cubic spline consists of 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. See the cross reference at
% the end for a spline library and other stuff.
% Cubing can be a real 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. And a pair of second influence points
% x2 and y2.
%
% The first influence point determines the SLOPE and the ENTHUASIASM with
% which the initial point is exited. The second influence point determines
% the SLOPE and the ENTHUASIASM with which the final point is entered. The
% other names for the enthuasiasm are the TENSION or then VELOCITY.
%
% Simple algebra lets you get 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 can solve these 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 us. See #190 BEZLNGTH.PS.
%
% Intermediate points along a spline can be found several ways. These were
% shown in #662 BEZCHORD.PS. One route is to use FLATTENPATH to break up
% the curve into straight line segments. A second is to use CONSTANT T to
% break the curve up into equal changes in t. And a third crude method
% was shown to break the curve into CONSTANT S (or constant length) segments.
%
% I wasn't happy with the latter since it was a 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. A given x can have one y value associated with it,
% three y values, or an infinite number of 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 an exact and 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 first worry about
% finding the 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 fairly close for curves that aren't
% bent very much. And a useful guess for ALL spline curves.
%
% Let's use a much 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 will have 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 code and an 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, corresponding 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 root to your cubic spline.
%
% 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 new 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 to 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 looping demo, we ask for the 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, you have three possible y values for every x.
% Figuring out 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 previous 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...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% RELATED GENIE DOWNLOADS:
%
%
% #670 - SERPPATH.PS Compiled SERPPATH.GPS runtime demo
% #669 - SERPPATH.GPS Glyphs along a serpentine path (!)
% #662 - BEZCHORD.PS Subdividing a cubic spline
% #637 - MIDPOINT.PS Bezier curve recursive midpoint rule.
% #632 - BEZSINE.PS Sine wave with eight Bezier curves.
% #630 - BEZIER7.PS Connect two Bezier curves.
% #626 - SPLINCAT.PS Prints full cubic spline catalog.
% #605 - BEZDIST.PS Distance from point to Bezier curve
% #618 - BEZ4PTS.PS Cubic spline thru four data points.
% #616 - HACK62.PS Cubic spline fundamentals, more.
% #601 - BEZIERU.PS A better way to do a Bezier curve.
% #590 - BEZIER.PS Bezier curve fast evaluation.
% #588 - FUZZYBEZ.PS Bezier curve from fuzzy data (!)
% #202 - TWIXTBEZ.FNT Flagwaver ffont & Bezier ct Utility
% #190 - BEZLNGTH.PS Bezier curve length utility
% #161 - SLEEZOID.PS Avuncular Sleezoid Utility
% #079 - MEOWWRRR.PS PostScript puss de resistance
%
% #537 - IMAGESIZ.TXT PS image file size reduction tricks
% #535 - RESBAROM.PS PostScript resolution barometer
% #533 - CSCRNFIX.TXT Fixing level 2 currentscreen blowups
% #528 - FUSONFIX.TXT SX/CX fusion roller maintenance tips
% #511 - NUTS9.PS PostScript for Hardware Hackers
% #500 - OPTICOMP.PS Compiled PS optimization study
% #488 - SXGAMMAC.PS Canon SX Gamma correction (!)
% #487 - LASGRES.TXT Desc: PhotoGrade resolution: How high?
% #469 - HACK55.TXT Digital halftones & PhotoGrade
% #459 - SPEEDUP.PS PostScript Speedup Secrets (!)
% #451 - LASGCAL.PS LaserWriter F/G calibration secrets
% #435 - TONERTRX.PS Terrific Toner Techniques (!)
% #389 - RTCTONER.TXT Cheryle White RECHARGER RTC 911108
% #388 - LASGNOTE.TXT LaserWriter G prliminary snoop
% #345 - RTCTHOMP.TXT PS laser printer repair RTC 911004
% #337 - SECRTEMP.PS PS Startup Secrets for PC users
% #264 - PGRAYSIT.MAC Unlimited gray scales in PageMaker
% #239 - CUSTMDOT.TXT Total PS Halftone Screen Control
% #231 - GRAFGRAY.PS Gobs of Great Guru Graphing Grays
% #180 - SCRNTEST.PS ScreenTest greyscale utility
% #179 - SCRNDOC.PS Docs for SCREENTEST utility
% #175 - DUPLXTST.PS Duplex paper/toner evaluator
% #144 - SPOTDOTS.GPS 300 DPI Spotfunction Patterns
% #141 - HMNYGRAY.PS How many PostScript grays?
% #128 - PINSCRTS.TXT PS Insider Secrets
% #67 - BLACKFLS.PS Black flasher
% #57 - RUBBERGR.PS Fine gray rubbergrid
%
% #517 - GONZO15.PTL Gonzo utility powertools upgrade
% #220 - GONZO13A.TXT Documentation for the Gonzo PTL
%
% #478-482 PSRT directory, cross reference & planner
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Write, call or GEnie [SYNERGETICS] email for your free PostScript
% Insider's Secrets and Hardware Hacker Insider Info brochures.
% Contact Don Lancaster's SYNERGETICS for reprint availability.
% Full consulting services available based on the concepts shown above.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FREE VOICE HELPLINE AND ADDITIONAL INFO: (602) 428-4073
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%