task02.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a64b3251",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d599a24b",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_views = 101\n",
    "n_features = 215\n",
    "data = np.loadtxt('data/data_matrix.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7485db17",
   "metadata": {},
   "source": [
    "### 1. Create the centroid subtracted matrix D from feature coordinates\n",
    "- Given m images and n features\n",
    "- $m$ = n_views = 101 and $n$ = n_features = 215\n",
    "- D is measurement matrix of size $2m \\times n$ which is $202\\times 215$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20c285ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "centroids = np.mean(data, axis=1).reshape((-1, 1, 2))\n",
    "x = np.swapaxes(data.reshape((n_views, 2, n_features)), 1, 2)\n",
    "x_hat = x - centroids\n",
    "D = np.swapaxes(x_hat, 1, 2).reshape(n_views * 2, n_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a0432ec",
   "metadata": {},
   "source": [
    "### 2. Compute SVD of D and enforce rank constraint \n",
    " - $D = UWV^{T}$\n",
    " - Create $U_{3}$ and $V_{3}$ by taking the first columns of $U$ and $V$\n",
    " - Create $W_{3}$ by taking the upper left $3\\times3$ block of $W$\n",
    " - Create motion ($M$) and shape ($S$) matrices:\n",
    " $$M = U_{3}W_{3}^{\\frac{1}{2}} \\qquad S = W_{3}^{\\frac{1}{2}}V_{3}^{T}$$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a656cac5",
   "metadata": {},
   "outputs": [],
   "source": [
    "U, W, V = la.svd(D)\n",
    "\n",
    "U = U[:, 0:3]\n",
    "V = V[0:3, :]\n",
    "W = np.diag(np.squeeze(W)[0:3])\n",
    "\n",
    "W_sqrt = np.sqrt(W)\n",
    "\n",
    "M = U @ W_sqrt\n",
    "\n",
    "S = W_sqrt @ V\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a631b04a",
   "metadata": {},
   "outputs": [],
   "source": [
    "M.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b7a5324",
   "metadata": {},
   "source": [
    "### 3. Eliminate affine ambiguity by enforcing orthonormality constraint \n",
    "If the motion matrix for an image is $M=\\begin{bmatrix} m_{1}^{T} \\\\ m_{2}^{T} \\end{bmatrix}$ then we solve for $L=CC^{T}$ which is symmetric and has 6 variables.\n",
    "\n",
    "$$\n",
    "m_{1}L\\; m_{2}^{T} = \n",
    "\\begin{bmatrix}\n",
    "m_{11} & m_{12} & m_{13} \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "l_{11} & l_{12} & l_{13} \\\\\n",
    "l_{12} & l_{22} & l_{23} \\\\\n",
    "l_{13} & l_{23} & l_{33} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "m_{11} \\\\\n",
    "m_{12} \\\\\n",
    "m_{13} \n",
    "\\end{bmatrix} = 1\n",
    "$$\n",
    "can be rewritten as :\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "m_{11}^{2} & 2m_{11}m_{12} & 2m_{11}m_{13} & m_{12}^{2} & 2m_{12}m_{13} & m_{13}^{2}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "l_{11} \\\\\n",
    "l_{12} \\\\\n",
    "l_{13} \\\\\n",
    "l_{22} \\\\\n",
    "l_{23} \\\\\n",
    "l_{33}\n",
    "\\end{bmatrix} = 1\n",
    "$$\n",
    "The same can be done for: \n",
    "$$\n",
    "m_{2}L\\; m_{2}^{T} = 1 \\\\\n",
    "m_{1}L\\; m_{2}^{T} = 0 \\\\\n",
    "$$\n",
    "\n",
    "- Finnaly $C$ can be recovered from $L$ by SVD or Cholesky\n",
    "decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bba8797",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.zeros((3 * n_views, 6))\n",
    "B = np.zeros((3 * n_views, 1))\n",
    "\n",
    "for i in range(n_views):\n",
    "    line_1 = i * 3\n",
    "    line_2 = i * 3 + 1\n",
    "    line_3 = i * 3 + 2\n",
    "\n",
    "    mi1 = M[i * 2]\n",
    "    mi2 = M[i * 2 + 1]\n",
    "\n",
    "    m11, m12, m13 = mi1\n",
    "    m21, m22, m23 = mi2\n",
    "\n",
    "    A[line_1] = np.array([m11 ** 2,\n",
    "                          2 * m11 * m12,\n",
    "                          2 * m11 * m13,\n",
    "                          m12 ** 2,\n",
    "                          2 * m12 * m13,\n",
    "                          m13 ** 2])\n",
    "    A[line_2] = np.array([m21 ** 2,\n",
    "                          2 * m21 * m22,\n",
    "                          2 * m21 * m23,\n",
    "                          m22 ** 2,\n",
    "                          2 * m22 * m23,\n",
    "                          m23 ** 2])\n",
    "    A[line_3] = np.array([m11 * m21,\n",
    "                          m12 * m21 + m11 * m22,\n",
    "                          m13 * m21 + m11 * m23,\n",
    "                          m12 * m22,\n",
    "                          m13 * m22 + m12 * m23,\n",
    "                          m13 * m23])\n",
    "\n",
    "    B[line_1] = 1\n",
    "    B[line_2] = 1\n",
    "    B[line_3] = 0\n",
    "\n",
    "X = la.lstsq(A, B, rcond=-1)[0]\n",
    "\n",
    "l11, l12, l13, l22, l23, l33 = np.squeeze(X)\n",
    "\n",
    "L = np.array([\n",
    "    [l11, l12, l13],\n",
    "    [l12, l22, l23],\n",
    "    [l13, l23, l33]\n",
    "])\n",
    "\n",
    "C = la.cholesky(L)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "080e7044",
   "metadata": {},
   "source": [
    "### 4. Update M and S\n",
    "\n",
    "- $M^{\\prime}=MC$ and $S^{\\prime}=C^{-1}S$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "235a6ca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "M_prime = M @ C\n",
    "S_prime = la.inv(C) @ S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22b56d81",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}