Direct Linear Transform (DLT)

DLT.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c6db7401-d049-4bbe-a656-d442de2e2bde",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "3ab2d64e-fbf5-49d1-bf12-d0ad7934bc92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dlt_rms_error(image_coords, spatial_coords, projector_est):\n",
    "    H = np.array([[0, -1], [1, 0]])\n",
    "    err = image_coords @ H @ projector_est @ spatial_coords.T\n",
    "    rmse = np.mean(err ** 2) ** (1/2)\n",
    "    return rmse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "643b8b3c-18aa-49a1-a9c3-7e47f6eac515",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dlt(image_coords, spatial_coords, threshold=0.1):\n",
    "    \"\"\"\n",
    "    Calculate the Direct Linear Transform (DLT).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    image_coords: array_like, Nx2\n",
    "        The 2D image coordinates.\n",
    "    spatial_coords: array_like, Nx3\n",
    "        The 3D spatial coordinates.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    projector: np.ndarray, 2x3\n",
    "        See notes.\n",
    "    \n",
    "    Notes\n",
    "    -----\n",
    "    Estimate the 2x3 projector P such that\n",
    "    \\\\alpha_{k} x_{k} = \\\\mathbf{P} y_{k}\n",
    "    where x are the 2D image coordinates, y are the 3D spatial coordinates, and alpha are unknown pesky scaling factors.\n",
    "    \n",
    "    https://en.wikipedia.org/wiki/Direct_linear_transformation#Example\n",
    "    \"\"\"\n",
    "    B = np.array([image_coords[:, 1] * spatial_coords[:, 0],\n",
    "                  -image_coords[:, 0] * spatial_coords[:, 0],\n",
    "                  image_coords[:, 1] * spatial_coords[:, 1],\n",
    "                  -image_coords[:, 0] * spatial_coords[:, 1],\n",
    "                  image_coords[:, 1] * spatial_coords[:, 2],\n",
    "                  -image_coords[:, 0] * spatial_coords[:, 2]]).T\n",
    "    u, s, vh = np.linalg.svd(B)\n",
    "    if np.abs(s[-1]) > threshold:\n",
    "        message = f'The projection is very noisy! The singular value should be close to 0, but got ~{s[-1]:.2g} instead. Check your data!'\n",
    "        warnings.warn(message, stacklevel=2)\n",
    "    a = vh[-1]\n",
    "    P = np.array([a[::2], a[1::2]])\n",
    "    rmse = dlt_rms_error(x, y, P)\n",
    "    if rmse > threshold:\n",
    "        message = f'The projection is very noisy! The Root-Mean-Square Error of a skew-symmetric projector operating on the data should be close to 0, but got ~{rmse:.2g} instead. Check your data!'\n",
    "        warnings.warn(message, stacklevel=2)\n",
    "    return P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "30b9da1e-80eb-40d5-8d6f-1e286377561d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.27564811, -0.52803998,  0.28600877],\n",
       "       [ 0.28151325, -0.58984647,  0.36908308]])"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y = np.random.random([11, 2]), np.random.random([11, 3]) # project 11 random points in the unit cube to random points on an imaging plane\n",
    "P = dlt(x, y)\n",
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92f741ce-1994-47af-8755-211cd1b89b40",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}