Re: Quick And Simple Bezier Curve Drawing or other
- Posted by begin Jul 02, 2018
- 1160 views
here another spline i use for timelines.
namespace QSPLINE ifdef BITS32 then constant FAIL = #3FFFFFFF elsedef constant FAIL = #3FFFFFFFFFFFFFFF end ifdef sequence SetUpSplineDerived = {} atom LOWEST_X = FAIL --Routine returns the table of 2nd derivatives needed for the cubic spline --interpolation routine spline(). x[*] and y[*] are the input vectors. --It is assumed that the x values are ordered but not necessarily equally spaced. --That is, xin[0] < xin[1] < xin[2] < ... xin[n-1]. public function SetUpSpline(sequence x, sequence y) atom tmp, dx, dx0, dx1, h integer n = length(x) if n != length(y) then SetUpSplineDerived = {} return false end if if n <= 1 then SetUpSplineDerived = {} return false end if integer im1, ip1, n1 = n - 1 sequence g = repeat(0.0, n) sequence temp = repeat(0.0, n1) temp = repeat(0.0, n1) SetUpSplineDerived = {} g[1] = 0.0 temp[1] = 0.0 for i = 2 to n1 do ip1 = i+1 im1 = i-1 dx = x[i] - x[im1] dx0 = x[ip1] - x[i] dx1 = x[ip1] - x[im1] if dx <= 0.0 or dx0 <= 0.0 or dx1 <= 0.0 then SetUpSplineDerived = {} return false end if tmp = dx/dx1 h = 2.0 + tmp*g[im1] g[i] = (tmp - 1)/h temp[i] = (y[ip1]-y[i])/dx0 - (y[i] -y[im1])/dx temp[i] = (temp[i]*6.0/dx1 - tmp*temp[im1])/h end for g[n1] = 0.0 for j = n-1 to 1 by -1 do g[j] = temp[j] + g[j] * g[j+1] end for LOWEST_X = x[1] SetUpSplineDerived = g return true end function --This routine returns the cubic spline interpolation at the point x0. --g[*] is the table of 2nd derivatives needed for the cubic spline formula --(see the spline0() function). --x[*] and y[*] are the input vectors of x-y data. --All of the vectors are of length n. --It is assumed that the x values are ordered but not necessarily equally spaced. --That is, xin[0] < xin[1] < xin[2] < ... xin[n-1]. public function Spliner(atom x0, sequence x, sequence y) if SetUpSplineDerived={} then return FAIL end if atom tmp, dx, xi, xi1, xi1d, xid, y0 integer i, top, bot, n = length(x) bot = 0 top = n while (top-bot) > 1 do -- find location in table through bisection i = shift_bits((top+bot), 1) if x[i] <= x0 then bot = i else top = i end if end while dx = x[top]-x[bot] xi1 = x[top]-x0 xi = x0-x[bot] xi1d = xi1/dx xid = xi/dx tmp = (SetUpSplineDerived[bot]*(xi1d*xi1-dx)*xi1 + SetUpSplineDerived[top]*(xid*xi-dx)*xi) / 6.0 y0 = tmp + y[bot]*xi1d + y[top]*xid return y0 end function public function spline_from_to(sequence x, sequence y, atom from_x1, atom to_x2, atom step_by=1) if SetUpSplineDerived={} then return {} end if if FAIL >= from_x1 and to_x2 <= from_x1 then return {} end if if step_by <= 0 then return {} end if sequence output={} atom counter = from_x1 while counter <= to_x2 do output &= Spliner(counter,x,y) -- interpolate with spline counter += step_by end while return output end function /* procedure splinetest() integer i, k, n=7, nout=19 atom delta, delta3, truth atom x0, y0, twopi_ = PI*2 sequence x, y, deriv y = repeat(0.0, n) x = repeat(0.0, n) delta = twopi_/(n-1) delta3 = delta/3 -- define unequally spaced x values x[1]=delta*0 x[2]=delta/3 x[3]=delta x[4]=delta*2.33333333 x[5]=delta*3 x[6]=delta*4.33333333 x[7]=delta*6 -- generate one cycle of a sinusoid; with nonequally spaced samples for i = 1 to n do y[i]=sin(x[i]) end for printf(1," interpolated original truth:\n") printf(1," x y k y[k] sin(x)\n\n") k=1 bool worked = SetUpSpline(x, y) -- find 2nd derivative of data for i = 0 to nout-1 do truth = sin(twopi_*i/(nout-1)) x0 = i*delta3 y0 = Spliner(x0, x, y) -- interpolate with spline if abs(x0-x[k]) < 10e-6 then printf(1,"%9.6f ", x0) printf(1,"%9.6f ", y0) printf(1,"%2d ", k) printf(1,"%9.6f ", y[k]) printf(1,"%9.6f\n", truth) k += 1 else printf(1,"%9.6f ", x0) printf(1,"%9.6f ", y0) printf(1,"%9.6f\n", truth) end if end for end procedure splinetest() maybe_any_key() */ /* sequence x={1,10,20,30,40,50} sequence y={0, 100, 50, 100, 50, 10} bool worked = SetUpSpline(x, y) ?spline_from_to(x, y, 1, 50, 1) maybe_any_key() */