Помогите найти библиотеку (Построение B-сплайнов)

← →
Dimon_St (2004-09-27 16:43) [0]

Здравствуйте. Нужна библиотека (на Delphi) для вычисления коэффицентов квадратичного и кубического B-сплайна.

← →
Рамиль © (2004-09-27 17:11) [1]

Кубический.
unit Splains;

interface
uses SysUtils, Windows, Messages, Dialogs;

type
TPolynom = class(TObject)
private
FKof: array[0..3] of extended;
function GetKof(index: integer): extended;
procedure SetKof(index: integer; const Value: extended);
public
constructor Create;
function Value(x: extended): extended;
property Kof[index: integer]: extended read GetKof write SetKof;
end;

TSplain = class(TObject)
private
FSplain: array of TPolynom;
Fx: array of extended ;
public
constructor Create(x, f: array of extended; xprob: extended);
function Value(x: extended): extended;
end;

implementation
uses Progonka;

{ TSplain }

constructor TSplain.Create(x, f: array of extended; Xprob: extended);
var
i, n: integer;
lam, mu: extended;
a, b, c, g, m, h: array of extended;
Polynom: TPolynom;
begin
n:=Length(x);
SetLength(Fx, n);
for i:=0 to n-1 do
Fx[i]:=x[i];
SetLength(m, N);
SetLength(a, N);
SetLength(b, N);
SetLength(c, N);
SetLength(g, N);
SetLength(h, N-1);
SetLength(FSplain, n-1);
for i:=1 to n-2 do
begin
a[i]:=2;
h[i]:=x[i+1]-x[i];
c[i]:=h[i]/(x[i+1]-x[i-1]);
b[i]:=1-c[i];
g[i]:=3*(c[i]*(f[i]-f[i-1])/(x[i]-x[i-1])+
b[i]*(f[i+1]-f[i])/(x[i+1]-x[i]));
end;
h[0]:=x[1]-x[0];
c[0]:=0; b[0]:=0; c[n-1]:=0; b[n-1]:=0;
lam:=(x[1]-x[0])/(x[2]-x[0]);
mu:=1-lam;
a[0]:=1; a[n-1]:=1;
g[0]:=(1+mu)*(f[1]-f[0])/(x[1]-x[0])-mu*(f[2]-f[1])/(x[2]-x[1]);
lam:=(x[n-1]-x[n-2])/(x[n-1]-x[n-3]);
g[n-1]:=-lam*(f[n-2]-f[n-3])/(x[n-2]-x[n-3])+(1+lam)*(f[n-1]-f[n-2])/(x[n-1]-x[n-2]);
SolutionProgonka(a, b, c, g, m);
for i:=0 to n-2 do
begin
Polynom:=TPolynom.Create;
Polynom.Kof[3]:= (2*f[i]-2*f[i+1]+h[i]*m[i]+h[i]*m[i+1])/((x[i+1]-x[i])*
(x[i+1]-x[i])*(x[i+1]-x[i]));
Polynom.Kof[2]:= (-3*f[i]*(x[i+1]+x[i])+3*f[i+1]*(x[i+1]+x[i])-
h[i]*m[i]*(2*x[i+1]+x[i])-h[i]*m[i+1]*(2*x[i]+
x[i+1]))/((x[i+1]-x[i])*(x[i+1]-x[i])*(x[i+1]-x[i]));
Polynom.Kof[1]:= (6*x[i]*x[i+1]*f[i]-6*x[i]*x[i+1]*f[i+1]+
h[i]*m[i]*(x[i+1]*x[i+1]+2*x[i]*x[i+1])+
h[i]*m[i+1]*(2*x[i]*x[i+1]+x[i]*x[i]))/((x[i+1]-
x[i])*(x[i+1]-x[i])*(x[i+1]-x[i]));
Polynom.Kof[0]:= ((x[i+1]*x[i+1]*x[i+1]-3*x[i]*x[i+1]*x[i+1])*f[i]+f[i+1]*
(3*x[i]*x[i]*x[i+1]-x[i]*x[i]*x[i])-h[i]*m[i]*x[i]*
x[i+1]*x[i+1]-h[i]*m[i+1]*x[i]*x[i]*x[i+1])/((x[i+1]-x[i])
*(x[i+1]-x[i])*(x[i+1]-x[i]));
FSplain[i]:=Polynom;
end;
Finalize(a);
Finalize(b);
Finalize(c);
Finalize(g);
Finalize(h);
Finalize(m);
end;

function TSplain.Value(x: extended): extended;
var i: integer;
begin
if (x<Fx[0]) or (x>Fx[High(Fx)]) then
ShowMessage("Идиот!!! Х вышел за пределы по которым ты строил сплайн.");
for i:=1 to high(Fx) do if x <= Fx[i] then break;
result:=FSplain[i-1].Value(x);
end;

{ TPolynom }

