UNIVERZA V LJUBLJANI
FAKULTETA ZA MATEMATIKO IN FIZIKO
FAKULTETA ZA RAČUNALNIŠTVO IN INFORMATIKO Računalništvo in matematika – 2. stopnja
Mark Baltič
SHEPARDOVA INTERPOLACIJA
Magistrsko delo
Mentor: izr. prof. dr. Emil Žagar
Ljubljana, 2021
Zahvala
Zahvaljujem se mentorju za natančno branje in dobre popravke, tako vsebinske kot tudi slogovne. Delo je po zaslugi teh vsekakor lažje berljivo. Naslednja zahvala gre moji družini za vso podporo, pa tudi da so mi pustili prostor in čas, ki je bil v času študija še kako pomemben. Zadnja zahvala gre moji partnerki Maši, za spremljanje na skoraj desetletje dolgi poti in pa za ogromno potrpežljivosti.
Kazalo
Program dela vii
1 Uvod 1
1.1 Notacija . . . 2
2 Splošno o interpolacijah 2 2.1 Enodimenzionalna interpolacija . . . 2
2.1.1 Interpolacija z najbližjim sosedom v eni dimenziji . . . 2
2.1.2 Linearna interpolacija v eni dimenziji . . . 3
2.1.3 Polinomska interpolacija v eni dimenziji . . . 3
2.2 Dvodimenzionalna interpolacija . . . 3
2.2.1 Interpolacija najbližjega soseda v dveh dimenzijah . . . 4
2.2.2 Linearna interpolacija v dveh dimenzijah . . . 5
2.2.3 Polinomska interpolacija v dveh dimenzijah . . . 7
3 Shepardova interpolacija 8 3.1 Klasični Shepardov operator . . . 9
3.1.1 Izbira parametra µ . . . 13
3.2 Taylorjeva razširitev . . . 18
3.2.1 Aproksimacija odvodov . . . 20
3.3 Lokalna Shepardova interpolacija . . . 24
3.3.1 Modificirani Shepardov operator . . . 24
3.3.2 k-Shepardov operator . . . 25
3.4 Triangulacijski Shepardov operator . . . 25
4 Preizkus na podatkih 33 4.1 Izbira znanih točk . . . 33
4.2 Algoritmi . . . 34
4.2.1 Algoritmi v preizkusu . . . 35
4.3 Gladka funkcijapeaks . . . 36
4.3.1 Časovna zahtevnost algoritmov . . . 36
4.3.2 Slike rešitev interpolacije gladke funkcije . . . 37
4.3.3 Napake rešitev interpolacije gladke funkcije . . . 44
4.4 Šmarna gora . . . 45
4.4.1 Umestitev podatkov v kontekst . . . 45
4.4.2 Priprava podatkov . . . 46
4.4.3 Slike interpolacij Šmarne gore . . . 46
4.4.4 Napake rešitev interpolacije Šmarne gore . . . 53
5 Zaključek 53
Literatura 55
Stvarno kazalo 57
Program dela
V magistrskem delu obravnavajte Shepardovo interpolacijo in njene izpeljanke. Opi- šite njene lastnosti in implementirajte algoritme. Delovanje preizkusite na primerih analitično podanih funkcij, pa tudi na realnih geodetskih podatkih.
Osnovna literatura
[2] F. Dell’Accio in F. Di Tommaso, Scattered data interpolation by shepard’s like methods: classical results and recent advances, Dolomites Research Notes on Approximation9(Special_Issue) (2016)
Podpis mentorja:
Shepardova interpolacija Povzetek
V magistrskem delu si ogledamo splošne pristope k problemu interpolacije točk, nato pa se osredotočimo na Shepardove interpolacije razpršenih točk v dveh dimen- zijah. Interpolant poskusimo izboljšati z uporabo lokalnosti, odvodov in triangula- cije območja. Predstavljene metode preizkusimo na analitični funkciji in na realnih geodetskih podatkih Šmarne gore pri Ljubljani.
English translation of the title Abstract
In this Master’s Thesis we study general approaches of the point interpolation pro- blem. Then we focus on the Shepard interpolation method on the scattered data points in two dimensions. We try to improve the interpolant by using locality, de- rivatives and by adding a triangulation to data points. First we test the presented methods with an analytical function, and then we move on to real geodetic data obtained from the mountain Šmarna gora, which is located near Ljubljana.
Math. Subj. Class. (2020): oznake kot 65D05, 65D15, 65D17http://www.ams.
org/msc/msc2020.html
Ključne besede: interpolacija, Shepard, razpršeni podatki, inverz razdalij, opera- tor
Keywords: interpolation, Shepard, scattered data, inverse distance, operator
1 Uvod
Ideja o interpolaciji točk se porodi naravno, morda že med sprehodom v naravi po zanimivem terenu. Ali bi lahko podali le nekaj točk na terenu in nato s pomočjo matematičnih modelov rekonstruirali njegovo obliko na računalniku? Kaj pa če smo zaradi težavnosti terena omejeni in lahko dostopamo le do nekaterih točk? Drugo vprašanje nas pripelje do temeljnega problema tega magistrskega dela – ali lahko iz razpršenih točk konstruiramo funkcijo, ki bo lokalno interpolirala izmerjene točke, globalno pa se ji približala. Ko govorimo o razpršenih točkah, imamo v mislih, da v točkah ne najdemo nobene strukture. Ne moremo jih smiselno postaviti na mrežo ali kakšno drugačno ureditev. Ravno zaradi tega je naloga še toliko bolj zanimiva, saj se ne moremo zanašati na urejenost točk – njihova pozicija je lahko naključna ali pa omejena z izvornim problemom.
(a) Regularne točke na mreži. (b) Razpršene točke.
Slika 1: Različni strukturi podatkov.
V nasprotnem primeru se interpolacija v višjih dimenzijah naravno razširi iz metod v eni dimenziji. Seveda lahko točkam v domeni dodelimo urejenost, npr. s triangulacijo, vendar je rešitev problema odvisna od izbire ureditve. Na kaj bi se torej lahko oprli pri računanju vrednosti neznanih točk, če ni urejenosti? Rešitev bomo poizkušali poiskati z upoštevanjem razdalj do vseh znanih točk in njihovih vrednosti. Na to idejo je v drugi polovici 20. stoletja prišel Donald Shepard1, ko je prvič predstavil Shepardovo interpolacijo.
Magistrsko delo sloni na [2], velik poudarek pa je tudi na implementaciji raziska- nih metod.
V drugem poglavju bomo definirali problem interpolacije na razpršenih točkah, predstavili nekaj enostavnih primerov interpolacij in nakazali izziv pri razširjanju metod iz ene v dve dimenziji.
V tretjem poglavju bomo raziskali Shepardovo interpolacijo ter rešili problem ničelnega gradienta s pomočjo razširitve na Taylorjev polinom. Predstavili bomo
1Ameriški matematik Donald Shepard je leta 1968 objavil znani članek [14] o dvodimenzional- nem interpolantu razpršenih podatkov, ki ga sedaj imenujemo Shepardov interpolant ali Shepardov operator.
tudi lokalno izpeljavo Shepardovega operatorja, nazadnje pa operator še razširili na triangulacije.
Četrto poglavje je namenjeno preizkusu na realnih podatkih in praktičnemu vre- dnotenju Shepardovih metod. V ta namen smo v programskem jeziku Matlab im- plementirali raziskane metode, zato bomo v tem poglavju namenili nekaj besed tudi implementaciji.
1.1 Notacija
Pri problemu interpolacije ponavadi nastopa množica točk iz območja Ω ⊂ Rd, za katere poznamo vrednosti funkcije f: Ω → R. Te točke bomo imenovali znane točke, označevali pa jih bomo z x1, . . . , xn v eni dimenziji ali x1, . . . ,xn v večih dimenzijah. Praktični primeri se bodo nanašali na dve dimenziji, zato bomo včasih točke pisali tudi po koordinatahxi = (xi, yi). Ko dobimo interpolant iz znanih točk, ga želimo vrednotiti v točkah, v katerih vrednosti funkcije še ne poznamo. Te točke bomo imenovali neznane točke, označevali jih bomo brez indeksov, torej kot x v eni dimenziji ali x v višjih dimenzijah (v primeru dveh dimenzij bo x = (x, y)). S simbolom ∥x∥ bomo vedno označevali evklidsko razdaljo.
2 Splošno o interpolacijah
Za začetek si oglejmo nekaj pristopov k problemu interpolacije. Opisali bomo tri enostavne metode v eni dimenziji, nato pa si pogledali razširitve teh metod v dve dimenziji.
2.1 Enodimenzionalna interpolacija
Za lažjo predstavo se lotimo interpolacije točk v eni dimenziji. Imamo n+ 1paroma različnih točk x0, . . . , xn ∈ X ⊂ R in pripadajoče funkcijske vrednosti funkcije f, ki jih označimo z fi = f(xi). Iščemo interpolacijsko funkcijo g, ki se bo v znanih točkah ujemala s funkcijo f, torej g(xi) = f(xi) za i = 0, . . . , n. Ko želimo izraču- nati vrednost funkcije f v neznani točki x, to vrednost aproksimiramo z vrednostjo interpolanta, tako da izračunamo g(x).
2.1.1 Interpolacija z najbližjim sosedom v eni dimenziji
Najenostavnejši pristop je interpolacija z najbližjim sosedom. Vsaki točki x ∈ R poiščemo najbližjo točko xj, kjer je j = arg minxj∈X|x− xj|. Nato za vrednost interpolanta v iskani točki x določimo kar vrednost v najbližji točki, torej g(x) = f(xj). Za enoličnost rešitve definiramo še vrednosti v točkah, ki so aritmetična sredina dveh zaporednih točk iz naraščajoče urejene množice X. To storimo tako, da za točko x = xi+x2i+1 določimo vrednost g(x) = f(xi)+f2(xi+1). Dobljena rešitev v točkah množice X ni zvezna. Primer take interpolacije si lahko ogledamo na sliki 2a.
2.1.2 Linearna interpolacija v eni dimenziji
Če znane točke uredimo naraščajoče, lahko med zaporednima točkama določimo linearen interpolant. Za točko x, za katero velja xi ≤ x≤ xi+1, določimo vrednost interpolanta tako, da izračunamo vrednost na daljici, ki jo določata točkixi inxi+1 po formuli g(x) = xfi+1−fi
i+1−xix+ fixi+1x −fi+1xi
i+1−xi . Dobimo lomljenko, ki je zvezna, a ni odvedljiva v znanih točkah. Primer take interpolacije si lahko ogledamo na sliki 2b.
2.1.3 Polinomska interpolacija v eni dimenziji
Če želimo gladek interpolant, raje poizkusimo s polinomsko aproksimacijo. V eni dimenziji nam pomaga izrek o enoličnem obstoju takega polinoma [13, izrek 10.1, str. 262]
Izrek 2.1. Za paroma različne točke x0, . . . , xn in vrednosti y0, . . . , yn obstaja na- tanko en polinompn stopnjen ali manj, za katerega veljapn(xi) = yi zai= 0, . . . , n.
Dokaz. Obstoj polinoma pnpokažemo s konstrukcijo. Za vsaki= 0, . . . , nvpeljimo Lagrangeev bazni polinom
ln,i(x) =
n
∏
k=0k̸=i
x−xk xi−xk.
Iskani polinom dobimo kot linearno kombinacijo Lagrangeevih baznih polinomov in vrednosti v znanih točkah, torej
pn(x) =
n
∑
k=0
ynln,k(x).
Hitro lahko preverimo, da velja pn(xi) =yi za vsaki= 0, . . . , n.
Dokažimo še enoličnost rešitve. Naj obstaja še en polinom rn stopnje manjše ali enaken, za katerega veljarn(xi) = yizai= 0, . . . , n. Oglejmo si polinomq =pn−rn. Dobljeni polinom je tudi stopnje manjše ali enake n. Velja tudi q(xi) = 0 za i = 0, . . . , n. Posledica osnovnega izreka algebre nam pove, da ima vsak nekon- stantni polinom stopnje n ali manj s kompleksnimi koeficient natankon ničel štetih z večkratnostmi. V našem primeru ima polinom q kar n + 1 različnih ničel, kar implicira q= 0. Velja torejpn =rn in enoličnost je dokazana.
Dobljen interpolant je polinom in torej gladka funkcija naR. Primer take inter- polacije si lahko ogledamo na sliki 2c.
2.2 Dvodimenzionalna interpolacija
Posplošitev metod v dve dimenziji ni najlažji zalogaj. V vsakem izmed primeru se bomo soočili z neprikladnostmi, naj bo to urejenost podatkov na mrežo, nezveznost ali pa triangulacija območja.
Sedaj imamo opravka z n+ 1 paroma različnimi točkami x0, . . . ,xn ∈ X ⊂ R2 in vrednostmi funkcije večih spremenljivkf v teh točkah. Interpolant je dvodimen- zionalen objekt v trodimenzionalnem prostoru.
(a) Interpolacija z najbljižjim sosedom na
petih točkah. (b) Linearna interpolacija na petih točkah.
(c) Polinomska interpolacija na petih točkah.
Slika 2: Enodimenzionalne interpolacije.
Opomba: Metode se da posplošiti tudi v višje dimenzije Rd, kjer je d >2, vendar se bomo osredotočili na dimenzijo 2, kjer smo tudi preizkušali interpolacijo na realnih podatkih.
2.2.1 Interpolacija najbližjega soseda v dveh dimenzijah Za vsako znano točko xi definiramo Voronojevo celico kot
Vi ={x∈Ω| ∥x−xi∥ ≤ ∥x−xj∥ zaj ̸=i}.
Množica Voronojevih celic za vse znane točke nam predstavlja particijo ravnine, ki jo imenujemo Voronojev diagram [17]. Pri interpolaciji najbližjega soseda je za točko x ∈ Vi vrednost funkcije enaka kar f(x) = fi. V kolikor točka meji na več Voronojevih celic, ji dodelimo povprečno vrednost funkcijskih vrednosti znanih točk, ki so inducirale te Voronojeve celice. Metoda je učinkovita za razpršene podatke,
Slika 3: Voronojev diagram na petindvajsetih točkah.
a je žal nezvezna na robovih Voronojevih celic. Primer take interpolacije si lahko ogledamo na sliki 4.
(a) Voronojev diagram.
(b) Interpolacija z najbližjim sosedom.
Slika 4: Dvodimenzionalna interpolacija z najbližjim sosedom na petih točkah. Pro- jekcija na xy ravnino nam pokaže Voronojev diagram.
2.2.2 Linearna interpolacija v dveh dimenzijah
Najprej si oglejmo primer, ko zaostrimo pogoj na urejenost znanih točk tako, da te ležijo na mreži. Sedaj lahko izvedemo bilinearno interpolacijo. Poznamo vrednosti funkcije f za točke na regularni mreži ((xi, yi))ni=0, kjer je xi = x0 +ihx in yi = x0 +ihy. Konstanti hx in hy predstavljata razdaljo med zaporednima točkama v abscisni in ordinatni smeri. Naj bo točka (x, y) znotraj pravokotnika, ki ga tvorijo točke(xi, yi),(xi+1, yi),(xi, yi+1),(xi+1, yi+1). Tedaj lahko vrednost interpolanta g v točki(x, y) izračunamo po formuli
g(x, y) =[
1−x x]
[ f(xi, yi) f(xi, yi+1) f(xi+1, yi) f(xi+1, yi+1)
] [1−y y
] .
Interpolant je linearen na daljicah med zaporednima točkama mreže, v notranjih točkah mreže pa je kvadratična funkcija dveh spremenljivk – od tod tudi ime bili- nearna interpolacija. Primer take interpolacije si lahko ogledamo na sliki sliki 5. Če
(a) Bilinearen interpolant nad štirimi sose- dnjimi točkami mreže.
(b) Bilinearna interpolacija funkcije f(x, y) = sinxcosy.
Slika 5: Primer bilinearnega interpolanta.
točke ne ležijo na regularni mreži, potem lahko dodamo strukturo domeni s trian- gulacijo. Pogosto zahtevamo, da je triangulacija regularna, kar pomeni, da se dva trikotnika sekata samo v oglišču ali vzdolž cele stranice. Za točko x znotraj triko- tnika, ki ga določajo znane točke xi,xj,xk, določimo interpolacijsko vrednost kot vrednost na ravnini, ki poteka skozi znane točke. To lahko izračunamo po formuli g(x, y) = ax+by+c, kjer so konstante a, b, c rešitve sistema linearnih enačb
⎡
⎣
xi yi 1 xj yj 1 xk yk 1
⎤
⎦
⎡
⎣ a b c
⎤
⎦=
⎡
⎣
f(xi, yi) f(xj, yj) f(xk, yk)
⎤
⎦.
Območje lahko trianguliramo na različne načine, interpolant pa je odvisen od tri- angulacije. Priljubljena je Delaunayeva triangulacija, ki maksimizira minimalni notranji kot v trikotnikih triangulacije. Oglejmo si primer na šestih točkah z dvema različnima triangulacijama (slika 6). Štirim točkam na ogliščih enotskega kvadrata dodelimo vednost 0, notranjima točkama pa vrednost 1. Interpoliramo območje z linearno interpolacijo. Na sliki 7 vidimo rezultat. Najbolj očitna razlika v tem pri- meru je ta, da sta pri drugi triangulaciji dve oglišči kvadrata povezani. Daljica, ki ju povezuje, je tako na višini0. Rezultat tega sta dva vrhova, ki sta pri Delaunayevi triangulaciji povezana z daljico na višini 1 (na sliki izgleda kot pobočje gore). Do- bljen interpolant prav tako ni odvedljiv na robovih triangulacije. Kot zanimivost še navedemo povezavo med Voronojevimi diagrami in Delaunayevo triangulacijo. Kon- struiramo Voronojev diagram za množico točk (xi)ni=0. Če si Voronojevi celici Vi in Vj delita stranico, potem z daljico povežemo točki xi in xj. Dobljena struktura je ravno Delaunayeva triangulacija (slika 8).
(a) Delaunayeva triangulacija. (b) Primer splošne triangulacije.
Slika 6: Dve triangulaciji šestih točk.
(a) Linearen interpolant nad Delaunayevo triangulacijo.
(b) Linearen interpolant nad triangulacijo s slike 6b.
Slika 7: Linearen interpolant nad dvema različnima triangulacijama šestih točk.
2.2.3 Polinomska interpolacija v dveh dimenzijah
V eni dimenziji smo lahko zadovoljni z izrekom o obstoju interpolacijskega polinoma, v višjih dimenzijah pa žal ni tako enostavno. Povzemimo ugotovitve iz [15, razdelek 1.1].
Definicija 2.2 (Haarov prostor). Naj območje Ω⊂Rd vsebuje vsajn točk. Linea- ren prostor zveznih funkcijV ⊂C(Ω) se imenuje Haarov prostor natanko tedaj, ko za vsako množico paroma različnih točk{x1, . . . ,xn} ⊆Ωin pripadajoče funkcijske vrednosti {f1, . . . , fn} ⊆ R obstaja natanko ena funkcija p ∈ V, za katero velja p(xi) = fi za vsak i= 1, . . . , n.
V razdelku 2.1.3 smo videli, da v eni dimenziji obstaja enoličen polinom stopnje vsajn, ki interpoliran+ 1paroma različnih točk. Prostor polinomov stopnje vsajn nad realnimi števili, označimo ga s Pn(R), je torejn-dimenzionalen Haarov prostor za vsako območje Ω⊆R, ki vsebuje vsaj n paroma različnih točk.
Sedaj si oglejmo izrek o eksistenci Haarovih prostorov v višjih dimenzijah.
Slika 8: Povezava med Voronojevim diagramom in Delaunayevo triangulacijo.
Izrek 2.3 (Mairhuber-Curtis). Naj bo d≥2 in naj območje Ω⊆Rd vsebuje vsaj n točk in eno notranjo točko. Potem na območju Ω ne obstaja noben Haarov prostor dimenzije n ≥2.
Po izreku ne obstaja enolična polinomska interpolacija v višjih dimenzijah, razen v nekaterih enostavnih primerih. Do neskončne množice rešitev pridemo že, če v trodimenzionalnem prostoru postavimo vse znane točke na premico p. Med slikami polinomov dveh spremenljivk, ki interpolirajo te točke, so vse ravnine v prostoru, na katerih leži premica p.
To nas napeljuje na spoznanje, da iščemo nove načine interpolacije razpršenih točk v večdimenzionalnih prostorih. Eno metodo si bomo podrobneje ogledali v naslednjem poglavju.
3 Shepardova interpolacija
Brezmrežne interpolacijske metode na razpršenih podatkih so v numerični matema- tiki aktualen problem. V zadnjem času so postale še posebej popularne metode z radialnimi baznimi funkcijami, o katerih si lahko več preberete v [15]. Naš inter- polant sodi k metodam, ki uporabljajo uteži z inverzi razdalj do znanih točk (ang.
inverse distance weighting). Najprej bomo izpeljali klasično verzijo interpolanta, nato pa ga poizkušali izboljšati.
Vseskozi naj bo X ={x1, . . . ,xn} množica n paroma različnih točk iz območja Ω⊂R2. Vsaki točkixi,i= 1, . . . , n, pripišemo vrednost funkcijef(xi) = fi. Oblika interpolanta s naj bo linearna kombinacija vrednosti funkcije v znanih točkah in baznih funkcij ali uteži w, s(x) = ∑n
i=1w(x)fi. Želimo, da je operator funkcija, ki interpolira znane točke, je zvezna in diferenciabilna.
3.1 Klasični Shepardov operator
Definicija 3.1. Za pozitiven parameter µ definiramo klasični Shepardov operator kot
Sµ[f](x) = {∑n
i=1Aµ,i(x)fi, zax∈/ X, fi, zax=xi ∈X, pri čemer je
Aµ,i(x) = ∥x−xi∥−µ
∑n
k=1∥x−xk∥−µ, i= 1, . . . , n, (3.1) kjer je ∥·∥ evklidska norma.
Intuicija za izbiro uteži izvira iz lastnosti neurejenosti znanih točk, saj se lahko tako zanašamo le na razdalje do le teh. Ker si želimo, da bi bližnje točke imele večji vpliv na vrednost funkcije, uporabimo inverz razdalje. Parameterµnam daje več svobode pri analiziranju in uporabi metode. Imenovalec v (3.1) služi kot normalizacija vseh uteži in, kot bomo videli v naslednjih ugotovitvah, prinese nekaj lepih lastnosti interpolantu.
Izpeljimo nekaj trditev o Shepardovem operatorju. Več analitičnih značilnosti tega operatorja je na voljo v [7] in [14].
Trditev 3.2. Shepardov operator je zvezna funkcija spremenljivke x.
Dokaz. V točkah x ∈/ X zveznost sledi iz zveznosti uteži (3.1), saj je operator linearna kombinacija zveznih funkcij. Preveriti moramo še zveznost v točkah množice X. Za vsak indeks i = 1, . . . , n, se vrednost operatorja Sµ[f](x) približuje fi, čim
∥x−xi∥ →0. Osredotočimo se na točkox iz okolice točke xi in si oglejmo, kaj se zgodi z utežmi, ko se razdalja med točkama manjša. Za utež Aµ,j(x), j ̸=i, velja
lim
∥x−xi∥→0Aµ,j(x) = lim
∥x−xi∥→0
∥x−xj∥−µ
n
∑
k=1
∥x−xk∥−µ
= lim
∥x−xi∥→0
∥xi−xj∥−µ
∥x−xi∥−µ+
n
∑
k=1k̸=i
∥xi−xk∥−µ
= 0,
saj je števec ulomka omejen, v imenovalcu pa so omejeni vsi členi razen člena∥x− xi∥−µ, ki pa gre v limiti proti neskončno. Imamo torej omejen števec, imenovalec pa raste čez vse meje, kar v limiti pošlje celoten ulomek proti 0.
Izračunajmo še limito za utež Aµ,i(x), tako da števec in imenovalec pomnožimo s členom ∥x−xi∥µ, torej
lim
∥x−xi∥→0Aµ,i(x) = lim
∥x−xi∥→0
∥x−xi∥−µ· ∥x−xi∥µ ( n
∑
k=1
∥x−xk∥−µ )
· ∥x−xi∥µ
= lim
∥x−xi∥→0
1
n
∑
k=1
∥x−xi∥µ
∥x−xk∥µ
= lim
∥x−xi∥→0
1 1 +
n
∑
k=1k̸=i
∥x−xi∥µ
∥xi−xk∥µ
= 1,
kjer zadnjo enakost dobimo iz dejstva, da je i-ti člen v imenovalcu enak1, v ostalih členih pa je imenovalec omejen, števec pa gre proti 0.
Od tod vidimo, da velja
∥x−xlimi∥→0Sµ[f](x) = lim
∥x−xi∥→0 n
∑
i=1
Aµ,i(x)fi =fi, kar pokaže, da je operator zvezen v vseh točkah domene Ω.
Lema 3.3. Za vsako točko x ∈ Ω so uteži Aµ,i(x), i = 1, . . . , n, nenegativne, kardinalne in tvorijo particijo enote.
Dokaz. Uteži so nenegativne, saj jih dobimo s potenciranji, seštevanji in deljenji razdalj, ki so nenegativne funkcije, prej naštete operacije pa objektom ne spremenijo predznaka. Kardinalnost tukaj pomeni, da se glede na znane točke iz množice X uteži vedejo kot Kroneckerjeva delta [16], torej da velja
Aµ,j(xi) =
{1, i=j, 0, i̸=j.
To lastnost smo pokazali v dokazu trditve 3.2. Kratek izračun nam tudi pokaže, da uteži tvorijo particijo enote
n
∑
i=i
Aµ,i(x) =
n
∑
i=i
∥x−xi∥−µ
n
∑
k=1
∥x−xk∥−µ
=
n
∑
i=i
∥x−xi∥−µ
n
∑
k=1
∥x−xk∥−µ
= 1.
Na sliki 9 vidimo sliko utežiA2,i, prirejeno znani točkixi. Opazimo, da v ostalih znanih točkah xj,j ̸=i, zavzame vrednost nič, najvišje vrednosti pa zavzame ravno v okolici točke xi.
Posledica 3.4. Shepardov operator je stabilen. Natančneje,
i=1,...,nmin fi ≤Sµ[f](x)≤ max
i=1,...,nfi, za vsako točko x∈Ω.
Slika 9: Slika klasične utežiA2,i nad znano točkoxi. Ostale znane točke so označene z rdečo barvo.
Dokaz. Naj bo m = mini=1,...,nfi in M = maxi=1,...,nfi. Ker so vse uteži Aµ,i(x) nenegativne, velja Aµ,i(x)fi ≤ Aµ,i(x)M in Aµ,i(x)fi ≥ Aµ,i(x)m za vsak i = 1, . . . , n. Računajmo
Sµ[f](x) =
n
∑
i=1
Aµ,i(x)fi ≤
n
∑
i=1
Aµ,i(x)M =M
n
∑
i=1
Aµ,i(x) =M,
kjer smo pri zadnji enakosti upoštevali, da uteži tvorijo particijo enote. Analogno dokažemo še spodnjo mejo, kot
Sµ[f](x) =
n
∑
i=1
Aµ,i(x)fi ≥
n
∑
i=1
Aµ,i(x)m=m
n
∑
i=1
Aµ,i(x) = m.
Trditev 3.5. Vrednosti Shepardovega operatorja Sµ[f](x) se približujejo povprečju vrednosti znanih točk 1n∑n
i=1fi, čim gre min
i=1,...,n∥x−xi∥ → ∞.
Dokaz. Označimo d:= mini=1,...,n∥x−xi∥ in si oglejmo, kaj se zgodi z i-to utežjo, ko pošljemo d → ∞. Vsi členi ∥x−xi∥µ se začnejo obnašati kot ∥x∥µ, kar nas privede do
d→∞lim Aµ,i(x) = lim
d→∞
∥x−xi∥−µ
n
∑
k=1
∥x−xk∥−µ
= lim
d→∞
∥x∥−µ
n
∑
k=1
∥x∥−µ
= lim
d→∞
∥x∥−µ n· ∥x∥−µ = 1
n.
(a) Enodimenzionalen primer.
Povprečje vrednosti točk znaša2.99.
(b) Dvodimenzionalen primer.
Povprečje vrednosti točk znaša 0.01.
Slika 10: Obnašanje Shepardovega operatorja (s parametrom µ = 2) v oddaljenih točkah. Na prvi sliki so interpolirane 4 točke z vrednostmi funkcije f(x) = x2, na drugi sliki pa 20 točk z vrednostmi funkcije f(x, y) =xe−x2−y2. V obeh primerih se vrednosti interpolanta približujejo povprečju znanih vrednosti, čim se oddaljujemo od množice X. Opazimo lahko tudi, da se interpolant z oddaljevanjem izravna.
Sedaj uporabimo to na Shepardovem operatorju, torej
d→∞lim Sµ[f](x) = lim
d→∞
n
∑
i=1
Aµ,i(x)fi =
n
∑
i=1
d→∞lim Aµ,i(x)fi
=
n
∑
i=1
1
nfi = 1 n
n
∑
i=1
fi.
Trditev 3.6. Shepardov operator ima algebraično stopnjo natančnosti enako 0.
Dokaz. Pokazati moramo, da Shepardov operator rekonstruira polinome stopnje 0, torej konstantne funkcije. Naj bo funkcijaf(x) = c, kjer jecpoljubno realno število, konstantna. Za vrednosti znanih točk torej velja fi =c zai = 1, . . . , n. Trdimo, da Shepardov operator vrne vrednost c za vsako točkox∈Ω. Računajmo:
Sµ[f](x) =
n
∑
i=1
Aµ,i(x)fi =
n
∑
i=1
Aµ,i(x)c=c
n
∑
i=1
Aµ,i(x) =c,
kjer pri zadnji enakosti upoštevamo dejstvo iz leme 3.3, da uteži tvorijo particijo enote.
Shepardov interpolant torej poizkuša interpolirati območje v bližini znanih točk iz množice X, ko se oddaljujemo od te množice, pa se interpolant izravna.
V nadaljevanju si bomo ogledali, kako se Shepardov interpolant obnaša s spre- minjanjem parametra µ.
(a) Enodimenzionalen primer. (b) Dvodimenzionalen primer.
Slika 11: Stabilnost Shepardovega operatorja (s parametrom µ= 2). Na prvi sliki so interpolirane 4 točke z vrednostmi funkcijef(x) =x2, na drugi sliki pa 20 točk z vrednostmi funkcije f(x, y) = xe−x2−y2. V obeh primerih se vrednosti interpolanta gibljejo med najnižjo in najvišjo vrednostjo znanih točk.
3.1.1 Izbira parametra µ
V svojem izvornem članku [14] se je Shepard osredotočal na primerµ= 2, saj v ute- žeh nastopa evklidska razdalja. S kvadriranjem se tako znebimo korena in postanejo uteži racionalne funkcije ter računsko bolj prikladne. Če upoštevamo, da imenovalec v izrazu za uteži služi normalizaciji, lahko števec primerjamo z »gravitacijskim zako- nom«, kjer je gravitacijska sila obratno sorazmerna s kvadratom razdalje F ∝r−2. Če si dopustimo parameter spreminjati, dobimo različne oblike rešitve. Z višanjem parametra namreč postaja rešitev v okolici znanih točk vedno bolj ploščata. O tem se prepričamo s sledečo trditvijo.
Trditev 3.7. Naj bo Sµ[f] Shepardov operator. Če pošljemo parameter µ → ∞, potem operator konvergira (po točkah) k odsekoma konstantni funkciji, ki sovpada z interpolacijo najbližjega soseda (razdelek 2.2.1).
Dokaz. Brez škode za splošnost lahko privzamemo, da je točkax∈Ωnajbližje znani točki x1, torej je |x−x1| <|x−xj| za j = 2, . . . , n. Izračunajmo, kaj se v limiti zgodi z utežjo Aµ,j(x)za j >1:
µ→∞lim Aµ,j(x) = lim
µ→∞
∥x−xj∥−µ
n
∑
k=1
∥x−xk∥−µ
· ∥x−x1∥µ
∥x−x1∥µ
= lim
µ→∞
(∥x−x1∥
∥x−xj∥ )µ
1 +
n
∑
k=2
(∥x−x1∥
∥x−xk∥
)µ = 0.
(3.2)
V zadnji enakosti upoštevamo, da je limµ→∞
(∥x−x
1∥
∥x−xj∥
)µ
= 0 za vsak j > 1, saj je
vrednost osnove potence manjša od 1. Oglejmo si še primer za utež Aµ,1(x):
µ→∞lim Aµ,1(x) = lim
µ→∞
∥x−x1∥−µ
n
∑
k=1
∥x−xk∥−µ
· ∥x−x1∥µ
∥x−x1∥µ
= lim
µ→∞
1 1 +
n
∑
k=2
(∥x−x1∥
∥x−xk∥
)µ = 1.
Pri tem ponovimo razmislek prejšnjega izračuna. Upoštevajmo izračunane limite pri izračunu vrednosti Shepardovega operatorja v točki x, da dobimo
µ→∞lim Sµ[f](x) = lim
µ→∞
n
∑
i=1
Aµ,i(x)fi =f1. (3.3) Shepardov operator v limiti vrne vrednost najbližje točke. Oglejmo si še poseben primer, ko točka območja nima le enega najbližjega soseda. Naj bo 1 ≤ d ≤ n in naj velja |x−x1| = |x−x2| = · · · =|x−xd| < |x−xj| za j = d+ 1, . . . , n. Za vsak l = 1. . . d, je limita enaka
µ→∞lim Aµ,l(x) = lim
µ→∞
∥x−xl∥−µ
n
∑
k=1
∥x−xk∥−µ
· ∥x−xl∥µ
∥x−xl∥µ
= lim
µ→∞
1 d+
n
∑
k=d+1
(∥x−xl∥
∥x−xk∥
)µ = 1 d.
Analogno kot v (3.2) lahko pokažemo, da zaj =d+1, . . . , n, veljalimµ→∞Aµ,j(x) = 0. Rezultata uporabimo pri izračunu vrednosti operatorja in dobimo
µ→∞lim Sµ[f](x) = lim
µ→∞
n
∑
i=1
Aµ,i(x)fi =
d
∑
i=1
1 dfi = 1
d
d
∑
i=1
fi, (3.4) kar je ravno povprečje vrednosti najbližjih sosedov. Iz (3.3) in (3.4) vidimo, da v limiti µ → ∞ Shepardov interpolant sovpada z definicijo interpolacije najbližjega soseda.
Na sliki 12 je vidno obnašanje interpolanta z višanjem parametra µ. Pri viso- kem parametru µ = 20, rezultat že spominja na rešitev interpolacije z najbližjim sosedom. Izbiro parametra µ pa lahko še bolj sprostimo, tako da vsaki znani točki xi dodelimo parameter µi, kot so to storili v [7]. S tem lahko reguliramo okolice posameznih znanih točk, saj z izbiro višje vrednosti parametra območje okrog točk bolj sploščimo. Globlje v izbiro primernih vrednosti µi se nismo spuščali, saj bi za dober algoritem potrebovali dodatno znanje, kot na primer določanje ploščatosti reliefa iz znanih točk in uporabo teh znanj za določanje smiselnih parametrov. Kot
(a)µ= 2 (b)µ= 3
(c)µ= 5 (d)µ= 20
Slika 12: Spreminjanje klasičnega Shepardovega operatorja z večanjem parametra µ.
zanimivost na sliki 13 pokažemo interpolant, pri katerem dodelimo višjo vrednost parametra µi najvišji in najnižji točki.
Slika 13: Na levi sliki imajo vse uteži ekponent enak µ = 2. Na desni sliki pa smo utežema, ki pripadata najvišji in najnižji točki, dodelili ekponent 6. Najvišja in najnižja točka sta oznaženi z rdečo barvo.
Do sedaj smo ugotovili, da je Shepardov operator zvezna funkcija za vsak po- zitiven parameter µ. Oglejmo si še, kako je z odvedljivostjo. Pomagali si bomo z alternativno reprezentacijo Shepardovega operatorja, ki sloni na drugačnem zapisu uteži.
Lema 3.8. Za vsako točko x∈Ω in znane točke xi ∈X za i= 1, . . . , n, naj bo
A˜µ,i(x) =
n
∏
j=1 j̸=i
∥x−xj∥µ
n
∑
k=1 n
∏
j=1 j̸=k
∥x−xj∥µ
. (3.5)
Potem je A˜µ,i(x) = Aµ,i(x).
Dokaz. V lemi 3.3 smo videli, da lahko uteži Aµ,i(x) zvezno razširimo na celotno domensko območje Ω, če za točko xj ∈X upoštevamo Aµ,i(xj) = δij. To sovpada z vrednostmi uteži A˜µ,i(x). Sedaj pokažimo še, da enakost velja tudi za točke x∈ Ω\X:
Aµ,i(x) = ∥x−xi∥−µ
n
∑
k=1
∥x−xk∥−µ
·
n
∏
j=1
∥x−xj∥µ
n
∏
j=1
∥x−xj∥µ
=
n
∏
j=1 j̸=i
∥x−xj∥µ
n
∑
k=1 n
∏
j=1 j̸=k
∥x−xj∥µ
= ˜Aµ,i(x).
Uteži A˜µ,i so prikladne, saj so zvezno definirane povsod in nimajo nobenih sin- gularnosti. S pomočjo teh uteži lahko zapišemo Shepardov operator enostavneje.
Posledica 3.9. Shepardov operator lahko zapišemo tudi v obliki Sµ[f](x) =
n
∑
i=1
A˜µ,i(x)fi za x∈Ω,
kjer so uteži A˜µ,i(x) dane z (3.5).
Ta oblika nam bo pomagala pri izračunih odvodov.
Izrek 3.10. Naj bodo X = {x1, . . . ,xn} znane točke, µ > 0 parameter in naj bo Sµ[f] Shepardov operator definiran na teh znanih točkah. Potem veljajo sledeče trditve:
1. Če je parameter µ ≤ 1, potem Shepardov operator v splošnem ni odvedljiv v znanih točkah X.
2. Če je parameter µ > 1, potem sta smerna odvoda Shepardovega operatorja v znanih točkah enaka 0.
V ostalih točkah x∈Ω\X je Shepardov operator odvedljiv.
Dokaz. Najprej bomo pokazali drugi del, ko jeµ >1, in se osredotočili na dimenzijo 2. Uporabili bomo zapis za Shepardov operator z utežmi (3.5). Za lažje zapise definiramo ri(x) = ∥x−xi∥ = √
(x−xi)2+ (y−yi)2. Utež A˜µ,i(x) lahko tako zapišemo kot
A˜µ,i(x) =
n
∏
j=1 j̸=i
rj(x)µ
n
∑
k=1 n
∏
j=1 j̸=k
rj(x)µ .
Izračunajmo parcialni odvod funkcijeriµ v smeri x
∂riµ
∂x (x) = µri(x)µ−11
2ri(x)−12(x−xi) =µri(x)µ−2(x−xi).
Preverimo, kaj se dogaja z odvodom, ko se točka x približuje točki xi (kar lahko ponazorimo z ri(x)→0), torej
lim
ri(x)→0
∂riµ
∂x(x) = 0, (3.6)
saj je µ >1. Sedaj definiramoBk(x) = ∏
j̸=krj(x). Hitro ugotovimo, da velja Bk(xi) = 0, zak ̸=i. (3.7) Izračunamo še odvod
∂Bk
∂x (x) = ∑
m̸=k
∂rmµ
∂x (x) ∏
l̸=m,l̸=k
rl(x)µ. Podobno kot zgoraj opazujemo odvod v bližini točke xi:
lim
ri(x)→0
∂Bk
∂x (x) = 0 zak ̸=i, (3.8)
saj v primeru, ko je m =i, uporabimo (3.6). V ostalih členih vsote pa v produktu nastopa faktor ri(x)µ. Uporabimo zapis za Shepardov operator z utežmi A˜µ,i(x) iz posledice 3.9 in ga prepišimo z novimi oznakami, torej
Sµ[f](x) =
n
∑
i=1
Bi(x)fi n
∑
i=k
Bk(x) .
Odvajajmo operator po pravilu za odvod ulomka
∂Sµ[f]
∂x (x) = [ ( n
∑
i=1
∂Bi
∂x (x)fi
) ( n
∑
i=k
Bk(x) )
− ( n
∑
i=1
Bi(x)fi ) ( n
∑
i=k
∂Bk(x)
∂x
) ]/( n
∑
i=k
Bk(x) )2
.
(3.9)
Izračunajmo, kaj se zgodi z izrazom (3.9) v limiti, ko pošljemo točko x proti točki xi. Če upoštevamo že znana rezultata (3.7) in (3.8) ugotovimo, da je veliko členov enakih 0. V naslednjem računu ohranimo le neničelne in dobimo
lim
ri(x)→0
∂Sµ[f]
∂x (x) = [∂Bi
∂x (x)
⏐
⏐
⏐
⏐
⏐x=xi
fiBi(xi)−Bi(xi)∂Bi
∂x (x)
⏐
⏐
⏐
⏐
⏐x=xi
fi ]/
Bi(xi)2 = 0.
Iz izraza (3.9) tudi vidimo, da je Shepardov operator odvedljiv v vseh neznanih točkah x∈Ω\X, saj je imenovalec vedno različen od 0.
Za parcialni odvod v smeri y pa opazimo, da velja ∂r
µ i
∂x(x) =µri(x)µ−2(y−yi).
Postopamo analogno kot v smeri x in pridemo do podobnega rezultata.
Za prvi del izreka, ko jeµ≤1, pa opazimo, da limita v izrazu (3.6) v splošnem ne obstaja. Iz tega sledi, da tudi odvod Shepardovega operatorja v splošnem ne obstaja.
(a)µ= 0.5 (b)µ= 1
Slika 14: Odvedljivost Shepardovega operatorja pri parametrihµ≤1. Interpoliranih je 5 točk z vrednostmi funkcije f(x) = x2.
Na sliki 14 vidimo, da v znanih točkah množiceX odvodi res ne obstajajo. S slik razberemo še več – pri parametrih µ≤ 1, se interpolant ne prilega izvirni funkciji, ker pa interpolira dane točke, v okolici teh nastanejo špice, zaradi katerih odvod tam ni definiran. Pri parametru µ = 1 pa levi in desni odvod v znanih točkah lahko obstajata, vendar sta različno predznačena [14]. Slika 15 nam potrdi, da je interpolant pri parametrih µ≥1 povsod odvedljiv.
3.2 Taylorjeva razširitev
V primeru, ko vzamemo parameter µ ≥ 1, nam lokalni ekstremi ali prevoji, ki so posledica ničelnih odvodov v znanih točkah, niso všeč, saj v splošnem to ni lastnost znanih točk. Interpolant bi želeli izboljšati z dodajanjem odvodov funkcij.
Ti v praktičnih primerih seveda niso vedno na voljo, zato naletimo še na problem računanja odvodov v razpršenih podatkih. Kako postopati v tem primeru, si bomo
(a)µ= 1.5 (b) µ= 2
Slika 15: Odvedljivost Shepardovega operatorja pri parametrihµ >1. Interpoliranih je 5 točk z vrednostmi funkcijef(x) =x2.
ogledali kasneje, zaenkrat privzemimo, da so odvodi ∂f∂x(xi),∂f∂y(xi) za i = 1, . . . , n znani, in po zgledu [5] poizkusimo izboljšati interpolant.
Definicija 3.11. Za funkcijo f in znano točko xi definiramo Taylorjev polinom dveh spremenljivk stopnje r kot
Tr[f,xi](x) =
r
∑
|k|=0
Dkf(xi)
k! (x−xi)k,
kjer jek = (k1, k2) multiindeks in|k|= k1+k2. Sedaj lahko za vsako točko x ∈Ω definiramo Shepard-Taylorjev operator kot
STr[f](x) =
n
∑
i=1
Aµ,i(x)Tr[f,xi](x). (3.10) Pri r = 0 se (3.10) ujema s klasičnim Shepardovim operatorjem. Ta je rekon- struiral konstantne polinome, z dodajanjem členov Taylorjevega polinoma pa lahko algebraično stopnjo interpolanta zvišamo (več o tem si lahko bralec prebere v [2]).
Omejili se bomo le na odvode prve stopnje, torejr = 1, zato lahko (3.10) prepišemo v obliko
ST1[f](x) =
n
∑
i=1
Aµ,i(x)(fi+ ∂f
∂x(xi)(x−xi) + ∂f
∂y(xi)(y−yi)).
Z dodajanjem vrednosti odvodov smo dvignili algebraično stopnjo za 1, iz česar sledi da Taylor-Shepardov operator rekonstruira polinome prve stopnje.
Oglejmo si dva primera, kjer imamo na voljo podatke o vrednostih funkcij in prvih odvodov za razpršene podatke.
Primer 3.12. Pričnemo z enodimenzionalnim primerom. Vzemimo 10 naključ- nih točk z intervala [−5,5] in jim dodelimo vrednosti funkcije f(x) = xe−x2 ter
(a) Interpolacija s klasičnim Shepardovim operatorjem.
(b) Interpolacija s Shepard-Taylorjevim ope- ratorjem ST1[f].
Slika 16: Preizkus Shepard-Taylorjevega operatorja na enodimenzionalnem primeru.
odvodaf′(x) = (1−2x2)e−x2. Preizkusimo klasični Shepardov operator in Shepard- Taylorjev operator pri parametru µ= 2.
Na sliki 16 vidimo, da smo s Shepard-Taylorjevim operatorjem rešili problem ničelnih odvodov v znanih točkah, prav tako se interpolant bolje prilega originalni funkciji f.
Primer 3.13. Primer 3.12 razširimo v dve dimenziji. Izvedimo interpolacijo na 50 naključnih točkah, ki jim pripišemo vrednosti funkcije in parcialnih odvodov
f(x, y) = xe−x2−y2, ∂f
∂x(x, y) = (1−2x2)e−x2−y2, ∂f
∂y(x, y) =−2xye−x2−y2. Rezultati so prikazani na sliki 17.
Interpolant se je zopet približal originalni funkciji. To nas navdaja z optimiz- mom, da je razširitev s Taylorjevim polinomom prava smer k izboljšanju interpo- lacije. Zato si poglejmo, kako bi lahko aproksimirali odvode le iz danih razpršenih točk.
3.2.1 Aproksimacija odvodov
Odvode lahko izračunamo s pomočjo aproksimacijskega polinoma na znanih točkah.
V eni dimenziji je to enostavno, saj nam izrek 2.1 zagotavlja obstoj takega polinoma.
V tem primeru aproksimiramo točke s polinomom, za vrednosti odvodov znanih točk pa vzamemo kar vrednosti odvodov aproksimanta v teh točkah. Oglejmo si enodimenzionalen primer.
Primer 3.14. Enodimenzionalne podatke lahko najprej interpoliramo s pomočjo Lagrangeevega polinoma , ki smo ga definirali v dokazu izreka 2.1. Nato pa odvajamo dobljeni polinom
p′n(x) =
n
∑
k=0
ynln,k′ (x), (3.11)
(a) Originalna funkcijaf(x, y) =xe−x2−y2. (b) Interpolacija s klasičnim Shepardovim operatorjem.
(c) Interpolacija s Shepard-Taylorjevim operatorjemST1[f].
Slika 17: Preizkus Shepard-Taylorjevega operatorja na dvodimenzionalnem primeru.
kjer se odvodi Lagrangeevih baznih polinomov izražajo kot
l′n,i(x) =
⎛
⎜
⎝
n
∏
k=0k̸=i
1 xi−xk
⎞
⎟
⎠
⎛
⎜
⎜
⎜
⎜
⎝
n
∑
k=0k̸=i n
∏
j=0 j̸=i j̸=k
x−xj
⎞
⎟
⎟
⎟
⎟
⎠ .
V (3.11) vstavimo znane točkexiin vrednosti uporabimo kot aproksimacije odvodov.
Rezultat za funkcijo f(x) =xe−x2 je prikazan na sliki 18.
(a) Lagrangeev polinom, ki interpolira znane točke.
(b) Taylor-Sheparov interpolant z uporabo odvodov Lagrangeevega polinoma.
Slika 18: Primer uporabe Taylor-Shepardovega operatorja z aproksimiranimi odvodi.
Če se premaknemo dimenzijo višje, pa se zatakne že pri problemu polinomske interpolacije, o čemer smo govorili v izreku 2.3. V [7] so se problema lotili z lokalnimi aproksimacijami ravnin v znanih točkah. Za vsako točko xi in tri njene najbližje sosede xj,xk,xl poiščemo ravninoaix+biy+ci, ki se jim najbolje prilega. Dobimo sistem enačb
aixi +biyi+ci =fi aixj +biyj+ci =fj aixk+biyk+ci =fk aixl+biyl+ci =fl, ki ga kompaktno lahko zapišemo v matrični obliki
⎡
⎢
⎢
⎣
xi yi 1 xj yj 1 xk yk 1 xl yl 1
⎤
⎥
⎥
⎦
⎡
⎣ ai bi ci
⎤
⎦=
⎡
⎢
⎢
⎣ fi fj
fk fl
⎤
⎥
⎥
⎦
(3.12)
Ker v sistemu (3.12) nastopa več enačb kot neznank, je to predoločen sistem.
Po [13, poglavje 5] problem rešimo po metodi najmanjših kvadratov. Metodo lahko
(a) Ravnina, ki se najbolje prilega štirim toč- kam.
(b) Ravnina, ki se najbolje prilega petim toč- kam.
Slika 19: Rešitev predoločenega sistema za aproksimacijo ravnine k danim točkam.
uporabimo tudi za konstrukcijo ravnine, ki se najbolje prilega petim točkam (znani točki in štirim najbližjim sosedom) tako, da v sistem dodamo še peto enačbo. Taki ravnini sta narisani na sliki 19. S pomočjo eksplicitnega izraza aproksimacijske ravnine aix+biy+c, lahko v točki xi določimo parcialna odvoda ∂f
∂x(xi) = ai in
∂f
∂y(xi) = bi. Kako se ravnina prilega grafu v izbrani točki, si lahko ogledamo na sliki 20.
Tako določimo aproksimacije odvodov vsem znanim točkam izXin nove informacije uporabimo za Taylor-Shepardov interpolant. Kako se obnese v primerjavi s klasičnim operatorjem, si lahko ogledamo na sliki 21.
Klasični operator smo s pomočjo odvodov že izboljšali. V naslednjem koraku pa nas zanima, ali bi lahko zanjšali vpliv bolj oddaljenih znanih točk na vrednosti funkcije, hkrati pa tudi naredili metodo hitrejšo.
(a) Ravnina, ki se najbolje prilega izbrani točki in najbližjim trem sosedom.
(b) Ravnina, ki se najbolje prilega izbrani točki in najbližjim štirim sosedom.
Slika 20: Določitev aproksimacijske ravnine na razpršene točke z vrednostmi funkcije f(x, y) =xe−x2−y2 in graf funkcije. Opomba: na sliki smo ravnino togo premaknili, da poteka skozi točko, katere parcialne odvode računamo iz aproksimacijske ravnine.
(a) Originalna funkcija. (b) Klasični Shepardov operator.
(c) Shepard-Taylorjev operator, odvodi iz- računani s pomočjo aproksimacijske ravnine na treh najbližjih sosedih.
(d) Shepard-Taylorjev operator, odvodi iz- računani s pomočjo aproksimacijske ravnine na štirih najbližjih sosedih.
Slika 21: Primerjava klasičnega Shepardovega operatorja s Shepard-Taylorjevim ope- ratorjem, pri čemer smo odvode v znanih točkah aproksimirali z ravninami. Apro- ksimirali smo funkcijo f(x, y) =xe−x2−y2 na 50 točkah.
3.3 Lokalna Shepardova interpolacija
Pri velikem številunznanih točk je lahko klasična Shepardova interpolacija iz defini- cije 3.1 zamudna, saj moramo za vsako neznano točko xizračunati razdalje do vseh n znanih točk. Prav tako lahko bolj oddaljene točke preveč vplivajo na vrednost interpolanta v znani točki. Zato lahko uteži Shepardovega operatorja prilagodimo tako, da upoštevajo le nekaj bližjih znanih točk, kar so predlagali tudi v [6].
Problema se lahko lotimo na dva načina. Prva možnost je, da fiksiramo radij vpliva, druga možnost pa, da fiksiramo število najbližjih sosedov, ki vplivajo na izbrano točko. V obeh primerih je metoda učinkovita, če za iskanje najbližjih sosedov uporabimo k-d drevo [15].
3.3.1 Modificirani Shepardov operator
Prva možnost je, da izberemo konstanten radijRin gledamo zaprte krogleB(x, R) = {xˆ ∈Ω| ∥x−x∥ ≤ˆ R}. Nato upoštevamo le razdalje do znanih točk, ki se nahajajo
v krogliB(x, R), torej točke iz množiceX∩B(x, R). Uteži tako ustrezno popravimo in jih zapišemo kot
W˜µ,i(x) = Wµ,i(x)
n
∑
j=1
Wµ,j(x) ,
kjer je
Wµ,i(x) =
( 1
∥x−xi∥ − 1 R
)µ +
.
V izrazu (3.13) predstavlja(·)+pozitivni del funkcije. Uteži točk zunaj radija vpliva so tako enake 0 in ne vplivajo na vrednosti točke. Dobljen interpolant je zvezen, paziti moramo le, da za vsako točkoxvelja|X∩B(x, R)| ≥1. Iskanje spodnje meje za takšen radij bi bilo zamudno, zato se raje posvetimo drugi metodi lokalizacije.
3.3.2 k-Shepardov operator
Fiksiramo število k ∈ N in nato za vsako neznano točko x izračunamo razdalje le do k najbližjih sosedov iz množice X. Ko uporabljamo razdalje do najbližjih k sosedov, za vsako točko x določimo minimalen radij Rk(x), za katerega velja
|B(x, Rk(x))∩X|=k. Nato lahko popravimo uteži in jih zapišemo kot W˜µ,ik (x) = Wµ,ik (x)
n
∑
j=1
Wµ,jk (x) ,
kjer je
Wµ,ik (x) =
( 1
∥x−xi∥− 1 Rk(x)
)µ +
. (3.13)
Z regulacijo števila sosedov k vplivamo tudi na časovno zahtevnost algoritma, saj se zmanjša število operacij, ko računamo le vrednosti neničelnih uteži. Interpolant tudi ni nujno zvezen, v skrajnem primeru pri k = 1 pa predstavlja interpolacijo z najbližjim sosedom 2.2.1. Primeri interpolacij pri različnih parametrihk so na sliki 22.
3.4 Triangulacijski Shepardov operator
Kljub temu, da smo v razdelku 2.2.2 pokazali, da triangulacija lahko vpliva na končno rešitev interpolacije, veliko metod temelji na njej. Uporabljajo lokalne po- linomske interpolante, ki bazirajo na vrednostih oglišč trikotnikov. Pri preučevanju interpolacijskih metod na razpršenih podatkih je leta 1982 F. Little [11] prilagodil Shepardov operator za triangulirana območja. Pri praktičnem preizkusu metode smo se omejili na Delaunayevo triangulacijo, več o primernih izbirah triangulacije in stopnji aproksimacije pa si lahko bralec ogleda v [3]. Ker bomo sedaj obravnavali interpolacije nad trikotniki si oglejmo baricentrične koordinate, o katerih je veliko napisano v [8].
(a)k= 1. (b)k= 8.
(c) k= 20. (d)k= 40.
Slika 22: Primerjava k-Shepardovih operatorjev na 50 točkah.
Definicija 3.15. Naj točke (x1, y1),(x2, y2),(x2, y2) določajo koordinate trikotnika T. Točki x = (x, y) lahko glede na trikotnik T določimo baricentrične koordinate λ1(x), λ2(x), λ3(x), ki jih izračunamo po formulah
λ1(x) = (y2−y3)(x−x3) + (x3−x2)(y−y3) (y2−y3)(x1−x3) + (x3−x2)(y1−y3), λ2(x) = (y3−y1)(x−x3) + (x1−x3)(y−y3)
(y2−y3)(x1−x3) + (x3−x2)(y1−y3), λ3(x) = 1−λ1(x, y)−λ3(x, y).
(3.14)
Za baricentrične koordinate velja sledeče. Naj bo x točka v trikotniku T in jo povežimo z oglišči trikotnika. Nastanejo trije novi trikotniki T1 = [x,x2,x3], T2 = [x,x1,x3], T3 = [x,x1,x2]. Baricentrične koordinate ponazarjajo razmerje ploščin malih trikotnikov proti prvotnemu trikotniku T, λi(x) = pl(Ti)
pl(T) za i = 1,2,3. Vi- dimo, da so izrazi za baricentrične koordinate (3.14) definirani za vsa realna števila x, y ∈R in jih tako posplošimo tudi na točke zunaj trikotnika.
Definicija 3.16. Naj bo T = {T1, . . . , Tm} triangulacija točk iz množice X, kjer je Tj = [xj1,xj2,xj3],xj1,xj2,xj3 ∈ X. Triangulacijske Shepardove bazne funkcije ali triangulacijske uteži, porojene s triangulacijo T, so definirane kot normaliziran produkt inverzov razdalj do oglišč trikotnikov, kar zapišemo kot
Bµ,j(x) =
3
∏
k=1
∥x−xjk∥−µ
m
∑
k=1 3
∏
l=1
∥x−xkl∥−µ
, j = 1, . . . , m, µ >0. (3.15)
Nadalje označimo linearen polinomski interpolant na trikotniku Tj, j = 1, . . . , m,
