Журнал о программированнии на языках Blitz3D, BlitzPlus, BlitzMax

Кубическая интерполяция

Материал из Blitz Et Cetera

Перейти к: навигация, поиск

Методы интерполяции

Интерполяцией называется процесс нахождения функции (в классическом случае - типа y = f(x)), график которой проходит через заданные точки (все они должны иметь различные координаты x). При этом, используют формулы исходя из целей и ресурсов: одни методы позволяют получить идеально плавную кривую, но требуют большого количества вычислений, которое значительно возрастает при увеличении количества точек, другие работают значительно быстрее, но результат более груб. Один из комплексных методов - метод Лагранжа. Пусть дан набор из n точек (xi, yi). Интерполирующая функция Лагранжа выглядит так:

Как видно, если подставить в нее одно из значений xi, то все слагаемые формулы будут равны 0, кроме одного, в котором все множители сократятся и останется yi. Вот программа, демонстрирующая метод Лагранжа (кстати, виден еще один недостаток этого метода - график слишком "колеблется"):

SeedRnd MilliSecs()


Const i=10
Dim c#(i,1)

Graphics 800,600,16

Color 255,0,0
For j=1 To i
 ;Проверка на совпадение x-координат
 Repeat
  ;Шаг, чтобы точки отстояли друг от друга
  x#=5*Rand(1,159)
  ;Цикл по всем предыдущим точкам
  For k=1 To j-1
   ;Если координаты совпадают - цикл прерывается
   If c#(k,0)=x# Then Exit
  Next
  ;Если цикл не был прерван, то k=(j-1)+1=j, иначе задаем новое значение x
  If k=j Then Exit
 Forever
 c#(j,0)=x#
 c#(j,1)=Rnd(200,400)
 Oval c#(j,0)-5,c#(j,1)-5,11,11
Next

Color 255,255,255
For x#=0 To 799
 ;Вычисление значения формулы
 y#=0
 For j=1 To i
  m#=c#(j,1)
  For k=1 To i
   If k<>j Then m#=m#*(x#-c#(k,0))/(c#(j,0)-c#(k,0))
  Next
  y#=y#+m#
 Next

 ;Корректное построение части графика (иначе будут разрывы)
 If x#>0 Then
  y1#=y#
  If y1#<-3 Then y1#=-3 ElseIf y1#>602 Then y1#=602
  If y2#<-3 Then y2#=-3 ElseIf y2#>602 Then y2#=602
  If y2#<y1# Then z#=y1#:y1#=y2#:y2#=z#
  For yy#=y1# To y2#
   Oval x#-2,yy#-2,5,5
  Next
 End If
 y2#=y#
Next

WaitKey

Чтобы упростить вычисления, ось x разбивают на отрезки так, чтобы каждый из них был ограничен x-координатами двух точек и не содержал в себе x-координат других точек. Затем для каждого отрезка находят интерполирующую функцию. Простейший вариант показан на рисунке:

Полученная функция обладает минимумом требований: она соединяет заданные точки и только. Этого бывает достаточно, но, как видите, такая функция обладает явным недостатком - острыми углами. В этом случае, интерполирующая функция - линейная (вида y = kx + b),

Трехмерная интерполирующая поверхность задается по схожим принципам. Хороший пример такой интерполяции - terrain в Blitz. Задаются только высоты точек на пересечении сетки, высоты всех остальных точек вычисляются по формуле (каждый квадрат этой сетки представлен двумя треугольниками) - в итоге из набора точек мы получаем поверхность.

Кубическая интерполяция

Но, вернемся к двумерной функции: для того, чтобы избавиться от острых углов, смежные интерполирующие функции должны иметь в общих точках одинаковую производную (она определяет касательную к функции в этой точке). В итоге получается система из четырех уравнений:

  • f(x1) = y1
  • f(x2) = y2
  • f'(x1) = y1'
  • f'(x2) = y2'