constructor TPolynom.Create;
var
i: integer;
begin
for i:=0 to 3 do FKof[i]:=0;
end;

function TPolynom.GetKof(index: integer): extended;
begin
result:=FKof[index];
end;

procedure TPolynom.SetKof(index: integer; const Value: extended);
begin
FKof[index]:=Value;
end;

function TPolynom.Value(x: extended): extended;
var
i: integer;
stepen: extended;
begin
result:=FKof[0];
stepen:=1;
for i:=1 to 3 do
begin
stepen:=stepen*x;
result:=result+FKof[i]*stepen;
end;
end;

end.

← →
Рамиль © (2004-09-27 17:11) [2]

unit Progonka; interface uses SysUtils, Windows, Messages; procedure SolutionProgonka(a, b, c, g:array of extended; var z: array of extended); implementation procedure SolutionProgonka(a, b, c, g: array of extended; var z: array of extended); var n, i, j: integer; u, v: array of extended; begin n:=Length(Z); SetLength(u, n+1); SetLength(v, n+1); v[0]:=0; u[0]:=0; for i:=1 to n do begin v[i]:=-b[i-1]/(a[i-1]+c[i-1]*v[i-1]); u[i]:=(g[i-1]-c[i-1]*u[i-1])/(a[i-1]+c[i-1]*v[i-1]); end; z[n-1]:=u[n]; for i:=n-1 downto 1 do begin z[i-1]:=v[i]*z[i]+u[i]; end; end; end.

← →
Рамиль © (2004-09-27 17:13) [3]

За качество не ручаюсь, отрывок из не моей курсовой на третьем курсе.

← →
MBo © (2004-09-27 17:15) [4]

>Рамиль
Это у тебя не B-сплайны, а обычные (полиномиальный базис).

>Dimon_St
http://astronomy.swin.edu.au/~pbourke/curves/spline/
да и сам поисковиком-то поработай...

← →
REA (2004-09-27 17:18) [5]

Кубический, за качество не ручаюсь - творчески поиздевался над кодом из инета:
(*************************************************************************
Ïðîöåäóðà Spline3BuildTable ñëóæèò äëÿ ïîñòðîéêè òàáëèöû êîýôôèöèåíòîâ
êóáè÷åñêîãî ñïëàéíà ïî çàäàííûì òî÷êàì è ãðàíè÷íûì óñëîâèÿì, íàêëàäûâàåìûì
íà ïðîèçâîäíûå.

Â äàëüíåéøåì ïîñòðåííàÿ òàáëèöà èñïîëüçóåòñÿ ôóíêöèåé Spline3Interpolate.

Ïàðàìåòðû:
N - ÷èñëî òî÷åê ìèíóñ 1.
DiffN - òèï ãðàíè÷íîãî óñëîâèÿ. "1" ñîîòâåòñòâóåò ãðàíè÷íûì óñëîâèÿì
íàêëàäûâàåìûì íà ïåðâûå ïðîèçâîäíûå, "2" - íà âòîðûå.
xs - ìàññèâ àáñöèññ îïîðíûõ òî÷åê ñ íîìåðàìè îò 0 äî N.
ys - ìàññèâ îðäèíàò îïîðíûõ òî÷åê ñ íîìåðàìè îò 0 äî N.
BoundL - ëåâîå ãðàíè÷íîå óñëîâèå. Åñëè DiffN ðàâíî 1, òî ïåðâàÿ ïðîèç-
âîäíàÿ íà ëåâîé ãðàíèöå ðàâíà BoundL, èíà÷å BoundL ðàâíà
âòîðàÿ.
BoundR - àíàëîãè÷íî BoundL

Âûõîäíûå çíà÷åíèÿ:
ctbl - â ýòîò ìàññèâ ïîìåùàåòñÿ òàáëèöà êîýôôèöèåíòîâ ñïëàéíà.
Ýòà òàáëèöà èñïîëüçóåòñÿ ôóíêöèåé Spline3Interpolate.
*************************************************************************)
procedure Spline3BuildTable(DiffN: Integer; x, y: TDoubleArray;
BoundL, BoundR: Double; var ctbl: TSplineArray);
var
C: Boolean;
N, G: Integer;
Tmp: Double;
nxm1: Integer;
I, J: Integer;
DX, DXJ, DYJ: Double;
DXJP1, DYJP1: Double;
DXP, YPPA, YPPB: Double;
PJ, b1, b2, b3, b4: Double;
begin
N := Length(X)-1;
//Body
g := (n+1) div 2;
repeat
i := g;
repeat
j := i-g;
c := True;
repeat
if x[j]<=x[j+g] then
c := False
else
begin
Tmp := x[j];
x[j] := x[j+g];
x[j+g] := Tmp;
Tmp := y[j];
y[j] := y[j+g];
y[j+g] := Tmp;
end;
j := j-1;
until not ((j>=0) and C);
i := i+1;
until not (i<=n);
g := g div 2;
until not (g>0);
For i := 0 To 4 Do
If High(ctbl[i])<>N Then
SetLength(ctbl[i], N+1);
N := N+1;
if DiffN=1 then
begin
b1 := 1;
b2 := 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-BoundL);
b3 := 1;
b4 := 6/(x[n-1]-x[n-2])*(BoundR-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
end else
begin
b1 := 0;
b2 := 2*BoundL;
b3 := 0;
b4 := 2*BoundR;
end;
nxm1 := n-1;
if n>=2 then
begin
if n>2 then
begin
dxj := x[1]-x[0];
dyj := y[1]-y[0];
j := 2;
while j<=nxm1 do
begin
dxjp1 := x[j]-x[j-1];
dyjp1 := y[j]-y[j-1];
dxp := dxj+dxjp1;
ctbl[1,j-1] := dxjp1/dxp;
ctbl[2,j-1] := 1-ctbl[1,j-1];
ctbl[3,j-1] := 6*(dyjp1/dxjp1-dyj/dxj)/dxp;
dxj := dxjp1;
dyj := dyjp1;
j := j+1;
end;
end;
ctbl[1,0] := -b1/2;
ctbl[2,0] := b2/2;
if n<>2 then
begin
j := 2;
while j<=nxm1 do
begin
pj := ctbl[2,j-1]*ctbl[1,j-2]+2;
ctbl[1,j-1] := -ctbl[1,j-1]/pj;
ctbl[2,j-1] := (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,J-2])/pj;
j := j+1;
end;
end;
yppb := (b4-b3*ctbl[2,nxm1-1])/(b3*ctbl[1,nxm1-1]+2);
i := 1;
while i<=nxm1 do
begin
j := n-i;
yppa := ctbl[1,j-1]*yppb+ctbl[2,j-1];
dx := x[j]-x[j-1];
ctbl[3,j-1] := (yppb-yppa)/dx/6;
ctbl[2,j-1] := yppa/2;
ctbl[1,j-1] := (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]*dx)*dx;
yppb := yppa;
i := i+1;
end;
i:=1;
while i<=n do
begin
ctbl[0,i-1] := y[i-1];
ctbl[4,i-1] := x[i-1];
Inc(i);
end;
end;
end;

