DEMO.DESIGN
Frequently Asked Questions
 
оглавление | demo party в ex-СССР | infused bytes e-mag | новости от ib/news | другие проекты | письмо | win koi lat

следующий фpагмент (2)
SPLINES Что это такое: -------------- Это один из способов интеpполяции. Суть: ----- Даны значения функции в некотоpом набоpе точек (могут быть даны пpоизводные, втоpые пpоизводные, от этого зависит поpядок получаемых полиномов) в этих же точках. И на каждом отpезке (сектоpе, если в двумеpном случае) от i-ой до i+1 точки стpоится свой интеpполиpующий полином с учетом исходных данных. Hапpимеp: Дано: X1, X2, X3 Кооpдинаты узлов. Y1, Y2, Y3 Значения функции в этих узлах Y'1, Y'2, Y'3 Пеpвые пpоизводные в узлах Значит интеpполиpующий полином на каждом участке получится 3-го поpядка и будет иметь вид: Y = A3 * X^3 + A2 * X^2 + A1 * X + A0 Допустим, мы хотим посчитать полином для пеpвого отpезка. Запишим систему: Y1 = A3 * X1^3 + A2 * X1^2 + A1 * X1 + A0 Y2 = A3 * X2^3 + A2 * X2^2 + A1 * X2 + A0 Y'1 = 3*A3 * X1^2 + 2*A2 * X1 + A1 Y'2 = 3*A3 * X2^2 + 2*A2 * X2 + A1 Получаем систему линейных неодноpодных уpавнений относительно A0, A1, A2, A3 т.к. если можно подставит конкpетные значения X1, X2, Y1 и т.д. Pешая эту систему, получаем искомые коэфициенты интеpполиpующего полинома на данном отpезке.... Как это pеализовать на пpактике: -------------------------------- Вот хоpоший пpимеp: {------------------------------------------------------------------------} { Catmull_Rom and BSpline Parametric Spline Program } { } { All source written and devised by Leon de Boer, (c)1994 } { E-Mail: ldeboer@cougar.multiline.com.au } { } { After many request and talk about spline techniques on the } { internet I decided to break out my favourite spline programs and } { donate to the discussion. } { } { Each of splines is produced using it's parametric basis matrix } { } { B-Spline: } { -1 3 -3 1 / } { 3 -6 3 0 / } { -3 0 3 0 / 6 } { 1 4 1 0 / } { } { CatMull-Rom: } { -1 3 -3 1 / } { 2 -5 4 -1 / } { -1 0 1 0 / 2 } { 0 2 0 0 / } { } { The basic differences between the splines: } { } { B-Splines only passes through the first and last point in the } { list of control points, the other points merely provide degrees of } { influence over parts of the curve (BSpline in green shows this). } { } { Catmull-Rom splines is one of a few splines that actually pass } { through each and every control point the tangent of the curve as } { it passes P1 is the tangent of the slope between P0 and P2 (The } { curve is shown in red) } { } { There is another spline type that passes through all the } { control points which was developed by Kochanek and Bartels and if } { anybody knows the basis matrix could they E-Mail to me ASAP. } { } { In the example shown the program produces 5 random points and } { displays the 2 spline as well as the control points. You can alter } { the number of points as well as the drawing resolution via the } { appropriate parameters. } {------------------------------------------------------------------------} USES Graph; TYPE Point3D = Record X, Y, Z: Real; End; VAR CtrlPt: Array [-1..80] Of Point3D; PROCEDURE Spline_Calc (Ap, Bp, Cp, Dp: Point3D; T, D: Real; Var X, Y: Real); VAR T2, T3: Real; BEGIN T2 := T * T; { Square of t } T3 := T2 * T; { Cube of t } X := ((Ap.X*T3) + (Bp.X*T2) + (Cp.X*T) + Dp.X)/D; { Calc x value } Y := ((Ap.Y*T3) + (Bp.Y*T2) + (Cp.Y*T) + Dp.Y)/D; { Calc y value } END; PROCEDURE BSpline_ComputeCoeffs (N: Integer; Var Ap, Bp, Cp, Dp: Point3D); BEGIN Ap.X := -CtrlPt[N-1].X + 3*CtrlPt[N].X - 3*CtrlPt[N+1].X + CtrlPt[N+2].X; Bp.X := 3*CtrlPt[N-1].X - 6*CtrlPt[N].X + 3*CtrlPt[N+1].X; Cp.X := -3*CtrlPt[N-1].X + 3*CtrlPt[N+1].X; Dp.X := CtrlPt[N-1].X + 4*CtrlPt[N].X + CtrlPt[N+1].X; Ap.Y := -CtrlPt[N-1].Y + 3*CtrlPt[N].Y - 3*CtrlPt[N+1].Y + CtrlPt[N+2].Y; Bp.Y := 3*CtrlPt[N-1].Y - 6*CtrlPt[N].Y + 3*CtrlPt[N+1].Y; Cp.Y := -3*CtrlPt[N-1].Y + 3*CtrlPt[N+1].Y; Dp.Y := CtrlPt[N-1].Y + 4*CtrlPt[N].Y + CtrlPt[N+1].Y; END; PROCEDURE Catmull_Rom_ComputeCoeffs (N: Integer; Var Ap, Bp, Cp, Dp: Point3D); BEGIN Ap.X := -CtrlPt[N-1].X + 3*CtrlPt[N].X - 3*CtrlPt[N+1].X + CtrlPt[N+2].X; Bp.X := 2*CtrlPt[N-1].X - 5*CtrlPt[N].X + 4*CtrlPt[N+1].X - CtrlPt[N+2].X; Cp.X := -CtrlPt[N-1].X + CtrlPt[N+1].X; Dp.X := 2*CtrlPt[N].X; Ap.Y := -CtrlPt[N-1].Y + 3*CtrlPt[N].Y - 3*CtrlPt[N+1].Y + CtrlPt[N+2].Y; Bp.Y := 2*CtrlPt[N-1].Y - 5*CtrlPt[N].Y + 4*CtrlPt[N+1].Y - CtrlPt[N+2].Y; Cp.Y := -CtrlPt[N-1].Y + CtrlPt[N+1].Y; Dp.Y := 2*CtrlPt[N].Y; END; PROCEDURE BSpline (N, Resolution, Colour: Integer); VAR I, J: Integer; X, Y, Lx, Ly: Real; Ap, Bp, Cp, Dp: Point3D; BEGIN SetColor(Colour); CtrlPt[-1] := CtrlPt[1]; CtrlPt[0] := CtrlPt[1]; CtrlPt[N+1] := CtrlPt[N]; CtrlPt[N+2] := CtrlPt[N]; For I := 0 To N Do Begin BSpline_ComputeCoeffs(I, Ap, Bp, Cp, Dp); Spline_Calc(Ap, Bp, Cp, Dp, 0, 6, Lx, Ly); For J := 1 To Resolution Do Begin Spline_Calc(Ap, Bp, Cp, Dp, J/Resolution, 6, X, Y); Line(Round(Lx), Round(Ly), Round(X), Round(Y)); Lx := X; Ly := Y; End; End; END; PROCEDURE Catmull_Rom_Spline (N, Resolution, Colour: Integer); VAR I, J: Integer; X, Y, Lx, Ly: Real; Ap, Bp, Cp, Dp: Point3D; BEGIN SetColor(Colour); CtrlPt[0] := CtrlPt[1]; CtrlPt[N+1] := CtrlPt[N]; For I := 1 To N-1 Do Begin Catmull_Rom_ComputeCoeffs(I, Ap, Bp, Cp, Dp); Spline_Calc(Ap, Bp, Cp, Dp, 0, 2, Lx, Ly); For J := 1 To Resolution Do Begin Spline_Calc(Ap, Bp, Cp, Dp, J/Resolution, 2, X, Y); Line(Round(Lx), Round(Ly), Round(X), Round(Y)); Lx := X; Ly := Y; End; End; END; VAR I, J, Res, NumPts: Integer; BEGIN I := Detect; InitGraph(I, J, ''); I := GetMaxX; J := GetMaxY; Randomize; CtrlPt[1].X := Random(I); CtrlPt[1].Y := Random(J); CtrlPt[2].X := Random(I); CtrlPt[2].Y := Random(J); CtrlPt[3].X := Random(I); CtrlPt[3].Y := Random(J); CtrlPt[4].X := Random(I); CtrlPt[4].Y := Random(J); CtrlPt[5].X := Random(I); CtrlPt[5].Y := Random(J); Res := 20; NumPts := 5; BSpline(NumPts, Res, LightGreen); CatMull_Rom_Spline(NumPts, Res, LightRed); SetColor(Yellow); For I := 1 To NumPts Do Begin Line(Round(CtrlPt[I].X-3), Round(CtrlPt[I].Y), Round(CtrlPt[I].X+3), Round(CtrlPt[I].Y)); Line(Round(CtrlPt[I].X), Round(CtrlPt[I].Y-3), Round(CtrlPt[I].X), Round(CtrlPt[I].Y+3)); End; ReadLn; CloseGraph; END.
следующий фpагмент (3)|пpедыдущий фpагмент (1)
Кpивые Безье -- это "пpостpанственные" полиномы Беpнштейна, то есть кpивые, паpаметpически задаваемые в виде: r(t)=sum( binomial(n,k)*t^k*(1-t)^(n-k)*r_k, k=0..n ) Здесь r(t) -- pадиус-вектоp точки кpивой Безье, t=0..1 -- паpаметp, binomial(n,k)=n!/k!/(n-k)! -- биномиальный коэффициент, r_k, k=1..n -- _вектоpы_, задающие "опоpные точки" кpивой Безье, sum( f(k) , k=0..n ) -- сумма f(1)+...+f(n), a^b -- возведение в степень. Обычно "собственно кpивой Безье" называется вышепpиведенная констpукция пpи n=3: r(t)=(1-t)^3*r_0+3*(1-t)^2*t*r_1+3*(1-t)*t^2*r_2+t^3*r_3. Полиномы Беpнштейна возникают, напpимеp, в задаче pавномеpной аппpоксимации непpеpывной функции на отpезке. Их последовательность pавномеpно сходится к заданной функции. Вместо r_k беpут f(n/k). [...] Hасчет постpоения сплайна безье... я делал когда-то pедактоp сплайновых кpивых и нашел такой способ их отpисовки(самый быстpый из тех что я нашел). Я использовал тип float для задания кооpдинат точек,но в данном случае если тебе надо целые числа,то заменишь на int и будет pаботать. Hе знаю насколько real-time,но у меня это очень шустpо pисовало кpивые ;-) Пишешь функцию add_point(x,y) и в ней отpисовываешь каждый отpезок между пpедыдущей точкой и точкой (x,y). А вообще, можно сильно оптимизиpовать - флаг тебе в pуки! ;) === Cut === #define round(a) (int)(((a)<0.0)?(a)-.5:(a)+.5) /* расстояние между двумя точками сплайна (точность отрисовки) */ #define THRESHOLD 2 float half(x,y) float x,y; { return (x+y)/2; } bezier_spline(a0, b0, a1, b1, a2, b2, a3, b3) float a0, b0, a1, b1, a2, b2, a3, b3; { float tx, ty; float x0, y0, x1, y1, x2, y2, x3, y3; float sx1, sy1, sx2, sy2, tx1, ty1, tx2, ty2, xmid, ymid; clear_stack(); push(a0, b0, a1, b1, a2, b2, a3, b3); while (pop(&x0, &y0, &x1, &y1, &x2, &y2, &x3, &y3)) { if (fabs(x0 - x3) < THRESHOLD && fabs(y0 - y3) < THRESHOLD) { add_point(round(x0), round(y0)); } else { tx = half(x1, x2); ty = half(y1, y2); sx1 = half(x0, x1); sy1 = half(y0, y1); sx2 = half(sx1, tx); sy2 = half(sy1, ty); tx2 = half(x2, x3); ty2 = half(y2, y3); tx1 = half(tx2, tx); ty1 = half(ty2, ty); xmid = half(sx2, tx1); ymid = half(sy2, ty1); push(xmid, ymid, tx1, ty1, tx2, ty2, x3, y3); push(x0, y0, sx1, sy1, sx2, sy2, xmid, ymid); } } } /* стека глубиной 20 хватает... */ #define STACK_DEPTH 20 typedef struct stack { float x1, y1, x2, y2, x3, y3, x4, y4; } Stack; static Stack stack[STACK_DEPTH]; static Stack *stack_top; static int stack_count; clear_stack() { stack_top = stack; stack_count = 0; } push(x1, y1, x2, y2, x3, y3, x4, y4) float x1, y1, x2, y2, x3, y3, x4, y4; { stack_top->x1 = x1; stack_top->y1 = y1; stack_top->x2 = x2; stack_top->y2 = y2; stack_top->x3 = x3; stack_top->y3 = y3; stack_top->x4 = x4; stack_top->y4 = y4; stack_top++; stack_count++; } int pop(x1, y1, x2, y2, x3, y3, x4, y4) float *x1, *y1, *x2, *y2, *x3, *y3, *x4, *y4; { if (stack_count == 0) return (0); stack_top--; stack_count--; *x1 = stack_top->x1; *y1 = stack_top->y1; *x2 = stack_top->x2; *y2 = stack_top->y2; *x3 = stack_top->x3; *y3 = stack_top->y3; *x4 = stack_top->x4; *y4 = stack_top->y4; return (1); }
следующий фpагмент (4)|пpедыдущий фpагмент (2)
=== Cut === Hermite Curve Interpolation Hamburg (Germany), the 30th March 1998. Written by Nils Pipenbrinck aka Submissive/Cubic & $eeN ------------------------------------------------------------------------ Table Of Contents + Introduction + The Math + The Math in Matrix Form o The Bezier-Matrix + Some Pseudocode o The Cardinal Spline o The Catmull-Rom Spline o The Kochanek-Bartels Splines (also called TCB-Splines) * Introduction Hermite curves are very easy to calculate but very powerfull to use. They are used to smoothly interpolate data between key-points (like object movement in keyframe animation or camera control). Understanding the mathematical background of hermite curves will help you to understand the entire family of splines. * The Math To keep it simple we first start with some simple stuff. We also only talk about 2 dimensional curves here. If you need a 3d curve just do the same with the z-coordinate what you did with y or x. Hermite curves work in in any dimension. To calculate a hermite curve you need the following vectors: o P1: the startpoint of the curve o T1: the tangent (e.g. direction and speed) how the curve lefts the startpoint o P2: he endpoint of the curve o T2: the tangent (e.g. direction and speed) how the curves enters the endpoint These 4 vectors are simply multiplied with 4 hermite basis functions and summed together. h1(s) = 2s^3 - 3s^2 + 1 h2(s) = -2s^3 + 3s^2 h3(s) = s^3 - 2s^2 + s h4(s) = s^3 - s^2 Take a closer look at functions h1 and h2: o h1 starts at 1 and goes slowly to 0. o h2 starts at 0 and goes slowly to 1. If you now multiply the startpoint with h1 and the endpoint with h2. Let s go from 0 to 1 to interpolate between start and endpoint. h3 and h4 are applied to the tangents in the same manner. They take care, that the curve bends into the desired direction at the start and endpoint. * The Math in Matrix Form All this stuff can be expessed with some vector and matrix algebra. I think the matrix-form is much better to understand. Vector S: The interpolation-point and it's powers up to 3: Vector C: The parameters of our hermite curve: Matrix h: The matrix form of the 4 hermite polyonials: | s^3 | | P1 | | 2 -2 1 1 | S = | s^2 | C = | P2 | h = | -3 3 -2 -1 | | s^1 | | T1 | | 0 0 1 0 | | 1 | | T2 | | 1 0 0 0 | To calculate a point on the curve you build the Vector S, multiply it with the matrix h and the again multiply with C. P = S * h * C A little side-note: Bezier-Curves This matrix-form is valid for all cubic polynomial curves.. The only thing that changes is the polynomial matrix.. For example if you want to draw a Bezier curve instead of hermites you might use this matrix: | -1 3 -3 1 | b = | 3 -6 3 0 | | -3 3 0 0 | | 1 0 0 0 | * Some Pseudocode Sure, this C-style pseudo-code won't compile.. C doesn't come with a power function, and unless you wrote yourself a vector-class any compiler would generate hundrets of erros and make you feel like an idiot. moveto (P1); // move pen to startpoint for (int t=0; t < steps; t++) { float s = (float)t / (float)steps; // scale s to go from 0 to 1 float h1 = 2s^3 - 3s^2 + 1; // calculate basis function 1 float h2 = -2s^3 + 3s^2; // calculate basis function 2 float h3 = s^3 - 2*s^2 + s; // calculate basis function 3 float h4 = s^3 - s^2; // calculate basis function 4 vector p = h1*P1 + // multiply and sum all funtions h2*P2 + // together to build the interpolated h3*T1 + // point along the curve. h4*T2; lineto (p) // draw to calculated point on the curve } o Cardinal splines Cardinal splines are just a subset of the hermite curves. They don't need the tangent points because they will be calculated out of the control points. We'll loose some of the flexibility of the hermite curves, but as a tradeoff the curves will be much easier to use. The formula for the tangents for cardinal splines is: T = a * ( P - P ) i i+1 i-1 a is a constant which affects the tightness of the curve. Write yourself a program and play around with it. ( a should be between 0 and 1, but this is not a must). o Catmull-Rom splines The Catmull-Rom spline again is just a subset of the cardinal splines. You only have to define a as 0.5, and you can draw and interpolate Catmull-Rom splines T = 0.5 * ( P - P ) i i+1 i-1 Catmull-Rom splines are great if you have some data-points and just want to smoothly interpolate between them. o The Kochanek-Bartels Splines (also called TCB-Splines) The kb-splines (mostly known from Autodesks 3d-Studio, Max and Newteks Lightwave) are nothing more than hermite curves and a hand full of formulas to calculate the tangents. These kind of curves have been introduced by D. Kochanek and R. Bartels in 1984 to give animators more control for keyframe animation. They introduced three control-values for each keyframe point: + Tension: How sharply does the curve bend? + Continuity: How hard is the change of speed and direction? + Bias: How is the direction of the curve as it passes through the keypoint The "incomming" Tangent equation: (1-t)*(1-c)*(1+b) TS = ----------------- * ( P - P ) i 2 i i-1 (1-t)*(1+c)*(1-b) + ----------------- * ( P - P ) 2 i+1 i The "outgoing" Tangent equation: (1-t)*(1+c)*(1+b) TD = ----------------- * ( P - P ) i 2 i i-1 (1-t)*(1-c)*(1-b) + ----------------- * ( P - P ) 2 i+1 i When you want to interpolate the curve you should use this vector: | P(i) | C = | P(i+1) | | TD(i) | | TS(i+1) | You might notice, that you always need the previous and following point if you want to calculate the curve. This might be a problem when you try to calculate keyframe data from lightwave or 3D-Studio. [...] Speed Control If you write yourself a keyframe-interpolation code and put it into a program you'll notice one problem: Unless you have your keyframes in fixed intervals you will have a sudden change of speed and direction whenever you pass a keyframe-point. This can be avoided if you take the number of key-positions (frames) between two keyframes into account: N is the number of frames (seconds, whatever) between two keypoints. 2 * N i-1 TD = TD * --------------- adjustment of outgoing tangent i i N + N i-1 i 2 * N i TS = TS * --------------- adjustment of incomming tangent i i N + N i-1 i ------------------------------------------------------------------------ === Cut ===
следующий фpагмент (5)|пpедыдущий фpагмент (3)
- Algorithms.. (2:5030/84) ------------------------------------- RU.ALGORITHMS - Msg : 971 of 971 From : Dmitry Ovcharov 2:5030/498.4 05 Apr 97 01:25:00 To : Alexey Drozdov 09 Apr 97 20:53:06 Subj : сплайн-интерполяция -------------------------------------------------------------------------------- Здравствуй, Alexey! Среда, 02 Апреля 1997 года, 01:42, Alexey Drozdov написал(а) to All: AD> Привет All! AD> Повтоpяю вопpос: AD> Какой есть хоpоший метод двумеpной сплайн-интеpполяции? AD> Я pеализовал с помощью кубических сплайнов, но больно негладко. Хоpоший метод двумеpной сплайн интеpполяции -- сглаживающие сплайны. Могу пpивести пpоцедуpу на Delphi. Писал уже достаточно давно пpичем на клиппеpе :), в Delphi только индексы испpавил, сами фоpмулы не помню. Там пятидиагональная система, мне лень было искать оптимальный способ ее pешения типа пpогонки, поэтому использовал обычный метод Гаусса (pазмеpы сглаживаемого невелики). Hу и классы всякие - диагpаммы, матpицы в виде обычных массивов можно записать. p -- паpаметp сглаживания. === Cut === procedure TDOSmoothDiagram.CopyFrom( Source : TDODiagram ); function min( a,b : Integer ) : Integer; begin if a<b then Result:= a else Result:= b; end; function max( a,b : Integer ) : Integer; begin if a>b then Result:= a else Result:= b; end; var H,B,f : TMatr; M,i,j,k : Integer; t : Double; begin M:= Source.Count-2; { Количество исходных точек минус 2 } if M<1 then { Сглаживать нечего. } begin inherited CopyFrom(Source); Exit; end; H:= TMatr.Create(3,M); { Вспомогательные матpицы. } B:= TMatr.Create(M,M); { Здесь можно соптимизиpовать. } f:= TMatr.Create(1,M); for i:= 1 to M do { Стpоим систему. } begin H[1,i]:= 1/(Source[i].X-Source[i-1].X); H[2,i]:= -1/(Source[i].X-Source[i-1].X)-1/ (Source[i+1].X-Source[i].X); H[3,i]:= 1/(Source[i+1].X-Source[i].X); end; for i:= 1 to M do begin B[i,i]:= (Source[i+1].X-Source[i-1].X)/3+ (sqr(H[1,i])+sqr(H[2,i])+sqr(H[3,i]))/p; if i<M then begin B[i+1,i]:= (Source[i+1].X-Source[i].X)/6+(H[1,i]+H[3,i])*H[2,i]/p; B[i,i+1]:= B[i+1,i]; if i<M-1 then begin B[i+2,i]:= H[1,i]*H[3,i]/p; B[i,i+2]:= B[i+2,i]; end; end; for j:= i+3 to M do begin B[i,j]:= 0; B[j,i]:= 0; end; f[1,i]:= H[1,i]*Source[i-1].Y+ H[2,i]*Source[i ].Y+ H[3,i]*Source[i+1].Y; end; for k:= 1 to M do { pешаем систему, пpямой ход. } begin f[1,k]:= f[1,k]/B[k,k]; for j:= min(k+2,M) downto k do B[k,j]:= B[k,j]/B[k,k]; for i:= k+1 to min(k+2,M) do begin f[1,i]:= f[1,i]-B[i,k]*f[1,k]; for j:= min(k+2,M) downto k do B[i,j]:= B[i,j]-B[i,k]*B[k,j]; end; end; { обpатный ход. } f[1,M-1]:= f[1,M-1]-B[M-1,M]*f[1,M]; for k:= M-2 downto 1 do f[1,k]:= f[1,k] - B[k,k+1]*f[1,k+1] - B[k,k+2]*f[1,k+2]; { Тепеpь в f -- втоpые пpоизводные сплайна в узлах, обычно обозн. как M } Self.Clear; for i:= 1 to Source.Count do begin t:= Source[i-1].Y; for j:= max(1,i-2) to min(i,M) do t:= t - H[1+i-j,j]*f[1,j]/p; { новое, сглаженное значение в узле } AddPoint(Source[i-1].X,t,Source[i-1].Marker); end; { А тепеpь X,t -- сглаженные узлы, дальше пpименяешь обычную фоpмулу, с использованием втоpых пpоизводных. } H.Free; B.Free; f.Free; end; === Cut === AD> Alexey С уважением, Дмитpий.
следующий фpагмент (6)|пpедыдущий фpагмент (4)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 16 of 61 From : jjh119@jester.usask.ca 2:5030/144.99 21 Jun 94 00:16:08 To : All 26 Jun 94 02:06:56 Subj : SOLVED: Curve through 3 points -------------------------------------------------------------------------------- Thanks again to all who helped out with my problem of drawing a curve through three points. The solution as many pointed out was a parametric Lagrange polynomial. My code for this function is given below: void DrawCurve(CDC* pDC, int x0, int y0, int x1, int y1, int x2, int y2) { long t; long t0 = 0; long t2 = 1024; long r0, r1, r2; int nsteps = 64; CPoint pts[65]; int i = 0; double d1 = sqrt(((double)x1-x0)*(x1-x0)+((double)y1-y0)*(y1-y0)); double d2 = sqrt(((double)x1-x2)*(x1-x2)+((double)y1-y2)*(y1-y2)); long t1 = __max(1L,(long)((d1*1024.0)/(d1+d2))); long tinc = 1024L / nsteps; t1 = __min(t1, 1023L); for (t = 0; t <= 1024; t += tinc) { r0 = (((t-t0)*(t-t1))<<10)/((t2-t0)*(t2-t1)); r1 = (((t-t1)*(t-t2))<<10)/((t0-t1)*(t0-t2)); r2 = (((t-t0)*(t-t2))<<10)/((t1-t0)*(t1-t2)); pts[i].x = (int)((r0*x2 + r1*x0 + r2*x1) >> 10; pts[i++].y = (int)((r0*y2 + r1*y0 + r2*y1)) >> 10; } pDC->Polyline(pts, nsteps); } I wanted to use integer arithmetic as much as possible. The only other thing of note is that t1 (the independent variable value at which the curve passes through x1,y1) is set to the proportion of the total distance between x0,y0 -> x1,y1 -> x2,y2. This allows the 'curve' to approximate a straight line when x1,y1 -> x0,y0 or x1,y1 -> x2, y2. Once again, thanks a lot for all your help James
следующий фpагмент (7)|пpедыдущий фpагмент (5)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 18 of 52 From : lmason@tartarus.uwa.edu.au 2:5030/144.99 30 Jun 94 13:58:56 To : All 07 Jul 94 02:49:14 Subj : Re: Drawing curve from a set of given points. -------------------------------------------------------------------------------- Chauncy Liu (cliu@osf.org) wrote: : Does anyone know how draw a curve line when a set of given points. : I use the spline but it did not go throught these given points. Does : anyone have an example I can look at? The easiest way I can think of would be to use a Catmull-Rom spline. In matrix form : (for a point q on the curve at t) q(t) = 0.5 * [ t^3 t^2 t 1 ] * [ -1 3 -3 1 ] * [ p1 ] [ 2 -5 4 -1 ] [ p2 ] [ -1 0 1 0 ] [ p3 ] [ 0 2 0 0 ] [ p4 ] Where p0,p1,p2,p3 are _points_ along the curve, (note however that the curve drawn actually only passes through points p2 and p3) (sorry my ascii matrices aren't too good :) (that was a constant by a 1x4 matrix * 4x4 matrix * 4x1 matrix) This becomes : q(t) = 0.5 * ((-p1 + 3*p2 -3*p3 + p4)*t*t*t + (2*p1 -5*p2 + 4*p3 - p4)*t*t + (-p1+p3)*t + 2*p2) For values of t between 0 and 1 the curve passes through P2 at t=0 and it passes through P3 at t=1. Another nice property of these curves is that the tangent vector at a point P is parallel to the line joining P's two surrounding points. To do more than two points just step through the array of points using the previous point, the current point and the next two points as the four points for the spline. For each of these segments draw a curve for 0<t<1. This curve will be between the current point and the next point. This method is not particularly fast though. A fast method for display of this class of splines is given in a SIGGRAPH 89 Paper, "A Recursive Evaluation Algorithm for a Class of Catmull-Rom Splines," by Barry, P. and R.Goldsman. Hope this helps, Llew Mason lmason@lethe.uwa.edu.au
следующий фpагмент (8)|пpедыдущий фpагмент (6)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 35 of 73 From : shoemake@thumper 2:5030/144.99 03 Sep 94 07:43:22 To : All 08 Sep 93 02:00:16 Subj : (1) Re: spline interp. of QUATERNIONS -------------------------------------------------------------------------------- In article <patrick.778511828@coral> patrick@coral.cs.jcu.edu.au (Pat Collins) writes: >Hi all, > I was wondering if you could direct me to where I could find the source >code for a spherical cubic spline interpolation routine for quaternions. >Watt^2 makes reference to a routine called squad() but it is not listed in >the book. Furthermore, is the extra calculation worth the effort? Or should >I just stay with spherical linear interpolation (slerp) ? The material on squad is taken from my Siggraph course notes, without citing that fact as I recall. Without going into the ethical (not to mention legal) implications of this practice, we can see its practical disadvantages, eh? These notes are formally unpublished, and implicitly copyright by their authors. They may also contain errors of varying seriousness, since they have not been reviewed (maybe not even by their authors). Nevertheless, they can be a valuable source of information. Still want to try squad? Then from the horse's mouth, here's an implementation: /* Spherical quadrangle curve from qu0 to qu1, for t IN [0..1], with * tangents manipulated by quadrangle points quq0 and quq1. All * quaternions are assumed to be of unit norm. */ Quat Qt_Squad(Quat qu0, Quat qu1, Quat quq0, Quat quq1, float t) { Quat qut; Quat qu01 = Qt_Slerp(qu0, qu1, t); Quat quq01 = Qt_Slerp(quq0, quq1, t); qut = Qt_Slerp(qu01, quq01, 2.0*t*(1.0 - t)); return (qut); } You've probably got this already, but for completeness here's slerp. /* Spherical linear interpolation of unit quaternions. Takes qu0 to qu1 * as t goes from 0 to 1. */ Quat Qt_Slerp(Quat qu0, Quat qu1, float t) { Quat qut; float omega, cosOmega, sinOmega, qu0Part, qu1Part; static float epsilon = 1.0e-5; cosOmega = Qt_Dot(qu0, qu1); if ((1.0 + cosOmega) > epsilon) { /* Usual case */ if ((1.0 - cosOmega) > epsilon) { /* Usual case */ omega = acos(cosOmega); sinOmega = sin(omega); qu0Part = sin((1.0 - t)*omega) / sinOmega; qu1Part = sin(t*omega) / sinOmega; } else { /* Ends very close */ qu0Part = 1.0 - t; qu1Part = t; } } else { /* Ends nearly opposite */ qu1.x = -qu0.y; qu1.y = qu0.x; qu1.z = -qu0.w; qu1.w = qu0.z; qu0Part = sin((0.5 - t) * PI); qu1Part = sin(t * PI); } qut = Qt_Add(Qt_Scale(qu0, qu0Part), Qt_Scale(qu1, qu1Part)); return (qut); } You might also want to check out John Schlag's Gem VIII.4 in the second volume of Graphics Gems. I would probably not use squad for an animation system these days, but slerp gives unacceptable discontinuities in the angular velocity. The original quaternion curve techniques are in my Siggraph 85 paper, "Animating Rotation with Quaternion Curves", which is also in my course notes: "Math for Siggraph", 1989 and 1991. Since then there have been a number of alternatives published by other authors. -- Ken
следующий фpагмент (9)|пpедыдущий фpагмент (7)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 19 of 89 From : gokturk@seas.gwu.edu 2:5030/315 03 Nov 95 04:47:58 To : All 08 Nov 95 02:27:42 Subj : Re: Spline Interpolation -------------------------------------------------------------------------------- X-RealName: Mehmet Gokturk In article <478b33$p9m@duke.cs.duke.edu>, geha@duke.cs.duke.edu (Ashish Gehani) wrote: | Given a grid of N x N points with intensity values in the range 0 to M, | I need to compute a spline so that I may find values of points in between the | grid points. Here are my spline codes for inbetweening. Note that take 4 points for each step and advance by one.. hope this helps ////////////////////////////////////////////////// SetCardinalMatrix ///////// // FUNC : SetCardinalMatrix // SUMMARY : Adapted from the paper by Smith, Spline Tutorial // Notes, Tech Memo No 77 Lucasfilm Ltd. May 84 // INPUT : a, tension value // *m, Matrix4 pointer where the cardinality // matrix is to be set up. // OUTPUT : returns 0 if no error // *m, matrix is constructed...(no allocation is done!) // int SetCardinalMatrix(float a,Matrix4 *m) { m->set(0,1,2.-a); m->set(0,2,a-2.); m->set(1,0,2.*a); m->set(1,1,a-3.); m->set(1,2,3.-2.*a); m->set(3,1,1.); m->set(0,3,a); m->set(2,2,a); m->set(0,0,-a); m->set(1,3,-a); m->set(2,0,-a); m->set(2,1,0.); m->set(2,3,0.); m->set(3,0,0.); m->set(3,2,0.); m->set(3,3,0.); return 0; } ////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////// SetBetaMAtrix /////////// // FUNC : SetBetaMatrix // SUMMARY : Adapted from the paper by Smith, Spline Tutorial // Notes, Tech Memo No 77 Lucasfilm Ltd. May 84 // INPUT : b0 bias // b1 tension // *m Matrix4 pointer to be set up // OUTPUT : returns 0 if no error // *m matrix constructed (no allocation done!) // int SetBetaMatrix(float b0,float b1,Matrix4 *m) { int i,j; float d,b2,b3; b2 = b0*b0; b3 = b0*b2; m->set(0,0,-2.*b3); m->set(0,1,2.*(b1+b3+b2+b0)); m->set(0,2,-2.*(b1+b2+b0+1.)); m->set(1,0,6.*b3); m->set(1,1,-3.*(b1+2.*b3+2.*b2)); m->set(1,2,3.*(b1+2.*b2)); m->set(2,0,-6.*b3); m->set(2,1,6.*(b3-b0)); m->set(2,2,6.*b0); m->set(3,0,2.*b3); m->set(3,1,b1+4.*(b2+b0)); m->set(0,3,2.); m->set(3,2,2.); m->set(1,3,0.); m->set(2,3,0.); m->set(3,3,0.); d = 1./(b1+2.*b3+4.*b2+4.*b0+2.); for(i=0;i<4;i++) for(j=0;j<4;j++) m->set(i,j,(*m)(i,j)*d); return 0; } ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// // FUNC : GetSplineValue // SUMMARY : Calculates the spline value from 4 consequtive // floating numbers and 0<u<1 . // INPUT : a,b,c,d four input numbers in order // u, the important alpha value.. varies from 0 to 1 // m, the precalculated spline matrix.. // OUTPUT : the spline value calculated from above 3 // float GetSplineValue(float a,float b,float c,float d,float u,Matrix4 *m) { float p0,p1,p2,p3; p0 = (*m)(0,0)*a + (*m)(0,1)*b + (*m)(0,2)*c + (*m)(0,3)*d; p1 = (*m)(1,0)*a + (*m)(1,1)*b + (*m)(1,2)*c + (*m)(1,3)*d; p2 = (*m)(2,0)*a + (*m)(2,1)*b + (*m)(2,2)*c + (*m)(2,3)*d; p3 = (*m)(3,0)*a + (*m)(3,1)*b + (*m)(3,2)*c + (*m)(3,3)*d; return(p3 + u*(p2+u*(p1+u*p0)) ); } ////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////// GetSplinePoint /////////////// // FUNC : GetSplinePoint // SUMMARY : Calculates the spline vector from 4 consequtive // vector values and 0<u<1. Unlike the getspline point // this time the matrix should be precalculated.(faster) // point values are not necessarily points or angles.. // INPUT : p1,p2,p3,p4 pointers to be calculated // u the important alpha value // m, the precalculated Spline matrix.see GetSplineValue // OUTPUT : returns a calculated intermediate point // Point GetSplinePoint(Point *p1,Point *p2,Point *p3,Point *p4,float u,Matrix4 * m) { Point k0,k1,k2,k3; k0 = (*m)(0,0)*(*p1) + (*m)(0,1)*(*p2) + (*m)(0,2)*(*p3) + (*m)(0,3)*(*p 4); k1 = (*m)(1,0)*(*p1) + (*m)(1,1)*(*p2) + (*m)(1,2)*(*p3) + (*m)(1,3)*(*p 4); k2 = (*m)(2,0)*(*p1) + (*m)(2,1)*(*p2) + (*m)(2,2)*(*p3) + (*m)(2,3)*(*p 4); k3 = (*m)(3,0)*(*p1) + (*m)(3,1)*(*p2) + (*m)(3,2)*(*p3) + (*m)(3,3)*(*p 4); return(k3 + u*(k2 + u*(k1 + u*k0))); } ////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////// GetSplineRQuat /////////////// // FUNC : GetSplineQuat // SUMMARY : Calculates the spline quaternion from 4 consequtive // quaternion values and 0<u<1. Unlike the getspline point // this time the matrix should be precalculated.(faster) // INPUT : q1,q2,q3,q4 quaternions to be calculated // u the important alpha value // m, the precalculated Spline matrix.see GetSplineValue // OUTPUT : returns a calculated intermediate quaternion // Quaternion GetSplineQuat(Quaternion *q1,Quaternion *q2,Quaternion *q3, Quaternion *q4,float u,Matrix4 *m) { Point t1; float t2; t1 = GetSplinePoint(&(q1->v),&(q2->v),&(q3->v),&(q4->v), u, m); t2 = GetSplineValue((q1->w),(q2->w),(q3->w),(q4->w), u, m); return(Quaternion(t1,t2)); } ////////////////////////////////////////////////////////////////////////////// end for now...originally continues.. -- # Mehmet Gokturk The George Washington University EECS # # gokturk@seas.gwu.edu (703) 893-9002 # # http://www.seas.gwu.edu/student/gokturk/. #
следующий фpагмент (10)|пpедыдущий фpагмент (8)
- [56] Available echoes... (2:5030/84) ------------------------- RU.ALGORITHMS - Msg : 37 of 43 From : Alex S Aganichev 2:5020/604.19 17 Jul 96 22:29:00 To : Andrey Litvinenko Subj : Re: help! -------------------------------------------------------------------------------- Greetings Andrey! Not so long ago (Tuesday July 16 1996) Andrey Litvinenko wrote to All: AL> Есть график работы прибора, построенный по точкам с довольно большой AL> дискретой. AL> Требуется нарисовать график с плавными переходами (т.е. сплайнами) AL> Если у кого-нибудь есть алгоритм для рисования таких графиков - AL> киньте please :)) Вот исходничек сплайн-аппpоксиматоpа под Borland C (извиняюсь за отсутствие комментаpиев и фpазочки в printf'ах, но, IMHO, pазобpаться немного можно): === Cut === /* * Splain.C * * Copyright (C) 1995, Alex. S. Aganichev * */ #include <stdio.h> #include <math.h> #include <graphics.h> #define MAX 300 double XMIN, XMAX, XGR, YMIN, YMAX, YGR; #define LEFT 20 #define TOP 20 #define RIGHT 20 #define BOTTOM 20 #define STEP 1 #define RADIUS 2 FILE *in1, /* Массив точек по х */ *in2; /* Массив точек по y */ int GraphDriver, GraphMode, MaxX, MaxY, MaxColors, XZ, YZ, XOFFS, YOFFS, n, k, j, i; double x [ MAX ], y [ MAX ], l [ MAX ], m [ MAX ], r [ MAX ], s [ MAX ], d, e, f, g, h, STEPX, STEPY, p; int main ( int argc, char **argv ) { if ( argc < 3 ) { printf ( "Splain\n", argv [ 0 ] ); printf ( "Usage: %s in1 in2\n", argv [ 0 ] ); return -1; } in1 = fopen ( argv [ 1 ], "rt" ); if ( in1 == NULL ) { fprintf ( stderr, "Could not open file %s\n", argv [ 1 ] ); return -1; } in2 = fopen ( argv [ 2 ], "rt" ); if ( in2 == NULL ) { fprintf ( stderr, "Could not open file %s\n", argv [ 2 ] ); return -1; } i = 0; for ( ;; ) { if ( fscanf ( in1, "%lf", x + i ) == EOF ) break; printf ( "%lf\t", x[i] ); if ( fscanf ( in2, "%lf", y + i ) == EOF ) break; printf ( "%lf\n", y[i] ); if ( i == MAX - 1 ) { fprintf ( stderr, "Warning: max (%d) number of dots reached.\n" ); fprintf ( stderr, "Press a key...\n" ); getch (); break; } ++ i; } fclose ( in2 ); fclose ( in1 ); if ( i < 3 ) { fprintf ( stderr, "Not enough data to splain.\n" ); return -1; } n = i; /* Соpтиpуем точки по возpастанию x методом Шелла */ for ( k = n >> 1; k; k >>= 1 ) for ( i = k; i < n; i ++ ) for ( j = i - k; ( j >= 0 ) && ( x [ j ] > x [ j + k ] ); j -= k ){ e = x [ j ]; x [ j ] = x [ j + k ]; x [ j + k ] = e; e = y [ j ]; y [ j ] = y [ j + k ]; y [ j + k ] = e; } XMIN = x [ 0 ]; XMAX = x [ n - 1 ]; YMIN = 1.7e308; YMAX = -1.7e308; for ( i = 0; i < n; i ++ ) { YMAX = ( y [ i ] > YMAX ) ? y [ i ] : YMAX; YMIN = ( y [ i ] < YMIN ) ? y [ i ] : YMIN; } printf ( "XMIN=%lf\tXMAX=%lf\tYMIN=%lf\tYMAX=%lf\n" "Would you like to change anything [Y/N]?\n", XMIN, XMAX, YMIN, YMAX ); switch ( getch () ) { case 'Y': case 'y': printf ( "XMIN=" ); scanf ( "%lf", &XMIN ); printf ( "XMAX=" ); scanf ( "%lf", &XMAX ); printf ( "YMIN=" ); scanf ( "%lf", &YMIN ); printf ( "YMAX=" ); scanf ( "%lf", &YMAX ); break; default: break; } XGR = XMAX - XMIN; YGR = YMAX - YMIN; if ( ( XGR <= 0.0 ) || ( YGR <= 0.0 ) ) { fprintf ( stderr, "Such parameters are not allowed !!!\n" ); return -99; } printf ( "STEPX=" ); scanf ( "%lf", &STEPX ); printf ( "STEPY=" ); scanf ( "%lf", &STEPY ); GraphDriver = DETECT; initgraph ( &GraphDriver, &GraphMode, "" ); if ( graphresult () != grOk ) { fprintf ( stderr, "Could not initialize graphics: check availability of driver.\n" ); return -1; } MaxColors = getmaxcolor () + 1; MaxX = getmaxx (); MaxY = getmaxy (); d = x [ 1 ] - x [ 0 ]; e = ( y [ 1 ] - y [ 0 ] ) / d; l [ 0 ] = 0; s [ 0 ] = 0; for ( i = 1; i < n; i ++ ) { h = d; d = x [ i + 1 ] - x [ i ]; f = e; e = ( y [ i + 1 ] - y [ i ] ) / d; l [ i ] = d / ( d + h ); r [ i ] = 1 - l [ i ]; s [ i ] = 6 * ( e - f ) / ( h + d ); } for ( i = 1; i < n; i ++ ) { p = 1 / ( r [ i ] * l [ i - 1 ] + 2 ); l [ i ] = - l [ i ] * p; s [ i ] = ( s [ i ] - r [ i ] * s [ i - 1 ] ) * p; } m [ n - 1 ] = 0; m [ n - 2 ] = l [ n - 2 ] = s [ n - 2 ]; for ( i = n - 3; i > 0; i -- ) { l [ i ] = l [ i ] * l [ i + 1 ] + s [ i ]; m [ i ] = l [ i ]; } setviewport ( LEFT, TOP + 16, MaxX - RIGHT - 16, MaxY - BOTTOM, 1 ); MaxX -= LEFT + RIGHT + 16; MaxY -= TOP + BOTTOM + 16; STEPX *= MaxX / XGR; STEPY *= MaxY / YGR; for ( k = 0; k < MaxX; k += STEP ) { i = 0; f = ( (double) k * XGR / (double) MaxX ) + XMIN; if ( f > x [ n - 1 ] ) { d = x [ n - 1 ] - x [ n - 2 ]; e = d * m [ n - 2 ] / 6.0 + ( y [ n - 1 ] - y [ n - 2 ] ) / d; e = e * ( f - x [ n - 1 ] ) + y [ n - 1 ]; } else if ( f <= x [ 0 ] ) { d = x [ 1 ] - x [ 0 ]; e = d * m [ 1 ] / 6.0 - ( y [ 1 ] - y [ 0 ] ) / d; e = y [ 0 ] - e * ( f - x [ 0 ] ); } else { while ( f > x [ i ] ) ++ i; d = x [ i ] - x [ i - 1 ]; h = f - x [ i - 1 ]; g = x [ i ] - f; p = d * d / 6.0; e = ( m [ i - 1 ] * g * g * g + m [ i ] * h * h * h ) / 6.0 / d; e += ( ( y [ i - 1 ] - m [ i - 1 ] * p ) * g + ( y [ i ] - m [ i ] * p ) * h ) / d; } putpixel ( k, MaxY - (int) ( ( e - YMIN ) * (double) MaxY / YGR ), EGA_WHITE ); } for ( i = 0; i < n; i ++ ) { circle ( (int) ( ( x [ i ] - XMIN ) * (double) MaxX / XGR ), MaxY - ( ( y [ i ] - YMIN ) * (double) MaxY / YGR ), RADIUS ); } MaxX = getmaxx (); MaxY = getmaxy (); setviewport ( 0, 0, MaxX, MaxY, 1 ); XZ = LEFT - XMIN * ( MaxX - LEFT - RIGHT - 16 ) / XGR; if ( XZ < LEFT ) XZ = LEFT; else if ( XZ > MaxX - RIGHT ) XZ = MaxX - RIGHT; YZ = MaxY - BOTTOM + YMIN * ( MaxY - TOP - BOTTOM - 16 ) / YGR; if ( YZ < TOP ) YZ = TOP; else if ( YZ > MaxY - BOTTOM ) YZ = MaxY - BOTTOM; line ( XZ, TOP, XZ, MaxY - BOTTOM ); line ( XZ - 2, TOP + 6, XZ, TOP ); line ( XZ + 2, TOP + 6, XZ, TOP ); line ( LEFT, YZ, MaxX - RIGHT, YZ ); line ( MaxX - RIGHT - 6, YZ - 2, MaxX - RIGHT, YZ ); line ( MaxX - RIGHT - 6, YZ + 2, MaxX - RIGHT, YZ ); XOFFS = ( XZ - LEFT ) % (int) STEPX; k = ( MaxX - LEFT - RIGHT - XOFFS - 16 ) / STEPX; for ( i = 0; i <= k + 1; i ++ ) line ( LEFT + XOFFS + i * STEPX, YZ - 2, LEFT + XOFFS + i * STEPX, YZ + 2 ); YOFFS = ( MaxY - BOTTOM - YZ ) % (int) STEPY; k = ( MaxY - TOP - BOTTOM - YOFFS - 16 ) / STEPY; for ( i = 0; i <= k + 1; i ++ ) line ( XZ - 2, MaxY - BOTTOM - YOFFS - i * STEPY, XZ + 2, MaxY - BOTTOM - YOFFS - i * STEPY ); getch (); closegraph (); return 0; } === Cut ===
следующий фpагмент (11)|пpедыдущий фpагмент (9)
- Computer Graphics (2:5030/84) --------------------------------- SU.GRAPHICS - Msg : 2668 of 2682 From : Dmitry Apanovitch 2:5020/902.6 19 Nov 97 22:25:00 To : All 23 Nov 97 01:15:16 Subj : [PROG] Reduce, обработка векторной графики ... ------------------------------------------------------------------------------- Мое почтение All! Здесь предстaвлен aлгоритм интерполяции пaрaметрической функции полиномaми 3-ей степени, т.е. f(t) = a*t^3 + b*t^2 + c*t + d. Причем, количество описывaющих функцию полиномов, в пределaх допустимой погрешности, - _минимaльно_. В CorelDraw Х.Х aнaлогичные действия выполняет оперaция сокрaщения узлов кривой - Reduce, только горaздо медленней и менее точно. Дaнный aлгоритм используется при экспорте рaзличных векторных изобрaжений в фaйлы формaтa Adobe Illustrator ( *.ai ), т.е. рaсчитывaются пaрaметры кривых Безье, минимaльным количеством которых можно описaть это векторное изобрaжение. {============================================================================} const NN = 1000; { Kоличество итерaций, требуемых для построения исходной функции } DltR = 0.001;{ Haибольшее допустимое отклонение угловых коэффициентов нa концaх куб. полиномa и в соответствующих точкaх исходной функции. Здесь зaдaется точность. Чем меньше число тем выше точность полученного результaтa} type { В этом клaссе описывaется куб. полином. [ KП ] } TReduceLink = class { Этот метод рaсчитывaет знaчение X и Y для пaрaметрa ii, 0<=ii<=1 } procedure GetLinkXY(ii: extended; var xx,yy: extended); { Этот метод рaсчитывaет по 4 точкaм знaчения коэффициентов KП: ax,ay - при 3-ей степени, bx,by - при 2-ей степени, cx,cy - при 1-ей степени, dx,dy - при 0-ей степени } procedure SetLinkXY(x0,y0,x1,y1,x2,y2,x3,y3: extended); private fax,fbx,fcx, FX fay,fby,fcy, FY: extended; public property AX: extended read fax write fax; property AY: extended read fay write fay; property BX: extended read fbx write fbx; property BY: extended read fby write fby; property CX: extended read fcx write fcx; property CY: extended read fcy write fcy; property DX: extended read fx write fx ; property DY: extended read fy write fy ; end; var ТТ: longint; { Kоличество KП описывaющих исходную функцию. } { Мaссив KП описывaющих исходную функцию. } Links: array [0..1000] of TReduceLink; { TReduceLink =========================================================} procedure TReduceLink.GetLinkXY(ii: extended; var xx,yy: extended); begin t2:=sqr(ii); {квaдрaт} t3:=ii*t2 ; {куб } xx:=fax*t3 + fbx*t2 + fcx*ii + FX; yy:=fay*t3 + fby*t2 + fcy*ii + FY; end; {процедурa рaсчетa коэффициентов кубического полиномa по 4 точкaм, где x.-aргументы} procedure CalcABCD_3D(vv0,vv1,vv2,vv3, x0,x1,x2,x3: extended; var a,b,c,d: extended); var M: array [0..3, 0..4] of extended; ii,jj,tt: TIndex; tmp: extended; begin M[0,0]:=1; M[0,1]:=x0; M[0,2]:=x0*x0; M[0,3]:=x0*x0*x0; M[0,4]:=vv0; M[1,0]:=1; M[1,1]:=x1; M[1,2]:=x1*x1; M[1,3]:=x1*x1*x1; M[1,4]:=vv1; M[2,0]:=1; M[2,1]:=x2; M[2,2]:=x2*x2; M[2,3]:=x2*x2*x2; M[2,4]:=vv2; M[3,0]:=1; M[3,1]:=x3; M[3,2]:=x3*x3; M[3,3]:=x3*x3*x3; M[3,4]:=vv3; ii:=0; repeat tmp:=M[ii,ii]; if tmp<>0 then begin {если эл-т глaвной диaгонaли<>0 то} for jj:=0 to 4 do {получaем 1 нa глaвной диaгонaли} M[ii,jj]:=M[ii,jj]/tmp; for tt:=0 to 3 do begin {вычитaем строку} tmp:=M[tt,ii]; if tt<>ii then for jj:=0 to 4 do M[tt,jj]:=M[tt,jj] - M[ii,jj]*tmp; end{for tt}; inc(ii); end{if tmp} else inc(ii); until ii=4; a:=M[3,4];{коэфф. при 3 степени} b:=M[2,4]; c:=M[1,4]; d:=M[0,4];{свободный член} end; procedure TReduceLink.SetLinkXY(x0,y0,x1,y1,x2,y2,x3,y3: extended); begin CalcABCD_3D( x0,x1,x2,x3, 0, 1/3, 2/3, 1, fax, fbx, fcx, ffx); CalcABCD_3D( y0,y1,y2,y3, 0, 1/3, 2/3, 1, fay, fby, fcy, ffy); FX:=ffx; FY:=ffy; end;
следующий фpагмент (12)|пpедыдущий фpагмент (10)
{ Дaннaя процедурa считaет знaчения функции xx и yy по aргументу ii. В дaнном случaе процедурa является aбстрaктной, т.к. функция не определенa. Если функция будет определенa нa интервaле [i0..i1], то необходимо в пределaх этой процедуры провести преобрaзовaние aргументa, примерно тaк: jj:=i0 + (i1-i0)*ii/NN } procedure GetFuncXY(ii: extended; var xx,yy: extended); { 0 <= ii =< NN } begin xx := FX(ii); yy := FY(ii); end; { Процедурa рaсчетa мaссивa KП. Результaтом ее рaботы является мaссив KП, минимaльным количеством которых можно описaть исходную функцию, в пределaх допустимой погрешности - DltR } procedure ReduceFunc; var I0, I1: longint; dd, xx,yy, Lconst, Lcr, x0,y0, x1,y1, x2,y2, x3,y3: extended; tmp: boolean; begin { I0 - нaчaло учaсткa исх. функции, I1 - конец учaсткa исх. функции. Им соотв. знaчения aргументa KП 0 и 1 соотв. } I0:=0; I1:=3; TT:=0; {Kоличество KП в мaссиве} while I0<NN do begin dd:=(I1-I0)/3; GetFuncXY(I0 , x0,y0); { знaчение функции в точке I0 } GetFuncXY(I0+0.0001, xx,yy); { знaчение функции в точке I0 + мaлое прирaщение, т.е. численно считaем производную. } Lconst:=ArcTan2(yy-y0, xx-x0);{ вычисляем уголовой коэфф. в I0 } repeat dd:=(I1-I0)/3; GetFuncXY(I0 +1*dd, x1,y1); GetFuncXY(I0 +2*dd, x2,y2); GetFuncXY(I0 +3*dd, x3,y3); ax:=9*( x3 - 3*x2 + 3*x1 - x0 )/2; { вычисляем коэфф. KП } bx:=9*( x2 - 2*x1 + x0 - 2*ax/9)/2; cx:=3*( x1 - x0 - ax/27 - bx/9 ); dx:=x0; ay:=9*( y3 - 3*y2 + 3*y1 - y0 )/2; by:=9*( y2 - 2*y1 + y0 - 2*ay/9 )/2; cy:=3*( y1 - y0 - ay/27 - by/9 ); dy:=y0; Lcr:=ArcTan2( cy, cx ); { отношение 2х производных в t=0 } GetFuncXY(I1-0.0001 , xx,yy); Lconst1:=ArcTan2(y3-yy, x3-xx); Lcr1:=ArcTan2( ay*(1-0.9999999*0.9999999*0.9999999) + by*(1-0.9999999*0.9999999) + cy*(1-0.9999999), ax*(1-0.9999999*0.9999999*0.9999999) + bx*(1-0.9999999*0.9999999) + cx*(1-0.9999999)); { ^ отношение 2х производных в t=1 } tmp:=( abs(Lconst1-Lcr1)<DltR ) and ( abs(Lconst-Lcr)<DltR ); { ^ если отклонение уголового коэфф. в KП(0)<DltR и KП(1)<DltR, то дaнный KП описывaет исх. функцию нa учaстке [I0..I1] и можно провести следующую итерaцию по уточнению коэфф. этого KП } I1:=I1+1; until not tmp or (I1>=NN+1); { выход из циклa если KП описaл учaсток функ. [I0..I1] и не подходит для описaния учaсткa [I0..I1+1], или если вся функция просчитaнa } if not tmp and not (I1>=NN+1) then I1:=I1-2 else I1:=I1-1; dd:=(I1-I0)/3; { делим учaсток функ. определяемой KП нa 3 чaсти } { берем соотв. знaчения функ. нa этом учaстке } GetFuncXY(I0 +1*dd, x1,y1); GetFuncXY(I0 +2*dd, x2,y2); GetFuncXY(I0 +3*dd, x3,y3); { создaем KП от этих знaчений, и зaносим его в мaссив ... } rl:=TReduceLink.Create; rl.SetLinkXY(x0,y0, x1,y1, x2,y2, x3,y3); Links[TT]:=rl; inc(TT); I0:=I1; I1:=I0+3; if I1>NN then I1:=NN; end {while I0<tt1...}; end; {===========================================================================} Приношу извинения зa достaточно кривое изложение сути aлгоритмa... Дaнный текст является чaстью весьмa рaботоспособной грaфической библиотеки, рaботaющей в облaсти мaкетировaния ценных бумaг (фоновые сетки, розетки, рaмки и т.д.). Библиотекa включaет в себя более 50 рaзличных объектов, преднaзнaченных для рaботы с векторной грaфикой и позволяет создaвaть прогрaммы обрaботки векторной грaфики подобные CorelDraw. С возможностью детaлизaции _векторных_ объектов до 10Е-11 степени. Если кого интерисуют подробности или есть вопросы или пожелaния и предложения то либо в эху, либо мылом нa Origin, либо по тел. (095)422-4957 Дмитрий. P.S. Можно зaлить чaсть библиотеки в эху (если кого интерисует )... Зa сим рaзрешите отклaняться, Митяй.
следующий фpагмент (13)|пpедыдущий фpагмент (11)
A simple spline based animation path Part 1: Given a set of keyframe positions, V0 ... Vn, construct a set of degree 3 Bezier curve segments that smoothly interpolates between them. A degree 3 Bezier curve requires 4 control points, the two endpoints of the curve, P1 and P4, and two others, P2 and P3. For the curve segment joining V[a] to V[a+1], P1 ... P4 are given by: P1 = V[a] P4 = V[a+1] D = length( V[a+1] - V[a] ) d1 = length( V[a+1] - V[a-1] ) d2 = length( V[a+2] - V[a] ) U1 = V[a+1] - V[a-1] U2 = V[a+2] - V[a] k = a constant < 1/2 P2 = P1 + ( D / (k*d1) ) * U1 P3 = P4 - ( D / (k*d2) ) * U2 Using this formulation, the curve segment will interpolate V[a] and V[a+1] between t=0 and t=1. Furthermore, any 2 adjoining curve segments in the path will be G1 continuous, but not necessarily G2 continuous. Note that the path is undefined before V[1] and after V[n-1]. This means that no V[0]...V[1] or V[n-1]...V[n] curve segments will occur as in the animation.
следующий фpагмент (14)|пpедыдущий фpагмент (12)
> : > I am writing a keyframer and I have implememnted LINEAR interpolations > : > (slerps) for rotations, but I am having trouble introducing continuity. > : > : I want to use the two quaternions and the tangents to that quaternions > : in quaternion-space to interpolate similar to Hermite-Cubic-Splines in > : 3d. But how to do that interpolation? > > Yeah, count me in too. I also want to do true quaternion interpolation > using Hermite curves (tens, cont, bias). I guess we can count in most of > the other programmers too that are extracting data from 3D Studio .3ds files. You're so right :-)) There was a code snippet on Autodesk's server some time ago: --------------8<--------------- How To Calculate KXP Spline Derivatives A Tips & Techniques document for use with the IPAS Toolkit Document # 746 Introduction: Several IPAS developers have requested information regarding how the derivatives for 3DS spline interpolation are computed for KXP development. This document contains a fairly uncommented example that will demonstrate 3D Studio's implementation of spline interpolation. Example 3D Studio Key Interpolation Code The following is an example of 3D Studio's key interpolation code for KXP IPAS developers. Note that position keys are handled just like any other track type except for rotations. Those are handled as detailed in the CompAB() function below. NOTE: This is unsupported code; use at your own risk. static void CompElementDeriv( float pp, float p, float pn, float *ds, float *dd,float ksm, float ksp, float kdm, float kdp ) { float delm, delp; delm = p - pp; delp = pn - p; *ds = ksm*delm + ksp*delp; *dd = kdm*delm + kdp*delp; } /*------------------------------------------------------- This computes the derivative at key, as a weighted average of the linear slopes into and out of key, the weights being determined by the tension and continuity parameters. Actually two derivatives are computed at key: "ds" is the "source derivative", or "arriving derivative" "dd" is the "destination derivative" or "departing derivative" ----------------------------------------------------------*/ static void CompDeriv( PosKey *keyp, PosKey *key, PosKey *keyn ) { int i; /* Full TCB computation */ float tm,cm,cp,bm,bp,tmcm,tmcp,ksm,ksp,kdm,kdp,delm,delp,c; float dt,fp,fn; /* fp,fn apply speed correction when continuity is 0.0 */ dt = .5 * (float)( keyn->frame - keyp->frame ); fp = ( (float)( key->frame - keyp->frame ) ) / dt; fn = ( (float)( keyn->frame - key->frame ) ) / dt; c = fabs( key->cont ); fp = fp + c - c * fp; fn = fn + c - c * fn; cm = 1.0 - key->cont; tm = 0.5 * ( 1.0 - key->tens ); cp = 2.0 - cm; bm = 1.0 - key->bias; bp = 2.0 - bm; tmcm = tm*cm; tmcp = tm*cp; ksm = tmcm*bp*fp; ksp = tmcp*bm*fp; kdm = tmcp*bp*fn; kdp = tmcm*bm*fn; for( i = X; i <= Z; i++ ) { CompElementDeriv( keyp->pos[i], key->pos[i], keyn->pos[i], &key->ds[i], &key->dd[i], ksm, ksp, kdm, kdp ); } } /* ----------------------------------------------------------- Compute the "a" and "b" terms at key "cur", which determine the incoming and outgoing tangents (in quaternion space ) -----------------------------------------------------------*/ static int CompAB( RotKey *prev, RotKey *cur, RotKey *next ) { int i; float qprev[4],qnext[4],q[4],qzero[4]; float qp[4],qm[4],qa[4],qb[4],qae[4],qbe[4]; float tm,cm,cp,bm,bp,tmcm,tmcp,ksm,ksp,kdm,kdp,c; float dt,fp,fn; if( prev != NULL ) { if( cur->angle > TWOPI-.00001 ) { COPY_POINT3( q, cur->axis ); q[3] = 0; qlog( q,qm ); } else { qcopy( qprev, prev->quat ); if( qdot( qprev, cur->quat ) < 0 ) qnegate( qprev ); qlndif( qprev, cur->quat, qm ); } } if( next != NULL ) { if( next->angle > TWOPI-.00001 ) { COPY_POINT3( q, next->axis ); q[3] = 0; qlog( q, qp ); } else { qcopy( qnext, next->quat ); if( qdot( qnext, cur->quat ) < 0 ) qnegate( qnext ); qlndif( cur->quat, qnext, qp ); } } if( prev == NULL ) qcopy( qm, qp ); if( next == NULL ) qcopy( qp, qm ); fp = fn = 1.0; cm = 1.0 - cur->cont; if( prev && next ) { dt = 0.5 * (float)(next->frame - prev->frame ); fp = ((float)(cur->frame - prev->frame))/dt; fn = ((float)(next->frame - cur->frame))/dt; c = fabs( cur->cont ); fp = fp + c - c * fp; fn = fn + c - c * fn; } tm = .5*(1.0 - cur->tens); cp = 2.0 - cm; bm = 1.0 - cur->bias; bp = 2.0 - bm; tmcm = tm * cm; tmcp = tm * cp; ksm = 1.0 - tmcm * bp * fp; ksp = -tmcp * bm * fp; kdm = tmcp * bp * fn; kdp = tmcm * bm * fn - 1.0; for( i = 0; i < 4; i++ ) { qa[i] = .5 * ( kdm * qm[i] + kdp * qp[i] ); qb[i] = .5 * ( ksm * qm[i] + ksp * qp[i] ); } qexp( qa, qae ); qexp( qb, qbe ); qmul( cur->quat, qae, cur->a ); qmul( cur->quat, qbe, cur->b ); return TRUE; } /* Subject: Re: Question re derivs.c (fwd) It looks like deriv.c didn't contain the ease function -- just the geometric derivatives. In that case, here's the ease function: */ /* Ease: remap parameter between two keys to apply eases Call this with a = easeFrom of key[n], b = easeTo of key[n+1], and u = (t-t[n]/(t[n+1]-t[n]) -------------------------------------*/ local float Ease(float u, float a, float b) { float k; float s = a+b; if (s==0.0) return(u); if (s>1.0) { a = a/s; b = b/s; } k = 1.0/(2.0-a-b); if (u<a) return((k/a)*u*u); else if (u<1.0-b) return( k*(2*u - a)); else { u = 1.0-u; return(1.0-(k/b)*u*u); } } --------------8<--------------- Please let me know when you got it working ;-) Bert -- bert@isg.cs.uni-magdeburg.de http://www.cs.uni-magdeburg.de/~bert/
следующий фpагмент (15)|пpедыдущий фpагмент (13)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 39 of 44 From : s795238@dutiwy.twi.tudelft.nl 2:5030/315 02 Dec 95 19:14:40 To : All 04 Dec 95 06:03:14 Subj : Re: Drawing Curves. -------------------------------------------------------------------------------- X-RealName: Menno Victor van der Star aofranco@newstand.syr.edu (00Soul) wrote: >I Need to draw a "best-fit" curve from a collection of 2-d points in an >array. Anyone know where I could find an algorithm to do such a thing? There are many ways to draw curves through series of points. Here are four procedures (in pascal) to draw curves : Curve : Curve through 3 points CubicBezierCurve : Curve through 4 points BSpline : Curve through any number of points Catmull_Rom_Spline : Curve trhough any number of points The procedures use a line procedure which should draw a line from (x1,y1) to (x2,y2). The type point is defined as : Point = Record x, y : Integer; End; Here are the procedures : Procedure Curve (x1, y1, x2, y2, x3, y3 : Integer; Segments : Word); { Draw a curve from (x1,y1) through (x2,y2) to (x3,y3) divided in 'Segments' segments } Var lsteps, ex, ey, fx, fy : LongInt; t1, t2 : Integer; Begin x2:=(x2 SHL 1)-((x1+x3) SHR 1); y2:=(y2 SHL 1)-((y1+y3) SHR 1); lsteps:=Segments; If (lsteps<2) then lsteps:=2; If (lsteps>128) then lsteps:=128; { Clamp value to avoid overcalculation } ex:=(LongInt (x2-x1) SHL 17) DIV lsteps; ey:=(LongInt (y2-y1) SHL 17) DIV lsteps; fx:=(LongInt (x3-(2*x2)+x1) SHL 16) DIV (lsteps*lsteps); fy:=(LongInt (y3-(2*y2)+y1) SHL 16) DIV (lsteps*lsteps); Dec (lsteps); While lsteps>0 Do Begin t1:=x3; t2:=y3; x3:=(((fx*lsteps+ex)*lsteps) SHR 16)+x1; y3:=(((fy*lsteps+ey)*lsteps) SHR 16)+y1; Line (t1,t2,x3,y3); Dec (lsteps); End; Line (x3,y3,x1,y1); End; Procedure CubicBezierCurve (x1, y1, x2, y2, x3, y3, x4, y4 : Integer; Segments : Word); { Draw a cubic bezier-curve using the basis functions directly } Var tx1, tx2, tx3, ty1, ty2, ty3, mu, mu2, mu3, mudelta : Real; xstart, ystart, xend, yend, n : Integer; Begin If (Segments<1) then Exit; If Segments>128 then Segments:=128; { Clamp value to avoid overcalculation } mudelta:=1/Segments; mu:=0; tx1:=-x1+3*x2-3*x3+x4; ty1:=-y1+3*y2-3*y3+y4; tx2:=3*x1-6*x2+3*x3; ty2:=3*y1-6*y2+3*y3; tx3:=-3*x1+3*x2; ty3:=-3*y1+3*y2; xstart:=x1; ystart:=y1; mu:=mu+mudelta; For n:=1 to Segments Do Begin mu2:=mu*mu; mu3:=mu2*mu; xend:=Round (mu3*tx1+mu2*tx2+mu*tx3+x1); yend:=Round (mu3*ty1+mu2*ty2+mu*ty3+y1); Line (xstart, ystart, xend, yend); mu:=mu+mudelta; xstart:=xend; ystart:=yend; End; End; Procedure BSpline (NumPoints : Word; Var Points : Array Of Point; Segments : Word); { Draw a BSpline approximating a curve defined by the array of points. } { Beware! A B-Spline generaly does not pass through the points defining } { it ! } Function Calculate (mu : Real; p0, p1, p2, p3 : Integer) : Integer; Var mu2, mu3 : Real; Begin mu2:=mu*mu; mu3:=mu2*mu; Calculate:=Round ((1/6)*(mu3*(-p0+3*p1-3*p2+p3)+ mu2*(3*p0-6*p1+3*p2)+ mu *(-3*p0+3*p2)+(p0+4*p1+p2))); End; Var mu, mudelta : Real; x1, y1, x2, y2, n, h : Integer; Begin If (NumPoints<4) Or (NumPoints>16383) then Exit; mudelta:=1/Segments; For n:=3 to NumPoints-1 Do Begin mu:=0; x1:=Calculate (mu,Points[n-3].x,Points[n-2].x,Points[n-1].x,Points[n].x); y1:=Calculate (mu,Points[n-3].y,Points[n-2].y,Points[n-1].y,Points[n].y); mu:=mu+mudelta; For h:=1 to Segments Do Begin x2:=Calculate (mu,Points[n-3].x,Points[n-2].x,Points[n-1].x,Points[n].x); y2:=Calculate (mu,Points[n-3].y,Points[n-2].y,Points[n-1].y,Points[n].y); Line (x1, y1, x2, y2); mu:=mu+mudelta; x1:=x2; y1:=y2; End; End; End; Procedure Catmull_Rom_Spline (NumPoints : Word; Var Points : Array Of Point; Segments : Word); { Draw a spline approximating a curve defined by the array of points. } { In contrast to the BSpline this curve will pass through the points } { defining is except the first and the last point. The curve will only } { pass through the first and the last point if these points are given } { twice after eachother, like this : } { Array of points : } { } { First point defined twice Last point defined twice } { |-----| |----------| } { (0,0),(0,0),(100,100),....,(150,100),(200,200),(200,200) } { the curve defined by these points will pass through all the points. } Function Calculate (mu : Real; p0, p1, p2, p3 : Integer) : Integer; Var mu2, mu3 : Real; Begin mu2:=mu*mu; mu3:=mu2*mu; Calculate:=Round ((1/2)*(mu3*(-p0+3*p1-3*p2+p3)+ mu2*(2*p0-5*p1+4*p2-p3)+ mu *(-p0+p2)+(2*p1))); End; Var mu, mudelta : Real; x1, y1, x2, y2, n, h : Integer; Begin If (NumPoints<4) Or (NumPoints>16383) then Exit; mudelta:=1/Segments; For n:=3 to NumPoints-1 Do Begin mu:=0; x1:=Calculate (mu,Points[n-3].x,Points[n-2].x,Points[n-1].x,Points[n].x); y1:=Calculate (mu,Points[n-3].y,Points[n-2].y,Points[n-1].y,Points[n].y); mu:=mu+mudelta; For h:=1 to Segments Do Begin x2:=Calculate (mu,Points[n-3].x,Points[n-2].x,Points[n-1].x,Points[n].x); y2:=Calculate (mu,Points[n-3].y,Points[n-2].y,Points[n-1].y,Points[n].y); Line (x1, y1, x2, y2); mu:=mu+mudelta; x1:=x2; y1:=y2; End; End; End; I hope this helps. Ciao, Menno ** Just because you're paranoid don't mean they're not after you... **
следующий фpагмент (16)|пpедыдущий фpагмент (14)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 92 of 138 From : Optix@sv.span.com 2:5030/144.99 19 Aug 94 17:16:00 To : All 21 Sep 93 01:28:14 Subj : (1) Correct Integer B-Spline Alg -------------------------------------------------------------------------------- Kerchink ..... BLAM! Yes, yours truly in an attempt to get the algorithm out quickly made a couple of tiny errors, which will cause the algorithm to fail. So below is a proper posting of the algorithm, complete with a couple of optimisations (many thanks to Gavin Cawley for the ideas). void IntegerBezier(int x1,int y1,int x2,int y2,int x3,int y3,int x4,int y4) { int x12,y12,x23,y23,x34,y34,x123,y123,x234,y234,x1234,y1234; // if the distance between points (x1,y1) and (x4,y4) is < 16 // // plot the line (x1,y1,x4,y4) if (abs(x4-x1)+abs(y4-y1) < 16) { Line(x1,y1,x4,y4); // Draw line from (x1,y1) -> (x4,y4) return; // Return up recursion tree } // Calculate all the mid-points of the lines // x12 = (x1+x2)/2; // The compiler will probably turn // y12 = (y1+y2)/2; // these into shifts, since they // x23 = (x2+x3)/2; // are integer variables // y23 = (y2+y3)/2; // (if the C compiler's decent) // x34 = (x3+x4)/2; // Are there any decent C compiler's // y34 = (y3+y4)/2; // left on the PC??? // x123 = (x12+x23)/2; // Best C compiler I ever used was // y123 = (y12+y23)/2; // Acorn C for RISC OS (for ARM 600) // x234 = (x23+x34)/2; // Best RISC processor available // y234 = (y23+y34)/2; // (Ooops, here come the flames!) // x1234 = (x123+x234)/2; y1234 = (y123+y234)/2; IntegerBezier(x1,x2,x12,y12,x123,y123,x1234,y1234); IntegerBezier(x1234,y1234,x234,y234,x34,y34,x4,y4); } The optimisations and correction are all on the distance checking if statement. Before I had a minus in the middle rather than a plus. Idiot. Secondly the routine used to calculate the distance using the good old distance formula which involved square roots (yipee!). Now the routine uses the 'manhattan distance' which works just fine. One last point the 16 in the if statement is the 'resoultion' of the curve. The greater this number the more jagged the curve, but the quicker it is generated. A preview kind of thing really. Don't make this value to small or you're code might crash if the recursion get too deep. Sorry 'bout the mistakes but I hope somebody finds some use for this. I haven't. Is this worthy of Graphics Gems V ????? |-) Chris Killpack =---------------------------------.-----------------------------------------= | Optix/Chris M Killpack | Sound & Vision BBS [+44] (0)1932 252323 | | InterNet : optix@sv.span.com | The UK's best Internet BBS! | =---------------------------------.-----------------------------------------=
следующий фpагмент (17)|пpедыдущий фpагмент (15)
>I have 3d control points for cubic Bezier curves. I need to draw these >curves. One complication is that a perspective transformation and a clip >needs to be done before drawing. Anyone have any pointers to some PD code >that will do this or papers describing ways to do it? Thanks, You can transform the control points with the perspective transformation onto 2D before rendering. The result will be the same as transforming each of the digitized bezier points in 3d. A very simple way to generate bezier curves is by recursion using subdivision Not the most efficient method but easy to write. p1------q1-------P2 | / \ | | r0--r1---r2 | |/ \ | q0 q2 | | | | | | P0 p3 For each set of control points p(0-3), after subdivision you get 2 sets of control points p0, q0, r0, r1 and r1, r2, q2, p3 where q0 = (p1+p0)/2 q1 = (p1+p2)/2 q2 = (p2+p3)/2 r0 = (q0+q1)/2 r2 = (q1+q2)/2 r1 = (r1+r2)/2 You stop the subdivision by either a fix number of recursion steps or when the polygon p0p1p2p3 is supfficiently flat to be approximated by a straight line. Other bezier drawing methods include forward differencing, and by evaluating the bezier polynomial (which can be made efficient by precalculating weights at fixed parametric points ) Anson Tsao Internet:atsao@tkk.win.net TKK Enterprises, Inc. CIS UserID: 76167,2273 Ontario, Canada (416) 338-9103
следующий фpагмент (18)|пpедыдущий фpагмент (16)
- [26] Usenet echoes (21:200/1) -------------------------------- COMP.GRAPHICS - Msg : 52 of 62 From : cschofie@novell.com 2:5030/144.99 09 Apr 94 18:51:00 To : All Subj : Re: Algorithm for generating Bezier curves in C or C++? -------------------------------------------------------------------------------- In article <2n9fi9$o37@portal.gmu.edu>, K@mason1.gmu.edu (Puru Subedi) wrote: > > Hi Netters, > > I am looking for a algorithm for Bezier curves in C or C++. can > any one give me some pointers. If you have any suggestions, please > E-mail to following address. Thank you in advance. > Here's a routine that I use: void DrawBezier (EndPoint_1, HandlePoint_1, EndPoint_2, HandlePoint_2) Rect *EndPoint_1; Rect *HandlePoint_1; Rect *EndPoint_2; Rect *HandlePoint_2; { #define BEZIER_STEPS 50 #define BEZIER_DIV pow (BEZIER_STEPS, 3) long AX, AY, BX, BY, CX, CY, DX, DY; long T, T2, M, M2, X, Y, PrevX, PrevY; AX = EndPoint_1->left + 2; AY = EndPoint_1->top + 2; BX = (HandlePoint_1->left + 2) * 3; BY = (HandlePoint_1->top + 2) * 3; CX = EndPoint_2->left + 2; CY = EndPoint_2->top + 2; DX = (HandlePoint_2->left + 2) * 3; DY = (HandlePoint_2->top + 2) * 3; MoveTo (AX, AY); PrevX = 0; PrevY = 0; for (T = 0; T < BEZIER_STEPS; T++) { T++; M = (BEZIER_STEPS - T); T2 = T * T; M2 = M * M; X = (M2 * (M * AX + T * BX) + T2 * (M * DX + T * CX)) / BEZIER_DIV; Y = (M2 * (M * AY + T * BY) + T2 * (M * DY + T * CY)) / BEZIER_DIV; if (X != PrevX || Y != PrevY) { LineTo (X, Y); PrevX = X; PrevY = Y; } } } Hope it helps. cschofie@novell.com
следующий фpагмент (19)|пpедыдущий фpагмент (17)
- [97] Computer Graphics (2:5030/84) ----------------------------- SU.GRAPHICS - Msg : 15 of 18 From : Alexander Kharitchev 2:5020/487.21 05 Aug 95 01:30:00 To : D.J. Wolf Subj : Фоpмyла кpивой Безье. -------------------------------------------------------------------------------- Hello D.J.! Monday July 31 1995 11:48, D.J. Wolf wrote to All: DJW> А не подcкажет ли yважаемый All фоpмyлy кpивой Безье? DJW> Hе дайте yмеpеть =) Функция на паскакале для кривой безье с двумя контрольнами точками ******************************************************************************* { K -- precision } Procedure DrawSpline(x1,y1,x2,y2,x3,y3,x4,y4:real;k:Word); Var rx,ry : real; n:word; Begin moveto(x1,y1); for n:=0 to k do Begin rx := {0} 1 * 1 * (1-n/k)*(1- n/k)*(1-n/k) * x1 + {1} 3 * n/k * (1-n/k)*(1- n/k) * x2 + {2} 3 * n*n/k/k * (1-n/k) * x3 + {3} 1 * n*n*n/k/k/k * 1 * x4; ry := {0} 1 * 1 * (1-n/k)*(1- n/k)*(1-n/k) * y1 + {1} 3 * n/k * (1-n/k)*(1- n/k) * y2 + {2} 3 * n*n/k/k * (1-n/k) * y3 + {3} 1 * n*n*n/k/k/k * 1 * y4; lineto(rx,ry); End; End; ******************************************************************************* Если тебе нужна именно формула -- попробуй вывести или покупай умные книжки в магазине "техническая книга" Alexander
следующий фpагмент (20)|пpедыдущий фpагмент (18)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 100 of 121 From : kenhill@indirect.com 2:5030/315 21 Jan 96 21:15:24 To : All 23 Jan 96 01:40:12 Subj : Re: How to calculate bezier curves? HELP! -------------------------------------------------------------------------------- kasaari@walrus.megabaud.fi (Kari Saarilahti) wrote: > Does someone know how to calculate bezier curves i.e. the mathematical > theory behind them? Any help appreciated. How to calculate Bezier curves and the mathematical theory behind them are two very different things. For theory, I'll refer you to Gerald Farin's book "Curves and Surfaces for Computer-Aided Geometric Design." Calculating points along a Bezier is easy. Here is the DeCasteau (Spelling?) algorithm in C: /* Structure of a 2d point */ typedef struct PNT2_tag { double x,y; } PNT2; typedef struct BEZIER_tag { int nPts; PNT2 pt[10]; } BEZIER; /* Compute the bezier point at t (0.0 <= t <= 1.0) */ void BezierPt(BEZIER *bez,double t,PNT2 *pt) { static PNT2 bezPts[10]; register i,j; memmove(bezPts,bez->pt,bez->nPts*sizeof(PNT2)); for (i = bez->nPts; i > 1; i--) for (j = 0; j < i; j++) { bezPts[j].x = (1 - t)*bezPts[j].x + t*bezPts[j+1].x; bezPts[j].y = (1 - t)*bezPts[j].y + t*bezPts[j+1].y; } (*pt) = bezPts[0]; } The algorithm should be clear from this code, the theory won't be. Also, you might want to get a copy of the FAQ for this group. It will save you time, and avoid needless repetition. Regards, Ken ------------------------------------------ kenhill@indirect.com kenhill@anvil5k.com Don't forget to pick up a copy of my new book: "Magneto-hydrodynamic Dynamo Theory--for Dummies"
следующий фpагмент (21)|пpедыдущий фpагмент (19)
- Usenet echoes (21:200/1) -------------------------- COMP.GRAPHICS.ALGORITHMS - Msg : 74 of 175 From : kenhill@indirect.com 2:5030/315 07 Nov 95 09:01:56 To : All 11 Nov 95 06:24:08 Subj : Re: Bezier points for sinewave? -------------------------------------------------------------------------------- X-RealName: Kenneth J. Hill bschor@vms.cis.pitt.edu (Bob Schor) wrote: >>How about 6 pts so we could aproximate an odd function with an odd >>polynomial? I suppose then I would use Reme's algorithm to compute a >>minimax polynomial fit of degree 5, and then solve for the Bezier >>points. If this is a "real world" (tm) problem, I have a small Reme >>program I wrote under Mathematica you are welcome to have. >>-Ken Hill >Ken, > I am interested in this algorithm (I've never heard of "Reme"), and would >be grateful for a reference or two (standard textbooks would be fine) for both >this algorithm and for "how to" figure out the Bezier control points given >some criterion on the waveform. > With some basic understanding of what is going on (after I look up the >references), I'd then would appreciate seeing your algorithm (I learn better >if I have some idea of what's going on]. >Bob Schor Sorry it took so long to respond, I've been away. Here is a routine to generate a Bezier frame for a Sine wave on [0,Pi]. It uses a least squares fit to a Chebyshev polynomial, and then computes Bezier points. The Bezier points are given as: 2 Pi 4 Pi 6 Pi Out[3]= {{0, 0.00659638}, {----, 1.13698}, {----, 3.23746}, {----, -3.23746}, 5 5 5 8 Pi > {----, -1.13698}, {2 Pi, -0.00659638}} 5 The error is about .007. I would round the first and last y value of the control points down to 0 so you have nice behaviour, even though this will increase the maximum error slightly. You asked for references: Any book on numerical analysis will discuss orthogonal polynomials and Chebyshev polynomials in particular. A good reference on Beziers is Farin's _Curves_and_Surfaces_for_CAGD_ cited below. As far as Reme's algorithm goes, I wrote my program based on a lecture given by Dr. Jackiewicz at Arizona State. Therefore I don't have any ready references, but I'm sure you could telnet to your university library and find some. Unfortunately, I seem to only have a hardcopy of my Reme program. If someone *really* needs it, I'll fax him or her a copy and that person can do the re-typing. In exchange I'd like the ascii file e-mailed back to me so I have a executable copy. (I still can't figure out how I lost my original :( ) (* Start of mathematica code *) (* Bernstein Basis function *) Bernstein[n_, i_, t_] := Binomial[n, i]*If[i==0,1,t^i]*If[i==n,1,(1 - t)^(n - i)]; (* Coefficients of Chebyshev polynomial *) c[n_,F_,a_,b_] := Table[If[i == 0,1,2]/Pi * NIntegrate[F[((b + a) + x(b - a))/2] * ChebyshevT[i,x] / Sqrt[1 - x^2], {x, -1+$MachineEpsilon, 1-$MachineEpsilon}],{i,0,n}]; (* M1 converts from the Chebyshev Basis to the monomial basis (t^i) *) M1[n_] := Table[D[ChebyshevT[i,t],{t,j}]/(j!),{i,0,n},{j,0,n}]/.t->0; (* M2 converts from the Bernstein Basis to the monomial basis *) (* See Farin, _Curves_and_Surfaces_for_CAGD_, 3rd Ed. page 59 *) M2[n_] := Table[(-1)^(j-i) Binomial[n,j] Binomial[j,i],{i,0,n},{j,0,n}]; (* Compute the Chebyshev coefficients and change to monomial basis *) a1 = N[c[5,Sin,0,2 Pi]].M1[5]; (* a1 is the function transformed into the range [-1,1], so transform to the Bernstein range [0,1] *) a2 = CoefficientList[Sum[a1[[i+1]] * (2 t -1)^i,{i,0,5}],t]; (* The 6 Bezier points-- *) (* The first vector {0,1,0,0,0,0} provides for the trivial parameterization x(t)=2 Pi t in the monomial (t^i) basis *) b = Transpose[{{0,2 Pi,0,0,0,0},a2}.Inverse[M2[5]]]; (* Plot the Bezier *) ParametricPlot[Sum[b[[i+1]] Bernstein[5,i,t],{i,0,5}],{t,0,1}] kenhill@indirect.com kenhill@anvil5k.com "Point of view is worth 80 IQ points" -Alan Kay
следующий фpагмент (22)|пpедыдущий фpагмент (20)
From: Steve Talent See "Curves and Surfaces for CAGD" by Gerald Farin. The deCasteljau algoritim for evaluating blossoms may be used to subdivide Bezier curves. This code will work with rational Bezier curves -- set w=1 for ordinary Bezier curves. Here's some code: void subdivide( int degree, double *c, /* input control polygon */ double t, /* parameter value of split */ double *c1, /* "left" control polygon */ double *c2 /* "right" control polygon */ ) { int i, j; double tb[4]; /* use blossom to find interior points */ for ( i = 0; i <= degree; i++ ) { for ( j = 0; j < degree-i; j++ ) tb[j] = 0.0; for ( j = degree-i; j < degree; j++ ) tb[j] = t; decasblos(degree,c,tb,&c1[4*i]); } for ( i = 0; i <= degree; i++ ) { for ( j = 0; j < degree-i; j++ ) tb[j] = t; for ( j = degree-i; j < degree; j++ ) tb[j] = 1.0; decasblos(degree,c,tb,&c2[4*i]); } } void decasblos( int degree, double *c, /* input control polygon */ double *t, double *p /* output point */ ) /* Modified version of Farin's code. This one expects [wx, wy, wz, w] to be passed in and returns result in same format in *p. This is an extension of decas.c that handles blossoms. */ { int r,i; int j; double t1,t0; double coeffa[20*4]; /* an auxiliary array. Change dim. if too small*/ /* save input */ for(i=0; i<=degree; i++) { for ( j = 0; j < 4; j++ ) { coeffa[i*4+j]=c[i*4+j]; } } for (r=1; r<= degree; r++) { t0 = t[r-1]; t1 = 1.0 - t0; for (i=0; i<= degree - r; i++) { for ( j = 0; j < 4; j++ ) { coeffa[i*4+j]= t1* coeffa[i*4+j] + t0 * coeffa[(i+1)*4+j] ; } } } for ( j = 0; j < 4; j++ ) { p[j] = coeffa[j]; } }

Всего 21 фpагмент(а/ов) |пpедыдущий фpагмент (21)

Если вы хотите дополнить FAQ - пожалуйста пишите.

design/collection/some content by Frog,
DEMO DESIGN FAQ (C) Realm Of Illusion 1994-2000,
При перепечатке материалов этой страницы пожалуйста ссылайтесь на источник: "DEMO.DESIGN FAQ, http://www.enlight.ru/demo/faq".