Puntwolken.ipynb

Kernel: Python 3 (Ubuntu Linux)

Interpolatie in meer dimensies

We zoeken een benadering van een gegeven functie $f : \mathbb{R}^d \rightarrow \mathbb{R}$ en datapunten $y_i = f(\mathbf{x}_i)$ voor $i=0, 1, \ldots m$

$$p(\mathbf{x}) = \sum_{i=1}^{m} c_i \phi_i(\mathbf{x}).$$

We hebben eerder al gekeken naar interpolatie van functies met behulp van polynomen in het geval dat $d = 1$, waar de $\phi_i$ polynomen van graad $\leq n$ waren. Als we precies $n + 1$ punten hebben (dus $m = n + 1$) kunnen we de coefficienten $c_i$ verkrijgen door een stelsel vergelijkingen op te lossen.

Als we dit willen uitbreiden naar hogere dimensies, moeten we geschikte basisfuncties vinden. Op basis van onze ervaring met lineaire interpolatie, kunnen we de interpolant in twee dimensies weergeven als

$$p(x,y) = c_0 + c_1x + c_2y + c_3xy.$$

We noemen dit bilineaire interpolatie.

Er zijn echter ook nog andere mogelijkheden voor de basisfuncties. Een daarvan zijn de zogenaamde radiale basisfuncties

$$\phi_i = \rho(\|x - x_i\|_2),$$

met

$$\rho(r) = e^{-(a r)^2},$$

of

$$\rho(r) = (1 + (a r)^2)^{-1},$$

waarbij $a$ een te kiezen schalingsparameter is. Ook hier kunnen we een stelsel van $m$ vergelijkingen met $m$ onbekenden opstellen om de coefficiënten $c$ te vinden zodat $f(x_i) = p(x_i)$.

Puntwolken

We kunnen de interpolatie ook gebruiken om oppervlakken in drie dimensies te construeren uit een puntwolk. We hebben in dat geval $m$ punten $\mathbf{x}_i$ op het oppervlak een 3D lichaam. Zulke metingen komen vaak uit een laserscan. Om het oorspronkelijke lichaam te reconstrueren voor visualisatie gaan we als volgt te werk. Allereerst hebben op ieder punt ook de normaalvectoren $\mathbf{n}_i$ nodig. We zoeken nu een functie $p(\mathbf{x})$ zodat

$$p(\mathbf{x}_i) = 0,$$$$p(\mathbf{x}_i + \delta\mathbf{n}_i) = 1,$$$$p(\mathbf{x}_i - \delta\mathbf{n}_i) = -1,$$

waar $\delta$ een nader te kiezen parameter is, en definieren vervolgens het gereconstrueerde lichaam als de contour waar $p(\mathbf{x}) = 0$. We zoeken een functie $f$ van de vorm

$$p(\mathbf{x}) = \sum_{i=1}^{3m} c_i \phi(\|\mathbf{x} - \widetilde{\mathbf{x}}_i\|_2),$$

waar $\phi(r)$ een radiale basisfunctie is en $\{\widetilde{\mathbf{x}}_j\}$ voor $j=1, 2, \ldots, 3m$ de set punten $\{\mathbf{x}_i, \mathbf{x}_i + \delta\mathbf{n}_i, \mathbf{x}_i - \delta\mathbf{n}_i\}$ voor $i=1,2,\ldots,m$ is.

We kunnen nu een stelsel van $3m$ vergelijkingen met $3m$ onbekenden $A\mathbf{c} = \mathbf{b}$ opstellen om de coefficienten $\mathbf{c}$ te vinden.

Vragen voor het verslag

Basis

• Voor een regelmatig rooster is kunnen we bilineare interpolatie zien als lineare interpolatie, toegepast in de twee richtingen. Geef in dit geval een uitdrukking voor de interpolatiefout en druk de interpolant uit in termen van Lagrange polynomen.

• Breid voorgaande uit naar interpolatie van willekeurige graad in 2 dimensies.

