jupyter/registration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Point Cloud Registration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "In order to lower the uncertainty of calculations done with `py4dgeo`, point cloud registration should be performed to tightly align the point clouds. `py4dgeo` supports this in two ways:\n",
    "\n",
    "* It allows you to apply arbitrary affine transformations to epochs and stores the transformation parameters as part of the epoch. You can use the tooling of your choice to calculate these.\n",
    "* It provides a number of registration algorithms that allow you to calculate the transformation directly in `py4dgeo`.\n",
    "\n",
    "This notebook show cases both ways of usage.\n",
    "\n",
    "Note: Be aware that an initial transformation is required that roughly aligns the two point clouds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import py4dgeo\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "We load two epochs and register them to each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "reference_epoch = py4dgeo.read_from_xyz(\"plane_horizontal_t1.xyz\")\n",
    "epoch = py4dgeo.read_from_xyz(\"plane_horizontal_t2.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "Registration algorithms use a simple function call interface. It returns an affine transformation as a `3x4` matrix which contains the Rotation matrix $\\mathbf{R}$ and the translation vector $\\mathbf{t}$:\n",
    "\n",
    "$$\n",
    "\\mathbf{T} = \\left(\\begin{array}{ccc}r_{11}&r_{12}&r_{13}&t_1\\\\r_{21}&r_{22}&r_{23}&t_2\\\\r_{31}&r_{32}&r_{33}&t_3\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "When calculating and applying such transformations in `py4dgeo`, it is always possible to specify the *reduction point* $\\mathbf{x_0}$. This allows shifting coordinates towards the origin before applying rotations - leading to a numerically robust operation. The algorithm for a transformation is:\n",
    "\n",
    "$$\n",
    "\\mathbf{Tx} = \\left(\\mathbf{R}(\\mathbf{x}-\\mathbf{x_0})\\right) + \\mathbf{t} + \\mathbf{x_0}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## Registration API"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "We demonstrate the use of registration algorithms within `py4dgeo` using our implementation of the standard Iterative Closest Point (ICP) algorithm. Registration algorithms are provided as simple Python functions taking two epochs and the reduction point as arguments. Some algorithms might require additional, algorithm-specific parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "trafo = py4dgeo.iterative_closest_point(\n",
    "    reference_epoch, epoch, reduction_point=np.array([0, 0, 0])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "The returned `trafo` object stores both the computed affine transformation, as well as the reduction point. It is important to store these togethers, as only the correct combination of these two quantities will yield the correct transformation result. The transformation `trafo` can now be applied to `epoch` by passing it to `epoch.transform`. Note that, by design, this operation is under explicit user control and not done implicitly while calculating the transformation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch.transform(trafo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "The `Epoch` class records applied transformations and makes them available through the `transformation` property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch.transformation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "Finally, we can export the transformed epoch to inspect the registration quality e.g. in CloudCompare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch.save(\"plane_horizontal_t2_transformed.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "To learn more, about parameters of our ICP implementations, you can have a look at the documentation of `py4dgeo.iterative_closest_point` or `py4dgeo.point_to_plane_icp`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "?py4dgeo.point_to_plane_icp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "## Using existing transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "If you have already computed a transformation using the tool of your choice, you can easily import it into `py4dgeo` by loading it as `numpy` arrays and pass it directly to `Epoch.transform`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch.transform(\n",
    "    rotation=np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),\n",
    "    translation=np.array([0, 0, 0]),\n",
    "    reduction_point=np.array([0, 0, 0]),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "It is both possible to define `rotation` and `translation` individually or to provide a `3x4` matrix to the parameter `affine_transformation`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "## ICP algorithm with stable areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "`py4dgeo` does not only implement standard ICP algorithms, but also an ICP algorithm that is better suited for geomorphological applications. It identifies stable areas in the point cloud and restricts the application of the ICP algorithm to these areas. It is based on the article \"Supervoxel-based targetless registration and identification of stable areas for deformed point clouds\" by Yuhui Yang and Volker Schwieger."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "The algorithm requires normals. You can either have `py4dgeo` calculate these normals or load precomputed ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "reference_epoch = py4dgeo.read_from_las(\n",
    "    \"plane_horizontal_t1_w_normals.laz\",\n",
    "    normal_columns=[\"NormalX\", \"NormalY\", \"NormalZ\"],\n",
    ")\n",
    "epoch = py4dgeo.read_from_las(\n",
    "    \"plane_horizontal_t2_w_normals.laz\",\n",
    "    normal_columns=[\"NormalX\", \"NormalY\", \"NormalZ\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25",
   "metadata": {},
   "source": [
    "As seen with [the standard ICP algorithm](#Registration-API), the API for the registration algorithm is a simple function call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "trafo = py4dgeo.icp_with_stable_areas(\n",
    "    reference_epoch,\n",
    "    epoch,\n",
    "    initial_distance_threshold=10,\n",
    "    level_of_detection=10,\n",
    "    reference_supervoxel_resolution=0.2,\n",
    "    supervoxel_resolution=0.2,\n",
    "    min_svp_num=1,\n",
    "    reduction_point=np.array([0, 0, 0]),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27",
   "metadata": {},
   "source": [
    "Apart from the usual arguments `reference_epoch`, `epoch` and `reduction_point`, the algorithm allows the following customizations. To achieve best performance, you should familiarize yourself with these parameters and tune them for your application:\n",
    "\n",
    "* `initial_distance_threshold`: The upper boundary of the distance threshold in the iteration. It can be (1) an empirical value manually set by the user according to the approximate accuracy of coarse registration, or (2) calculated by the mean and standard of the nearest neighbor distances of all points.\n",
    "* `level_of_detection`: The lower boundary (minimum) of the distance threshold in the iteration. It can be  (1) an empirical value manually set by the user according to the approximate uncertainty of laser scanning measurements in different scanning configurations and scenarios (e.g., 1 cm for TLS point clouds in short distance and 4 cm in long distance, 8 cm for ALS point clouds, etc.), or (2) calculated by estimating the standard deviation from local modeling (e.g., using the level of detection in M3C2 or M3C2-EP calculations).\n",
    "* `reference_supervoxel_resolution`: The approximate size of generated supervoxels for the reference epoch. It can be (1) an empirical value manually set by the user according to different surface geometries and scanning distance (e.g., 2-10 cm for indoor scenes, 1-3 m for landslide surface), or (2) calculated by 10-20 times the average point spacing (original resolution of point clouds). In both cases, the number of points in each supervoxel should be at least 10 (i.e., min_svp_num = 10).\n",
    "* `supervoxel_resolution`:  The same as `reference_supervoxel_resolution` but for a different epoch.\n",
    "* `min_svp_num`: Minimum number of points for supervoxels to be taken into account in further calculations."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}