← →
REA (2004-09-27 17:20) [6]

(*************************************************************************
Âû÷èñëåíèå çíà÷åíèÿ ñïëàéíà â òî÷êå ïî òàáëèöå êîýôôèöèåíòîâ

Ôóíêöèÿ Spline3Interpolate ïî ïîñòðîåííîé ïðîöåäóðîé Spline3BuildTable
òàáëèöå êîýôôèöèåíòîâ âû÷èñëÿåò çíà÷åíèå èíòåðïîëèðóåìîé ôóíêöèè â
çàäàííîé òî÷êå.

Ïàðàìåòðû:
N - ÷èñëî òî÷åê ìèíóñ 1
C - òàáëèöà êîýôôèöèåíòîâ, ïîñòðîåííàÿ ïðîöåäóðîé
Spline3BuildTable
X - òî÷êà, â êîòîðîé âåäåì ðàñ÷åò

Ðåçóëüòàò:
çíà÷åíèå êóáè÷åñêîãî ñïëàéíà, çàäàííîãî òàáëèöåé C, â òî÷êå X
*************************************************************************)
function Spline3Interpolate(N: Integer; const c: TSplineArray; X: Double):Double;
var
I: Integer;
L: Integer;
Half: Integer;
First: Integer;
Middle: Integer;
begin
L := N;
First := 0;
while L>0 do
begin
Half := L div 2;
Middle := First+Half;
if C[4,Middle]<X then
begin
First := Middle+1;
L := L-Half-1;
end else L := Half;
end;
I := First-1;
if I<0 then I := 0;
Result := c[0,I]+(X-c[4,i])*(C[1,I]+(X-c[4,i])*(C[2,I]+C[3,i]*(X-c[4,i])));
end;

Типы забыл еще:
TDoubleArray = Array Of Double;
TSplineArray = Array[0..4] Of TDoubleArray;

← →
Рамиль © (2004-09-27 17:23) [7]

MBo © (27.09.04 17:15) [4]
Сам спалайнами не занимался.. Странно только, вроде автор курсовой делал B-сплайны :o)

← →
Dimon_St (2004-09-27 17:33) [8]

Спасибо.

← →
Dimon_St (2004-09-27 18:21) [9]

А вообще, есть ли готовые Дельфийские библиотеки по оптимизации, численным методам,дифф.у-ния и т.д. Я не могу поверить, что нет готовых математических библиотек для Delphi.

Правильно, не верь.
Все есть:))

← →
Dimon_St (2004-09-27 18:28) [11]

Вот только где?

Помогите найти библиотеку (Построение B-сплайнов) Найти похожие ветки