Metody interpolacji

Copyright (c) 2007 Krystian Komisarek

www.spider.dathox.com

Wstęp

Jest to prosty i krótki artykuł poświęcony temat interpolacji liniowej, kubicznej, splinów. No to ‘fajnie jak w kombajnie’ wiemy teraz wszystko:). Omówię matematykę związaną z krzywymi interpolacji oraz przedstawię przykładowy kod. W artykule pojawi się kilka ciekawych optymalizacji i usprawnień do standardowego podejścia. Oczywiście wszystko będzie tak wytłumaczone, że nie jest potrzebna żadna wiedza tajemna ani pochodzenie z rodziny Einsteina, aby to zrozumieć. Wymagania optymalne aby czytać artykuł to znajomość języka polskiego i podstawy matematyki na poziomie liceum. Jak zwykle postaram się odpowiadać na najważniejsze pytanie "dlaczego?". Na koniec artykułu dla wytrwałych zajmę się interpolacją wektorów i obrotów zapisanych kwaternionami.

Co to interpolacja oraz do czego służy?

Zanim przejdę do sedna sprawy musze odpowiedzieć  na pewne ważne pytanie nurtujące nie jedną osobę, która przypadkowo mogła wybrać mój artykuł nie wiedząc tak naprawdę, o czym będzie i czeka z niecierpliwością, aż wytłumaczę, czym tak naprawdę jest interpolacja. Synonim do słowa „interpolacja” to „przejście od do w czasie”. Tak wiem dość rozbudowany i dziwny synonim, ale krótszego nie wymyślę, dlatego korzysta się z uproszczonej nazwy tego mechanizmu :).

Jeśli nadal nic nie jest jasne to wytłumaczę w inny sposób. Wyobraźmy sobie 2 punkty A i B w dowolnej przestrzeni R1 (dla tych, którzy nie wiedzą, co mam namyśli to R oznacza wymiar w liczbach rzeczywistych np. R3 = 3D), R2, R3, a jak ktoś potrafi to nawet w R4 lol. Teraz wyobraźmy sobie, że mija czas w przedziale t <0,1> (ci, co sobie wyobrażali R4 pewnie mają teraz problem :). Teraz opiszmy problem, który ten artykuł ma za zadanie rozwiązać. Jak przedostać się z punktu A do B w czasie t <0,1>? Dokładniej. W jaki sposób określić pozycję punktu pośredniego pomiędzy A i B w czasie t gdzie t = 0 odpowiada punktowi A, a t = 1 to B? Odpowiedź prosta Interpolując te punkty po czasie t :). P(t) to funkcja która zwróci nam wyszukiwany wynik.

Pozostaje do rozwikłania jeszcze jedna zagadka, po co nam to ? Interpolacja może nam się przydać w wielu przypadkach pogranicza nas tylko wyobraźnia. Konkretniej mam na myśli:

• Ruch kamery - ustalamy linie punktów kamery oraz linie punktów oka. Interpolujemy.
• Animacja kości - kości zapisujemy najczęściej jako pozycja i rotacja klatek kluczowych.
• Skinning dwie kluczowe klatki można interpolować podobnie.
• Stretching wagi wierzchołków mają wpływ na rozciąganie siatki.
• Płynny ruch gracza po polach w planszowej grze.
• Generowanie brył obrotowych o dowolnych kształtach
• Płynna zmiana koloru - gradienty.

To tylko wybrane przypadki w których ja korzystam z interpolacji, tak naprawdę jest o wiele więcej dlatego zachęcam do skrupulatnego czytania artykułu.

Postać ogólna funkcji interpolującej

Mógłbym pominąć ten krok wypisać wzory na kolejne rodzaje interpolacji, które są dość proste, ale niewiele z nich da się wywnioskować bezpośrednio. Czytelnicy artykułu uznaliby, że mówię prawdę i tym samym używaliby ich jak czarnej magii. Nie zrobię tak jednak, ponieważ nie było to celem mojego artykułu a raczej rozwikłanie zagadki matematycznej. Aby cokolwiek zrobić trzeba poznać odrobinę teorii związanej z zagadnieniem w tym przypadku musimy dowiedzieć się jak wygląda ogólna postać funkcji interpolacji, a wtedy stanie się jasność. Podam teraz szablonowy wzór interpolacji:

Postać ogólna funkcji interpolującej.

P(t) = A*(1-g(t)) +B*(g(t))

P(t) - Szukany punkt interpolowany w czasie t.
g(t) - Funkcja interpolująca (wygładzająca)
A, B - Punkty, pomiędzy którymi interpolujemy

Dlaczego tak jest? Wypada wytłumaczyć, aby to zrobić przekształćmy wzór do odrobinę innej formy:

Geometryczna forma wzoru na inteprolację.

P(t) = (B-A)*g(t) + A

P(t) - Szukany punkt interpolowany w czasie t.
B-A - wektor przesunięcia z punktu A do B
g(t) - Funkcja interpolująca, skaluje nam wektor przesunięcia.
A - Przenosi wynik interpolacji o wektor A.

