{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [ "meta", "draft" ] }, "source": [ "# Euclidean geometry" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "tags": [ "hide" ] }, "outputs": [], "source": [ "# Import directives\n", "\n", "#%pylab notebook\n", "%pylab inline\n", "\n", "pylab.rcParams['figure.figsize'] = (6, 6)\n", "\n", "#import warnings\n", "#warnings.filterwarnings('ignore')\n", "\n", "import math\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import axes3d\n", "\n", "from ipywidgets import interact" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2D transformations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def plot2d(x, fmt=\"ok\"):\n", " plt.axis('equal')\n", "\n", " plt.axis([-5, 5, -5, 5])\n", " plt.xticks(np.arange(-5, 5, 1))\n", " plt.yticks(np.arange(-5, 5, 1))\n", "\n", " plt.axhline(y=0, color='k')\n", " plt.axvline(x=0, color='k')\n", "\n", " plt.plot(x[:,0], x[:,1], fmt)\n", "\n", " plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rotation around the origin" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "R(\\theta) =\n", "\\begin{pmatrix}\n", " \\cos \\theta & -\\sin \\theta \\\\\n", " \\sin \\theta & \\cos \\theta\n", "\\end{pmatrix}\n", "$$\n", "\n", "\n", "\n", "$$\n", "\\begin{pmatrix}\n", " A'_{x} \\\\\n", " A'_{y}\n", "\\end{pmatrix}\n", "=\n", "\\begin{pmatrix}\n", " \\cos \\theta & -\\sin \\theta \\\\\n", " \\sin \\theta &\\cos \\theta \\\\\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", " A_{x} \\\\\n", " A_{y}\n", "\\end{pmatrix}\n", "$$\n", "\n", "See: https://fr.wikipedia.org/wiki/Matrice_de_rotation" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define initial points\n", "A = np.array([[0., 0.],\n", " [1., 0.],\n", " [1., 1.],\n", " [0., 0.]])\n", "\n", "# Define the rotation angle\n", "theta = np.radians(30)\n", "\n", "# Define the rotation matrix\n", "R = np.array([[np.cos(theta), -np.sin(theta)],\n", " [np.sin(theta), np.cos(theta)]])\n", "\n", "# Rotate points\n", "Aprime = np.dot(R, A.T).T\n", "\n", "# Print and plot\n", "print(A)\n", "print(Aprime)\n", "\n", "plot2d(A, fmt=\"-ok\")\n", "plot2d(Aprime, fmt=\"-or\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3D transformations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def plot3d(x, axis=None, fmt=\"ok\"):\n", " if axis is None:\n", " fig = plt.figure()\n", " axis = axes3d.Axes3D(fig)\n", " \n", " axis.scatter(x[:,0], x[:,1], x[:,2], fmt)\n", " axis.plot(x[:,0], x[:,1], x[:,2], fmt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rotation around the x axis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "R_{\\mathbf {x}}(\\theta) =\n", "\\begin{pmatrix}\n", " 1& 0 & 0 \\\\\n", " 0& \\cos \\theta & -\\sin \\theta \\\\\n", " 0& \\sin \\theta & \\cos \\theta\n", "\\end{pmatrix}\n", "$$\n", "\n", "See: https://fr.wikipedia.org/wiki/Matrice_de_rotation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rotation around the y axis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "R_{\\mathbf {y}}(\\theta )\n", "=\n", "\\begin{pmatrix}\n", " \\cos \\theta & 0 & \\sin \\theta \\\\\n", " 0 & 1 & 0 \\\\\n", " -\\sin \\theta & 0 & \\cos \\theta\n", "\\end{pmatrix}\n", "$$\n", "\n", "See: https://fr.wikipedia.org/wiki/Matrice_de_rotation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rotation around the z axis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "R_{\\mathbf {z}}(\\theta )\n", "=\n", "\\begin{pmatrix}\n", " \\cos \\theta & -\\sin \\theta & 0 \\\\\n", " \\sin \\theta & \\cos \\theta & 0 \\\\\n", " 0 & 0 & 1\n", "\\end{pmatrix}\n", "$$\n", "\n", "See: https://fr.wikipedia.org/wiki/Matrice_de_rotation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define initial points\n", "A = np.array([[0., 0., 0.],\n", " [1., 0., 0.],\n", " [0., 0., 0.],\n", " [0., 1., 0.],\n", " [0., 0., 0.],\n", " [0., 0., 1.]])\n", "\n", "# Define the rotation angle\n", "theta = np.radians(90)\n", "\n", "# Define the rotation matrices\n", "Rx = np.array([[1., 0., 0.],\n", " [0., np.cos(theta), -np.sin(theta)],\n", " [0., np.sin(theta), np.cos(theta)]])\n", "\n", "Ry = np.array([[np.cos(theta), 0., np.sin(theta)],\n", " [0., 1., 0. ],\n", " [-np.sin(theta), 0., np.cos(theta)]])\n", "\n", "Rz = np.array([[np.cos(theta), -np.sin(theta), 0.],\n", " [np.sin(theta), np.cos(theta), 0.],\n", " [0., 0., 1.]])\n", "\n", "# Rotate points\n", "Ax = np.dot(Rx, A.T).T\n", "Ay = np.dot(Ry, A.T).T\n", "Az = np.dot(Rz, A.T).T\n", "\n", "# Plot\n", "fig = plt.figure()\n", "ax = axes3d.Axes3D(fig)\n", "plot3d(A, axis=ax, fmt=\"-ok\")\n", "plot3d(Ax, axis=ax, fmt=\":or\")\n", "plot3d(Ay, axis=ax, fmt=\":og\")\n", "plot3d(Az, axis=ax, fmt=\":ob\")\n", "ax.text(1, 0, 0, \"x\", color=\"r\")\n", "ax.text(0, 1, 0, \"y\", color=\"g\")\n", "ax.text(0, 0, 1, \"z\", color=\"b\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rotation around a given axis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rotation around the axis defined by the unit vector $\\overrightarrow{u} \\begin{pmatrix} u_{x} \\\\ u_{y} \\\\ u_{z} \\end{pmatrix}$ (i.e. with $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}=1$) by an angle $\\theta$.\n", "\n", "$$\n", "R = \n", "\\begin{pmatrix}\n", " u_{x}^{2}(1-c)+c & u_{x}u_{y}(1-c)-u_{z}s & u_{x}u_{z}(1-c)+u_{y}s \\\\\n", " u_{x}u_{y}(1-c)+u_{z}s & u_{y}^{2}(1-c)+c & u_{y}u_{z}(1-c)-u_{x}s \\\\\n", " u_{x}u_{z}(1-c)-u_{y}s & u_{y}u_{z}(1-c)+u_{x}s & u_{z}^{2}(1-c)+c\n", "\\end{pmatrix}\n", "$$\n", "\n", "where $c = \\cos(\\theta)$ and $s = \\sin(\\theta)$\n", "\n", "See https://fr.wikipedia.org/wiki/Matrice_de_rotation" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "scrolled": false }, "outputs": [], "source": [ "# Define initial points\n", "A = np.array([[0., 0., 0.],\n", " [1., 0., 0.],\n", " [0., 0., 0.],\n", " [0., 1., 0.],\n", " [0., 0., 0.],\n", " [0., 0., 1.]])\n", "\n", "# Define the rotation angle\n", "theta = np.radians(10)\n", "\n", "u = np.array([1., 1., 0.])\n", "ux, uy, uz = u[0], u[1], u[2]\n", "\n", "c = np.cos(theta)\n", "s = np.sin(theta)\n", "\n", "# Define the rotation matrices\n", "R = np.array([[ux**2 * (1-c) + c, ux*uy * (1-c) - uz*s, ux*uz * (1-c) + uy*s],\n", " [ux*uy * (1-c) + uz*s, ux**2 * (1-c) + c, uy*uz * (1-c) - ux*s],\n", " [ux*uz * (1-c) - uy*s, uy*uz * (1-c) + ux*s, uz**2 * (1-c) + c]])\n", "\n", "# Rotate points\n", "Ar = np.dot(R, A.T).T\n", "\n", "# Plot\n", "fig = plt.figure()\n", "ax = axes3d.Axes3D(fig)\n", "plot3d(A, axis=ax, fmt=\"-ok\")\n", "plot3d(np.array([np.zeros(3), u]), axis=ax, fmt=\"--ok\")\n", "plot3d(Ar, axis=ax, fmt=\":or\")\n", "ax.text(1, 0, 0, \"x\", color=\"k\")\n", "ax.text(0, 1, 0, \"y\", color=\"k\")\n", "ax.text(0, 0, 1, \"z\", color=\"k\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Projections" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Project 2D points on a line" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Line equation" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "@interact(a=(-5., 5., 0.1), b=(-5., 5., 0.1), c=(-5., 5., 0.1))\n", "def plot(a, b, c):\n", " plt.axis('equal')\n", " plt.axis([-5, 5, -5, 5])\n", " plt.xticks(np.arange(-5,5,1))\n", " plt.yticks(np.arange(-5,5,1))\n", " plt.axhline(y=0, color='k')\n", " plt.axvline(x=0, color='k')\n", "\n", " x = np.array([-10., 10.])\n", " f = lambda x: a/(-b) * x + c/(-b)\n", " \n", " try:\n", " plt.plot(x, f(x))\n", " except ZeroDivisionError:\n", " print(\"b should not be equal to 0\")\n", "\n", " plt.grid()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "#### Distance from a point to a line" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Line defined by an equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of a line in the plane given by the equation $ax + by + c = 0$, where $a$, $b$ and $c$ are real constants with $a$ and $b$ not both zero, the distance from the line to a point $(x_0,y_0)$ is\n", "$$\\operatorname{distance}(ax+by+c=0, (x_0, y_0)) = \\frac{|ax_0+by_0+c|}{\\sqrt{a^2+b^2}}.$$\n", "\n", "The point on this line which is closest to $(x_0,y_0)$ has coordinates:\n", "$$x = \\frac{b(bx_0 - ay_0)-ac}{a^2 + b^2}$$\n", "and\n", "$$y = \\frac{a(-bx_0 + ay_0) - bc}{a^2+b^2}.$$\n", "\n", "For mor information, see https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Setup the plot\n", "\n", "def plot(a, b, c, p, p2):\n", " plt.axis('equal')\n", " plt.axis([-5, 5, -5, 5])\n", " plt.xticks(np.arange(-5,5,1))\n", " plt.yticks(np.arange(-5,5,1))\n", " plt.axhline(y=0, color='k')\n", " plt.axvline(x=0, color='k')\n", "\n", " x = np.array([-10., 10.])\n", " f = lambda x: a/(-b) * x + c/(-b)\n", " \n", " plt.plot(x, f(x))\n", " plt.scatter(*p)\n", "\n", " # Plot the projection point\n", "\n", " plt.scatter(*p2)\n", " plt.plot((p2[0], p[0]), (p2[1], p[1]))\n", " #plt.arrow(*p2, *p) # TODO: doesn't work...\n", "\n", " plt.grid()\n", "\n", "# Define the distance and projection functions\n", "\n", "def distance(a, b, c, p):\n", " d1 = (a*p[0] + b*p[1] + c)\n", " d2 = math.sqrt(math.pow(a, 2.) + math.pow(b, 2.))\n", " d = abs(d1)/d2\n", " return d\n", "\n", "def projection(a, b, c, p):\n", " p2 = ((b*(b*p[0] - a*p[1]) - a*c)/(math.pow(a,2.)+math.pow(b,2.)),\n", " (a*(-b*p[0] + a*p[1]) - b*c)/(math.pow(a,2.)+math.pow(b,2.)))\n", " return p2\n", "\n", "# Define the line and the point\n", "\n", "a = 2.\n", "b = 1.\n", "c = -2.\n", "\n", "p = (-4., 2.)\n", "\n", "# Compute the distance and the projection point on the line\n", "\n", "d = distance(a, b, c, p)\n", "p2 = projection(a, b, c, p)\n", "\n", "print(\"Distance:\", d)\n", "print(\"Projection point:\", p2)\n", "\n", "# Plot the line and the point\n", "\n", "plot(a, b, c, p, p2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Line defined by two points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the line passes through two points $P_1 = (x_1, y_1)$ and $P_2 = (x_2, y_2)$ then the distance of $(x_0, y_0)$ from the line is:\n", "\n", "$$\\operatorname{distance}(P_1, P_2, (x_0, y_0)) = \\frac{|(y_2-y_1)x_0-(x_2-x_1)y_0+x_2 y_1-y_2 x_1|}{\\sqrt{(y_2-y_1)^2+(x_2-x_1)^2}}.$$\n", "\n", "The denominator of this expression is the distance between $P_1$ and $P_2$.\n", "The numerator is twice the area of the triangle with its vertices at the three points, $(x_0, y_0)$, $P_1$ and $P_2$.\n", "\n", "The expression is equivalent to $h=\\frac{2A}{b}$, which can be obtained by rearranging the standard formula for the area of a triangle: $A = \\frac{1}{2} bh$, where $b$ is the length of a side, and $h$ is the perpendicular height from the opposite vertex." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# TODO..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Line defined by a point and an angle" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def angle_point_to_equation(angle_degree, p):\n", " angle_radian = math.radians(angle_degree)\n", " \n", " a = math.tan(angle_radian)\n", " b = -1\n", " c = -math.tan(angle_radian) * p[0] + p[1]\n", " \n", " return a, b, c\n", "\n", "angle_degree = 30\n", "p0 = (3, 2)\n", "\n", "a, b, c = angle_point_to_equation(angle_degree, p0)\n", "\n", "p = (-4., 2.)\n", "\n", "# Compute the distance and the projection point on the line\n", "\n", "d = distance(a, b, c, p)\n", "p2 = projection(a, b, c, p)\n", "\n", "print(\"Distance:\", d)\n", "print(\"Projection point:\", p2)\n", "\n", "# Plot the line and the point\n", "\n", "plot(a, b, c, p, p2)\n", "plt.scatter(*p0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Project 3D points on a plane without perspective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "TODO..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Project 3D points on a plane with perspective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following variables are defined to describe this transformation:\n", "\n", "* $\\mathbf{a}_{x,y,z}$: the 3D position of a point A that is to be projected.\n", "* $\\mathbf{c}_{x,y,z}$: the 3D position of a point C representing the camera.\n", "* $\\mathbf{\\theta}_{x,y,z}$: the orientation of the camera (represented by Tait–Bryan angles).\n", "* $\\mathbf{e}_{x,y,z}$: the viewer's position **relative to the display surface which goes through point C** representing the camera (i.e. $\\mathbf{e}$ doesn't change with the position of the camera C and typically $\\mathbf{e} = (0, 0, -1)^T$ for a view angle of $2*\\cos^{-1}(1/\\sqrt{2})$ radians when the projection is made on the x/y plane).\n", "\n", "Which results in:\n", "\n", "* $\\mathbf {b} _{x,y}$: the 2D projection of a $\\mathbf{a}$\n", "\n", "$$\n", "\\begin{bmatrix}\n", " \\mathbf {d} _{x} \\\\\n", " \\mathbf {d} _{y} \\\\\n", " \\mathbf {d} _{z} \\\\\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", " 1 & 0 & 0 \\\\\n", " 0 & \\cos(\\mathbf{\\theta}_{x}) & \\sin(\\mathbf{\\theta}_{x}) \\\\\n", " 0 & -\\sin(\\mathbf{\\theta}_{x}) & \\cos(\\mathbf{\\theta}_{x}) \\\\\n", "\\end{bmatrix}\n", "~\n", "\\begin{bmatrix}\n", " \\cos(\\mathbf{\\theta}_{y}) & 0 & -\\sin(\\mathbf{\\theta}_{y}) \\\\\n", " 0 & 1 & 0 \\\\\n", " \\sin(\\mathbf{\\theta}_{y}) & 0 & \\cos(\\mathbf{\\theta}_{y}) \\\\\n", "\\end{bmatrix}\n", "~\n", "\\begin{bmatrix}\n", " \\cos(\\mathbf{\\theta}_{z}) & \\sin(\\mathbf{\\theta}_{z}) & 0 \\\\\n", " -\\sin(\\mathbf{\\theta}_{z}) & \\cos(\\mathbf{\\theta}_{z}) & 0 \\\\\n", " 0 & 0 & 1 \\\\\n", "\\end{bmatrix}\n", "~\n", "\\left(\n", "\\begin{bmatrix}\n", " \\mathbf{a}_{x} \\\\\n", " \\mathbf{a}_{y} \\\\\n", " \\mathbf{a}_{z} \\\\\n", "\\end{bmatrix}\n", "-\n", "\\begin{bmatrix}\n", " \\mathbf{c}_{x} \\\\\n", " \\mathbf{c}_{y} \\\\\n", " \\mathbf{c}_{z} \\\\\n", "\\end{bmatrix}\n", "\\right)\n", "$$\n", "\n", "The transformed point $\\boldsymbol d$ of $\\boldsymbol d$ can then be projected onto the 2D plane using the formula (here, **x/y is used as the projection plane**; literature also may use x/z):\n", "\n", "$$\n", "\\begin{array}{lcl}\n", " \\mathbf{b}_x & = & \\frac{\\mathbf{e}_z}{\\mathbf{d}_z} \\mathbf{d}_x - \\mathbf{e}_x \\\\\n", " \\mathbf{b}_y & = & \\frac{\\mathbf{e}_z}{\\mathbf{d}_z} \\mathbf{d}_y - \\mathbf{e}_y \\\\\n", "\\end{array}\n", "$$\n", "\n", "Or, in matrix form using homogeneous coordinates, the system\n", "\n", "$$\n", "\\begin{bmatrix}\n", " \\mathbf{f}_x \\\\\n", " \\mathbf{f}_y \\\\\n", " \\mathbf{f}_z \\\\\n", " \\mathbf{f}_w\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", " 1 & 0 & -\\frac{\\mathbf{e}_x}{\\mathbf{e}_z} & 0 \\\\\n", " 0 & 1 & -\\frac{\\mathbf{e}_y}{\\mathbf{e}_z} & 0 \\\\\n", " 0 & 0 & 1 & 0 \\\\\n", " 0 & 0 & \\frac{1}{\\mathbf{e}_z} & 0\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", " \\mathbf{d}_x \\\\\n", " \\mathbf{d}_y \\\\\n", " \\mathbf{d}_z \\\\\n", " 1\n", "\\end{bmatrix}\n", "$$\n", "\n", "in conjunction with an argument using similar triangles, leads to division by the homogeneous coordinate, giving\n", "\n", "$$\n", "\\begin{array}{lcl}\n", " \\mathbf{b}_x & = & \\mathbf{f}_x / \\mathbf{f}_w \\\\\n", " \\mathbf{b}_y & = & \\mathbf{f}_y / \\mathbf{f}_w \\\\\n", "\\end{array}\n", "$$\n", "\n", "The distance of the viewer from the display surface, $\\mathbf{e}_{z}$, directly relates to the field of view, where $\\alpha = 2 \\cdot \\tan^{-1}(\\frac{1}{\\mathbf{e}_{z}})$ is the viewed angle.\n", "\n", "Note: here we assume that you map the points $(-1,-1)$ and $(1,1)$ to the corners of your viewing surface.\n", "\n", "See: https://en.wikipedia.org/wiki/3D_projection#Perspective_projection" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "scrolled": false }, "outputs": [], "source": [ "# Define initial points to project\n", "a = np.array([0., 1., 2.])\n", "\n", "# Define camera's position\n", "c = np.array([0., 0., 0.])\n", "\n", "# Define viewer's position\n", "e = np.array([0., 0., -1.])\n", "\n", "# Define the orientation of the camera\n", "theta = np.array([np.radians(0),\n", " np.radians(0),\n", " np.radians(0)])\n", "\n", "theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]\n", "\n", "# Define the rotation matrices\n", "Rx = np.array([[1., 0., 0.],\n", " [0., np.cos(theta_x), np.sin(theta_x)],\n", " [0., -np.sin(theta_x), np.cos(theta_x)]])\n", "\n", "Ry = np.array([[np.cos(theta_y), 0., -np.sin(theta_y)],\n", " [0., 1., 0. ],\n", " [np.sin(theta_y), 0., np.cos(theta_y)]])\n", "\n", "Rz = np.array([[np.cos(theta_z), np.sin(theta_z), 0.],\n", " [-np.sin(theta_z), np.cos(theta_z), 0.],\n", " [0., 0., 1.]])\n", "\n", "d = np.dot(Rx, Ry)\n", "d = np.dot(d, Rz)\n", "d = np.dot(d, a-c)\n", "\n", "## TODO: which version is correct ? The one above or the one below ?\n", "#d = a - c\n", "#d = np.dot(Rz, d)\n", "#d = np.dot(Ry, d)\n", "#d = np.dot(Rx, d)\n", "\n", "print(\"d:\", d)\n", "\n", "b = np.array([e[2]/d[2] * d[0] - e[0],\n", " e[2]/d[2] * d[1] - e[1]])\n", "\n", "print(\"b:\", b)\n", "\n", "# Alternative to compute b\n", "Rf = np.array([[1., 0., -e[0]/e[2], 0.],\n", " [0., 1., -e[1]/e[2], 0.],\n", " [0., 0., 1., 0.],\n", " [0., 0., 1./e[2], 0.]])\n", "\n", "f = np.dot(Rf, np.concatenate([d, np.ones(1)]))\n", "\n", "b = np.array([f[0]/f[3],\n", " f[1]/f[3]])\n", "\n", "print(\"b:\", b)\n", "\n", "plot2d(np.array([b, b]), \"ok\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Multiple points version" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "@interact(theta_x=(-90., 90., 1.), theta_y=(-90., 90., 1.), theta_z=(-90., 90., 1.))\n", "def projection(theta_x, theta_y, theta_z):\n", " \n", " # Define initial points to project\n", " A = np.array([[-1., 0., 1.],\n", " [ 1., 0., 1.],\n", " [-1., 0., 2.],\n", " [ 1., 0., 2.],\n", " [-1., 0., 5.],\n", " [ 1., 0., 5.],\n", " [-1., 0., 15.],\n", " [ 1., 0., 15.]])\n", "\n", " # Define camera's position\n", " c = np.array([0., -2., 0.])\n", " C = np.tile(c, (A.shape[0], 1))\n", "\n", " # Define viewer's position\n", " e = np.array([0., 0., -1.])\n", "\n", " # Define the orientation of the camera\n", " theta = np.radians(np.array([theta_x, theta_y, theta_z]))\n", "\n", " theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]\n", "\n", " # Define the rotation matrices\n", " Rx = np.array([[1., 0., 0.],\n", " [0., np.cos(theta_x), np.sin(theta_x)],\n", " [0., -np.sin(theta_x), np.cos(theta_x)]])\n", "\n", " Ry = np.array([[np.cos(theta_y), 0., -np.sin(theta_y)],\n", " [0., 1., 0. ],\n", " [np.sin(theta_y), 0., np.cos(theta_y)]])\n", "\n", " Rz = np.array([[np.cos(theta_z), np.sin(theta_z), 0.],\n", " [-np.sin(theta_z), np.cos(theta_z), 0.],\n", " [0., 0., 1.]])\n", "\n", " d = np.dot(Rx, Ry)\n", " d = np.dot(d, Rz)\n", " d = np.dot(d, (A-C).T)\n", "\n", " ## TODO: which version is correct ? The one above or the one below ?\n", " #d = a - c\n", " #d = np.dot(Rz, d)\n", " #d = np.dot(Ry, d)\n", " #d = np.dot(Rx, d)\n", "\n", " print(\"d:\", d)\n", "\n", " b = np.array([e[2]/d[2] * d[0] - e[0],\n", " e[2]/d[2] * d[1] - e[1]])\n", "\n", " print(\"b:\", b)\n", "\n", " # Alternative to compute b\n", " Rf = np.array([[1., 0., -e[0]/e[2], 0.],\n", " [0., 1., -e[1]/e[2], 0.],\n", " [0., 0., 1., 0.],\n", " [0., 0., 1./e[2], 0.]])\n", "\n", " # Add a line of ones\n", " d = np.vstack([d, np.ones(d.shape[1])])\n", "\n", " f = np.dot(Rf, d)\n", "\n", " b = np.array([f[0]/f[3],\n", " f[1]/f[3]])\n", "\n", " print(\"b:\", b)\n", "\n", " plot2d(b.T, \"ok\")\n", " plot2d(b.T, \"-k\")" ] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Tags", "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python3" }, "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.5.3" } }, "nbformat": 4, "nbformat_minor": 1 }