%! % % ************************************************************************ % ************************************************************************ % % BEZIER CURVE THROUGH FUZZY DATA % % ************************************************************************ % 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 % % SUMMARY: Guru Don Lancaster gives us a tutorial and heavy duty set of % utilities that allow you to convert noisy data points into % smooth cubic splines. This is a preliminary release. % % This is especially useful in capturing fonts from bitmaps, % plotting engineering curves, or smoothing scanned artwork. % % Yes, the key constraint problems of preserving the input % or angle and staying accurate on the end points are specifically % handled. % % 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 the modified code AS A STANDARD ASCII TEXTFILE % using a unique filename. % % (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. % %%%%%%%%%%%% (A) INSTALL AND RUN GONZO UTILITIES %%%%%%%%%%%%%%%%%% % The demo uses my gonzo utilities from % http://www.tinaja.com/psutils/gonzo20.ps (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities gonzo begin % activate the gonzo utilities ps.util.1 begin % printerror nuisance begin %%%%%%%%%%%% (B) INTRO AND THEORY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % How do you take a bunch of noisy dots and run a smooth curve through % them? This problem crops up over and over again when designing typogrophy % or doing engineering charts and graphs. % Even trickier, you can't just tie a bunch of cubic splines together. More % often than not, you want to add CONSTRAINTS. For instance, adjacent splines % should have matching endpoints. For best appearance, the entry and exit % angles should match. And ideally, the rate of change of curvature should % also match at the joints. % One VERY handy tool: given four or more noisy data points that are within % the range of what a single cubic spline can handle, find a very good % Bezier curve that precisely goes through the end points, while smoothly % coming reasonably near all of the others. AND starts off at a precisely % selected angle, letting you smoothly continue from your previous spline. % This key tool would solve the general problem of taking a noisy % character shape or an engineering graph, and drawing a smooth and efficient % group-of-splines PostScript curve through everything. % One simple older tool set involved the CURVETRACE routines found in my % GONZO20.PS and in the older gonzo utilities. See our classic MEOWWRRR.PS % for a use example. The curvetracing routines take a series of data values % of x0-y0-ang0, x1-y1-ang1, up through ... xn-yn-angn. They then draw some % "fairly weak" splines through the selected points, using one spline per % data triad. Slopes and endpoints all match. % While very useful for cartoons and digitized signatures, my curvetracing % routines are intolerant of fuzzy or noisy data. Although they form a very % sparse data set, curvetracing also uses lots of weak splines, instead of % a few good ones. Angles also have to be specified for every point. % Now, the CORRECT way to handle the problem appears in CURVEFITTING WITH % PIECEWISE PARAMETIC CUBICS by Michael Plass and Maureen Stone (with a % little help by some guy named Warnock). This appeared in the July 1983 % issue of COMPUTER GRAPHICS on pages 229-239, and otherwise known as the % 1983 SIGGRAPH PROCEEDINGS. % The only little problem I had with this paper was doing the vector space % Newton-Raphson iterations of the Hermite basis functions. Seriously, it % was painfully obvious that bunches of time and effort would be required % to master the postgraduate math involved. Especially while translating % everything to convenient PostScript end user procs. % So, what I thought we would do here is some dumb and slow tools that may % let you do much of the same stuff in a more intuitive fashion. What follows % is open ended, intended as both tools you can use now to smooth fuzzy % data, and as intermediate steps to fancier and faster code that uses more % elegant math. % Naturally, once you smooth your data once, you steal the plans and just % use the final splines for all repeat copies of the same job. % We will attack three different constraint problems: Type A holds % the entry angle and the end points. Type B holds the exit angle and % the end points. Type C holds the end points only. % We'll use my gonzo utilities to simplify our study here. Lets start off % by doing a persistent download of GONZO15.PTL and then activating the % guru tools by saying, of all things... % Now let's pick a bunch of slightly noisy data points, carefully arranged % as [x0 y0 x1 y1 ... xn yn]. We want at least four points, but can have % as many as you need. The assumption is made that one spline can handle % this portion of the curve we are trying to fit. One spline can produce % a smoothly flowing curve that has at most one cusp, one loop, or two % inflection points. For brevity, we'll use whole numbers here at first... /rawdata [2 3 5 8 8 14 11 17 14 17 16 15 18 11] def % Next, let's draw a fancy gray grid... 100 100 10 setgrid 20 20 showgrid % do a simple Gonzo dot-plotter... /dot { currentpoint newpath 0.150 0 360 arc fill } def /plotdots {/mat exch def 0 2 mat length 1 sub {/posn exch def mat posn get mat posn 1 add get mt dot} for} def rawdata plotdots % and optionally look at them... % copypage % Let's assume our previous spline left us with a required entrance angle % of sixty degrees... /entang 60 def % The key to any successive approximation scheme is to start out with % something reasonable as our first guess. As a wild guess, and % remembering our constraints, let's try a rule of... % FIRST influence point at TWICE the distance to the second % data point ALONG THE ENTANG AXIS. % % SECOND influence point at TWICE the distance from the next % to the last data point IN THE DIRECTION OF THAT DATA POINT. % We'll keep things simple and obvious here for a while, rather than % going for speed or compact code... /mat {rawdata} def /delta1 {mat 2 get mat 0 get sub dup mul mat 3 get mat 1 get sub dup mul add sqrt} def /initsplineguess{ /x0a {mat 0 get} def /y0a {mat 1 get} def /x1a {delta1 entang cos mul 2 mul x0a add } def /y1a {delta1 entang sin mul 2 mul y0a add} def /x2a {mat dup length 4 sub get dup mat dup length 2 sub get sub 2 mul add} def /y2a {mat dup length 3 sub get dup mat dup length 1 sub get sub 2 mul add} def /x3a {mat dup length 2 sub get} def /y3a {mat dup length 1 sub get} def } def initsplineguess % Let's see what our first cut looks like... line1 /showspline {x0a y0a moveto x1a y1a x2a y2a x3a y3a curveto stroke} def showspline grestore % get off grid for now showpage % show initial guess on first page % Which ain't half bad for amateurs. %%%%%%%%%%%%%%% (C) POINT OPTIMIZATION %%%%%%%%%%%%%%%%%%% % Now, instead of ultra fancy math, % let's use some slow and intuitive successive approximations to get % us started understanding what's coming down. % The game plan from here: We find how close the curve comes to each one % of the intermediate points. Then we find the total sum of the errors % squared for the points. Then we adjust x1, y1, x2, and y2 incrementally % to minimize this error. In this case, x1 and Y1 are locked together by % entang. With some luck, things will rapidly converge for resonable % types of curves that are likely to occur. % But first a review: There are two forms for a Bezier curve. These are % the graph space and the equation space. % In graph space, a Bezier curve starts at x0,y0 and ends at x3 y3. % It leaves x0,y0 with a direction and an enthuasiasm set by the first % control point x1,y1 and enters x3,y3 with a direction and an % enthuasiasm set by the second control point x2,y2. % Here are the same Bezier formulas in the algebraic equation space: % % x(t) = A(t^3) + B(t^2) + C(t) + x0 % y(t) = D(t^3) + E(t^2) + F(t) + y0 % % Note that (t) varies from 0 to 1 from the % beginning to the end of your spline curve. % You can think of a Bezier spline as a three dimensional spiral in % x, y, and z space. As you go from the initial values of x and y % to the final values of x and y, t goes exactly from 0 to 1. It is % sometimes useful to think of t as TIME, starting at t=0 and ending % with t=1. % As the red book tells us, here is how to get from equation space % to graph space: % % x1 = x0 + C/3 y1 = y0 + F/3 % x2 = x1 + C/3 + B/3 y2 = y1 + F/3 + E/3 % x3 = x0 + C + B + A y3 = y0 + F + E + D % You can use simple algebra to get the far more useful inverse form. % Here is how to get from graph space back to equation space .. % % A = x3 - 3x2 + 3x1 + x0 D = y3 - 3y2 + 3 y1 + y0 % B = 3x2 - 6x1 + 3x0 E = 3y2 - 6y1 + 3y0 % C = 3x1 - 3x0 F = 3y1 - 3y0 % This grabs six graph values off the stack... /grabgraf {/y3 exch def /x3 exch def /y2 exch def /x2 exch def /y1 exch def /x1 exch def /y0 exch def /x0 exch def} def % Having predefined all four x and y values this finds x(t) given tt in % the range of 0 to 1 ... /xtt {x3 x2 3 mul sub x1 3 mul add x0 sub tt 3 exp mul x2 3 mul x1 6 mul neg add x0 3 mul add tt dup mul mul add x1 3 mul x0 3 mul neg add tt mul add x0 add} def % this finds y(t) given tt in the range of 0 to 1 ... /ytt {y3 y2 3 mul sub y1 3 mul add y0 sub tt 3 exp mul y2 3 mul y1 6 mul neg add y0 3 mul add tt dup mul mul add y1 3 mul y0 3 mul neg add tt mul add y0 add} def % this finds the magnitude of the error distance between a given tt and % points x, y on the stack /finderrdist {ytt sub exch xtt sub dup mul exch dup mul add sqrt} def % In general, tt does NOT travel uniformly from 0 to 1. It tends to % advance more violently on the "more bent" portions of the curve. % Let's guess a tt value for each data point, initially assuming an % equal tt spacing... /initttguess {/ttguess [ 0 mat length 2 div dup 1 exch div /inc exch def 2 sub cvi {dup inc add} repeat 1 ] def} def initttguess % Now, let's improve each tt value by incrementing it in distinc steps to % find the minimum distance. What we are doing here is sliding along the % t space until we are as close to each noisy point as we can get. % We will try half percent tt steps initially. This seems good enough for % now... /distinc .005 def % We first try a direction to see if the error gets larger or smaller. % Then we continue making the error smaller till we overshoot. After % an overshoot, we use the "parabola fit interpolation" of Appendix A % to put a better guess into the ttguess array. /improvettguess { /tt curtt def /curinc distinc def mat curpos 2 mul get /curx exch def mat curpos 2 mul 1 add get /cury exch def mark curx cury finderrdist /tt tt curinc add def curx cury finderrdist exch dup 2 index sub 3 1 roll 2 index 0 gt {pop}{exch pop exch neg exch /curinc curinc neg def /tt tt curinc add def} ifelse /tt tt curinc add def {curx cury finderrdist dup 2 index exch sub 3 -1 roll pop exch 1 index 0 gt {/tt tt curinc add def}{pop exit} ifelse} loop curinc 0 lt {exch /curinc curinc neg def /tt tt curinc 2 mul add def} if interpolate tt add distinc sub dup /tt exch def ttguess exch curpos exch put cleartomark} def % Here's the interpolation used (see Appendix A) ... /interpolate {neg exch div dup 1 exch sub exch 1 add dup 0 eq {0.00000001} if div distinc mul 2 div} def % this takes care of optomizing all the inside values, given the eight % graphspace Bezier points... /setnearest { grabgraf 1 1 ttguess length 2 sub {/curpos exch def ttguess curpos get /curtt exch def improvettguess} for} def x0a y0a x1a y1a x2a y2a x3a y3a setnearest % At this point, the ttguess array has values that accurate approximate % the nearest point on the x(t), y(t) curve to the points to be fitted % finderrors measures the error distances for all the intermediate points % and returns an array of those actual distances /finderrors {mark 1 1 ttguess length 2 sub {/posn exch def /tt ttguess posn get def mat posn 2 mul get mat posn 2 mul 1 add get finderrdist } for] /currenterrors exch def} def finderrors % Executing finderrors give us a value of .146415 for the first point % first pass, which "looks about right" and converges to the .146321 % actual distance found at higher resolution. Other values of .569515, % .846671, .789305, and .442002 also "look about right" on the graph. % To find the total sum-of-squares error, each error value is squared % and added... /findsoserror { 0 currenterrors {dup mul add} forall /currenterror exch def } def % findsoserror % Our current sum-of-squares error is 1.86033. % The key tools are now in place. We can find the error for any given % Bezier curve fit to our noisy data. Now the fun begins... % In general, there are eight variables in the Bezier graph space for % any given spline. Fixing the end points fixes four of these variables. % fixing the entang entry angle locks x1 and y1 to each other. We % thus have three variables left that we can adjust to try and get the % best fit we can. % The trial increment is the normalized percentage change we are going % to make for each test improvement... /trialinc .005 def % use half percent increment steps % This finds the length of the end-to-end vector. It will be handy to % normalize the incremental test steps... /chord { mat 0 get mat dup length 2 sub get sub dup mul mat 1 get mat dup length 1 sub get sub dup mul add sqrt} def % This combines the previous utilities into a single error finder /findcurrenterror {x0a y0a x1a y1a x2a y2a x3a y3a setnearest finderrors findsoserror currenterror} def % findcurrenterror % uncomment to activate % these three routines are used to optimize x2a. x2a is incremented once % to determine a direction. It is then incremented in the correct % direction until a minimum is found... /incx1a {/x1a x1a trialinc chord mul add def} def /decx1a {/x1a x1a trialinc chord mul sub def} def /optimizex1 {mark findcurrenterror incx1a findcurrenterror 2 copy gt {{incx1a findcurrenterror 2 copy lt{exit} if} loop decx1a } { exch decx1a {decx1a findcurrenterror 2 copy lt{exit} if} loop incx1a} ifelse cleartomark} def % optimizex1 % uncomment to activate % these three routines are used to optimize x2a. x2a is incremented once % to determine a direction. It is then incremented in the correct % direction until a minimum is found... /incx2a {/x2a x2a trialinc chord mul add def} def /decx2a {/x2a x2a trialinc chord mul sub def} def /optimizex2 {mark findcurrenterror incx2a findcurrenterror 2 copy gt {{incx2a findcurrenterror 2 copy lt{exit} if} loop decx2a } { exch decx2a {decx2a findcurrenterror 2 copy lt{exit} if} loop incx2a} ifelse cleartomark} def % optimizex2 % uncomment to activate % these three routines are used to optimize y1a. Y1a is incremented once % to determine a direction. It is then incremented in the correct % direction until a minimum is found... /incy1a {/y1a y1a trialinc chord mul add def} def /decy1a {/y1a y1a trialinc chord mul sub def} def /optimizey1 {mark findcurrenterror incy1a findcurrenterror 2 copy gt {{incy1a findcurrenterror 2 copy lt{exit} if} loop decy1a } { exch decy1a {decy1a findcurrenterror 2 copy lt{exit} if} loop incy1a} ifelse cleartomark} def % optimizey1 % uncomment to activate % these three routines are used to optimize y2a. Y2a is incremented once % to determine a direction. It is then incremented in the correct % direction until a minimum is found... /incy2a {/y2a y2a trialinc chord mul add def} def /decy2a {/y2a y2a trialinc chord mul sub def} def /optimizey2 {mark findcurrenterror incy2a findcurrenterror 2 copy gt {{incy2a findcurrenterror 2 copy lt{exit} if} loop decy2a } { exch decy2a {decy2a findcurrenterror 2 copy lt{exit} if} loop incy2a} ifelse cleartomark} def % optimizey2 % uncomment to activate % these three routines are used to optimize x1a and y1a, WHEN THEY ARE % CONSTRAINED TO LIE ALONG THE ENTANG AXIS. They get incremented once % to determine a direction. They are then incremented in the correct % direction until a minimum is found... /incx1ay1a {/x1a x1a trialinc chord mul entang cos mul add def /y1a y1a trialinc chord mul entang sin mul add def} def /decx1ay1a {/x1a x1a trialinc chord mul entang cos mul sub def /y1a y1a trialinc chord mul entang sin mul sub def} def /optimizex1y1 {mark findcurrenterror % put two errors on stack; what we incx1ay1a % have and what we could get findcurrenterror 2 copy gt % save these errors If inc is high { {incx1ay1a findcurrenterror 2 copy lt{exit} if} loop decx1ay1a } { exch decx1ay1a {decx1ay1a findcurrenterror 2 copy lt{exit} if} loop incx1ay1a } ifelse cleartomark } def % these three routines are used to optimize x2a and y2a, WHEN THEY ARE % CONSTRAINED TO LIE ALONG THE EXTANG AXIS. They get incremented once % to determine a direction. They are then incremented in the correct % direction until a minimum is found... /incx2ay2a {/x2a x2a trialinc chord mul entang cos mul add def /y2a y2a trialinc chord mul entang sin mul add def} def /decx2ay2a {/x2a x2a trialinc chord mul entang cos mul sub def /y2a y2a trialinc chord mul entang sin mul sub def} def /optimizex2y2 {mark findcurrenterror incx2ay2a findcurrenterror 2 copy gt {{incx2ay2a findcurrenterror 2 copy lt{exit} if} loop decx2ay2a} { exch decx2ay2a {decx2ay2a findcurrenterror 2 copy lt{exit} if} loop incx2ay2a } ifelse cleartomark} def /extang 0 def % default % this reports the currentspline to host for recording... /makestring {10 string cvs print ( ) print} def /reportsplinetohost {(\r\r) print flush x0a makestring y0a makestring (moveto\r) print flush x1a makestring y1a makestring x2a makestring y2a makestring x3a makestring y3a makestring (curveto\r\r) print flush (The total least squares error is ) print flush findcurrenterror makestring (.\r\r) print flush (The entry angle is ) print y1a y0a sub x1a x0a sub dup 0 eq {pop 0.000001} if atan 100 mul round 100 div dup 180 gt {360 sub} if makestring (degrees.\r\r\r) print flush (The exit angle is ) print y3a y2a sub x3a x2a sub dup 0 eq {pop 0.000001} if atan 100 mul round 100 div dup 180 gt {360 sub} if makestring (degrees.\r\r\r) print flush} def % reportsplinetohost % uncomment to use % this combines the optimization process... %%%%%%%%%(D) HOLD ENTRY ANGLE, ENTRY AND EXIT POINTS %%%%%%%%%%%%%% % "A" spline constraints are as follows: Entry angle specified; entry % and exit points are exact and fixed /findbestAspline {/mat exch def initsplineguess initttguess optimizey2 optimizex2 optimizex1y1} def /tryAagain {optimizey2 optimizex2 optimizex1y1} def %%%%%%%%%%%% (D) HOLD EXIT ANGLE, ENTRY AND EXIT POINTS %%%%%%%%%%%% % "B" spline constraints are as follows: Exit angle specified; entry % and exit points are exact and fixed /findbestBspline {/mat exch def initsplineguess initttguess optimizey1 optimizex1 optimizex2y2} def /tryBagain {optimizey1 optimizex1 optimizex2y2} def %%%%%%%%%%%% (E) HOLD ENTRY AND EXIT POINTS ONLY %%%%%%%%%%%%%%%%%% % "C" spline constraints are as follows: Entry and exit points % are exact and fixed; angles may vary /findbestCspline {/mat exch def initsplineguess initttguess optimizey1 optimizex2 optimizex1 optimizey2} def /tryCagain {optimizey1 optimizex2 optimizex1 optimizey2} def % and finally a control point showing utility... /showcontrolpoints{gsave x0a y0a moveto 1.5 dup scale dot grestore gsave x1a y1a moveto 2 dup scale dot grestore gsave x2a y2a moveto 2.5 dup scale dot grestore gsave x3a y3a moveto 3 dup scale dot grestore} def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % And that's where we are for now. The code currently takes around % five seconds per typical spline optimization on Distiller. The code % can obviously be dramatically simplified and improved. But hey, it % seems to work. % And I have run roughshod over traditional math here. There's no % guarantee that this will converge or will converge independently % of how the points are adjusted. % Obvious extensions are to design constraints for end points held % only and for exit angle fixed. % Let me know what you want to see further here. %%%%%%%%%%%%% (F) SOME EXAMPLES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (A) Our original problem... % Run the Gonzo utilities (already on in this example, so commented) % (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities % gonzo begin % activate the gonzo utilities % ps.util.1 begin % % printerror % nuisance begin 100 100 10 setgrid % start a gonzo grid 20 20 showgrid % and draw part of it line1 % a nice line width /entang 60 def % entry angle in degrees /distinc .005 def % distance increment as size fraction /trialinc .005 def % approximation increment as size fraction [2 3 5 8 8 14 11 17 14 17 16 15 18 11] % the data set, 4 points minimum dup plotdots % show the data points on the graph findbestAspline % find the best A constrained spline showspline % draw it reportsplinetohost % report key values to host grestore % get off grid showpage % show results % (B) Spline with inflection point and required exit angle of 90 degrees... % Run the Gonzo utilities (already on in this example, so commented) % (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities % gonzo begin % activate the gonzo utilities % ps.util.1 begin % % printerror % nuisance begin 100 100 10 setgrid % start a gonzo grid 20 20 showgrid % and draw part of it line1 % a nice line width /entang 70 def % entry angle in degrees /distinc .05 def % distance increment as size fraction /trialinc .05 def % approximation increment as size fraction [2 10 5 12 8 11 11 8 14 6 17 5 19 10] % the data set, 4 points minimum dup plotdots % show the data points on the graph /distinc 0.1 def /trialinc 0.1 def findbestAspline % find the best A constrained spline showspline /distinc 0.05 def /trialinc 0.05 def tryAagain % make a second pass showspline /distinc 0.005 def /trialinc 0.005 tryAagain % and a third showspline tryAagain % and a fourth showspline % draw it reportsplinetohost % report key values to host grestore % get off grid showpage % show results % Note: If the optimization run seems to take forever, drop the % distinc and trialinc values as shown here. In this case, three % passes were needed to converge on a good solution. If the % results are not close enough, this means it is unrealistic to % try and fit the data points with a single spline. Try splitting % the data in half and try again. % Here is the error versus the number of passes for this curve... % 1st trial ...... 2.80568 % 2nd trial ...... 1.92304 % 3rd trial ...... 0.88609 % 4th trial ...... 0.08672 % 5th trial ...... 0.08671 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (C) Spline with a loop ... % A cubic spline can have a cusp discontinuity or even a loop in it. % The algorithm used here only can work towards a nearby minimum, % which is not necessarily the best possible minimum. % The code only works properly IF YOUR FIRST GUESS IS A GOOD ONE. % If you want to do a loop, you have to use a good initial guess % that is loop specific. Like so... % Run the Gonzo utilities (already on in this example, so commented) % (C:\\windows\\desktop\\gonzo\\gonzo.ps) run % run the gonzo utilities % gonzo begin % activate the gonzo utilities % ps.util.1 begin % % printerror % nuisance begin save /loopsnap exch def % save because of a custom guess 100 100 10 setgrid % start a gonzo grid 20 20 showgrid % and draw part of it line1 % a nice line width /entang 30 def % entry angle in degrees /distinc .005 def % distance increment as size fraction /trialinc .005 def % approximation increment as size fraction % This special guess is data specific.... /initsplineguess{ /x0a 2 def /y0a 5 def /x1a 25 def /y1a 15 def /x2a -5 def /y2a 15 def /x3a 18 def /y3a 5 def } def initsplineguess showcontrolpoints showspline [2 5 5 7 8 8 12 12 13 14 12 17 10 18 8 17 7 14 8 12 12 8 15 7 18 5] % the data set, 4 points minimum dup plotdots % show the data points on the graph /distinc 0.05 def /trialinc 0.05 def findbestCspline % find the best C constrained spline (still at it...) print flush 4 {tryCagain} repeat showspline reportsplinetohost % report key values to host grestore % get off grid showpage % show results loopsnap restore % get rid of special guess % As you can see, it is far more accurate to split the loop up into at % least a pair of splines. It is unrealistic to try and do all but the % simplest of loops, especially if the input or exit angle must be held. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SUPPLEMENT A - INTERPOLATING ALONG A PARABOLA; % ============================================= % Interpolating along a parabola to find a minimum distance to an % off-curve point. % This scheme finds the minimum of a parabola given three points near % that minimum. It is exact if the parabola is exact with no higher order % terms. Most very short pieces of most cubic splines can be accurately % approximated by a perfect parabola. % Let y1, y2, and y3 be three points of equal x spacing d along a parabola % near its minimum. Let y2 be an x spacing k away from the origin. % Let the parabola be of form y = ax^2 + b. Then % y1 = a(k-d)^2 + b % y2 = ak^2 + b % y3 = a(k+d)^2 + b % take some differences... % y2 - y1 = ak^2 + b - ak^2 +2adk -ad^2 - b = a (2dk - d^2) % y3 - y2 = ak^2 + 2akd +2ad^2 +b -ak^2 - b = a (2dk + d^2) % and divide them, calling D the ratio of the y differences... % D = (y2 - y1)/(y3 - y2) = ad (2k - d) / ad (2k + d) = (2k - d)/(2k + d) % solve for k, the offset between the minimum and the y1 value... % 2kD + dD = 2k - d % 2kD - 2k = -dD - d % 2k(1 - D) = dD + d = d (D + 1) % k = d (1 - D) / 2(1 + D) % k is the distance from the x distance of y2 to the minimum value. % Curiously, the math is independent of the a and b scale and offset % values of the parabola itself. "If you've seen one parabola, you've % seen them all." % Assume that the stack holds (y3 - y2) on top of (y2 - y1) and that % inca is the x spacing. Adjustment from the x position of y2 will be... % /interpolate {neg exch div dup 1 exch sub exch 1 add % dup 0 eq {0.00000001} if div inca mul 2 div} def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%EOF