Naturalnie to właśnie funkcja g(t) wpływa na sposób interpolacji pomiędzy punktami może to być dowolna funkcja spełniająca następujące warunki:

• g(0) = 0
• g(1) = 1
• zbiór wartości funkcji na przedziale <0,1> należy do przedziału <0,1>.
• ciągła (nie jest konieczne jeśli nie chcemy otrzymać płynnego interpolowania).

Interpolacja liniowa (linear)

Najprostszy sposób interpolacji, jaki można sobie wyobrazić to przenoszenie punktu po linii prostej z równomiernym rozkładem. Wytłumaczę to w bardzo prymitywny sposób, jaki znam. Interpolacja nazywa się liniowa, ponieważ do funkcja interpolująca g(t) jest liniowa zapiszmy ją w następujący sposób wraz z warunkami:

g(t) = a*t + b
g(0) = 0 = a*0 + b = 0 => b = 0
g(1) = 1 =a*1 + b = 0 => a = 1

Mamy więc naszą szukaną funkcje interpolującą w pełnej okazałości. Co z pewnością zachwyca nas prostotą :). Wraz z wykresem.

g(t) = t

Interpolacja liniowa wyraża się więc wzorem:

P(t) = A*(1-t) +B*(t)

Tutaj wprowadzimy pewną modyfikację wzoru zauważmy, że korzystając z równania zapisanego geometrycznie obliczeń będzie mniej, ponieważ pozbywamy się jednego mnożenia na rzecz odejmowania wektorów (to taka mała optymalizacja).

Wzór na interpolacje liniową gdzie A B to punkty kontrolne t czas na przedziale <0,1>

P_linear(t) = (B-A)*t + A

Na koniec kod naszej funkcji zajmującej się interpolacją pisany w Delphi:

Liniowy sposób interpolacji jest bardzo prosty i szybki ma jednak jedną wadę interpolacja wykonuje się ciągłe z tą samą prędkością, tak więc przejście z punktu do drugiego punktu odbywa się ostrym startem i ostrym finiszem. W niektórych przypadkach jednak to podejście się bardzo przydaje.

Interpolacja kubiczna/sześcienna (cubic)

Drugi rodzaj interpolacji niweluje największy problem poprzedniego a więc ostry start i finisz. Tym razem startując z punktu 0 wraz ze wzrostem t punkt powoli przyspieszy a w momencie osiągnięcia t=0.5 będzie miał największą prędkość później znowu zwolni aż w punkcie 1 będzie miał prędkość 0. Tak jak poprzednim razem interpolacja otrzymała swoją nazwę od sposobu, w jakim powstała. Interpolacja  kubiczna - sześcienna powstała z wielomianu trzeciego stopnia. Sprawdźmy, więc i tym razem warunki początkowe oraz wzór ogólny wielomianu:

g(t) = a*t³ + b*t² + c*t + d

Po nałożeniu wymaganych warunków
g(0) = 0 = 0+ 0 + 0 + d => d = 0
g(1) = 1 = a + b + c

Z tych warunków wzoru nie wyprowadzimy. Przyjrzyjmy się wykresowi funkcji, jaką chcielibyśmy uzyskać. Wiemy także, że pochodna funkcji g(t) w punktach krańcowych musi się równać 0. Dla niewtajemniczonych w interpretacji geometrycznej wartość pochodnej w punkcie jest współczynnikiem kierunkowym stycznej do wykresu w tym punkcie. My chcemy, aby te współczynniki w punkcie 0 i 1 były równe 0 wtedy na krańcowych wartościach zmiana prędkości będzie najmniejsza.

Obliczmy więc pochodną:
g'(t) = 3*a*t² + 2*b*t + c
g'(0) = 0 = 3*a*0 + 2*b*0 + c => c = 0
g'(1) = 0 = 3*a + 2*b + 0 => a = -2*b/3

Wróćmy do:
g(1) = 1 = a + b + c = -2*b/3 + b = 1 => b = 3 i a = -2

Kolejny raz sukces. Mamy więc naszą szukaną funkcje interpolacyjną kubiczną .
g(t) = -2*t3 + 3*t2 = (3 - 2*t)*t2

Dla ciekawskich dodam, że wyprowadzona funkcja ma nawet swojego imiennika jest to jedna z czterech  podstawowych funkcji mieszania Hermita używanych w interpolacji. Czasem nazywanych też funkcjami wtapiania "blending functions". Mało tego drugą funkcję można uzyskać z zależności g2(t) := 1-g(t). Pozostałe dwie funkcje poznamy w dalszej części artykułu.

Interpolacja kubiczna wyraża się, więc wzorem:

P_cubic(t) = A*g(t) + B*(1-g(t))
P_cubic(t) = (B-A)*(3 - 2*t)*t² + A

Na koniec kod naszej funkcji zajmującej się interpolacją kubiczną:

