{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Frames, Points and Point Clouds\n",
    "\n",
    "`PointCloud` is a numba-jittable class, dedicated to large sets of points in Euclidean space. It is a subclass of the `Vector` class of the `Neumann` library, thus being an estension of NumPy's `ndarray` class equipped with a transformation mechanism and other goodies. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sigmaepsilon.mesh.space import PointCloud\n",
    "from sigmaepsilon.mesh.triang import triangulate\n",
    "\n",
    "coords, *_ = triangulate(size=(800, 600), shape=(3, 3))\n",
    "points = PointCloud(coords)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A remarkable feature is that sliced objects retain their indices in the original pointcloud:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 4, 5])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points[3:6].inds\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get the index of the closest point (according to the standard Euclidean metric)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.index_of_closest(points.center())\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or several indices of points being closest to an array of targets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2], dtype=int64)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.index_of_closest(points[:3])\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likewise, for the indices of the points being the furthest from one or more targets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8, 6, 6], dtype=int64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.index_of_furthest(points[:3])\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions like `closest` and `furthest` return either a `Point` or a `PointCloud`, depending on the number of targets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([400., 300.,   0.])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.closest(points.center())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sigmaepsilon.mesh.space.point.Point"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(points.closest(points.center()))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.closest(points.center()).id\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PointCloud([[  0.   0.   0.]\n",
       " [400.   0.   0.]\n",
       " [800.   0.   0.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.closest(points[:3])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PointCloud([[800. 600.   0.]\n",
       " [  0. 600.   0.]\n",
       " [  0. 600.   0.]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.furthest(points[:3])\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get the bounds calling the `bounds` method on an instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0., 800.],\n",
       "       [  0., 600.],\n",
       "       [  0.,   0.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.bounds()\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `PointCloud` class is an instance of the `Vector` class defined in `sigmaepsilon.math` and therefore is equipped with all the mechanisms that the library provides. For instance, to apply a 90 degree rotation about the Z axis and then to move the points along the X axis:       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-590.,   10.],\n",
       "       [   0.,  800.],\n",
       "       [   0.,    0.]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "(\n",
    "    points\n",
    "    .rotate(\"Space\", [0, 0, np.pi / 2], \"XYZ\")\n",
    "    .move(np.array([10.0, 0., 0.]))\n",
    "    .bounds()\n",
    ")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba JIT compilation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data and the indices are both accessible inside numba-jitted functions, even in `nopython` mode. See code block below shows the available attributes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0., 0., 1.]),\n",
       " 0.0,\n",
       " array([0., 0.]),\n",
       " Array([[0., 0., 1.],\n",
       "        [0., 0., 0.]]),\n",
       " array([1, 2]))"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from numba import jit\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def numba_nopython(arr):\n",
    "    return arr[0], arr[0, 0], arr.x, arr.data, arr.inds\n",
    "\n",
    "\n",
    "c = np.array([[0, 0, 0], [0, 0, 1.0], [0, 0, 0]])\n",
    "pc = PointCloud(c, inds=np.array([0, 1, 2]))\n",
    "numba_nopython(pc[1:])\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference frames"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sigmaepsilon.mesh` relies on the mechanism introduced in `sigmaepsilon.math` to account for different reference frames and it introduces a new `CartesianFrame` class to be used for geometrical applications. You can then use these frames when defining a pointcloud. If you define a pointcloud without explicity specifying a frame, it is assumed that the points are embedded in the ambient frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0., 800.],\n",
       "       [  0., 600.],\n",
       "       [  0.,   0.]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points = PointCloud(coords)\n",
    "points.bounds()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-6.0000000e+02,  4.8985872e-14],\n",
       "       [ 0.0000000e+00,  8.0000000e+02],\n",
       "       [ 0.0000000e+00,  0.0000000e+00]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sigmaepsilon.mesh.space import CartesianFrame\n",
    "\n",
    "ambient_frame = CartesianFrame(dim=3)\n",
    "definition_frame = ambient_frame.rotate(\n",
    "    \"Space\", [0, 0, np.pi / 2], \"XYZ\", inplace=False)\n",
    "\n",
    "points = PointCloud(coords, frame=definition_frame)\n",
    "points.bounds()\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When calling the `bounds` method on a pointcloud, we can specify a target reference frame. Of yourse, the bounds of the pointcloud in the definition frame is the same as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0., 800.],\n",
       "       [  0., 600.],\n",
       "       [  0.,   0.]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.bounds(definition_frame)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access the frame and the coordinates of a pointcloud through properties `frame` and `array`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6.123234e-17,  1.000000e+00,  0.000000e+00],\n",
       "       [-1.000000e+00,  6.123234e-17,  0.000000e+00],\n",
       "       [ 0.000000e+00,  0.000000e+00,  1.000000e+00]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.frame\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([[  0.,   0.,   0.],\n",
       "       [400.,   0.,   0.],\n",
       "       [800.,   0.,   0.],\n",
       "       [  0., 300.,   0.],\n",
       "       [400., 300.,   0.],\n",
       "       [800., 300.,   0.],\n",
       "       [  0., 600.,   0.],\n",
       "       [400., 600.,   0.],\n",
       "       [800., 600.,   0.]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.array\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to get the coorinates of a pointcloud in a target frame, use the `show` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([[ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00],\n",
       "       [ 2.4492936e-14, -4.0000000e+02,  0.0000000e+00],\n",
       "       [ 4.8985872e-14, -8.0000000e+02,  0.0000000e+00],\n",
       "       [ 3.0000000e+02,  1.8369702e-14,  0.0000000e+00],\n",
       "       [ 3.0000000e+02, -4.0000000e+02,  0.0000000e+00],\n",
       "       [ 3.0000000e+02, -8.0000000e+02,  0.0000000e+00],\n",
       "       [ 6.0000000e+02,  3.6739404e-14,  0.0000000e+00],\n",
       "       [ 6.0000000e+02, -4.0000000e+02,  0.0000000e+00],\n",
       "       [ 6.0000000e+02, -8.0000000e+02,  0.0000000e+00]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_frame = ambient_frame.rotate(\n",
    "    \"Space\", [0, 0, np.pi], \"XYZ\", inplace=False)\n",
    "points.show(target_frame)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important difference between the `CartesianFrame` class in `sigmaepsilon.mesh` and `Neumann` is that the former also supports the concept of an origo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([[2., 0., 0.],\n",
       "       [1., 1., 0.],\n",
       "       [1., 0., 1.]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ambient_frame = CartesianFrame(dim=3)\n",
    "\n",
    "origo = np.array([1.0, 0.0, 0.0])\n",
    "definition_frame = CartesianFrame(dim=3, origo=origo)\n",
    "\n",
    "coords = np.array([\n",
    "    [1.0, 0.0, 0.0],\n",
    "    [0.0, 1.0, 0.0],\n",
    "    [0.0, 0.0, 1.0]\n",
    "])\n",
    "points = PointCloud(coords, frame=definition_frame)\n",
    "points.show(ambient_frame)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 ('.sigeps': venv)",
   "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.8.10"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5facf25dadae24d0f6f3d9b821e9851478f51388ee31821a60476e833f1169c6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}