Night time photography calculations (sky coverage, star trail length)

2023-08-14-photography-calculations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Camera calculations, Pavel Peřina, 2023-08-14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import some libraries\n",
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are constants for Fujifilm X100V APS-C camera\n",
    "sensor_width  = 23.5e-3\n",
    "sensor_height = 15.6e-3\n",
    "focal_length  = 23.0e-3\n",
    "image_width   = 6240 \n",
    "image_height  = 4160"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Camera resolution:        26.0 MPix\n",
      "Pixel size:               3.77 µm\n",
      "Field of view horizontal: 54.1 °\n",
      "Field of view vertical:   37.5 °\n"
     ]
    }
   ],
   "source": [
    "# Pixel and sensor size\n",
    "pixel_size = sensor_width / image_width\n",
    "print(\"Camera resolution:        {:.3} MPix\".format(float(image_width * image_height) / 1.0e6))\n",
    "print(\"Pixel size:               {:.3} µm\".format(pixel_size * 1e6))\n",
    "\n",
    "# Field of view\n",
    "h_fov_deg = math.atan2(sensor_width  / 2, focal_length) * 2.0 / math.pi * 180.0\n",
    "v_fov_deg = math.atan2(sensor_height / 2, focal_length) * 2.0 / math.pi * 180.0\n",
    "print(\"Field of view horizontal: {:.3} °\".format(h_fov_deg))\n",
    "print(\"Field of view vertical:   {:.3} °\".format(v_fov_deg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Star trail length per second is 0.444 pixels\n",
      "Which means star moves by one pixel per 2.25 seconds\n"
     ]
    }
   ],
   "source": [
    "########################################\n",
    "# Star trails, Pavel Peřina 2028-08-14\n",
    "########################################\n",
    "\n",
    "# For focal length in meters and shutter speed in seconds, return trail length\n",
    "# on a sensor in meters. Do not use long shutter speed times longer than six\n",
    "# hours, star will get behind a camera\n",
    "def star_trail_length(focal_length: float, shutter_speed: float) -> float:\n",
    "    radians_per_second = 2.0 * math.pi / (24 * 60 * 60)\n",
    "    radians = radians_per_second * shutter_speed\n",
    "    # tan(radians) = trail_length / focal_length\n",
    "    trail_length: float = math.tan(radians) * focal_length\n",
    "    return trail_length\n",
    "\n",
    "star_trail_pixels_per_second = star_trail_length(focal_length, 1.0) / pixel_size\n",
    "star_trail_seconds_per_pixel = 1.0 / star_trail_pixels_per_second\n",
    "print(\"Star trail length per second is {:.3} pixels\".format(star_trail_pixels_per_second))\n",
    "print(\"Which means star moves by one pixel per {:.3} seconds\".format(star_trail_seconds_per_pixel))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sky coverage calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "####################################\n",
    "# Sky coverage (Monte-Carlo method)\n",
    "# Pavel Peřina, 2023-08-14\n",
    "####################################\n",
    "\n",
    "# Random point on a sphere\n",
    "# This function generates random point with coordinates from -1..1. \n",
    "# Checks if it's inside a unit sphere and repeat until condition is met\n",
    "# Then normalize point putting it onto sphere's surface\n",
    "def random_point_on_sphere():\n",
    "    while True:\n",
    "        x, y, z = np.random.uniform(-1, 1, 3)\n",
    "        mag_squared = x**2 + y**2 + z**2\n",
    "        if mag_squared <= 1:\n",
    "            mag = np.sqrt(mag_squared)\n",
    "            return np.array([x/mag, y/mag, z/mag, 1.0])\n",
    "\n",
    "# Assume we have camera looking in -z direction, y is up, x is right\n",
    "# This right-handed coordinate system is commonly used in 3D computer graphics\n",
    "# Returns camera normalized clipping planes (unit normal vector, distance)\n",
    "def camera_clipping_planes(sensorWidth: float, sensorHeight: float, focalLength: float):\n",
    "    hFovHalf = math.atan2(sensorWidth  / 2.0, focalLength)\n",
    "    vFovHalf = math.atan2(sensorHeight / 2.0, focalLength)\n",
    "    # With a lower FoV cos goes to 1.0 and sin goes to 0.0\n",
    "    result = [\n",
    "        np.array([+math.cos(hFovHalf), 0.0, -math.sin(hFovHalf), 0.0]), # left\n",
    "        np.array([-math.cos(hFovHalf), 0.0, -math.sin(hFovHalf), 0.0]), # right\n",
    "        np.array([0.0, +math.cos(vFovHalf), -math.sin(vFovHalf), 0.0]), # bottom\n",
    "        np.array([0.0, -math.cos(vFovHalf), -math.sin(vFovHalf), 0.0]), # top\n",
    "    ]\n",
    "    return result\n",
    "\n",
    "def inside_fov_ratio(sensorWidth: float, sensorHeight: float, focalLength: float, iteration_count: int) -> float:\n",
    "    clipping_planes = camera_clipping_planes(sensorWidth, sensorHeight, focalLength)\n",
    "    rejected_count: int = 0\n",
    "    for _ in range(iteration_count):\n",
    "        point = random_point_on_sphere()\n",
    "        for clipping_plane in clipping_planes:\n",
    "            signed_distance: float = np.dot(clipping_plane, point)\n",
    "            if (signed_distance < 0):\n",
    "                rejected_count += 1\n",
    "                break\n",
    "    accepted_count: int = iteration_count - rejected_count\n",
    "    return float(accepted_count) / iteration_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percents of sky captured: 9.16%\n"
     ]
    }
   ],
   "source": [
    "ratio = inside_fov_ratio(sensor_width, sensor_height, focal_length, 10000)\n",
    "# half of the sky is below horizon, convert ratio to percents\n",
    "print(\"Percents of sky captured: {}%\".format(ratio * 2.0 * 100.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "# Camera database\n",
    "# Pavel Peřina, 2023-08-16\n",
    "##########################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sqlite3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_db():\n",
    "    # Connect to an in-memory database\n",
    "    conn = sqlite3.connect(':memory:')\n",
    "    cursor = conn.cursor()\n",
    "    conn.row_factory = sqlite3.Row\n",
    "\n",
    "\n",
    "    cursor.execute('''\n",
    "    CREATE TABLE mounts (\n",
    "        mount_id TEXT PRIMARY KEY NOT NULL,\n",
    "        mount_name TEXT, \n",
    "        sensor_width REAL,\n",
    "        sensor_height REAL\n",
    "    )\n",
    "    ''')\n",
    "    cursor.executemany(\n",
    "        \"INSERT INTO mounts (mount_id, mount_name, sensor_width, sensor_height) VALUES (?, ?, ?, ?)\", \n",
    "        [('FujiX', 'APS-C', 23.5e-3, 15.6e-3),\n",
    "         ('X100V', 'APS-C', 23.5e-3, 15.6e-3),\n",
    "         ('M43', 'Micro 4/3', 17.3e-3, 13.0e-3),\n",
    "         ('FF', 'Full Frame', 36.0e-3, 24.0e-3),\n",
    "         ('1inch', '1\"', 12.8e-3, 9.6e-3)\n",
    "         ])\n",
    "\n",
    "    # Create the 'cameras' table\n",
    "    cursor.execute('''\n",
    "    CREATE TABLE cameras (\n",
    "        camera_id INTEGER PRIMARY KEY,\n",
    "        mount_id INTEGER,\n",
    "        camera_name TEXT,\n",
    "        image_width INTEGER,\n",
    "        image_height INTEGER\n",
    "    )''')\n",
    "    cursor.executemany(\n",
    "        \"INSERT INTO cameras (mount_id, camera_name, image_width, image_height) VALUES (?,?,?,?)\",\n",
    "         [('X100V', 'Fujifilm X100V', 6240, 4160),\n",
    "          ('M43', 'Panasonic G80', 4592, 3448),\n",
    "          ])\n",
    "\n",
    "    cursor.execute('''\n",
    "    CREATE TABLE lenses (\n",
    "        lens_id INTEGER PRIMARY KEY,\n",
    "        lens_name TEXT,\n",
    "        focal_length REAL,\n",
    "        mount_id TEXT NOT NULL\n",
    "    )''')\n",
    "    cursor.executemany(\n",
    "        \"INSERT INTO lenses (mount_id, lens_name, focal_length) VALUES (?,?,?)\",\n",
    "         [('X100V', 'X100V 23mm', 23.0e-3),\n",
    "          ('FujiX', 'FujiX 14mm', 14.0e-3),\n",
    "          ('FujiX', 'FujiX 18mm', 18.0e-3),\n",
    "          ('FujiX', 'FujiX 16mm', 16.0e-3),\n",
    "          ('FujiX', 'FujiX 23mm', 23.0e-3),\n",
    "          ('M43', 'Panasonic 20mm', 20.0e-3),\n",
    "          ('M43', 'Panasonic 14mm', 14.0e-3),\n",
    "          ('M43', 'Panasonic 25mm', 25.0e-3), \n",
    "          ('M43', 'Panasonic 8mm',   8.0e-3), \n",
    "          ('M43', 'Panasonic 14mm', 18.0e-3),\n",
    "          ('M43', 'Panasonic 12mm', 12.0e-3),\n",
    "          ('FF', 'Full Frame 14mm', 14.0e-3),\n",
    "          ('FF', 'Full Frame 18mm', 18.0e-3),\n",
    "          ('FF', 'Full Frame 24mm', 24.0e-3),\n",
    "          ('FF', 'Full Frame 28mm', 28.0e-3),\n",
    "          ('FF', 'Full Frame 35mm', 35.0e-3),\n",
    "          ('FF', 'Full Frame 50mm', 50.0e-3),\n",
    "          ])\n",
    "    \n",
    "    conn.commit()\n",
    "    return conn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "db_conn = create_db()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "| Camera | Resolution [MPix] | Pixel Size [um] | Focal Length[mm] | Seconds per Pixel [1/s] |\n",
      "|---|---:|--:|--:|--:|\n",
      "| Fujifilm X100V | 26.0 | 3.77 | 23.0 | 2.25\n",
      "| Panasonic G80 | 15.8 | 3.77 |  8.0 | 6.48\n",
      "| Panasonic G80 | 15.8 | 3.77 | 12.0 | 4.32\n",
      "| Panasonic G80 | 15.8 | 3.77 | 14.0 | 3.70\n",
      "| Panasonic G80 | 15.8 | 3.77 | 18.0 | 2.88\n",
      "| Panasonic G80 | 15.8 | 3.77 | 20.0 | 2.59\n",
      "| Panasonic G80 | 15.8 | 3.77 | 25.0 | 2.07\n"
     ]
    }
   ],
   "source": [
    "cursor = db_conn.cursor()\n",
    "cursor.execute(\"SELECT mount_name, sensor_width, sensor_height, camera_name, image_width, image_height, focal_length FROM cameras JOIN lenses on lenses.mount_id = cameras.mount_id JOIN mounts on mounts.mount_id = cameras.mount_id ORDER BY sensor_width DESC, focal_length\")\n",
    "results = cursor.fetchall()\n",
    "print(\"| Camera | Resolution [MPix] | Pixel Size [um] | Focal Length[mm] | Seconds per Pixel [1/s] |\")\n",
    "print(\"|---|---:|--:|--:|--:|\")\n",
    "for row in results:\n",
    "    focal_length: float = row[\"focal_length\"]\n",
    "    sensor_width: float = row[\"sensor_width\"]\n",
    "    sensor_height: float = row[\"sensor_height\"]\n",
    "    image_width: int = row[\"image_width\"]\n",
    "    image_height: int = row[\"image_height\"]\n",
    "    pixel_size = sensor_width / image_width\n",
    "    pixels = image_width * image_height\n",
    "    star_trail_pixels_per_second = star_trail_length(focal_length, 1.0) / pixel_size\n",
    "    star_trail_seconds_per_pixel = 1.0 / star_trail_pixels_per_second\n",
    "    print(\"| {} | {:.3} | {:.3} | {:4.1f} | {:4.2f}\".format(row[\"camera_name\"], float(pixels)/1.0e6, pixel_size*1e6, focal_length*1e3, star_trail_seconds_per_pixel))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "| Camera Format | Focal Length[mm] | HFov[°] | VFoV[°] | Sky Coverage[%] |\n",
      "|----|---:|--:|--:|--:|\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "| Full Frame | 14 | 104.3 | 81.2 | 34.0 |\n",
      "| Full Frame | 18 | 90.0 | 67.4 | 25.4 |\n",
      "| Full Frame | 24 | 73.7 | 53.1 | 17.3 |\n",
      "| Full Frame | 28 | 65.5 | 46.4 | 13.6 |\n",
      "| Full Frame | 35 | 54.4 | 37.8 |  9.5 |\n",
      "| Full Frame | 50 | 39.6 | 27.0 |  4.9 |\n",
      "| APS-C | 14 | 80.0 | 58.2 | 20.1 |\n",
      "| APS-C | 16 | 72.6 | 52.0 | 16.6 |\n",
      "| APS-C | 18 | 66.3 | 46.9 | 13.8 |\n",
      "| APS-C | 23 | 54.1 | 37.5 |  9.4 |\n",
      "| APS-C | 23 | 54.1 | 37.5 |  9.3 |\n",
      "| Micro 4/3 |  8 | 94.5 | 78.2 | 30.7 |\n",
      "| Micro 4/3 | 12 | 71.6 | 56.9 | 18.0 |\n",
      "| Micro 4/3 | 14 | 63.4 | 49.8 | 13.8 |\n",
      "| Micro 4/3 | 18 | 51.3 | 39.7 |  9.3 |\n",
      "| Micro 4/3 | 20 | 46.8 | 36.0 |  7.9 |\n",
      "| Micro 4/3 | 25 | 38.2 | 29.1 |  5.4 |\n"
     ]
    }
   ],
   "source": [
    "cursor = db_conn.cursor()\n",
    "cursor.execute(\"SELECT mount_name, sensor_width, sensor_height, focal_length FROM lenses JOIN mounts on mounts.mount_id = lenses.mount_id ORDER BY sensor_width DESC, focal_length\")\n",
    "results = cursor.fetchall()\n",
    "print(\"| Camera Format | Focal Length[mm] | HFov[°] | VFoV[°] | Sky Coverage[%] |\")\n",
    "print(\"|----|---:|--:|--:|--:|\")\n",
    "for row in results:\n",
    "    focal_length: float =row[\"focal_length\"]\n",
    "    sensor_width: float =row[\"sensor_width\"]\n",
    "    sensor_height: float = row[\"sensor_height\"]\n",
    "    hFovDeg = math.atan2(sensor_width  / 2, focal_length) * 2.0 / math.pi * 180.0\n",
    "    vFovDeg = math.atan2(sensor_height / 2, focal_length) * 2.0 / math.pi * 180.0\n",
    "    coverage = inside_fov_ratio(sensor_width, sensor_height, focal_length, 100000) * 2.0 * 100.0\n",
    "    print(\"| {} | {:2.0f} | {:4.1f} | {:4.1f} | {:4.1f} |\".format(row[\"mount_name\"], focal_length*1e3, hFovDeg, vFovDeg, coverage))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.11.4 (main, Jun  7 2023, 00:00:00) [GCC 13.1.1 20230511 (Red Hat 13.1.1-2)]\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print(sys.version)"
   ]
  }
 ],
 "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.11.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Cross link to blog post: https://www.pavelp.cz/posts/eng-photo-perseids2/