{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Point Based Registration\n", "\n", "Imagine you have collected points\n", "\n", "* 3 points, of anatomical locations, in a CT scan\n", "* 3 points, of the same anatomical locations, using a tracked pointer\n", "\n", "The first 3 points use image co-ordinates, the second set of 3 points use tracker coordinates. In order to use image data in the tracker coordinate system, we must find a transformation that will map between image and tracker coordinates." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Jupyter notebook sets the cwd to the folder containing the notebook.\n", "# So, you want to add the root of the project to the sys path, so modules load correctly.\n", "import sys\n", "sys.path.append(\"../../\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# All imports for this notebook\n", "import numpy as np\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "\n", "# Note that the scikit-surgery libraries provide point-based registration using Arun's method and matrix utilities.\n", "import sksurgerycore.algorithms.procrustes as pbr\n", "import sksurgerycore.transforms.matrix as mu" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAMeElEQVR4nO3dX6icd53H8fdnkxbTuthKDsUkZZOLkqUIS+Qg1YJI65LuKqYXi3RB6YrQm1WrSCTdm94WKqIXIoRaLVhapIZa3MVYWkX2pvSkKaRNDJZq2yStPSJRkUD/+N2LM3WTY5vknHlmJt+c9wvKmfllZp7vQ9N3nzzzTCZVhSSpn7+b9QCSpNUx4JLUlAGXpKYMuCQ1ZcAlqan109zYxo0ba+vWrdPcpCS1d+DAgd9V1dzy9akGfOvWrSwsLExzk5LUXpIX3m7dUyiS1JQBl6SmDLgkNWXAJakpAy5JTZ0z4EnuTfJqkmdOW3tvkkeT/Gr088rJjilJWu58jsC/B9y0bG0P8FhVXQM8Nro/EQ8fPM71dz3Otj3/zfV3Pc7DB49PalOS1Mo5A15VvwB+v2x5F3Df6PZ9wM0DzwUsxfuOfYc4fvIUBRw/eYo79h0y4pLE6s+BX1VVL49uvwJcNdA8Z7h7/1FOvf7mGWunXn+Tu/cfncTmJKmVsd/ErKVvhHjHb4VIcluShSQLi4uLK3rtEydPrWhdktaS1Qb8t0neBzD6+eo7PbCq9lbVfFXNz839zUf5z2rTFRtWtC5Ja8lqA/4IcOvo9q3Aj4YZ50y7d25nwyXrzljbcMk6du/cPonNSVIr5/zLrJI8AHwU2JjkGHAncBfwgySfA14APjWJ4W7esRlYOhd+4uQpNl2xgd07t/91XZLWskzzS43n5+fLv41QklYmyYGqml++7icxJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKbGCniSLyd5NskzSR5I8q6hBpMknd2qA55kM/BFYL6q3g+sA24ZajBJ0tmNewplPbAhyXrgMuDE+CNJks7HqgNeVceBrwEvAi8Df6iqny5/XJLbkiwkWVhcXFz9pJKkM4xzCuVKYBewDdgEXJ7k08sfV1V7q2q+qubn5uZWP6kk6QzjnEL5GPDrqlqsqteBfcCHhxlLknQu4wT8ReC6JJclCXAjcGSYsSRJ5zLOOfAngIeAp4BDo9faO9BckqRzWD/Ok6vqTuDOgWaRJK2An8SUpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTYwU8yRVJHkryyyRHknxoqMEkSWe3fsznfxP4SVX9W5JLgcsGmEmSdB5WHfAk7wE+AvwHQFW9Brw2zFiSpHMZ5xTKNmAR+G6Sg0nuSXL58gcluS3JQpKFxcXFMTYnSTrdOAFfD3wA+HZV7QD+DOxZ/qCq2ltV81U1Pzc3N8bmJEmnGyfgx4BjVfXE6P5DLAVdkjQFqw54Vb0CvJRk+2jpRuDwIFNJks5p3KtQvgDcP7oC5Xngs+OPJEk6H2MFvKqeBuYHmkWStAJ+ElOSmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqauyAJ1mX5GCSHw8xkCTp/AxxBH47cGSA15EkrcBYAU+yBfg4cM8w40iSzte4R+DfAL4K/OWdHpDktiQLSRYWFxfH3Jwk6S2rDniSTwCvVtWBsz2uqvZW1XxVzc/Nza12c5KkZcY5Ar8e+GSS3wAPAjck+f4gU0mSzmnVAa+qO6pqS1VtBW4BHq+qTw82mSTprLwOXJKaWj/Ei1TVz4GfD/FakqTz4xG4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaWnXAk1yd5GdJDid5NsntQw4mSTq79WM89w3gK1X1VJK/Bw4kebSqDg80myTpLFZ9BF5VL1fVU6PbfwKOAJuHGkySdHaDnANPshXYATzxNr92W5KFJAuLi4tDbE6SxAABT/Ju4IfAl6rqj8t/var2VtV8Vc3Pzc2NuzlJ0shYAU9yCUvxvr+q9g0zkiTpfIxzFUqA7wBHqurrw40kSTof4xyBXw98BrghydOjf/51oLkkSeew6ssIq+p/gQw4iyRpBfwkpiQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmhor4EluSnI0yXNJ9gw1lCTp3Nav9olJ1gHfAv4ZOAY8meSRqjo81HAADx88zt37j3Li5Ck2XbGB3Tu3c/OOzUNuQpImYtL9WnXAgQ8Cz1XV8wBJHgR2AYMF/OGDx7lj3yFOvf4mAMdPnuKOfYcAjLikC9o0+jXOKZTNwEun3T82WhvM3fuP/nXn33Lq9Te5e//RITcjSYObRr8m/iZmktuSLCRZWFxcXNFzT5w8taJ1SbpQTKNf4wT8OHD1afe3jNbOUFV7q2q+qubn5uZWtIFNV2xY0bokXSim0a9xAv4kcE2SbUkuBW4BHhlmrCW7d25nwyXrzljbcMk6du/cPuRmJGlw0+jXqt/ErKo3knwe2A+sA+6tqmcHm4z/P9HvVSiSuplGv1JVg73YuczPz9fCwsLUtidJF4MkB6pqfvm6n8SUpKYMuCQ1ZcAlqSkDLklNGXBJamqqV6EkWQReWOXTNwK/G3CcDtzntWGt7fNa218Yf5//oar+5pOQUw34OJIsvN1lNBcz93ltWGv7vNb2Fya3z55CkaSmDLgkNdUp4HtnPcAMuM9rw1rb57W2vzChfW5zDlySdKZOR+CSpNMYcElqqkXAk9yU5GiS55LsmfU8k5bk6iQ/S3I4ybNJbp/1TNOQZF2Sg0l+POtZpiHJFUkeSvLLJEeSfGjWM01aki+Pfk8/k+SBJO+a9UxDS3JvkleTPHPa2nuTPJrkV6OfVw6xrQs+4EnWAd8C/gW4Fvj3JNfOdqqJewP4SlVdC1wH/Oca2GeA24Ejsx5iir4J/KSq/hH4Jy7yfU+yGfgiMF9V72fpewRume1UE/E94KZla3uAx6rqGuCx0f2xXfABBz4IPFdVz1fVa8CDwK4ZzzRRVfVyVT01uv0nlv7Dvqi/xSLJFuDjwD2znmUakrwH+AjwHYCqeq2qTs52qqlYD2xIsh64DDgx43kGV1W/AH6/bHkXcN/o9n3AzUNsq0PANwMvnXb/GBd5zE6XZCuwA3hitpNM3DeArwJ/mfUgU7INWAS+OzptdE+Sy2c91CRV1XHga8CLwMvAH6rqp7OdamquqqqXR7dfAa4a4kU7BHzNSvJu4IfAl6rqj7OeZ1KSfAJ4taoOzHqWKVoPfAD4dlXtAP7MQH+svlCNzvvuYul/XpuAy5N8erZTTV8tXbs9yPXbHQJ+HLj6tPtbRmsXtSSXsBTv+6tq36znmbDrgU8m+Q1Lp8huSPL92Y40cceAY1X11p+sHmIp6BezjwG/rqrFqnod2Ad8eMYzTctvk7wPYPTz1SFetEPAnwSuSbItyaUsvenxyIxnmqgkYenc6JGq+vqs55m0qrqjqrZU1VaW/v0+XlUX9ZFZVb0CvJTkra8ovxE4PMORpuFF4Lokl41+j9/IRf7G7WkeAW4d3b4V+NEQL7rqb6Wflqp6I8nngf0svWt9b1U9O+OxJu164DPAoSRPj9b+q6r+Z4YzaXhfAO4fHZg8D3x2xvNMVFU9keQh4CmWrrQ6yEX4sfokDwAfBTYmOQbcCdwF/CDJ51j6K7U/Nci2/Ci9JPXU4RSKJOltGHBJasqAS1JTBlySmjLgktSUAZekpgy4JDX1f+CxB7Hz8mxFAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Define 3 points, as if they were in an image\n", "image_points = np.zeros((3,3))\n", "image_points[0][0] = 0\n", "image_points[0][1] = 0\n", "image_points[0][2] = 0\n", "image_points[1][0] = 10\n", "image_points[1][1] = 0\n", "image_points[1][2] = 0\n", "image_points[2][0] = 0\n", "image_points[2][1] = 10\n", "image_points[2][2] = 0\n", "\n", "# Draw them in 2D. Its a triangle.\n", "plt.scatter(image_points[:,0], image_points[:,1])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAL1klEQVR4nO3dT2icdR7H8c9nGxdTEat0rLbKpsiSPRR3lTmoB/dQtT3I1sMeFIW6Cr0piFSsgt4WoYIIglL80x5KL6X+uWgtXryoMFo11VqFVWvS1o6UetCArXz3kKmk00xmMvPMPPNN3q9LZn7zJM/3oe3bJ0+eMY4IAQDy+VPZAwAAukPAASApAg4ASRFwAEiKgANAUiOD3NnKlStjbGxskLsEgPQ+/vjjnyKi0rw+0ICPjY2pVqsNcpcAkJ7t7+da5xIKACRFwAEgKQIOAEkRcABIioADQFIEHACSGuhthN144+CUtu8/omOnp7V6xai2bhjXXTesKXssACjdUAf8jYNT2rZvQtNnfpckTZ2e1rZ9E5JExAEseUN9CWX7/iN/xPuc6TO/a/v+IyVNBADDY6gDfuz09ILWAWApGeqAr14xuqB1AFhKhjrgWzeMa/SiZeetjV60TFs3jJc0EQAMj6H+Iea5H1RyFwoAXGioAy7NRJxgA8CFhvoSCgCgNQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiqbcBtv2r7pO1Dc7z2qO2wvbI/4wEAWunkDHynpI3Ni7avlXSHpKMFzwQA6EDbgEfE+5JOzfHSc5IekxRFDwUAaK+ra+C2N0maiojPOth2i+2a7Vq9Xu9mdwCAOSw44LaXS3pC0lOdbB8ROyKiGhHVSqWy0N0BAFro5gz8OklrJX1m+ztJ10j6xPZVRQ4GAJjfgn8jT0RMSLry3PNGxKsR8VOBcwEA2ujkNsI9kj6QNG570vaD/R8LANBO2zPwiLinzetjhU0DAOgY78QEgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiKgANAUgQcAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSbQNu+1XbJ20fmrW23fZXtj+3/brtFf0dEwDQrJMz8J2SNjatHZC0LiKul/S1pG0FzwUAaKNtwCPifUmnmtbejYizjacfSrqmD7MBAOZRxDXwByS9XcDXAQAsQE8Bt/2kpLOSds+zzRbbNdu1er3ey+4AALN0HXDb90u6U9K9ERGttouIHRFRjYhqpVLpdncAgCYj3XyS7Y2SHpP0z4j4tdiRAACd6OQ2wj2SPpA0bnvS9oOSXpB0qaQDtj+1/VKf5wQANGl7Bh4R98yx/EofZgEALADvxASApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiKgANAUgQcAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSahtw26/aPmn70Ky1K2wfsP1N4+Pl/R0TANCskzPwnZI2Nq09Lum9iPirpPcazwEAA9Q24BHxvqRTTcubJO1qPN4l6a6C5wIAtNHtNfBVEXG88fiEpFWtNrS9xXbNdq1er3e5OwBAs55/iBkRISnmeX1HRFQjolqpVHrdHQCgoduA/2j7aklqfDxZ3EgAgE50G/C3JG1uPN4s6c1ixgEAdKqT2wj3SPpA0rjtSdsPSnpG0u22v5F0W+M5AGCARtptEBH3tHhpfcGzAAAWgHdiAkBSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiKgANAUgQcAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgqZ4CbvsR21/YPmR7j+2LixoMADC/rgNue42khyVVI2KdpGWS7i5qMADA/Hq9hDIiadT2iKTlko71PhIAoBNdBzwipiQ9K+mopOOSfo6Id4saDAAwv14uoVwuaZOktZJWS7rE9n1zbLfFds12rV6vdz8pAOA8vVxCuU3StxFRj4gzkvZJuqV5o4jYERHViKhWKpUedgcAmK2XgB+VdJPt5bYtab2kw8WMBQBop5dr4B9J2ivpE0kTja+1o6C5AABtjPTyyRHxtKSnC5oFALAAvBMTAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiKgANAUgQcAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASKqngNteYXuv7a9sH7Z9c1GDAQDmN9Lj5z8v6Z2I+LftP0taXsBMAIAOdB1w25dJulXS/ZIUEb9J+q2YsQAA7fRyCWWtpLqk12wftP2y7UuaN7K9xXbNdq1er/ewOwDAbL0EfETSjZJejIgbJP0i6fHmjSJiR0RUI6JaqVR62B0AYLZeAj4paTIiPmo836uZoAMABqDrgEfECUk/2B5vLK2X9GUhUwEA2ur1LpSHJO1u3IHyP0n/6X0kAEAnegp4RHwqqVrQLACABeCdmACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACTV6//Mqu/eODil7fuP6Njpaa1eMaqtG8Z11w1ryh4LANrqd7+GOuBvHJzStn0Tmj7zuyRp6vS0tu2bkCQiDmCoDaJfQ30JZfv+I38c/DnTZ37X9v1HSpoIADoziH4NdcCPnZ5e0DoADItB9GuoA756xeiC1gFgWAyiX0Md8K0bxjV60bLz1kYvWqatG8ZbfAYADIdB9Guof4h57kI/d6EAyGYQ/XJEFPbF2qlWq1Gr1Qa2PwBYDGx/HBEX/PrKob6EAgBojYADQFIEHACSIuAAkBQBB4CkBnoXiu26pO8HtsPirJT0U9lDDNhSO+aldrwSx5zJXyKi0rw40IBnZbs21y08i9lSO+aldrwSx7wYcAkFAJIi4ACQFAHvzI6yByjBUjvmpXa8EsecHtfAASApzsABICkCDgBJEfB52F5he6/tr2wftn1z2TP1m+1HbH9h+5DtPbYvLnumotl+1fZJ24dmrV1h+4DtbxofLy9zxqK1OObtjb/bn9t+3faKMmcs2lzHPOu1R22H7ZVlzFYUAj6/5yW9ExF/k/R3SYdLnqevbK+R9LCkakSsk7RM0t3lTtUXOyVtbFp7XNJ7EfFXSe81ni8mO3XhMR+QtC4irpf0taRtgx6qz3bqwmOW7Wsl3SHp6KAHKhoBb8H2ZZJulfSKJEXEbxFxutypBmJE0qjtEUnLJR0reZ7CRcT7kk41LW+StKvxeJekuwY6VJ/NdcwR8W5EnG08/VDSNQMfrI9a/DlL0nOSHpOU/g4OAt7aWkl1Sa/ZPmj7ZduXlD1UP0XElKRnNXNmclzSzxHxbrlTDcyqiDjeeHxC0qoyhynBA5LeLnuIfrO9SdJURHxW9ixFIOCtjUi6UdKLEXGDpF+0+L6tPk/juu8mzfzHa7WkS2zfV+5Ugxcz99amPzvrlO0nJZ2VtLvsWfrJ9nJJT0h6quxZikLAW5uUNBkRHzWe79VM0Bez2yR9GxH1iDgjaZ+kW0qeaVB+tH21JDU+nix5noGwfb+kOyXdG4v/TSHXaebk5DPb32nmktEntq8qdaoeEPAWIuKEpB9sn/sV0uslfVniSINwVNJNtpfbtmaOeVH/4HaWtyRtbjzeLOnNEmcZCNsbNXMt+F8R8WvZ8/RbRExExJURMRYRY5o5Sbux8W89JQI+v4ck7bb9uaR/SPpvyfP0VeO7jb2SPpE0oZm/H4vqrceSZHuPpA8kjduetP2gpGck3W77G818J/JMmTMWrcUxvyDpUkkHbH9q+6VShyxYi2NeVHgrPQAkxRk4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkNT/ARITwYZkXY+5AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Define 3 points, as if they were in tracker space\n", "tracker_points = np.zeros((3,3))\n", "tracker_points[0][0] = 5\n", "tracker_points[0][1] = 5\n", "tracker_points[0][2] = -1000\n", "tracker_points[1][0] = 15\n", "tracker_points[1][1] = 5\n", "tracker_points[1][2] = -1000\n", "tracker_points[2][0] = 5\n", "tracker_points[2][1] = 15\n", "tracker_points[2][2] = -1000\n", "\n", "# Draw them in 2D. Its a triangle, same point order, different location.\n", "plt.scatter(tracker_points[:,0], tracker_points[:,1])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# Compute Transformation from image to tracker\n", "R, t, FRE = pbr.orthogonal_procrustes(tracker_points, image_points)\n", "T = mu.construct_rigid_transformation(R,t)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1.00000000e+00 -4.26642159e-17 0.00000000e+00]\n", " [ 3.58404070e-17 1.00000000e+00 0.00000000e+00]\n", " [ 0.00000000e+00 0.00000000e+00 1.00000000e+00]]\n", "[[ 5.]\n", " [ 5.]\n", " [-1000.]]\n", "1.4503892858778862e-15\n", "[[ 1.00000000e+00 -4.26642159e-17 0.00000000e+00 5.00000000e+00]\n", " [ 3.58404070e-17 1.00000000e+00 0.00000000e+00 5.00000000e+00]\n", " [ 0.00000000e+00 0.00000000e+00 1.00000000e+00 -1.00000000e+03]\n", " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]\n" ] } ], "source": [ "print(R)\n", "print(t)\n", "print(FRE)\n", "print(T)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1.00000000e+00 -8.53284318e-17 0.00000000e+00 -2.66453526e-15]\n", " [ 7.16808141e-17 1.00000000e+00 0.00000000e+00 -2.66453526e-15]\n", " [ 0.00000000e+00 0.00000000e+00 1.00000000e+00 9.09494702e-13]\n", " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]\n" ] } ], "source": [ "# Construct inverse\n", "R2, t2, FRE2 = pbr.orthogonal_procrustes(image_points, tracker_points)\n", "T2 = mu.construct_rigid_transformation(R2,t2)\n", "print(np.matmul(T,T2))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.44015551225972854\n" ] } ], "source": [ "# Add noise to 1 point, recalculate, look at FRE\n", "image_points[0][0] = 1\n", "R, t, FRE = pbr.orthogonal_procrustes(tracker_points, image_points)\n", "print(FRE)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 5.58210085 5.06944883 -1000. ]\n", " [ 14.57914373 4.83875543 -1000. ]\n", " [ 4.83875543 15.09179574 -1000. ]]\n" ] } ], "source": [ "# Transform image points into tracker space\n", "transformed_image_points = np.transpose(np.matmul(R, np.transpose(image_points)) + t)\n", "print(transformed_image_points)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAUu0lEQVR4nO3df2zV9b3H8dcbqMNq+SF0uq1McHHI2lJsi5QxwI2JJXWVLZcwdknK1UnUObdlKWPX7MLcckciiUHnhDoUXEyXSPjRXDdT52gwU4hF6y8kMASkKlKq/Ci1k8L7/nHaM1pKe9pzeo6f8nwkJ+d8f5zzfX8Oh1c/53O+53PM3QUACM+gVBcAAOgbAhwAAkWAA0CgCHAACBQBDgCBGpLMg40ePdrHjh2bzEMCQPB27tx51N0zO69PaoCPHTtWtbW1yTwkAATPzA52tZ4hFAAIFAEOAIEiwAEgUEkdAwdwvtOnT6u+vl4tLS2pLgUpNnToUGVlZSktLS2m/QlwIMXq6+uVkZGhsWPHysxSXQ5SxN3V2Nio+vp6jRs3Lqb7MIQCpFhLS4tGjRpFeF/kzEyjRo3q1TsxAhz4DCC8IfX+dUCAA0Cgwgrwqirpnnsi1wAS4tixY/rDH/6QkMdavny5Vq5cmZDHQs/CCfCqKmnBAumRRyLXhDiQEBcK8NbW1qTWkezjDQThBHh1tdTcHLnd3CxVV2vYMMns35dhw1JbIhCipUuXat++fZo0aZImT56s6dOnq7S0VF/72tckSXPnzlVBQYGys7NVUVERvd+zzz6r/Px85eXladasWec97mOPPaY5c+bok08+0b59+1RcXKyCggJNnz5du3fvliQtWrRId955p6ZMmaIlS5Ykp8EDibsn7VJQUOB9tmWLe3q6uxS53rLFJT/vAoRm165dvdo/I6Pjaz4jI77j79+/37Ozs93dfevWrZ6enu7vvPNOdHtjY6O7uzc3N3t2drYfPXrUjxw54llZWdH92vdZtmyZP/DAA/7www97aWmpt7S0uLv7t771Ld+zZ4+7u2/fvt2/+c1vurt7WVmZl5SUeGtra3yNGEC6ej1IqvUuMjWc88BLS6XKykhPfPbsyDJwETp5svvleN1www0dzkN+6KGHtGnTJknSoUOHtHfvXjU0NGjGjBnR/a644oro/k8++aTGjBmjzZs3Ky0tTU1NTXrxxRc1b9686D7/+te/orfnzZunwYMHJ7YRF4lwAlyKhDbBDfSryy67LHq7pqZGf/vb3/TSSy8pPT1dN954Y4/nKefm5qquri76hZSzZ89qxIgRqqur6/F46J1wxsC7kJHR/TKAnmVkZOjkBbrxx48f18iRI5Wenq7du3dr+/btkqSioiJt27ZN+/fvlyR99NFH0ftcf/31WrNmjUpLS/X+++9r2LBhGjdunJ5++mlJkWHb1157rZ9bdXEIOsBPnOg4An7iRKorAvpfojsuo0aN0rRp05STk6Py8vIO24qLi9Xa2qoJEyZo6dKlKioqkiRlZmaqoqJC3/ve95SXl6f58+d3uN83vvENrVy5UiUlJTp69KieeuoprV27Vnl5ecrOztaWLVviKxqSJIuMjydHYWGh84MOQEdvv/22JkyYkOoy8BnR1evBzHa6e2HnfYPugQPAxYwAB4BAEeAAECgCHAACRYADQKAIcAAIFAEOXOQSOZ1sTxYsWKCJEyfqwQcfTMrxOqupqdEtt9xy3vra2lrde++9KagoPmF9lR5AwrUH+N13333ettbWVg0ZkpiYOHz4sF5++WX985//jPk+iTx+dwoLC1VYeN5p1p959MCBi9y508mWl5erpqYm5illL7/8ct13333Ky8tTUVGRPvzwQ0nS008/rZycHOXl5WnGjBmSpNmzZ+u9997TpEmT9MILL6iurk5FRUWaOHGivvvd7+rjjz+WJN1444366U9/qsLCQq1atUqLFi3SXXfdpaKiIl1zzTWqqanRbbfdpgkTJmjRokXRWqqrqzV16lTl5+dr3rx5ampqkhSZ9va6665Tfn6+Nm7c2OVzcG7PfPny5SorK9P06dN19dVXa+PGjVqyZIlyc3NVXFys06dPS5Luv/9+TZ48WTk5OVq8eLHavxT58ssva+LEidHnMycnR5J05swZlZeXa/LkyZo4caLWrFkT/z9eV1MU9tclrulkgQGqt9PJuntkeuUf/ShyHadzp5N1j31KWXd3SV5VVeXu7uXl5f6b3/zG3d1zcnK8vr7e3d0//vjjLo+Tm5vrNTU17u7+q1/9yn/yk5+4u/vMmTP9rrvuiu5XVlbm8+fP97Nnz/rmzZs9IyPDX3/9dT9z5ozn5+f7q6++6g0NDT59+nRvampyd/cVK1b4r3/9a//kk088KyvL9+zZ42fPnvV58+Z5SUnJec/B1q1bo+uXLVvm06ZN808//dTr6ur80ksv9b/85S/u7j537lzftGlTh+fE3X3hwoXR5yE7O9tffPFFd3f/xS9+EW3zmjVros9PS0uLFxQUdHiO2/VmOtkee+Bm9riZHTGzN7vY9nMzczMbHf+fEgAxScKvU3U1pWx7L7t9SllJuuSSS6I914KCAh04cECSNG3aNC1atEiPPfaYzpw5c97jHz9+XMeOHdPMmTMlSWVlZdq2bVt0e+e5Vb7zne/IzJSbm6srr7xSubm5GjRokLKzs3XgwAFt375du3bt0rRp0zRp0iStX79eBw8e1O7duzVu3Dhde+21MjMtXLgwpvbPmTNHaWlpys3N1ZkzZ1RcXCwpMtNiexu3bt2qKVOmKDc3V3//+9/11ltv6dixYzp58qSmTp0qSfrBD34Qfczq6mo9+eSTmjRpkqZMmaLGxsbo89hXsQwurZP0e0lPnrvSzMZImi3p3bgqANA7Xfw6VaKnWY51Stm0tLToL6kPHjw4+rNoq1ev1o4dO/TMM8+ooKBAO3fu7PPxJelzn/ucJGnQoEHR2+3Lra2tGjx4sG666SZVVlZ2uN+FprDtybnHO7eN7cdraWnR3XffrdraWo0ZM0bLly/vcZpdd9fDDz+sm2++uU81daXHHri7b5P0URebHpS0RFLyZsMCEPlBk/T0yO309MhyHLqbTla68JSy3dm3b5+mTJmi+++/X5mZmTp06FCH7cOHD9fIkSP1wgsvSJL+9Kc/RXvjfVFUVKR//OMf0Q9IT506pT179ui6667TgQMHtG/fPkk6L+D7qj2sR48eraamJm3YsEGSNGLECGVkZGjHjh2SpD//+c/R+9x888169NFHo2Poe/bs0alTp+Kqo08f75rZrZLec/fX2v8ydbPvYkmLJenLX/5yXw4H4FwJ/nWqc6eTnTNnjkpKSjpsLy4u1urVqzVhwgSNHz8+OqVsd8rLy7V37165u2bNmqW8vDwdPHiwwz7r16/XnXfeqebmZl1zzTV64okn+tyGzMxMrVu3TgsWLIj+2s9vf/tbffWrX1VFRYVKSkqUnp6u6dOnd/vHKlYjRozQHXfcoZycHF111VWaPHlydNvatWt1xx13aNCgQZo5c6aGDx8uSfrhD3+oAwcOKD8/X+6uzMxMbd68Oa46YppO1szGSvo/d88xs3RJWyXNdvfjZnZAUqG7H+3pcZhOFjgf08kOLE1NTbr88sslSStWrNAHH3ygVatWxXz/3kwn25ce+FckjZPU3vvOkvSKmd3g7of78HgAMGA888wz+t3vfqfW1lZdffXVWrduXb8dq9cB7u5vSPp8+3JveuAAMNDNnz//vLNo+ksspxFWSnpJ0ngzqzez2/u/LODiEstQJga+3r4OeuyBu/uCHraP7dURAXQwdOhQNTY2atSoUerppAAMXO6uxsZGDR06NOb7MBcKkGJZWVmqr69XQ0NDqktBig0dOlRZWVkx70+AAymWlpbW4VuPQKyYzAoAAkWAA0CgCHAACBQBDgCBIsABIFAEOAAEigAHgEAR4AAQKAIcAAJFgANAoAhwAAgUAQ4AgSLAASBQBDgABIoAB4BAEeAAECgCHAACRYADQKAIcAAIFAEOAIEiwAEgUAQ4AASKAAeAQBHgABAoAhwAAkWAA0CgCHAACFSPAW5mj5vZETN785x1D5jZbjN73cw2mdmI/i0TANBZLD3wdZKKO617TlKOu0+UtEfSLxNcFwCgBz0GuLtvk/RRp3XV7t7atrhdUlY/1AYA6EYixsBvk/TXC200s8VmVmtmtQ0NDQk4HABAijPAzew+Sa2SnrrQPu5e4e6F7l6YmZkZz+EAAOcY0tc7mtkiSbdImuXunrCKAAAx6VOAm1mxpCWSZrp7c2JLAgDEIpbTCCslvSRpvJnVm9ntkn4vKUPSc2ZWZ2ar+7lOAEAnPfbA3X1BF6vX9kMtAIBe4JuYABAoAhwAAkWAA0CgCHAACBQBDgCBIsABIFAEOAAEigAHgEAR4AAQKAIcAAJFgANAoAhwAAgUAQ4AgSLAASBQBDgABIoAB4BAEeAAECgCHAACRYADQKAIcAAIFAEOAIEiwAEgUAQ4AASKAAeAQBHgABAoAhwAAkWAA0CgCHAACBQBDgCB6jHAzexxMztiZm+es+4KM3vOzPa2XY/s3zIBAJ3F0gNfJ6m407qlkp5392slPd+2DABIoh4D3N23Sfqo0+pbJa1vu71e0twE1wUA6EFfx8CvdPcP2m4flnTlhXY0s8VmVmtmtQ0NDX08HACgs7g/xHR3l+TdbK9w90J3L8zMzIz3cACANn0N8A/N7AuS1HZ9JHElAQBi0dcAr5JU1na7TNKWxJQDAIhVLKcRVkp6SdJ4M6s3s9slrZB0k5ntlfTttmUAQBIN6WkHd19wgU2zElwLAKAX+CYmAASKAAeAQBHgABAoAhwAAkWAA0CgCHAACBQBDgCBIsABIFAEOAAEigAHgEAR4AAQKAIcAAJFgANAoAhwAAgUAQ4AgSLAASBQBDgABIoAB4BAEeAAECgCHAACRYADQKAIcAAIFAEOAIEiwAEgUAQ4AASKAAeAQBHgABAoAhwAAkWAA0Cg4gpwM/uZmb1lZm+aWaWZDU1UYQCA7vU5wM3sS5LulVTo7jmSBkv6fqIKAwB0L94hlCGSLjWzIZLSJb0ff0kAgFj0OcDd/T1JKyW9K+kDScfdvbrzfma22Mxqzay2oaGh75UCADqIZwhlpKRbJY2T9EVJl5nZws77uXuFuxe6e2FmZmbfKwUAdBDPEMq3Je139wZ3Py1po6SvJ6YsAEBP4gnwdyUVmVm6mZmkWZLeTkxZAICexDMGvkPSBkmvSHqj7bEqElQXAKAHQ+K5s7svk7QsQbUAAHqBb2ICQKAIcAAIFAEOAIEiwAEgUAQ4AASKAAeAQBHgABAoAhwAAkWAA0CgCHAACBQBDgCBIsABIFAEOAAEigAHgEAR4AAQKAIcAAJFgANAoAhwAAgUAQ4AgSLAASBQBDgABIoAB4BAEeAAECgCHAACRYADQKAIcAAIFAEOAIEiwAEgUAQ4AAQqrgA3sxFmtsHMdpvZ22Y2NVGFAQC6NyTO+6+S9Ky7/4eZXSIpPQE1AQBi0OcAN7PhkmZIWiRJ7v6ppE8TUxYAoCfxDKGMk9Qg6Qkze9XM/mhml3XeycwWm1mtmdU2NDTEcTgAwLniCfAhkvIlPeru10s6JWlp553cvcLdC929MDMzM47DAQDOFU+A10uqd/cdbcsbFAl0AEAS9DnA3f2wpENmNr5t1SxJuxJSFQCgR/GehfJjSU+1nYHyjqT/ir8kAEAs4gpwd6+TVJigWgAAvcA3MQEgUAQ4AASKAAeAQBHgABAoAhwAAkWAA0CgCHAACBQBDgCBIsABIFAEOAAEigCXpKoq6Z57ItcAEIigA3zYMMns35dhw/rwIFVV0oIF0iOPRK4JcQAJkpCM6kbQAX7yZPfLMamulpqbI7ebmyPLAJAACcmobgQd4Akxe7aU3vZbzOnpkWUACEC884GHr7RUqqyM9Lxnz44sA8CFVFV9ZvLC3D1pByssLPTa2tqEPd6wYR3fkmRkSCdOJOzhAaCj9s/Mmpsj79grK7sN8URllJntdPfzfnsh6CGUEyck939fCG8A/aqXn5n1d0YFHeAAkFSfsc/MGAMHgFh9xj4zI8ABoDdKS1Me3O0YQgGAQBHgABAoAhwAAkWAA0CgCHAACBQBDgCBIsABIFAEOAAEKqmTWZlZg6SDSTtg/xot6Wiqi0iii6m9F1NbJdobgqvdPbPzyqQG+EBiZrVdzQ42UF1M7b2Y2irR3pAxhAIAgSLAASBQBHjfVaS6gCS7mNp7MbVVor3BYgwcAAJFDxwAAkWAA0CgCPA+MLMRZrbBzHab2dtmNjXVNfUXM/uZmb1lZm+aWaWZDU11TYlkZo+b2REze/OcdVeY2XNmtrftemQqa0ykC7T3gbbX8utmtsnMRqSyxkTpqq3nbPu5mbmZjU5FbYlCgPfNKknPuvt1kvIkvZ3ievqFmX1J0r2SCt09R9JgSd9PbVUJt05Scad1SyU97+7XSnq+bXmgWKfz2/ucpBx3nyhpj6RfJruofrJO57dVZjZG0mxJ7ya7oEQjwHvJzIZLmiFprSS5+6fufiy1VfWrIZIuNbMhktIlvZ/iehLK3bdJ+qjT6lslrW+7vV7S3KQW1Y+6aq+7V7t7a9vidklZSS+sH1zg31aSHpS0RFLwZ3AQ4L03TlKDpCfM7FUz+6OZXZbqovqDu78naaUiPZUPJB139+rUVpUUV7r7B223D0u6MpXFJNltkv6a6iL6i5ndKuk9d38t1bUkAgHee0Mk5Ut61N2vl3RKA+stdlTb2O+tivzR+qKky8xsYWqrSi6PnGcbfE8tFmZ2n6RWSU+lupb+YGbpkv5b0v+kupZEIcB7r15SvbvvaFveoEigD0TflrTf3Rvc/bSkjZK+nuKakuFDM/uCJLVdH0lxPf3OzBZJukXSf/rA/XLIVxTpjLxmZgcUGSp6xcyuSmlVcSDAe8ndD0s6ZGbj21bNkrQrhSX1p3clFZlZupmZIm0dkB/YdlIlqaztdpmkLSmspd+ZWbEiY8Kl7t6c6nr6i7u/4e6fd/ex7j5Wkc5Yftv/6SAR4H3zY0lPmdnrkiZJ+t8U19Mv2t5lbJD0iqQ3FHm9DJivIUuSmVVKeknSeDOrN7PbJa2QdJOZ7VXkXciKVNaYSBdo7+8lZUh6zszqzGx1SotMkAu0dUDhq/QAECh64AAQKAIcAAJFgANAoAhwAAgUAQ4AgSLAASBQBDgABOr/ARVv5jNfpHcKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot tracker points and transformed image points\n", "fig = plt.figure()\n", "ax1 = fig.add_subplot(111)\n", "ax1.scatter(tracker_points[:,0], tracker_points[:,1], s=10, c='b', marker=\"s\", label='tracker')\n", "ax1.scatter(transformed_image_points[:,0], transformed_image_points[:,1], s=10, c='r', marker=\"o\", label='transformed image')\n", "plt.legend(loc='upper right');\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "mphy0026", "language": "python", "name": "mphy0026" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }