% ************************************************************************ % ************************************************************************ % % 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.PSL % 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: (928) 428-4073 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%