Zauważmy, że w podobny sposób możemy wyprowadzić dowolną funkcję interpolująca pod warunkiem ze znamy jej własności. Warto, więc dokładnie zrozumieć teorię by można dzięki temu stworzyć na przykład funkcję hybrydową do interpolacji liniowej i kubicznej gdzie punkt będzie powoli ruszał z miejsca i gwałtownie się zatrzymywał. Nie wyprowadzę tutaj wzoru, ale radzę poćwiczyć by zyskać trochę praktyki.

Interpolacja splejnów (Cubic Hermite Spline)

To już ostatnia metoda interpolacji, jaką mam zamiar zaprezentować. Jest odrobinę trudniejsza tym samym uważam, że jest najciekawsza. Podobnie jak poprzednio, tym razem postaramy się rozwinąć odrobinę metodę sześcienną. Przedtem mieliśmy ustaloną prędkość początkową 0 oraz końcową 0. Chcielibyśmy tym razem mieć wpływ na prędkość początkowa oraz końcową punktów, co da nam największą płynność w przechodzeniu pomiędzy kolejnymi punktami kontrolnymi. Nadal będziemy posługiwać się wielomianem trzeciego stopnia tyle, że pochodne będziemy mieli ustalone Wprowadzimy, więc:

Av - Oznacza prędkość początkową w punkcie A.
Bv - Oznacza prędkość końcową w punkcie B.

Zauważmy się że prędkość wypływa na nową pozycję dlatego postaramy się przesunąć wynik interpolacji kubicznej o odpowiednią wartość co prawdopodobnie;) da nam oczekiwany wynik. Szukamy takiego równania:

P_spline(t) = P_cubic(t) + P_vel(t)

Pvel(t) - będzie naszym przesunięciem zależnym od prędkości w punktach krańcowych w czasie t. Kolejny raz zastanówmy się, jakie własności ma nasza szukana funkcja, aby otrzymać to, czego chcemy. Prędkość punktów ma pewien wpływ na pozycje końcową. Przez wpływ mam namyśli przesuniecie pozycji punktu względem oryginalnej interpolacji kubicznej. Wpływ Av na pozycje punktu musi być zerowy w miejscu t = 0 szybko wzrastać i spadać powoli zmierzając do argumentu 1 jak na rysunku po prawo. Podobne własności ma Bv. Szukana funkcja powinna mieć takową postać:
P_vel(t) = Av*(m(t)) + Bv*(n(t))

Podobnie jak w poprzednim przypadku korzystając z wiadomości o stycznych i miejscach zerowych wielomianu możemy wyprowadzić wzory na m(t) i n(t) obliczenia jednak pominę i pozostawię do samodzielnego sprawdzenia się. Wynik naszych działań to kolejne funkcje Hermita:

g(t) = -2*t3 + 3*t2
h(t) = 2*t3 - 3*t2 +1
m(t) = t3 - 2t2 + t
n(t) = t3 - t

Na rysunku obok przedstawione są wykresy wszystkich funkcji mieszania.

W rezultacie po podstawieniu i paru przekształceniach optymalizacyjnych nasza funkcja interpolacji prędkości jest odrobinę nieczytelna, ale najważniejsze, że wiemy, dlaczego właśnie tak wygląda:

P_vel(t) = Av*(t*(t-1)² ) + Bv*(t2*(t-1)) = ((Av+Bv)*t - Av)*t*(t-1)
P_spline(t) = (B-A)*(3 - 2*t)*t² + ((Av+Bv)*t - Av)*t*(t-1)

Jest jeszcze druga wersja tego wzoru, którą warto używać w przypadku, gdy współczynniki funkcji Hermita mamy wcześniej obliczone wygląda ona następująco.

P_spline(t) = A*g(t) + B*(1-g(t)) + Av*m(t) + Bv*(n(t))

Obliczanie stycznej w punkcie (Spline Tangent)

Musimy sobie zdawać ze interpolacja nie ogranicza się tylko do dwóch punktów zyskuje na przydatności, gdy stworzymy n punktów kontrolnych, które w odpowiedniej kolejności chcemy przechodzić wraz z mijającym czasem. Pojawia się jeszcze jedna kwestia do rozwiązania, a mianowicie nie wiemy, w jaki sposób liczyć wektory styczne dla wierzchołków. W tym przypadku także z pomocą przychodzi nam matematyka, istnieje wiele metod określenia wartości stycznej. Najprostszy i moim zdaniem najlepszy sposób to skorzystać z wierzchołków pobliskich. Dlatego też my ograniczymy się do wpływu tylko dwóch sąsiadów poprzednika oraz następnika. Załóżmy, że nasza krzywa jest cykliczna tzn. zamknięta. P_k to wybrany punkt z krzywej. T_k (od Tangent) to nasza styczna w punkcie P_k. Potrzebny nam będzie jeszcze a - współczynnik dobrany doświadczalnie określający siłę wpływu sąsiadów na styczną (ja wybieram 0.7).

Wzór na styczną w punkcie kontronym:
T_k= a *(P_k+1-_k-1)