• Bekijk eens in 1D hoe de radiale basisfuncties het doen als je een gegeven functie $f(x)$ wilt benaderen (test bijvoorbeeld op $f(x) = x^k$ voor verschillende waarden van $k$). Gebruik in dat geval een regelmatig rooster van $m$ punten $x_i \in [0,1]$ en stel een stelsel vergelijkingen op zodat $f(x_i) = p(x_i)$, waar $p(x) = \sum_{i=0}^{m-1} c_i \rho(|x-x_i|)$.

• Test de precisie (voor welke orde polynoom de methode exact is) en orde (hoe snell de fout naar nul gaat als functie van $h = 1/(m-1)$) van de interpolatie numeriek door de methode.

• Beschrijf hoe de matrix $A$ eruit ziet en onderzoek de relevante eigenschappen als functie van $a$.

• Bespreek mogelijke voordelen van interpolatie met radiële basisfuncties ten opzichte van globale of stuksgewijs polynomiale interpolatie.

Extra

• Test de methode uit in 2D op de bijgeleverde puntwolkendataset, probeer goede waarden voor de paramaters $a$, $\delta$ te vinden

• Zou je ook polynomen als basisfuncties kunnen gebruiken? Waarom wel of niet?

In de onderstaande code gaan we een bi-lineaire interpolant opstellen met steunpunten $(0,0)$, $(0,1)$, $(1,0)$ en $(1,1)$ en functiewaarden $1$, $2$, $3$, $4$. We doen dit door de coefficienten van het polynoom $$p(x,y) = c_0 + c_1 x c_2 y + c_3 xy$$ te bepalen.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# lineare interpolatie in 2 dimensies

# punten
xi = [0,0,1,1]
yi = [0,1,0,1]

# functiewaarden
fi = [1,2,3,4]

# basisfuncties
p0 = lambda x,y : 1
p1 = lambda x,y : x
p2 = lambda x,y : y
p3 = lambda x,y : x*y

# stel vergelijkingen op
A = np.array([[p0(xi[0],yi[0]), p1(xi[0],yi[0]), p2(xi[0],yi[0]), p3(xi[0],yi[0])], [p0(xi[1],yi[1]), p1(xi[1],yi[1]), p2(xi[1],yi[1]), p3(xi[1],yi[1])], [p0(xi[2],yi[2]), p1(xi[2],yi[2]), p2(xi[2],yi[2]), p3(xi[2],yi[2])], [p0(xi[3],yi[3]), p1(xi[3],yi[3]), p2(xi[3],yi[3]), p3(xi[3],yi[3])]])

# los op
c = np.linalg.solve(A,fi)

# interpolant
p = lambda x,y : c[0]*p0(x,y) + c[1]*p1(x,y) + c[2]*p2(x,y) + c[3]*p3(x,y)

# grafiek
xx,yy = np.meshgrid(np.linspace(0,1,10),np.linspace(0,1,10))

fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot_surface(xx,yy,p(xx,yy))

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f9f4acbad30>

In de onderstaande code gaan we een interpolant opstellen met behulp van radiale basisfuncties

# steunpunten
m = 3
xi = np.linspace(0,1,m)

# te interpoleren functie
f = lambda x : x**2 - x

# basisfunctie
phi = lambda r : np.exp(-r**2)

# stel matrix op
A = np.zeros((m,m))

for i in range(m):
    for j in range(m):
        A[i,j] = phi(np.abs(xi[i] - xi[j]))

# coefficienten
c = np.linalg.solve(A,f(xi))

# interpolant
def p(x):
    y = 0*x
    for i in range(m):
        y = y + c[i]*phi(np.abs(x - xi[i]))
    return y

# plot
x = np.linspace(0,1,100)
plt.plot(x,f(x),x,p(x))

# fout
np.linalg.norm(f(x) - p(x))

0.086420625275182614

Hieronder vind je code om een 2d puntwolken dataset in te lezen

X = np.loadtxt('points_2D.dat',delimiter=',')
N = np.loadtxt('normals_2D.dat',delimiter=',')

plt.plot(X[:,0],X[:,1],'*')
plt.quiver(X[:,0],X[:,1],-N[:,0],-N[:,1])
plt.xlim(-1,1)
plt.ylim(-1,1)

(-1, 1)