Чтобы функция могла быть определена для любых возможных значений, она должна иметь четыре параметра (по числу уравнений). Простейший вариант - кубический многочлен: y = ax3 + bx2 + cx + d. Путем введения вспомогательной формулы x' = x + x1, переменной x3=x2-x1 и преобразований, получаем:

  • d = y1
  • c = y1'
  • b = (3y2 - y2'x3 - 2cx3 - 3d) / x32
  • a = (y2' - 2bx3 - c) / 3x32

Производные можно выбирать следующим образом: проведем прямую через две точки, смежные с данной; касательная к кривой в этой точке будет параллельна этой прямой и, соответственно, иметь тот же коэффициент, равный производной. Коэффициент прямой - это k в уравнении прямой y = kx + b. Для прямой, проходящей через две заданные точки он будет равен (y2 - y1) / (x2 - x1). Естественно, x2 <> x1.

Одно условие: точки должны быть отсортированы по возрастанию координаты x. А вот программа построения интерполяционной кривой:

SeedRnd MilliSecs()


Dim ptx#(81)
Dim pty#(81)

Graphics 800,600

;Создается цепь из точек
x#=0
y#=300
Color 255,0,0
Repeat
 q=q+1
 ptx#(q)=x#
 pty#(q)=y#
 Oval x#-4,y#-4,9,9
 x#=x#+Rnd(30,100)
 y#=y#+Rnd(-100,100)
Until x#>=800

Color 255,255,255
;Цикл по всем отрезкам (для крайней правой точки отрезка нет)
For n=1 To q-1
 d#=pty#(n)
 If n=1 Then
  ;Если берется начальный отрезок, то производная равна 0 (т.к. смежная точка
  ; слева, необходимая для определения коэффициента отсутствует
  c#=0
 Else
  ;Вычисление коэффициента, равного производной N1 по формуле
  c#=(pty#(n+1)-pty#(n-1))/(ptx#(n+1)-ptx#(n-1))
 End If
 ;Аналогично вычисляется производная N2
 If n=q Then
  dy2#=0
 Else
  dy2#=(pty#(n+2)-pty#(n))/(ptx#(n+2)-ptx#(n))
 End If
 ;Вычисление остальных коэффициентов многочлена
 x3#=ptx#(n+1)-ptx#(n)
 xx3#=x3#*x3#
 b#=(3*pty#(n+1)-dy2#*x3#-2*c#*x3#-3*d#)/xx3#
 a#=(dy2#-2*b#*x3#-c#)/(3*xx3#)
 ;Построение отрезка кривой
 For x#=0 To x3#
  xx#=x#*x#
  y#=a#*xx#*x#+b#*xx#+c#*x#+d#
  x1#=x#+ptx#(n)  
  If x1#>0 Then
   y1#=y#
   If y1#<-3 Then y1#=-3 ElseIf y1#>602 Then y1#=602
   If y2#<-3 Then y2#=-3 ElseIf y2#>602 Then y2#=602
   If y2#<y1# Then z#=y1#:y1#=y2#:y2#=z#
   For yy#=y1# To y2#
    Rect x1#-1,yy#-1,3,3
   Next
  End If
  y2#=y#
 Next
Next

WaitKey

Интерполяция применяется не только для построения визуальных кривых и изогнутых поверхностей, но и для описания траекторий, то есть создания функций, где ось OX определяет время. В программе, приведенной ниже, объект, изображенный в виде круга, должен достигнуть заданных точек в определенное время и вернуться в начало. Для этого, зададим две интерполирующие функции x = f1(t) и y = f2(t) по принципу, изложенному выше:

SeedRnd MilliSecs()


Const q=10

;Так как координаты x и y обрабатываются одинаково, создаются массивы, где
; хранятся эти координаты и коэффициенты для функций на текущем отрезке времени
Dim ptc#(q+2,1)
Dim a#(1)
Dim b#(1)
Dim c#(1)
Dim d#(1)
;Массив для времени, в которое объект должен посетить заданную точку
Dim tim(q+2)
;Массив координат объекта
Dim oc#(1)

Graphics 800,600
SetFont LoadFont ("Arial",16)

x#=0
y#=300
Color 255,0,0
For n=1 To q
 x#=Rnd(50,750)
 y#=Rnd(50,550)
 ptc#(n,0)=x#
 ptc#(n,1)=y#
 tim(n)=t
 Oval x#-4,y#-4,9,9
 Text x#,y#+4,.001*t+"s",True
 t=t+Rand(1000,3000)
Next

;Задаются значения параметров крайних точек для обеспечения цикличности
For nn=0 To 1
 ptc#(0,nn)=ptc#(q,nn)
 ptc#(q+1,nn)=ptc#(1,nn)
 ptc#(q+2,nn)=ptc#(2,nn)
Next
tim(q+1)=t
tim(q+2)=t+tim(2)
tim(0)=tim(q)-t

Color 255,255,255
;Сохраняем фон
i=CreateImage(800,600)
GrabImage i,0,0
SetBuffer BackBuffer()

;Cчетчик времени устанавливается за пределами интервала цикла, чтобы
; вычислить коэффициенты для начального отрезка времени, пройдя через условие
t=t+1
n=q+1
Repeat
 If t>tim(n+1) Then
  n=n+1
  ;Если номер узла вышел за пределы массива - возврат на узел 1, обнуление
  ; счетчика времени
  If n>q Then
   n=1
   ms=0
   tbeg=MilliSecs()
  End If
  For nn=0 To 1
   d#(nn)=ptc#(n,nn)
   c#(nn)=(ptc#(n+1,nn)-ptc#(n-1,nn))/(tim(n+1)-tim(n-1))
   dy2#=(ptc#(n+2,nn)-ptc#(n,nn))/(tim(n+2)-tim(n))
   x3#=tim(n+1)-tim(n)
   xx3#=x3#*x3#
   b#(nn)=(3*ptc#(n+1,nn)-dy2#*x3#-2*c#(nn)*x3#-3*d#(nn))/xx3#
   a#(nn)=(dy2#-2*b#(nn)*x3#-c#(nn))/(3*xx3#)
  Next
 End If

 ;Вычисление координат объекта
 For nn=0 To 1
  v#=t-tim(n)
  vv#=v#*v#
  oc#(nn)=a#(nn)*vv#*v#+b#(nn)*vv#+c#(nn)*v#+d#(nn)
 Next

 ;Отображение фона, объекта и текущего времени
 DrawBlock i,0,0
 Oval oc#(0)-9,oc#(1)-9,19,19
 Text 0,0,"Time:"+(.001*t)+"s"
 Flip

 t=MilliSecs()-tbeg

Until KeyHit(1)

Автор: Матвей Меркулов (E-mail: MattMerkulov_sobaka_gmail.com, ICQ: 392-274-050)

Источник — «http://blitzetc.blitzmax.ru/index.php/%D0%9A%D1%83%D0%B1%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B0%D1%8F_%D0%B8%D0%BD%D1%82%D0%B5%D1%80%D0%BF%D0%BE%D0%BB%D1%8F%D1%86%D0%B8%D1%8F»

Другие

Друзья