diff --git "a/programmierspa\303\237/Distances.ipynb" "b/programmierspa\303\237/Distances.ipynb"
new file mode 100644
index 0000000000000000000000000000000000000000..21135cce01e977fbed02684335dd30af57721631
--- /dev/null
+++ "b/programmierspa\303\237/Distances.ipynb"
@@ -0,0 +1,229 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let us write some functions to create points in the 2D plane. We start with a well-known configuration ..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "\n",
+ "def points_house():\n",
+ " return np.asarray([(0,0),(2,0),(0,2),(2,2),(1,3)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "... and visualise it in a plot:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "p = points_house()\n",
+ "\n",
+ "plt.subplots()[1].set_aspect(1)\n",
+ "plt.grid()\n",
+ "plt.scatter(*zip(*p))\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Well done. Let's add some more functions:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import random\n",
+ "\n",
+ "from sklearn.datasets import make_blobs\n",
+ "\n",
+ "# Randomly create count points within a width x height rectangle \n",
+ "def points_rectangle(count, width, height):\n",
+ " return np.asarray([(random.random()*width, random.random()*height) for i in range(count)])\n",
+ "\n",
+ "# Randomly create counts points within a circle \n",
+ "def points_circle(count):\n",
+ " X, y = make_blobs(count, centers=[[0,0]])\n",
+ " return X\n",
+ "\n",
+ "# Create evenly spaced points on a quadratic lattice of size k\n",
+ "def points_lattice(k):\n",
+ " return np.asarray([(i,j) for j in range(k) for i in range(k)])\n",
+ "\n",
+ "r = points_rectangle(10, 2, 3)\n",
+ "c = points_circle(100)\n",
+ "l = points_lattice(5)\n",
+ "\n",
+ "plt.subplots()[1].set_aspect(1)\n",
+ "plt.grid()\n",
+ "plt.scatter(*zip(*r), label=\"rectangle\")\n",
+ "plt.scatter(*zip(*c), label=\"circle\")\n",
+ "plt.scatter(*zip(*l), label=\"lattice\")\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "from scipy.spatial import distance_matrix\n",
+ "from scipy.linalg import svd\n",
+ "import icp\n",
+ "\n",
+ "def rmse(a, b):\n",
+ " return np.sqrt(((a - b) ** 2).mean())\n",
+ "\n",
+ "# apply rotation R and translation t to A\n",
+ "def apply(A, R, t):\n",
+ " result = np.empty_like(A)\n",
+ " for i in range(np.shape(A)[0]):\n",
+ " result[i] = R @ A[i] + t\n",
+ " return result\n",
+ "\n",
+ "# see https://math.stackexchange.com/questions/156161/\n",
+ "def get_m(d):\n",
+ " m = np.empty_like(d)\n",
+ " shp = np.shape(d)\n",
+ " for i in range(shp[0]):\n",
+ " for j in range(shp[1]):\n",
+ " m[i,j] = (d[1,j]**2 + d[i,1]**2 - d[i,j]**2) / 2\n",
+ " return m\n",
+ "\n",
+ "# cf. https://math.stackexchange.com/questions/156161/\n",
+ "def points_from_distances(d):\n",
+ " # create intermediate matrices\n",
+ " m = get_m(d)\n",
+ " # print(\"m =\", m, sep=\"\\n\")\n",
+ " # eigenvalue decomposition M = USU'\n",
+ " u, s, v = svd(m, full_matrices=True)\n",
+ " # print(\"s = \", s)\n",
+ " # re-estimate points\n",
+ " x = u @ np.sqrt(np.diag(s))\n",
+ "\n",
+ " # extract points\n",
+ " q = x[:,0:2]\n",
+ " print(\"q = \", q, sep=\"\\n\")\n",
+ "\n",
+ " return q\n",
+ "\n",
+ "def plot_points(p, q, qr):\n",
+ " fig, ax = plt.subplots()\n",
+ " plt.grid()\n",
+ " plt.scatter(*zip(*p), label=\"points\")\n",
+ " # plt.scatter(*zip(*q), label=\"restored points\")\n",
+ " plt.scatter(*zip(*qr), label=\"restored and rotated points\")\n",
+ " ax.legend()\n",
+ " ax.set_aspect(1)\n",
+ " plt.show()\n",
+ " \n",
+ "def noise_and_restore(p, scale):\n",
+ " print(\"p =\", p, sep=\"\\n\")\n",
+ " \n",
+ " # measure their distances\n",
+ " d = distance_matrix(p, p)\n",
+ " print(\"d =\", d, sep=\"\\n\")\n",
+ " \n",
+ " # add noise\n",
+ " noise = np.random.normal(0, scale, (np.shape(d)))\n",
+ " print(\"noise =\", noise, sep=\"\\n\")\n",
+ " d += noise\n",
+ " print(\"d =\", d, sep=\"\\n\")\n",
+ " d_rmse = rmse(noise, np.zeros_like(noise))\n",
+ " print(\"RMSE distances:\", d_rmse)\n",
+ "\n",
+ " # restore points\n",
+ " q = points_from_distances(d)\n",
+ "\n",
+ " # rotate points (https://github.com/ClayFlannigan/icp)\n",
+ " T, R, t = icp.best_fit_transform(q, p)\n",
+ " print(\"T =\", T, sep=\"\\n\")\n",
+ " qr = apply(q, R, t)\n",
+ " print(\"qr =\", qr, sep=\"\\n\")\n",
+ "\n",
+ " return q, qr\n",
+ " \n",
+ "if __name__ == '__main__':\n",
+ " # create random points\n",
+ " # p = points_rectangle(5, 2, 4)\n",
+ " # p = points_circle(100)\n",
+ " # p = points_house()\n",
+ " # p = points_fixed()\n",
+ " p = points_lattice(10)\n",
+ "\n",
+ " x = []\n",
+ " y = []\n",
+ " for e in np.arange(-1, -10, -1.0):\n",
+ " q, qr = noise_and_restore(p, 10**e)\n",
+ " \n",
+ " # compute distance between new and old points\n",
+ " delta = p - qr\n",
+ " print(\"p - qr =\", delta, sep=\"\\n\")\n",
+ "\n",
+ " # compute difference\n",
+ " qr_rmse = rmse(p, qr)\n",
+ " print(\"RMSE qr:\", qr_rmse)\n",
+ "\n",
+ " x.append(e)\n",
+ " y.append(qr_rmse)\n",
+ "\n",
+ " # plot RMSE\n",
+ " plt.plot(x, y)\n",
+ " plt.xlabel(\"Gaussian scale (?)\")\n",
+ " plt.ylabel(\"RMSE\")\n",
+ " plt.grid()\n",
+ " plt.show()\n",
+ " \n",
+ " # plot points\n",
+ " # plot_points(p, q, qr)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "language_info": {
+ "name": "python",
+ "pygments_lexer": "ipython3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}