{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "595bfddc-6636-4f76-8487-4aa2953366a7",
   "metadata": {},
   "source": [
    "This notebook analyzes a video capture from the laser camera where a far mirror was set up 2.35m across the room.  The laser diode and camera were both placed on the vibration isolated mount of the STM.  In front of the laser diode, a ø2mm apterture was placed to somewhat contain the beam spread.\n",
    "\n",
    "From the video, we recover the center position of the laser in each frame.  Then we can analyze the movement of the laser beam over time and extract frequency content.  This should give us an indication of the transmittance of the vibration isolation over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6a7b4902-23e3-4e4d-a525-1a2c4112be54",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import ndimage\n",
    "import time\n",
    "\n",
    "capture = cv2.VideoCapture(\"./laser-img/2025-09-07-011400.webm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96660ca6-8d5f-493d-8154-791a19c0d635",
   "metadata": {},
   "source": [
    "Fetch video capture metadata for processing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7e80a62f-aca8-4d2c-95df-a3e29c5f3e24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Capture Info: 3000 frames with dimension 1920×1080 at 30 Hz, totalling 100.00 seconds\n"
     ]
    }
   ],
   "source": [
    "frame_count = int(capture.get(cv2.CAP_PROP_FRAME_COUNT))\n",
    "frame_rate = int(capture.get(cv2.CAP_PROP_FPS))\n",
    "frame_time = 1 / frame_rate\n",
    "frame_width = int(capture.get(cv2.CAP_PROP_FRAME_WIDTH))\n",
    "frame_height = int(capture.get(cv2.CAP_PROP_FRAME_HEIGHT))\n",
    "\n",
    "# Overwrite frame_count for limiting analysis window\n",
    "frame_count = 3000\n",
    "\n",
    "print(f\"Capture Info: {frame_count} frames with dimension {frame_width}×{frame_height} at {frame_rate} Hz, totalling {frame_count * frame_time:.2f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d648a6-1fde-4ea5-bdc9-c93bcbaf6147",
   "metadata": {},
   "source": [
    "## Laser Beam Position Recovery\n",
    "The laser beam position is recovered using the `ndimage.center_of_mass()` function.  This function hinges on the background of the image being completely dark — otherwise there will be a positional bias towards the center of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "462111d4-7609-43f0-8e34-1c7dcfcc63c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_laser_beam(image):\n",
    "    image_monochrome = np.mean(image, axis=2)\n",
    "    intensity = np.mean(image)\n",
    "    y, x = ndimage.center_of_mass(image_monochrome)\n",
    "    return np.array((x, y, intensity))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502c7df6-13fc-4a0c-b720-cec6326701c2",
   "metadata": {},
   "source": [
    "## Video Capture Processing\n",
    "Process the capture frame by frame to recover the laser beam position in each frame.  For verification of the recovery function, a few shots from the capture are visualized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "30816936-7cbc-4d29-9afd-b4082a687e66",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Snapshot capture of frame      0 at   0.00s: Intensity 9.282\n",
      "Snapshot capture of frame    333 at  11.10s: Intensity 10.440\n",
      "Snapshot capture of frame    666 at  22.20s: Intensity 12.150\n",
      "Snapshot capture of frame    999 at  33.30s: Intensity 12.801\n",
      "Snapshot capture of frame   1332 at  44.40s: Intensity 12.652\n",
      "Snapshot capture of frame   1665 at  55.50s: Intensity 11.099\n",
      "Snapshot capture of frame   1998 at  66.60s: Intensity 11.897\n",
      "Snapshot capture of frame   2331 at  77.70s: Intensity 12.612\n",
      "Snapshot capture of frame   2664 at  88.80s: Intensity 12.257\n",
      "Snapshot capture of frame   2997 at  99.90s: Intensity 12.728\n",
      "Video processing took 257.32 seconds (0.086 seconds per frame).\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 65512 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 122408 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 119108 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 142996 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 138844 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 85344 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 143032 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 129736 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 141420 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 132348 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "processing_start = time.monotonic()\n",
    "\n",
    "laser_position = np.empty((frame_count, 3), np.dtype('float32'))\n",
    "\n",
    "current_frame = 0\n",
    "ret = True\n",
    "while current_frame < frame_count and ret:\n",
    "    ret, frame = capture.read()\n",
    "    if ret:\n",
    "        current_laser = find_laser_beam(frame)\n",
    "        current_frame_time = current_frame * frame_time\n",
    "\n",
    "        laser_position[current_frame] = current_laser\n",
    "\n",
    "        if current_frame % (frame_count // 9) == 0:\n",
    "            print(f\"Snapshot capture of frame {current_frame:6} at {current_frame_time:6.2f}s: Intensity {current_laser[2]:.3f}\")\n",
    "            plt.figure()\n",
    "            plt.imshow(frame)\n",
    "            plt.axvline(x=current_laser[0], color=\"red\")\n",
    "            plt.axhline(y=current_laser[1], color=\"red\")\n",
    "            plt.figtext(0.5, 0.05, f\"Frame {current_frame} at {current_frame_time:.2f} seconds; Intensity {current_laser[2]:.3f}\", ha=\"center\")\n",
    "    \n",
    "    current_frame += 1\n",
    "\n",
    "processing_time = time.monotonic() - processing_start\n",
    "print(f\"Video processing took {processing_time:.2f} seconds ({processing_time / frame_count:.3f} seconds per frame).\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1c58de-57b6-47d7-8e6d-d1f9046ad1ca",
   "metadata": {},
   "source": [
    "Derive data for the plots below from the raw analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "26fdc3f5-72a9-48e6-824f-790128a80772",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_axis = np.linspace(0, frame_count * frame_time, frame_count)\n",
    "\n",
    "laser_displacement_zeroed = laser_position[:, 0:2] - np.nanmean(laser_position[:, 0:2], axis=0)\n",
    "\n",
    "velocity_xy = np.diff(laser_position[:, 0:2])\n",
    "velocity_abs = np.nan_to_num(np.linalg.norm(velocity_xy, axis=1))\n",
    "\n",
    "max_intensity = np.max(laser_position[:, 2])\n",
    "laser_intensity = laser_position[:, 2] / max_intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72bba9f4-6158-415a-b853-44fe70cdec94",
   "metadata": {},
   "source": [
    "### Time-dependent Laser Position and Intensity plots\n",
    "The following figure shows the laser displacement, velocity, and intensity over the whole capture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1593299b-d463-4a64-80ae-4bed7c035dcc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 490528 characters]",
      "text/plain": [
       "<Figure size 1800x1600 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(4, 1, figsize=(18, 16))\n",
    "ax_displacement = ax[0:3]\n",
    "ax_intensity = ax[3]\n",
    "\n",
    "ax_displacement[0].plot(time_axis, laser_displacement_zeroed[:, 0] / frame_width, label=\"Y position\")\n",
    "ax_displacement[1].plot(time_axis, laser_displacement_zeroed[:, 1] / frame_width, label=\"X position\")\n",
    "ax_displacement[2].plot(time_axis, velocity_abs, label=\"Velocity\")\n",
    "for ax_d in ax_displacement:\n",
    "    ax_d.set_xlabel(\"time / seconds\")\n",
    "    ax_d.set_ylabel(\"displacment / sensor width\")\n",
    "    ax_d.legend()\n",
    "\n",
    "ax_intensity.plot(time_axis, laser_intensity, label=\"Intensity\")\n",
    "ax_intensity.set_xlabel(\"time / seconds\")\n",
    "ax_intensity.set_ylabel(\"relative intensity\")\n",
    "ax_intensity.set_ylim(0, 1)\n",
    "ax_intensity.legend()\n",
    "\n",
    "fig.text(0.5, 0.05, \"Laser X/Y, Velocity, and Intensity over time\", ha=\"center\")\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bc7871c6-caf4-4e16-9218-3c928c8290d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "\n",
    "frequency_x, power_spectral_density_x = signal.periodogram(np.nan_to_num(laser_position[:, 1]), frame_rate)\n",
    "frequency_y, power_spectral_density_y = signal.periodogram(np.nan_to_num(laser_position[:, 0]), frame_rate)\n",
    "\n",
    "amplitude_spectral_density_x = np.sqrt(power_spectral_density_x)\n",
    "amplitude_spectral_density_y = np.sqrt(power_spectral_density_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b59c5f-da6b-40b9-9509-b92521bc4a7a",
   "metadata": {},
   "source": [
    "## Spectral Analysis\n",
    "The following figure looks at the spectral density in the displacement along X and Y axis.  This gives us a good insight to see which vibration frequencies were transmitted through the isolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "ce7eeef9-d657-40ce-b199-2b4f09dc54be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 223832 characters]",
      "text/plain": [
       "<Figure size 1600x800 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2, 1, figsize=(16, 8))\n",
    "\n",
    "ax[1].axvline(x=2.822, color=\"red\")\n",
    "ax[1].annotate(\" theoretical spring resonance\", xy=(2.822, 1080 / 2 - 150))\n",
    "\n",
    "ax[0].loglog(frequency_x, amplitude_spectral_density_x, label=\"X axis\")\n",
    "ax[1].loglog(frequency_y, amplitude_spectral_density_y, label=\"Y axis\")\n",
    "\n",
    "for a in ax:\n",
    "    a.set_xlabel(\"frequency / Hz\")\n",
    "    a.set_ylabel(\"amplitude spectral density\")\n",
    "    a.set_xlim(0.5, 15)\n",
    "    a.set_ylim(1, 1080 / 2)\n",
    "    a.legend()\n",
    "    a.grid(axis=\"x\", which=\"both\")\n",
    "    a.grid(axis=\"y\", which=\"major\")\n",
    "\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bf08176-0009-4e40-bae2-4fafa5aa04eb",
   "metadata": {},
   "source": [
    "Of note is the highlighted _theoretical spring resonance_ which we had previously calculated using the known spring constant and mass of the system.  Seeing a very visible peak at this frequency along the Y axis (the axis of the springs) is very nice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "9ae71da9-d3f9-4e55-b382-98b850bce3a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.save(\"./laser-results/2025-09-07-011400-analyzed.npy\", laser_position)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "074aa31c-1966-401c-bd20-538efe31c716",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stm-env",
   "language": "python",
   "name": "stm-env"
  },
  "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.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
