tutorials/tutorial12/tutorial-12-cdmd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "A6kN-H5uItGh"
   },
   "source": [
    "# Foreground Modeling using Dynamic Mode Decomposition\n",
    "\n",
    "This tutorial is inspired by the paper [Compressed dynamic mode decomposition for background modeling](https://arxiv.org/abs/1512.04205) by Erichson et al.\n",
    "\n",
    "Author: [Josh Myers-Dean](https://joshmyersdean.github.io/)\n",
    "\n",
    "Packages needed: **PyDMD**, NumPy, pandas, opencv-python, Matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pKUxF2JhaRkO"
   },
   "source": [
    "## Download data\n",
    "\n",
    "First we will download the data, [SegTrackV2](https://web.engr.oregonstate.edu/~lif/SegTrack2/dataset.html) from Oregon State. This is an older binary segmentation dataset that offers a good test bed for this method. This cell could take some seconds to execute since it downloads and unzips a compressed file of size 200MB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "x9RD4MBQhP65",
    "outputId": "7c8b8b96-9c9c-4817-9081-bb951a301d24"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "print(\"Downloading data...\")\n",
    "if not os.path.exists(\"SegTrackv2\"):\n",
    "    !wget https://web.engr.oregonstate.edu/~lif/SegTrack2/SegTrackv2.zip >& /dev/null\n",
    "    !unzip -qq SegTrackv2.zip\n",
    "    print(\"Done!\")\n",
    "else:\n",
    "    print(\"Data already present\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "id": "FIFSwmDlgoGm",
    "outputId": "27417941-7c5c-493f-94ed-814ab7c2da95"
   },
   "outputs": [],
   "source": [
    "from pydmd import CDMD, DMD\n",
    "from pydmd.plotter import plot_eigs\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import style\n",
    "\n",
    "plt.gray()\n",
    "style.use(\"tableau-colorblind10\")\n",
    "from matplotlib import animation\n",
    "\n",
    "from IPython.display import HTML\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import cv2\n",
    "from typing import Tuple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nZE1yHh6fuvz"
   },
   "source": [
    "We peak at the data to examine the amount of frames per object video. Note that the videos with small number frames will most likely have poor results. It is recommended to use the 'Frog' and 'Worm' videos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 238
    },
    "id": "9N8cG-BqC9vN",
    "outputId": "655b293f-cfed-4e79-f496-a192e080dc77"
   },
   "outputs": [],
   "source": [
    "videos = sorted(os.listdir(\"SegTrackv2/JPEGImages\"))[1:]\n",
    "vid_metadata = {}\n",
    "for i in videos:\n",
    "    vid_metadata[i] = len(os.listdir(f\"SegTrackv2/JPEGImages/{i}\"))\n",
    "df = pd.DataFrame(vid_metadata.items(), columns=[\"Name\", \"Number of Frames\"])\n",
    "valid = [\n",
    "    \"birds_of_paradise\",\n",
    "    \"birdfall\",\n",
    "    \"frog\",\n",
    "    \"monkey\",\n",
    "    \"parachute\",\n",
    "    \"soldier\",\n",
    "    \"worm\",\n",
    "]\n",
    "df = df[df.Name.isin(valid)]\n",
    "display(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jcgBV8TiqJj8"
   },
   "outputs": [],
   "source": [
    "OBJ = \"frog\"  # Change this to desired object\n",
    "assert OBJ in df[\"Name\"].to_list(), \"Object not found in dataset\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Zof3WxAnf1Wy"
   },
   "source": [
    "## Methods needed for the tutorial\n",
    "Below we define (and comment) some methods which we're going to use later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RvQu_G-MrrOg"
   },
   "outputs": [],
   "source": [
    "def get_video_dmd(\n",
    "    object: str = \"frog\", noise: bool = False, noise_amt: float = 0.01\n",
    ") -> Tuple[np.ndarray, Tuple[int, int]]:\n",
    "    \"\"\"\n",
    "    Retreives a video in matrix format (i.e., each column is a frame)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    object :\n",
    "      name of video\n",
    "    noise :\n",
    "      boolean to use noise or not\n",
    "    noise_amt :\n",
    "      standard deviation for noise ~ N(0, noise_amt)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Tuple:\n",
    "      Matrix representation of the video\n",
    "      Original shape of frame\n",
    "\n",
    "    \"\"\"\n",
    "    imgs = []\n",
    "    jpeg_dir = \"SegTrackv2/JPEGImages/\"\n",
    "    shape = None\n",
    "    for i in sorted(os.listdir(f\"{jpeg_dir}{object}\")):\n",
    "        tmp = cv2.imread(f\"{jpeg_dir}{object}/{i}\", cv2.IMREAD_GRAYSCALE)\n",
    "        shape = tmp.shape\n",
    "        tmp = tmp.reshape(-1).astype(np.float32) / 255\n",
    "        imgs.append(tmp)\n",
    "    vid = np.vstack(imgs).T\n",
    "    if noise:\n",
    "        vid += np.random.normal(0, noise_amt, vid.shape)\n",
    "        vid = vid.clip(0, 1)\n",
    "    return np.vstack(imgs).T, shape\n",
    "\n",
    "\n",
    "def get_video(object: str = \"frog\") -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Retreives a video in tensor format (i.e., frames x heigh x width)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    object :\n",
    "      name of video\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Numpy array:\n",
    "      Tensor representation of the video\n",
    "\n",
    "    \"\"\"\n",
    "    imgs = []\n",
    "    jpeg_dir = \"SegTrackv2/JPEGImages/\"\n",
    "    for i in sorted(os.listdir(f\"{jpeg_dir}{object}\")):\n",
    "        tmp = cv2.imread(f\"{jpeg_dir}{object}/{i}\", cv2.IMREAD_GRAYSCALE)\n",
    "        tmp = tmp.astype(np.float32) / 255\n",
    "        imgs.append(tmp)\n",
    "    return np.asarray(imgs)\n",
    "\n",
    "\n",
    "def calc_iou(pred: np.ndarray, truth: np.ndarray, thresh: float = 0.1) -> float:\n",
    "    \"\"\"\n",
    "    Helper method to calculate IoU for single frame\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pred :\n",
    "      background subtracted video\n",
    "    truth :\n",
    "      segmentation ground truth\n",
    "    thresh :\n",
    "      cut off for deciding if a pixel is foreground\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float:\n",
    "      IoU of  a frame\n",
    "\n",
    "    \"\"\"\n",
    "    pred[pred < thresh] = 0\n",
    "    pred[pred >= thresh] = 1\n",
    "    intersection = np.logical_and(pred, truth).sum()\n",
    "    union = np.logical_or(pred, truth).sum()\n",
    "    return intersection / union if union > 0 else 0\n",
    "\n",
    "\n",
    "def calc_miou(pred, truth, thresh=0.1) -> float:\n",
    "    \"\"\"\n",
    "    Calculate average IoU for a video\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pred :\n",
    "      background subtracted video\n",
    "    truth :\n",
    "      segmentation ground truth\n",
    "    thresh :\n",
    "      cut off for deciding if a pixel is foreground\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float:\n",
    "      mIoU of  a video\n",
    "\n",
    "    \"\"\"\n",
    "    assert pred.shape == truth.shape, \"Pred and truth must be same shape\"\n",
    "    pred = pred.copy()\n",
    "    miou = 0\n",
    "    for i in range(pred.shape[0]):\n",
    "        iou = calc_iou(pred[i, :, :], truth[i, :, :], thresh=thresh)\n",
    "        miou += iou\n",
    "    return miou / pred.shape[0]\n",
    "\n",
    "\n",
    "def f1_score(y_true: np.ndarray, y_pred: np.ndarray, beta: int = 1) -> float:\n",
    "    \"\"\"\n",
    "    Calculate F1 score.\n",
    "    The original implmentation is written by Michal Haltuf on Kaggle.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    y_true :\n",
    "      segmentation ground truth\n",
    "\n",
    "    y_pred :\n",
    "      background subtracted video\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float:\n",
    "      f1 score of a video\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    tp = (y_true * y_pred).sum()\n",
    "    tn = ((1 - y_true) * (1 - y_pred)).sum()\n",
    "    fp = ((1 - y_true) * y_pred).sum()\n",
    "    fn = (y_true * (1 - y_pred)).sum()\n",
    "\n",
    "    epsilon = 1e-7\n",
    "\n",
    "    precision = tp / (tp + fp + epsilon)\n",
    "    recall = tp / (tp + fn + epsilon)\n",
    "\n",
    "    f1 = (1 + beta**2) * (precision * recall)\n",
    "    f1 /= beta**2 * precision + recall + epsilon\n",
    "    return f1\n",
    "\n",
    "\n",
    "def calc_f1(pred: np.ndarray, truth: np.ndarray, thresh: float = 0.1) -> float:\n",
    "    \"\"\"\n",
    "    Calculate f1 score for a video\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pred :\n",
    "      background subtracted video\n",
    "    truth :\n",
    "      segmentation ground truth\n",
    "    thresh :\n",
    "      cut off for deciding if a pixel is foreground\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float:\n",
    "      f1 score of a video\n",
    "\n",
    "    \"\"\"\n",
    "    assert pred.shape == truth.shape, \"Pred and truth must be same shape\"\n",
    "    pred = pred.copy()\n",
    "    truth = truth.copy()\n",
    "    pred[pred < thresh] = 0\n",
    "    pred[pred >= thresh] = 1\n",
    "    pred = pred.astype(np.uint8)\n",
    "    truth = truth.astype(np.uint8)\n",
    "    return f1_score(truth.flatten(), pred.flatten())\n",
    "\n",
    "\n",
    "def get_gt_video(object: str = \"frog\") -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Retreives a video of GT in tensor format (i.e., frames x heigh x width)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    object :\n",
    "      name of video\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Numpy array:\n",
    "      Tensor representation of the GT video\n",
    "\n",
    "    \"\"\"\n",
    "    imgs = []\n",
    "    jpeg_dir = \"SegTrackv2/GroundTruth/\"\n",
    "    gt_dir = \"SegTrackv2/JPEGImages/\"\n",
    "    jpeg_list = sorted(os.listdir(f\"{jpeg_dir}{object}\"))\n",
    "    gt_list = sorted(os.listdir(f\"{gt_dir}{object}\"))\n",
    "    valid = list(set(jpeg_list).intersection(set(gt_list)))\n",
    "    for i in sorted(valid):\n",
    "        tmp = f\"{gt_dir}/{i}\"\n",
    "        tmp = cv2.imread(f\"{jpeg_dir}{object}/{i}\", cv2.IMREAD_GRAYSCALE)\n",
    "        tmp = tmp.astype(np.float32) / 255\n",
    "        imgs.append(tmp)\n",
    "    return np.asarray(imgs)\n",
    "\n",
    "\n",
    "def play_video(object: str = \"frog\", interval: int = 10) -> None:\n",
    "    \"\"\"\n",
    "    Helper function to play original video\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    object :\n",
    "      name of video\n",
    "    interval :\n",
    "      delay in ms between frames\n",
    "\n",
    "    \"\"\"\n",
    "    video = get_video(object=object)\n",
    "\n",
    "    fig = plt.figure()\n",
    "    im = plt.imshow(video[0, :, :])\n",
    "\n",
    "    plt.close()  # this is required to not display the generated image\n",
    "\n",
    "    def init():\n",
    "        im.set_data(video[0, :, :])\n",
    "\n",
    "    def animate(i):\n",
    "        im.set_data(video[i, :, :])\n",
    "        return im\n",
    "\n",
    "    anim = animation.FuncAnimation(\n",
    "        fig, animate, init_func=init, frames=video.shape[0], interval=interval\n",
    "    )\n",
    "    return HTML(anim.to_html5_video())\n",
    "\n",
    "\n",
    "def play_gt_video(object: str = \"frog\", interval: int = 10):\n",
    "    \"\"\"\n",
    "    Helper function to play GT video\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    object :\n",
    "      name of video\n",
    "    interval :\n",
    "      delay in ms between frames\n",
    "\n",
    "    \"\"\"\n",
    "    video = get_gt_video(object=object)\n",
    "\n",
    "    fig = plt.figure()\n",
    "    im = plt.imshow(video[0, :, :])\n",
    "\n",
    "    plt.close()  # this is required to not display the generated image\n",
    "\n",
    "    def init():\n",
    "        im.set_data(video[0, :, :])\n",
    "\n",
    "    def animate(i):\n",
    "        im.set_data(video[i, :, :])\n",
    "        return im\n",
    "\n",
    "    anim = animation.FuncAnimation(\n",
    "        fig, animate, init_func=init, frames=video.shape[0], interval=interval\n",
    "    )\n",
    "    return HTML(anim.to_html5_video())\n",
    "\n",
    "\n",
    "def get_video_removed(video: np.ndarray, bg: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Helper function to subtract background from video\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    video :\n",
    "      original video\n",
    "    bg :\n",
    "      background model\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    np.ndarray :\n",
    "      predicted foreground\n",
    "\n",
    "    \"\"\"\n",
    "    for i in range(video.shape[0]):\n",
    "        video[i, :, :] -= bg\n",
    "    return video\n",
    "\n",
    "\n",
    "def play_video_removed(\n",
    "    bg: np.ndarray,\n",
    "    object: str = \"frog\",\n",
    "    mask: bool = False,\n",
    "    interval: int = 10,\n",
    "    thresh: float = 0.1,\n",
    ") -> None:\n",
    "    \"\"\"\n",
    "    Helper function to play foreground video\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    bg :\n",
    "      background model\n",
    "    object :\n",
    "      name of video\n",
    "    mask :\n",
    "      Show binary version\n",
    "    interval :\n",
    "      delay in ms between frames\n",
    "    thresh :\n",
    "      value to decide if pixel is foreground\n",
    "\n",
    "    \"\"\"\n",
    "    video = get_video(object=object)\n",
    "\n",
    "    fig = plt.figure()\n",
    "    subbed = video[0, :, :] - bg\n",
    "    if mask:\n",
    "        subbed[subbed > thresh] = 1\n",
    "        subbed[subbed <= thresh] = 0\n",
    "    else:\n",
    "        subbed = subbed.clip(0, 1)\n",
    "    im = plt.imshow(subbed)\n",
    "\n",
    "    plt.close()  # this is required to not display the generated image\n",
    "\n",
    "    def init():\n",
    "        subbed = video[0, :, :] - bg\n",
    "        if mask:\n",
    "            subbed[subbed > thresh] = 1\n",
    "            subbed[subbed <= thresh] = 0\n",
    "        else:\n",
    "            subbed = subbed.clip(0, 1)\n",
    "        im.set_data(subbed)\n",
    "\n",
    "    def animate(i):\n",
    "        subbed = video[i, :, :] - bg\n",
    "        if mask:\n",
    "            subbed[subbed > thresh] = 1\n",
    "            subbed[subbed <= thresh] = 0\n",
    "        else:\n",
    "            subbed = subbed.clip(0, 1)\n",
    "        im.set_data(subbed)\n",
    "        return im\n",
    "\n",
    "    anim = animation.FuncAnimation(\n",
    "        fig, animate, init_func=init, frames=video.shape[0], interval=interval\n",
    "    )\n",
    "    return HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wpV2jksDf6ny"
   },
   "source": [
    "## Let's peak at the data! \n",
    "\n",
    "The parameter `interval` is the delay (in ms) between frames. For shorter videos, use a higher interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 310
    },
    "id": "sbsAybOjtLW1",
    "outputId": "a65fd4e1-8977-4a18-c198-d416fd0dbd89"
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "play the video we will model,\n",
    "change interval based on video size\n",
    "\"\"\"\n",
    "play_video(OBJ, interval=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 310
    },
    "id": "krCa_vhE-2Oz",
    "outputId": "a0481ec6-8236-4e30-e568-e9413ceb4e1c"
   },
   "outputs": [],
   "source": [
    "# show gt video\n",
    "play_gt_video(OBJ, interval=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a DMD instance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oTD2LLy9f_qq"
   },
   "source": [
    "Time to fit our video to (c)DMD, add noise, tinker with SVD rank, etc. We can also visualize our background model! Since we average the first `K` modes, we can visualize how our model changes wrt `K`. \n",
    "\n",
    "Recall that each column of our matrix will be a frame in the video and unless we have a long video, the system will be overdetermined. Note that in the paper the authors use an optimization scheme to decide on the number of modes; we are going to choose it emprically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "GaRi8Tz3-m_W",
    "outputId": "fe7dfd7f-1c2e-4ed8-821b-b3096e05d3ca"
   },
   "outputs": [],
   "source": [
    "use_noise = False\n",
    "noise = 0.01\n",
    "video, shape = get_video_dmd(OBJ, use_noise, noise)  # get video\n",
    "print(f\"Condition number of video matrix is {np.linalg.cond(video): .3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "VX1O0IZLedpk"
   },
   "outputs": [],
   "source": [
    "comp = True  # use compressed DMD\n",
    "svd_rank = 0  # rank=0 will automatically detect rank, try other values next!\n",
    "optim = True  # Use optimized DMD\n",
    "compression = [\"linear\", \"sparse\", \"uniform\", \"sample\"]\n",
    "cmat = compression[2]  # compression matrix\n",
    "\n",
    "if comp:\n",
    "    dmd = CDMD(svd_rank=svd_rank, opt=optim, compression_matrix=cmat).fit(video)\n",
    "else:\n",
    "    dmd = DMD(svd_rank=svd_rank, opt=optim).fit(video)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "Ye4Zjqmd-gGE",
    "outputId": "16f21a9e-3755-44af-c8d6-a6c8bb513135"
   },
   "outputs": [],
   "source": [
    "modes = dmd.reconstructed_data.T.reshape(video.shape[1], shape[0], shape[1])\n",
    "\n",
    "# Try changing the value of K, amount of modes we use\n",
    "K = min(100, modes.shape[0])\n",
    "\n",
    "modes = np.abs(modes)  # deal with complex values\n",
    "bg = np.zeros_like(modes[0, :, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(5, 5, figsize=(16, 16))\n",
    "axes = axes.flatten()\n",
    "\n",
    "idx = 0\n",
    "for k in range(K):\n",
    "    bg += modes[k, :, :]\n",
    "    if k % 4 == 0:\n",
    "        if idx >= len(axes):\n",
    "            continue\n",
    "        axes[idx].axis(\"off\")\n",
    "        axes[idx].imshow(bg / (k + 1))\n",
    "        axes[idx].set_title(f\"K = {k}\")\n",
    "        idx += 1\n",
    "plt.suptitle(\"Background Model using varying amount of Modes\", y=0.92)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bg /= K\n",
    "plt.imshow(bg)\n",
    "plt.title(\"Background Model\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And an example frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 248
    },
    "id": "Qy6JF7s0xS3V",
    "outputId": "e0b82aa3-b119-4ebc-c3d8-d6a333686e1a"
   },
   "outputs": [],
   "source": [
    "tmp = get_video(OBJ)\n",
    "img = tmp[0, :, :]\n",
    "fg = (img - bg).clip(0, 1)\n",
    "plt.title(\"Foreground for Frame 0\")\n",
    "plt.imshow(fg)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lr--8mHUgOgp"
   },
   "source": [
    "## Results\n",
    "Time to get some quanitative (and qualitative) results! We will examine mean intersection over union (`mIoU`) and `F1`\n",
    "score, respectively. These metrics are a function of the threshold we choose for deciding foreground vs background. \n",
    "\n",
    "We can also visualize the eigenvalues, modes, and dynamics of our video computed from DMD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "j-V0DMjT57G0",
    "outputId": "4d2a8981-75c4-490d-e99a-d6e6c20c1de3"
   },
   "outputs": [],
   "source": [
    "video_removed = get_video_removed(tmp, bg).clip(0, 1)\n",
    "gt = get_gt_video(OBJ)\n",
    "\n",
    "video_removed.shape, gt.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute `mIoU` and `F1` as a function of threshold value, takes a bit to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds = np.linspace(0, 1, 10)\n",
    "\n",
    "mious = []\n",
    "f1s = []\n",
    "for thresh in thresholds:\n",
    "    mious.append(calc_miou(video_removed, gt, thresh=thresh))\n",
    "    f1s.append(calc_f1(video_removed, gt, thresh=thresh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "NdqBdO-b6My1",
    "outputId": "0448d8d2-a3ff-46d0-e140-7f4c845c2a05"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 5))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(mious, \"bo-\")\n",
    "plt.title(\"mIoU vs Threshold\")\n",
    "plt.ylabel(\"mIoU\")\n",
    "plt.xlabel(\"Threshold\")\n",
    "plt.xticks(range(len(thresholds)), np.round(thresholds, 2))\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(f1s, \"ro-\")\n",
    "plt.title(\"F1 Score vs Threshold\")\n",
    "plt.ylabel(\"F1\")\n",
    "plt.xlabel(\"Threshold\")\n",
    "plt.xticks(range(len(thresholds)), np.round(thresholds, 2))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot the video with the background removed. You can try playing with the threshold, keep in mind that it shoud remain inside $[0,1]$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 310
    },
    "id": "GcIHsnJFwZvh",
    "outputId": "94be630d-473b-45a7-900c-d000b904ceba"
   },
   "outputs": [],
   "source": [
    "# show binary output or not, if False thresh doesn't matter\n",
    "use_mask = False\n",
    "play_video_removed(bg, OBJ, mask=use_mask, thresh=0.03)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print the distances from the unit circle of the first 6 eigenvalues in `dmd.eigs`, and plot all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 611
    },
    "id": "IXNLX9gN0AY4",
    "outputId": "1526d755-7e21-4f16-f2f3-49ee0e3b21ea"
   },
   "outputs": [],
   "source": [
    "for idx, eig in enumerate(dmd.eigs[:6]):\n",
    "    print(\n",
    "        f\"Eigenvalue {eig}: distance from unit circle {np.abs(np.abs(eig)-1): .5f}\"\n",
    "    )\n",
    "\n",
    "plot_eigs(dmd, show_axes=True, show_unit_circle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "I1tVYZMxm8AZ"
   },
   "source": [
    "We also plot the first 6 modes and dynamics. The modes are hard to disentangle when SVD rank is larger than 3 but we can see the slow varying dynamic, which is our background mode!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 499
    },
    "id": "4JOSgXIx4A9w",
    "outputId": "19a563f7-f75d-4121-ed28-0765e8b708a8"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 8))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "for idx, mode in enumerate(dmd.modes.T[:6]):\n",
    "    plt.plot(mode.real, alpha=0.5, label=f\"Mode {idx}\")\n",
    "plt.title(\"Modes\")\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "for idx, dynamic in enumerate(dmd.dynamics[:6]):\n",
    "    plt.plot(dynamic.real, label=f\"Mode {idx}\")\n",
    "plt.title(\"Dynamics\")\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "tutorial-12-cdmd",
   "provenance": []
  },
  "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}