Кубические сплайны. Вычислительная математика
Язык: Borland Pascal 7.0
Курс: Вычислительная математика. Кубические сплайны
В статье реализован алгоритм работы с кубическими сплайнами и построение 3D модели.
[1] Постановка задачи
Пусть задана сетка значений точек поверхностей некоторой фигуры. Используя кубические сплайны, получить некоторые значения, которые впоследствии будут использоваться для построения данной фигуры в 3D-пространстве. Упор делается на вычисление сплайнов.
[2] Данные
Константы N=20 Создаем новые типы MyType = Extended Тип представления вещественных чисел Vector = array[0..N] of MyType Используется при вычислении сплайнов Mass = array[0..N,0..N+1] of MyType Используется при вычислении сплайнов Переменные p : array[1..472,1..3] of MyType; выходной массив точек Alpha : array[0..N+1] of MyType; для метода прогонки Betha : array[0..N+1] of MyType; для метода прогонки S1,S2 : Mass; для вычислении сплайнов Q1,Q2,ff1,ff2 : Vector; для вычислении сплайнов FileName : String; Название файла с сеткой MyCircle,MyCar : Boolean; Кроме того используется ряд других переменных, имеющих второстепеное значение
[3] Структура программы
Структура программы строится по следующей иерархической системе процедур и
функций :
1 Solve - Процедура вычисления коэффициентов для сплайнов
Для сплайна 1
вход : S1, ff1 выход : Q1
Для сплайна 2
вход : S2, ff2 выход : Q2
2 OutputData - Процедура получения выходного массива точек
вход : Sp1, Sp2 выход : p
3 MainDraw - Процедура визуализации фигуры
Вспомогательные функции и процедуры
S1,S2 - функции для первого и второго сплайнов соответственно
[4] Описание алгоритма
Для каждого сплайна Sp1 и Sp2 используется метод прогонки, который решает
систему ленточных уравнений (трехдиагональная система линейных уравнений).
S1*Q1=ff1 и S2*Q2=ff2. В результате находим вектора Q1 и Q2.
Далее, используя функции Sp1 и Sp2, получаем конкретные значения точек на
заданных поверхностях.
Получив сплайны, составляем выходной массив массив точек, который делится
на несколько областей [1..472]
1..36 - точки окружности, принадлежащей верхней поверхности
37..72 - точки окружности, принадлежащей нижней поверхности
73..172 - 100 точек на ближней кромке верхней поверхности
173..272 - 100 точек на дальней кромке верхней поверхности
273..372 - 100 точек на ближней кромке нижней поверхности
373..472 - 100 точек на дальней кромке нижней поверхности
Используя процедуру MainDraw, выводим на экран 3D-изображение фигуры.
[5] Описание управляющих клавиш
Выход из прогрыммы - Esc
Посмотреть статистику - F1
Увеличить рамеры и угол обзора финуры - 1
Уменьшить рамеры и угол обзора финуры - 2
Поворот по часовой стрелке - клавиша ->
Поворот против часовой стрелки - клавиша <-
Останов движения фигуры - клавиша <Курсор-Вниз>
Плоское изображение фигуры - 0
Увеличение размеров фигуры - 1
Уменьшение размеров фигуры - 2
Сдвиг фигуры вверх - клавиша <Курсор-Вверх>
Сдвиг фигуры вниз - клавиша <Курсор-Вниз>
Редактирование данных - F2
Показать или спрятать окружности - +<+>
[6] Список литературы
1 А.А.Самарский, Е.С.Николаев "Методы решения сеточных уравнений";
2 С.Б.Стечкин, Ю.Н.Субботин "Сплайны в вычислительной математике";
3 Курс лекций по вычислительной математике.
Приложение (Текст программы)
{$N+}
Program Kursovik; {Вычислительная математика : Кубические сплайны}
Uses Crt,Graph,Graph256;
Const N=20;
Type MyType = Extended;
Vector = array[0..N] of MyType;
Mass = array[0..N,0..N+1] of MyType;
Var
p : array[1..472,1..3] of MyType; {выходной массив точек}
Alpha : array[0..N+1] of MyType;
Betha : array[0..N+1] of MyType;
S1,S2 : Mass;
Q1,Q2,ff1,ff2 : Vector;
h,x : MyType;
a,b : MyType;
ii,j,k,i : integer;
FileName : String;{Название файла с сеткой}
Error : Boolean;
Dr_x,Dr_y: integer;
page : integer;
MyCircle,MyCar : Boolean;
MyDelay: integer; {Задержка между перерисовками}
t : Integer;
stroka : string;
fil : Text;
Setka : Vector;
H2,z :MyType;
Gd,Gm,Gx,Gy,x1,y1,x2,y2 : integer;
_x,_h,_y,_rx,_ry,_rr : MyType;
Size : Integer;
Delta : integer; {Поворот вправо или влево}
Key : Char;
Function IntToStr(qwe: integer): string;
var resstr: string;
begin
str(qwe,resstr);IntToStr:=resstr
end;
Function FloatToStr(qwe: MyType): string;
var resstr: string;
begin
str(qwe:3:3,resstr); FloatToStr:=resstr
end;
Procedure Introduction;
begin
textcolor(black); ClrScr; TextColor(LightRed); gotoxy(24,7);
Writeln('К У Р С О В А Я Р А Б О Т А'); TextColor(White);
gotoxy(16,9); writeln(' по дисциплине : "Вычислительная математика"');
TextColor(Yellow); gotoxy(16,11);
writeln(' по теме : "Построение 3D изображения фигуры,');gotoxy(25,12);
writeln('используя кубические сплайны"');TextColor(White);gotoxy(25,15);
writeln('Студент :',' ':10,'Шаров Е.Н.');gotoxy(25,16);
writeln('Группа : ',' ':10,'ПА-97'); gotoxy(25,17);
writeln('Преподаватель : Лобанов С.В.');TextColor(DarkGray);
gotoxy(32,23);write('Рыбинск, 2000')
end;
Procedure LoadString (X, Y: Integer; len, fcol, bcol: Byte; Var stroka: String);
Var tempstr,mid: string; c: char;
Begin
tempstr:=stroka;mid := stroka + '_';SetColor (fcol);OutTextXY (X, Y, mid);
Repeat
c := ReadKey;
If (c = #8) And (Length (mid) > 1) Then
Begin
SetColor (bcol);OutTextXY (X, Y, mid);
Delete (mid, Length (mid) - 1, 2);mid := mid + '_'
End;
If ( (Ord (c) > 31) And (Ord (c) < 240) ) And (Length (mid) < len + 1) Then
Begin
SetColor (bcol);OutTextXY (X, Y, mid);
Delete (mid, Length (mid), 1); mid := mid + c + '_'
End;
SetColor (fcol);OutTextXY (X, Y, mid)
Until (c = #13) or (c=#27);
SetColor (bcol);OutTextXY (X, Y, mid);Delete (mid, Length (mid), 1);
SetColor (fcol);OutTextXY (X, Y, mid);stroka := mid;
if c=#27 then begin Error:=True; Stroka:=TempStr; end
End;
Function Sp1(x:MyType):MyType;
var xx1,xx2 : MyType;
begin
xx1:=A+h; i:=0;
While x>xx1 do begin xx1:=xx1+h;Inc(i) end;
Sp1:=q1[i]*(xx1-x)*(xx1-x)*(xx1-x)/(6*h) +
q1[i+1]*(x-xx1+h)*(x-xx1+h)*(x-xx1+h)/(6*h) +
ff1[i]*(xx1-x)/h +ff1[i+1]*(x-xx1+h)/h -
q1[i]*h*(xx1-x)/6 -q1[i+1]*h*(x-xx1+h)/6
end;
Function Sp2(x:MyType):MyType;
var xx1,xx2 : MyType;
begin
xx1:=A+h;i:=0;
While x>xx1 do begin xx1:=xx1+h;Inc(i) end;
Sp2:=q2[i]*(xx1-x)*(xx1-x)*(xx1-x)/(6*h)+
q2[i+1]*(x-xx1+h)*(x-xx1+h)*(x-xx1+h)/(6*h) +
ff2[i]*(xx1-x)/h+ff2[i+1]*(x-xx1+h)/h-
q2[i]*h*(xx1-x)/6-q2[i+1]*h*(x-xx1+h)/6
end;
Procedure Solve;
begin
Error:=False; h:=(B-A)/N;x:=A;
{$I-}
Assign(fil,FileName);Reset(fil);
{$I+}
If ioresult <>0
then begin Error:=True; exit end;
readln(fil,stroka);
for i:=0 to N do
begin {Считывание верхней поверхности}
readln(fil,stroka); val(stroka,x,j); setka[i]:=x
end;
For i:=0 to N do begin ff1[i]:=setka[i]; x:=x+h end;
readln(fil,stroka);
for i:=0 to N do
begin {Считывание верхней поверхности}
readln(fil,stroka);val(stroka,x,j);setka[i]:=x
end;
For i:=0 to N do begin ff2[i]:=setka[i]; x:=x+h end;
Close(fil);
{Составляем матрицу 1}
For i:=0 to N do For j:=0 to N+1 do S1[i,j]:=0.0;
For i:=1 to N-1 do
begin
S1[i,i-1]:=h/6; S1[i,i]:=4*H/3; S1[i,i+1]:=h/6;
S1[i,N+1]:=(ff1[i+1]-2*ff1[i]+ff1[i-1])/h
end;
S1[0,0]:=1; S1[0,N+1]:=0; S1[N,N]:=1; S1[1,N+1]:=0;
{Составляем матрицу 2}
For i:=0 to N do For j:=0 to N+1 do S2[i,j]:=0.0;
for i:=0 to n do begin alpha[i]:=0; betha[i]:=0 end;
For i:=1 to N-1 do
begin
S2[i,i-1]:=h/6; S2[i,i]:=4*H/3; S2[i,i+1]:=h/6;
S2[i,N+1]:=(ff2[i+1]-2*ff2[i]+ff2[i-1])/h
end;
S2[0,0]:=1; S2[0,N+1]:=0; S2[N,N]:=1; S2[1,N+1]:=0;
Alpha[1]:=S1[0,1]/S1[0,0]; Betha[1]:=S1[0,N+1];
for i:= 1 to N-2 do Alpha[i+1]:=S1[i,i+1]/(S1[i,i]-S1[I,I-1]*Alpha[i]);
for i:= 1 to N do
Betha[i+1]:=(S1[i,i-1]*Betha[i]+S1[i,N+1])/(S1[i,i]-S1[I,I-1]*Alpha[i]);
Q1[n]:=Betha[N+1];
for i:= N-1 downto 0 do Q1[i]:=Alpha[i+1]*Q1[i+1]+Betha[i+1];
Alpha[1]:=S2[0,1]/S2[0,0]; Betha[1]:=S2[0,N+1];
for i:= 1 to N-2 do Alpha[i+1]:=S2[i,i+1]/(S2[i,i]-S2[I,I-1]*Alpha[i]);
for i:= 1 to N do
Betha[i+1]:=(S2[i,i-1]*Betha[i]+S2[i,N+1])/(S2[i,i]-S2[I,I-1]*Alpha[i]);
Q2[n]:=Betha[N+1];
for i:= N-1 downto 0 do Q2[i]:=Alpha[i+1]*Q2[i+1]+Betha[i+1];
end;
Procedure OutputData;
begin
{Нахождение точек заданной фигуры}
for t:=1 to 36 do {36 точек на окружность}
begin {Верхняя окружность}
p[t,1]:=_rr*cos((t-1)*10*pi/180)+_rx;
p[t,2]:=_rr*sin((t-1)*10*pi/180)+_ry;
p[t,3]:=Sp1(p[t,1]);
{Нижняя окружность}
p[t+36,1]:=_rr*cos((t-1)*10*pi/180)+_rx;
p[t+36,2]:=_rr*sin((t-1)*10*pi/180)+_ry;
p[t+36,3]:=Sp2(p[t+36,1])
end;
_h:=0;
for t:=73 to 172 do
begin
{100 точек на ближней кромке верхней поверхности}
p[t,1]:=_h;p[t,2]:=0;p[t,3]:=Sp1(p[t,1]);
{100 точек на дальней кромке верхней поверхности}
p[t+100,1]:=_h;p[t+100,2]:=_y;p[t+100,3]:=SP1(p[t+100,1]);
{100 точек на ближней кромке нижней поверхности}
p[t+200,1]:=_h;p[t+200,2]:=0;p[t+200,3]:=Sp2(p[t+200,1]);
{100 точек на дальней кромке нижней поверхности}
p[t+300,1]:=_h;p[t+300,2]:=_y;p[t+300,3]:=Sp2(p[t+300,1]);
_h:=_h+_x/100
end;
{Смещение}
for t:=1 to 472 do
begin
p[t,1]:=p[t,1]- b/2;
p[t,2]:=p[t,2]- _ry/2;
end;
end;
Procedure DrawXY;
var LocStop : boolean;
Procedure Dr;
begin
ClearDevice; SetColor(Blue);Line(0,Gy div 2,Gx,Gy div 2); Line(0,0,0,Gy+300);
{Первый сплайн}
SetColor(LightCyan);z:=0;x1:=0;y1:=Round(sp1(z)*5);
While (zGetMaxY+150 then Gy:=GetMaxY+150
else begin
setactivepage(not(page)); page:=not(page);
Dr; setvisualpage(page);
end
end
else
if key=#72 then
begin
Dec(Gy,10);
if Gy15 then
begin
Dr_x:=Dr_x-4;Dr_y:=Dr_y-1;setactivepage(not(page));
page:=not(page);Dr; setvisualpage(page)
end
end else LocStop:=True
end
end;
CloseGraph;Gd:=EGA;Gm:=EGAlo;InitGraph(Gd,Gm,'C:\BP\BGI');
key:=#1
end;
{В данной процедуре происходит рисование и управление фигурой}
Procedure MainDraw;
var
xc,yc : integer; {центры экрана}
sx,sy : MyType; {Коэффициенты для пропорций}
x1,x2,y1,y2 : MyType;
Angle : Integer; {Угол аксонометрии}
ssx,ssy : integer;
Procedure OneStep1(os: integer);
begin
x1:=sx*p[os,1]*cos(Angle*pi/180); y1:=sy*p[os,1]*sin(Angle*pi/180);
x1:=x1-sx*p[os,2]*sin(Angle*pi/180); y1:=y1+sy*p[os,2]*cos(Angle*pi/180);
y1:=y1+p[os,3];
end;
Procedure OneStep2(os: integer);
begin
x2:=sx*p[os,1]*cos(Angle*pi/180); y2:=sy*p[os,1]*sin(Angle*pi/180);
x2:=x2-sx*p[os,2]*sin(Angle*pi/180); y2:=y2+sy*p[os,2]*cos(Angle*pi/180);
y2:=y2+p[os,3];
end;
Procedure Step(st:integer);
begin
OneStep2(st);line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
x1:=x2; y1:=y2
end;
Procedure Draw;
begin
ClearDevice;
{Вывод детали на экран}
if MyCircle
then
begin
setcolor(11); OneStep1(36); for t:=1 to 36 do Step(t);
setcolor(14); OneStep1(72); for t:=37 to 72 do Step(t);
end;
setcolor(11); OneStep1(73);for ii:=73 to 172 do Step(ii);
OneStep1(173);for ii:=173 to 272 do Step(ii);
setcolor(14); OneStep1(273);for ii:=273 to 372 do Step(ii);
OneStep1(373);for ii:=373 to 472 do Step(ii);
setcolor(11);OneStep2(273);OneStep1(373);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
{for сar's window}
if MyCar
then
begin
OneStep2(97);OneStep1(197);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
OneStep2(104);OneStep1(204);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
OneStep2(164);OneStep1(264);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
OneStep2(143);OneStep1(243);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
SetColor(14); OneStep2(368);OneStep1(468);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
OneStep2(278);OneStep1(378);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
end;
OneStep2(272);OneStep1(172);
line(xc+round(x1),yc-round(y1),xc+round(x2),yc-round(y2));
setvisualpage(page);setactivepage(not(page));page:=not(page);
end;
Procedure Help;
begin
page:=not(page);SetActivePage(page);SetVisualPage(page);ClearDevice;
SetColor(LightCyan);Rectangle(0,0,GetMaxX,GetMaxY);
OutTextXY(200,10,'Cells (N): '+IntToStr(n));
OutTextXY(200,20,'X : '+FloatToStr(_x));
OutTextXY(200,30,'Y : '+FloatToStr(_y));
if MyCircle
then
begin
OutTextXY(200,50,'Xcircle : '+FloatToStr(_rx));
OutTextXY(200,60,'Ycircle : '+FloatToStr(_ry));
OutTextXY(200,70,'Radius circle : '+FloatToStr(_rr));
end;
OutTextXY(200,90,'File with Cells : '+FileName);
OutTextXY(200,100,'Interval : [ '+FloatToStr(a)+', '+FloatToStr(B)+']');
key:=readkey;if key=#0 then key:=readkey;key:=#0
end;
{Езменение переменных величин}
Procedure Options;
var
j1,j2,j3 : integer;
j4 : MyType;
Procedure MenuOptions;
begin
page:=not(page);
setactivepage(page);
setvisualpage(page);
cleardevice; SetColor(LightCyan);
OuttextXY(230,50, 'A : '+FloatToStr(a));
OuttextXY(230,60, 'B : '+FloatToStr(b));
OuttextXY(230,70, 'X : '+FloatToStr(_x));
OuttextXY(230,80, 'Y : '+FloatToStr(_y));
OuttextXY(230,90, 'Xcircle : '+FloatToStr(_rx));
OuttextXY(230,100,'Ycircle : '+FloatToStr(_ry));
OuttextXY(230,110,'Rcircle : '+FloatToStr(_rr));
OutTextXY(230,120,'File with Cells : '+FileName);
OutTextXY(240,180,'Esc to exit'); setvisualpage(page);
end;
begin
j1:=1; j2:=j1; MenuOptions;
While key<>#27 do
begin
SetColor(White); OutTextXY(160,j1*10+40,'-->');
key:=readkey;
if key=#0 then
begin
key:=readkey;
Case key of
#72 : begin {Нажата <стрелка вниз>}
j1:=j1-1;if j1=0 then j1:=8
end;
#80 : begin {Нажата <стрелка вниз>}
j1:=j1+1; if j1=9 then j1:=1
end
end
end
else
if key=#13 then {Редактирование данных}
begin
case j1 of
1:begin
stroka:=FloatToStr(a);Error:=False;
LoadString(262,50,5,10,0,stroka);
if not error then
begin
val(stroka,j4,j3);
if (j3=0) and (j4a) then begin b:=j4; Solve;OutPutData end
end;
MenuOptions
end;
3:begin
stroka:=FloatToStr(_x);Error:=False;
LoadString(262,70,5,10,0,stroka);
if not error then
begin
val(stroka,j4,j3);
if (j3=0) and (j4<17)and (j4>0) then
begin _x:=j4;Solve;OutPutData end
end;
MenuOptions
end;
4:begin
stroka:=FloatToStr(_y);Error:=False;
LoadString(262,80,5,10,0,stroka);
if not error then
begin
val(stroka,j4,j3);
if (j3=0) and (j4>0) then
begin _y:=j4;Solve;OutPutData end
end;
MenuOptions
end;
5:begin
stroka:=FloatToStr(_rx); Error:=False;
LoadString(310,90,5,10,0,stroka);
if not error then
begin
val(stroka,j4,j3);
if j3=0 then
begin _rx:=j4; Solve; OutPutData end
end;
MenuOptions
end;
6:begin
stroka:=FloatToStr(_ry);Error:=False;
LoadString(310,100,5,10,0,stroka);
if not error then
begin
val(stroka,j4,j3);
if j3=0 then
begin _ry:=j4;Solve;OutPutData end
end;
MenuOptions
end;
7:begin
stroka:=FloatToStr(_rr);Error:=False;
LoadString(310,110,5,10,0,stroka);
if not error then
begin
val(stroka,j4,j3);
if j3=0 then begin _rr:=j4;Solve;OutPutData end
end;
MenuOptions
end;
8:begin
stroka:=FileName;Error:=False;
LoadString(374,120,12,10,0,FileName);
If not Error then Solve;
if not error then OutPutData else FileName:=stroka;
MenuOptions
end
end
end;
SetColor(Black); OutTextXY(160,j2*10+40,'-->');j2:=j1
end;
key:=#1
end;
begin
gd:=Ega;gm:=EGALo;InitGraph(gd,gm,'C:\BP\BGI');
if grOk<>0
then
begin WriteLn('Ошибка подключения графического режима.');Halt(0) end;
xc:=GetMaxX div 2; yc:=GetMaxY-90;MyCircle:=True;MyCar:=False;
Dr_x:=45;Dr_y:=10; Size:=5; Angle:=35; sx:=14; sy:=2;setvisualpage(page);
setactivepage(not(page)); page:=not(page); Draw; key:=#0;Delta:=2;
While Key <>#27 do
begin
if not keypressed and (Delta<>0)
then
begin
Angle:=Angle+Delta;
if Angle<=0 then Angle:=360-Angle
else
if Angle>=360 then Angle:=Angle-360;
Draw; Delay(MyDelay);
end
else
begin
key:=readkey;
if key=#0
then
begin
key:=readkey;
if key=#59 then Help else
if key=#60 then Options else
if key=#73 then
begin
Dec(MyDelay,10); if MyDelay<10 then MyDelay:=10
end;
if key=#75 then Delta:=2 else
if key=#77 then Delta:=-2 else
if key=#80 then Delta:=0;
if key=#81 then
begin
Inc(MyDelay,10); if MyDelay>400 then MyDelay:=400
end
end
else
begin
if (Key='1') and (Size<8) then
begin
sx:=sx+2.5;sy:=sy+1.5;Draw;Inc(Size);
end;
if (Key='2') and (Size>4) then
begin
sx:=sx-2.5;sy:=sy-1.5;Draw;Dec(Size)
end;
if key='0'
then
begin
DrawXY;setvisualpage(page);setactivepage(not(page));
page:=not(page);Draw
end;
if key='+'then
begin MyCircle:=not(MyCircle); Draw end;
if key='*' then
begin MyCar:=not(MyCar);Draw end
end
end
end;
CloseGraph
end;
Begin
Introduction;
{Строим сетку значений}
FileName:='setka3.in';a:=0.0; b:=14.0; Solve; MyDelay:=100;
{Плокое изображение фигуры}
_x:=B; _y:=7; _rx:=3; _ry:=_y/2; _rr:=1.8; OutPutData;
ReadKey; MainDraw
End.
Файл с данными setka3.in
[Top] 0 2 2.5 3 3.5 4 10 9.85 9.5 9.2 8.4 8 7.5 7 6.3 11.5 13.3 14.3 2.05 1 0 [Bottom] 0 -5 -4.8 -6.8 -8.2 -8.2 -6.8 -4.8 -5 -5 -5 -5 -5 -4.8 -6.8 -8.2 -8.2 -6.8 -4.8 -5 0
Файл с данными setka1.in
[Top] 0 4 6.7 8.4 9.3 9.9 10 9.85 9.55 9.2 8.7 8.2 7.6 6.95 6.2 5.35 4.3 3.3 2.05 1 0 [Bottom] 0 1.5 2.7 3.4 3.9 4.1 4.2 4.3 4.26 4.2 4.1 3.9 3.63 3.22 2.8 2.3 1.7 1.2 0.65 0.25 0
Apache — это кросплатформаенное программное обеспечение, относящееся к классу http-серверов. Поддерживается множеством операционных систем: Windows, Linux, MacOS и т.д. Одним из ключевых факторов в вопросе использования данного web-сервера является гибкость настройки и надежность выполнения операций. Apache включает в себя множество дополнительных модулей, позволяющих работать с различными базами данных, контролировать аутентификацию пользователей и т.д. |
Интересные материалы на сайте:
|
|