Skip to content
Snippets Groups Projects
Commit 93921009 authored by Prof. Dr. Robert Jäschke's avatar Prof. Dr. Robert Jäschke
Browse files

moved to notebooks repo

parent 03619fd7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
programmierspaß/Distances.ipynb deleted 100644 → 0
+ 0
− 229
View file @ 03619fd7
%% Cell type:markdown id: tags:
Let us write some functions to create points in the 2D plane. We start with a well-known configuration ...
%% Cell type:code id: tags:
```
import numpy as np
def points_house():
return np.asarray([(0,0),(2,0),(0,2),(2,2),(1,3)])
```
%% Cell type:markdown id: tags:
... and visualise it in a plot:
%% Cell type:code id: tags:
```
import matplotlib.pyplot as plt
p = points_house()
plt.subplots()[1].set_aspect(1)
plt.grid()
plt.scatter(*zip(*p))
plt.show()
```
%% Cell type:markdown id: tags:
Well done. Let's add some more functions:
%% Cell type:code id: tags:
```
import random
from sklearn.datasets import make_blobs
# Randomly create count points within a width x height rectangle
def points_rectangle(count, width, height):
return np.asarray([(random.random()*width, random.random()*height) for i in range(count)])
# Randomly create counts points within a circle
def points_circle(count):
X, y = make_blobs(count, centers=[[0,0]])
return X
# Create evenly spaced points on a quadratic lattice of size k
def points_lattice(k):
return np.asarray([(i,j) for j in range(k) for i in range(k)])
r = points_rectangle(10, 2, 3)
c = points_circle(100)
l = points_lattice(5)
plt.subplots()[1].set_aspect(1)
plt.grid()
plt.scatter(*zip(*r), label="rectangle")
plt.scatter(*zip(*c), label="circle")
plt.scatter(*zip(*l), label="lattice")
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
```
```
%% Cell type:code id: tags:
```
from scipy.spatial import distance_matrix
from scipy.linalg import svd
import icp
def rmse(a, b):
return np.sqrt(((a - b) ** 2).mean())
# apply rotation R and translation t to A
def apply(A, R, t):
result = np.empty_like(A)
for i in range(np.shape(A)[0]):
result[i] = R @ A[i] + t
return result
# see https://math.stackexchange.com/questions/156161/
def get_m(d):
m = np.empty_like(d)
shp = np.shape(d)
for i in range(shp[0]):
for j in range(shp[1]):
m[i,j] = (d[1,j]**2 + d[i,1]**2 - d[i,j]**2) / 2
return m
# cf. https://math.stackexchange.com/questions/156161/
def points_from_distances(d):
# create intermediate matrices
m = get_m(d)
# print("m =", m, sep="\n")
# eigenvalue decomposition M = USU'
u, s, v = svd(m, full_matrices=True)
# print("s = ", s)
# re-estimate points
x = u @ np.sqrt(np.diag(s))
# extract points
q = x[:,0:2]
print("q = ", q, sep="\n")
return q
def plot_points(p, q, qr):
fig, ax = plt.subplots()
plt.grid()
plt.scatter(*zip(*p), label="points")
# plt.scatter(*zip(*q), label="restored points")
plt.scatter(*zip(*qr), label="restored and rotated points")
ax.legend()
ax.set_aspect(1)
plt.show()
def noise_and_restore(p, scale):
print("p =", p, sep="\n")
# measure their distances
d = distance_matrix(p, p)
print("d =", d, sep="\n")
# add noise
noise = np.random.normal(0, scale, (np.shape(d)))
print("noise =", noise, sep="\n")
d += noise
print("d =", d, sep="\n")
d_rmse = rmse(noise, np.zeros_like(noise))
print("RMSE distances:", d_rmse)
# restore points
q = points_from_distances(d)
# rotate points (https://github.com/ClayFlannigan/icp)
T, R, t = icp.best_fit_transform(q, p)
print("T =", T, sep="\n")
qr = apply(q, R, t)
print("qr =", qr, sep="\n")
return q, qr
if __name__ == '__main__':
# create random points
# p = points_rectangle(5, 2, 4)
# p = points_circle(100)
# p = points_house()
# p = points_fixed()
p = points_lattice(10)
x = []
y = []
for e in np.arange(-1, -10, -1.0):
q, qr = noise_and_restore(p, 10**e)
# compute distance between new and old points
delta = p - qr
print("p - qr =", delta, sep="\n")
# compute difference
qr_rmse = rmse(p, qr)
print("RMSE qr:", qr_rmse)
x.append(e)
y.append(qr_rmse)
# plot RMSE
plt.plot(x, y)
plt.xlabel("Gaussian scale (?)")
plt.ylabel("RMSE")
plt.grid()
plt.show()
# plot points
# plot_points(p, q, qr)
```
%% Cell type:code id: tags:
```
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Impressum | Datenschutzerklärung