3. Star mapping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 3: Star Mapping\n",
    "\n",
    "Prepared by the International Centre for Neuromorphic Systems (ICNS) at Western Sydney University, adapted from the short course presented at the AMOS Conference 2023.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "***What you will learn***\n",
    "In this tutorial, we will work with real-world event-based data to learn how to load and process it. \n",
    "\n",
    "By the end of this tutorial, you will...\n",
    "* Be able to load and view event-based star data\n",
    "* Warp event temporally to stabilise stars\n",
    "* Find the centres of the stars\n",
    "* Perform an astrometry fit on the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting Started\n",
    "\n",
    "This tutorial builds upon the previous tutorials and exercises but requires the installation of a few additional packages.\n",
    "\n",
    "First, we will need the fantastic [astropy](https://www.astropy.org/) library to work with star coordinate systems. We will also make use of a Python package called [astrometry](https://github.com/neuromorphicsystems/astrometry), which is a Python wrapper around Astrometry.net that accepts star positions and can be used to perform plate solving. These two packages can be installed using the following command:\n",
    "\n",
    "```sh\n",
    "python3 -m pip install astropy requests astrometry\n",
    "```\n",
    "\n",
    "We will also require a few additional packages for plotting and clustering. These can be installed using the following command:\n",
    "```sh\n",
    "python3 -m pip install scikit-image matplotlib-label-lines ipywidgets\n",
    "```\n",
    "\n",
    "**Note:** If you are using `Anaconda`, don't forget to use `conda` instead of pip!\n",
    "\n",
    "If you are running this tutorial on Google Colaboratory, simply run the Python code below to install the required libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install the required dependencies (if not already installed). This is required to install the required packages in Google Colaboratory.\n",
    "!pip3 install event_stream astropy requests astrometry scikit-image matplotlib-label-lines ipywidgets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Downloading the data\n",
    "\n",
    "If you are running this tutorial on Google Colab, you will need to download the required data file from Github. The command below will automatically download the correct file and place it in the same directory as this tutorial. If you are running the tutorial locally, you may already have this file locally and do not need to download it again. The data file can be found in the Github repository at https://github.com/neuromorphicsystems/tutorials/raw/main/data/star_cut.es.\n",
    "\n",
    "Note that this data file is approximtely 21 MB in size. It may take a few minutes to download."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "\n",
    "\n",
    "# Download the sample data from Github (55 MB)\n",
    "stars_path, _ = urllib.request.urlretrieve(\n",
    "    \"https://github.com/neuromorphicsystems/tutorials/raw/main/data/stars_cut.es\",\n",
    "    \"stars_cut.es\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to load the data. We for simplicity and ease-of-use, we will load all the data into a single variable in memory. Then we'll accumulate all events into one frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import event_stream\n",
    "import numpy\n",
    "from tonic import transforms\n",
    "\n",
    "\n",
    "with event_stream.Decoder(stars_path) as decoder:\n",
    "    width = decoder.width\n",
    "    height = decoder.height\n",
    "    events = numpy.concatenate([chunk for chunk in decoder])\n",
    "\n",
    "# Tonic requires the events to be a structured numpy array with a 'p' field\n",
    "events = events.view(dtype=events.dtype.descr[:3] + [(('on', 'p'), '?')])\n",
    "\n",
    "image_transform = transforms.ToImage(sensor_size=(decoder.width, decoder.height, 2))\n",
    "accumulated_frame = image_transform(events)\n",
    "# Flip the image vertically and sum the two channels\n",
    "accumulated_frame = numpy.flip(accumulated_frame.sum(0), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualising the event-based data\n",
    "\n",
    "Let's plot the accumulated frame to show a summary of all the data in the file. We can do that as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.colors\n",
    "import matplotlib.pyplot\n",
    "import matplotlib.transforms\n",
    "\n",
    "matplotlib.rcParams[\"figure.figsize\"] = [12, 8]\n",
    "\n",
    "matplotlib.pyplot.imshow(\n",
    "    accumulated_frame,\n",
    "    norm=matplotlib.colors.LogNorm(vmax=numpy.percentile(accumulated_frame, 99.9)),\n",
    "    cmap=\"magma\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Event Warping\n",
    "\n",
    "Now, we will warp the events based on the known sidereal motion of the stars. We can do that using the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx = -1.861  # px/s\n",
    "vy = 0.549  # px/s\n",
    "\n",
    "t0 = events[\"t\"][0]\n",
    "duration = events[\"t\"][-1] - t0\n",
    "\n",
    "warped_events = numpy.zeros(len(events), dtype=events.dtype)\n",
    "\n",
    "# warp the events (shear operation)\n",
    "warped_events[\"t\"] = events[\"t\"]\n",
    "if vx > 0:\n",
    "    warped_events[\"x\"] = events[\"x\"] + numpy.floor(\n",
    "        (vx / 1e6) * (duration - (events[\"t\"] - t0))\n",
    "    )\n",
    "else:\n",
    "    warped_events[\"x\"] = events[\"x\"] + numpy.floor((-vx / 1e6) * (events[\"t\"] - t0))\n",
    "if vy > 0:\n",
    "    warped_events[\"y\"] = events[\"y\"] + numpy.floor(\n",
    "        (vy / 1e6) * (duration - (events[\"t\"] - t0))\n",
    "    )\n",
    "else:\n",
    "    warped_events[\"y\"] = events[\"y\"] + numpy.floor((-vy / 1e6) * (events[\"t\"] - t0))\n",
    "warped_events[\"on\"] = events[\"on\"]\n",
    "\n",
    "# accumulate the warped events\n",
    "warped_accumulated_frame = numpy.ones(\n",
    "    (\n",
    "        height + int(numpy.floor(abs(vy / 1e6) * duration)),\n",
    "        width + int(numpy.floor(abs(vx / 1e6) * duration)),\n",
    "    )\n",
    ")\n",
    "numpy.add.at(\n",
    "    warped_accumulated_frame,\n",
    "    (warped_accumulated_frame.shape[0] - 1 - warped_events[\"y\"], warped_events[\"x\"]),\n",
    "    1,\n",
    ")\n",
    "\n",
    "# display the result\n",
    "matplotlib.pyplot.imshow(\n",
    "    warped_accumulated_frame,\n",
    "    norm=matplotlib.colors.LogNorm(),\n",
    "    cmap=\"magma\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets\n",
    "\n",
    "@ipywidgets.interact(vx = -1.861, vy = 0.549)\n",
    "def warp_events(vx, vy):\n",
    "    t0 = events[\"t\"][0]\n",
    "    duration = events[\"t\"][-1] - t0\n",
    "    warped_events = numpy.zeros(len(events), dtype=events.dtype)\n",
    "\n",
    "    # warp the events (shear operation)\n",
    "    warped_events[\"t\"] = events[\"t\"]\n",
    "    if vx > 0:\n",
    "        warped_events[\"x\"] = events[\"x\"] + numpy.floor(\n",
    "            (vx / 1e6) * (duration - (events[\"t\"] - t0))\n",
    "        )\n",
    "    else:\n",
    "        warped_events[\"x\"] = events[\"x\"] + numpy.floor((-vx / 1e6) * (events[\"t\"] - t0))\n",
    "    if vy > 0:\n",
    "        warped_events[\"y\"] = events[\"y\"] + numpy.floor(\n",
    "            (vy / 1e6) * (duration - (events[\"t\"] - t0))\n",
    "        )\n",
    "    else:\n",
    "        warped_events[\"y\"] = events[\"y\"] + numpy.floor((-vy / 1e6) * (events[\"t\"] - t0))\n",
    "    warped_events[\"on\"] = events[\"on\"]\n",
    "\n",
    "    # accumulate the warped events\n",
    "    warped_accumulated_frame = numpy.ones(\n",
    "        (\n",
    "            height + int(numpy.floor(abs(vy / 1e6) * duration)),\n",
    "            width + int(numpy.floor(abs(vx / 1e6) * duration)),\n",
    "        )\n",
    "    )\n",
    "    numpy.add.at(\n",
    "        warped_accumulated_frame,\n",
    "        (warped_accumulated_frame.shape[0] - 1 - warped_events[\"y\"], warped_events[\"x\"]),\n",
    "        1,\n",
    "    )\n",
    "\n",
    "    # display the result\n",
    "    matplotlib.pyplot.imshow(\n",
    "        warped_accumulated_frame,\n",
    "        norm=matplotlib.colors.LogNorm(),\n",
    "        cmap=\"magma\",\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering the image to remove noise\n",
    "\n",
    "The above image looks far better, but there are still some events that have left small trails in the data. We can clean this up slightly by applying a simple median filter to the image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import PIL.Image\n",
    "import PIL.ImageChops\n",
    "import PIL.ImageFilter\n",
    "\n",
    "image = PIL.Image.fromarray(warped_accumulated_frame)\n",
    "filtered_image = image.filter(PIL.ImageFilter.MedianFilter(5))\n",
    "filtered_accumulated_frame = numpy.array(filtered_image)\n",
    "\n",
    "matplotlib.pyplot.imshow(\n",
    "    filtered_accumulated_frame,\n",
    "    norm=matplotlib.colors.LogNorm(),\n",
    "    cmap=\"magma\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Labelling stars in the image\n",
    "\n",
    "To find the stars within this image, we will use a labeling technique provided by `scikit.image`. We will be using the [label](https://scikit-image.org/docs/stable/api/skimage.measure.html#skimage.measure.label) function, which looks for connected regions within the image and assigns them individual and unique labels. The label function considers pixels connected if they are adjacent and have the same value. \n",
    "\n",
    "In order to improve the labeling process, we will take the above image and make a binary image out of it. We can do this using the 99.5% percentile to remove a large portion of the noise in the field of view. We can use the code below to make the binary mask of our star map: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nonzero = filtered_accumulated_frame[filtered_accumulated_frame > 0.0]\n",
    "binary_mask = filtered_accumulated_frame > numpy.percentile(nonzero, 99.5)\n",
    "\n",
    "matplotlib.pyplot.imshow(\n",
    "    binary_mask,\n",
    "    cmap=\"gray\",\n",
    "    interpolation=None,\n",
    "    interpolation_stage=\"rgba\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below uses the binary mask to label the stars in the image. It generates a random colour map based on the number of stars identified in the recording."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage.measure\n",
    "\n",
    "labels_array, maximum_label = skimage.measure.label(binary_mask, connectivity=1, background=0, return_num=True)\n",
    "\n",
    "# manually generate a colormap to visualise the result of connected components analysis\n",
    "labels_colors = numpy.array([[0.0, 0.0, 0.0, 1.0]])\n",
    "while len(labels_colors) < maximum_label + 1:\n",
    "    new_colors = matplotlib.colormaps[\"Set2\"](numpy.linspace(1/16, 15/16, 8))[0:maximum_label + 1 - len(labels_colors)]\n",
    "    numpy.random.shuffle(new_colors)\n",
    "    labels_colors = numpy.concatenate([labels_colors, new_colors])\n",
    "\n",
    "matplotlib.pyplot.imshow(\n",
    "    labels_array + 0.5,\n",
    "    cmap=matplotlib.colors.ListedColormap(labels_colors),\n",
    "    interpolation=None,\n",
    "    interpolation_stage=\"rgba\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Labeling the raw events\n",
    "\n",
    "The code above detected stars using the accumulated frame of events generated from the event-based data. We therefore have labelled regions of the above image and now need to map that back to the raw events received from the camera. To do this, we will extend the definition of the event to include an additional value indicating the label to which the event is attached.\n",
    "\n",
    "We can do this on the warped events and apply these labels to the events directly. We could also reverse the event-warping and label the original events too. However, this isn't necessary for star mapping.\n",
    "\n",
    "Note that we added a default background label, and therefore any events that are not stars will be mapped to the background label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labelled_events = numpy.zeros(\n",
    "    len(warped_events),\n",
    "    dtype=[\n",
    "        (\"t\", \"<u8\"),\n",
    "        (\"x\", \"<u2\"),\n",
    "        (\"y\", \"<u2\"),\n",
    "        (\"on\", \"?\"),\n",
    "        (\"label\", \"<u4\"),\n",
    "    ],\n",
    ")\n",
    "\n",
    "labelled_events[\"t\"] = warped_events[\"t\"]\n",
    "labelled_events[\"x\"] = warped_events[\"x\"]\n",
    "labelled_events[\"y\"] = warped_events[\"y\"]\n",
    "labelled_events[\"on\"] = warped_events[\"on\"]\n",
    "labelled_events[\"label\"] = labels_array[warped_accumulated_frame.shape[0] - 1 - warped_events[\"y\"], warped_events[\"x\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the star centres\n",
    "\n",
    "Now that we have labelled events, we can then iterate through them and calculate the mean position of the events to get the star centre. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "centers = numpy.zeros((maximum_label + 1, 2))\n",
    "for label in range(1, maximum_label + 1):\n",
    "    centers[label] = (\n",
    "        numpy.mean(labelled_events[\"x\"][labelled_events[\"label\"] == label]),\n",
    "        numpy.mean(labelled_events[\"y\"][labelled_events[\"label\"] == label]),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the labelled stars\n",
    "\n",
    "The code below plots a weighted version of the star map and the plots circles based on the detected star circles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labelled_accumulated_frame = numpy.zeros(\n",
    "    (\n",
    "        3,\n",
    "        height + int(numpy.floor(abs(vy / 1e6) * duration)),\n",
    "        width + int(numpy.floor(abs(vx / 1e6) * duration)),\n",
    "    )\n",
    ")\n",
    "for label in range(1, maximum_label + 1):\n",
    "    source_events = labelled_events[labelled_events[\"label\"] == label]\n",
    "    source_alpha = numpy.zeros(\n",
    "        (labelled_accumulated_frame.shape[1], labelled_accumulated_frame.shape[2])\n",
    "    )\n",
    "    numpy.add.at(\n",
    "        source_alpha,\n",
    "        (\n",
    "            source_alpha.shape[0] - 1 - source_events[\"y\"],\n",
    "            source_events[\"x\"],\n",
    "        ),\n",
    "        1,\n",
    "    )\n",
    "    source_alpha = (source_alpha / source_alpha.max()) ** 0.25\n",
    "    source_accumulated_frame = numpy.zeros(labelled_accumulated_frame.shape)\n",
    "    source_accumulated_frame[0, :, :] = labels_colors[label][0]\n",
    "    source_accumulated_frame[1, :, :] = labels_colors[label][1]\n",
    "    source_accumulated_frame[2, :, :] = labels_colors[label][2]\n",
    "    labelled_accumulated_frame = numpy.multiply(\n",
    "        labelled_accumulated_frame, 1.0 - source_alpha\n",
    "    ) + numpy.multiply(source_accumulated_frame, source_alpha)\n",
    "labelled_accumulated_frame = numpy.transpose(labelled_accumulated_frame, (1, 2, 0))\n",
    "matplotlib.pyplot.imshow(labelled_accumulated_frame)\n",
    "matplotlib.pyplot.scatter(\n",
    "    x=centers[:, 0],\n",
    "    y=source_alpha.shape[0] - 1 - centers[:, 1],\n",
    "    marker=\"o\",\n",
    "    facecolor=\"none\",\n",
    "    edgecolors=\"#ffffff\",\n",
    "    linewidth=0.4,\n",
    "    s=100,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performing the astrometry fit\n",
    "\n",
    "**WARNING: The following code will download the astrometry caches files when first run. This will be approximately 650 Mb in size.**\n",
    "\n",
    "Now that we have a list of star centres, we can now use the `astrometry` library to perform an astrometric fit on the data. The `astrometry` library is build around `Astrometry.Net` and was specially written to perform plate-solving on lists of star positions, rather than images. It provides a high-level programmatic interface and automatically handles the downloading of the necessary index files.\n",
    "\n",
    "To perform the star mapping, we can use the following code to do the plate-solving:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astrometry\n",
    "import logging\n",
    "\n",
    "logging.getLogger().setLevel(logging.INFO)\n",
    "\n",
    "# see https://github.com/neuromorphicsystems/astrometry for documentation\n",
    "\n",
    "solver = astrometry.Solver(\n",
    "    astrometry.series_5200_heavy.index_files(\n",
    "        cache_directory=\"astrometry_cache\",\n",
    "        scales={6},\n",
    "    )\n",
    ")\n",
    "\n",
    "solution = solver.solve(\n",
    "    stars=centers,\n",
    "    size_hint=None,\n",
    "    position_hint=None,\n",
    "    solution_parameters=astrometry.SolutionParameters(),\n",
    ")\n",
    "\n",
    "assert solution.has_match()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the results\n",
    "\n",
    "Now that we have a fit, we can use `astropy` to convert from the RA/DEC coordinates back to pixel space. We can also retrieve the star magnitudes from the `astrometry` fit. we will use this to scale the star accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.wcs\n",
    "import labellines\n",
    "\n",
    "match = solution.best_match()\n",
    "wcs = astropy.wcs.WCS(match.wcs_fields)\n",
    "\n",
    "# convert stars RA/Dec to pixels\n",
    "stars = wcs.all_world2pix(\n",
    "    [[star.ra_deg, star.dec_deg] for star in match.stars],\n",
    "    0,\n",
    ")\n",
    "\n",
    "# retrieve stars magnitudes (requires series_5200_heavy)\n",
    "magnitudes = numpy.array([star.metadata[\"mag\"] for star in match.stars])\n",
    "scaled_magnitudes = 1.0 - (magnitudes - magnitudes.min()) / (\n",
    "    magnitudes.max() - magnitudes.min()\n",
    ")\n",
    "\n",
    "# draw the accumulated frame\n",
    "matplotlib.pyplot.imshow(labelled_accumulated_frame)\n",
    "\n",
    "# plot grid lines (one line per degree)\n",
    "grid_ra = numpy.arange(\n",
    "    numpy.floor(match.center_ra_deg - match.scale_arcsec_per_pixel / 3600 * width * 2),\n",
    "    numpy.ceil(match.center_ra_deg + match.scale_arcsec_per_pixel / 3600 * width * 2),\n",
    ")\n",
    "grid_dec = numpy.arange(\n",
    "    numpy.floor(match.center_dec_deg - match.scale_arcsec_per_pixel / 3600 * width * 2),\n",
    "    numpy.ceil(match.center_dec_deg + match.scale_arcsec_per_pixel / 3600 * width * 2),\n",
    ")\n",
    "for ra in grid_ra:\n",
    "    points = wcs.all_world2pix([[ra, dec] for dec in grid_dec], 0)\n",
    "    matplotlib.pyplot.plot(\n",
    "        points[:, 0],\n",
    "        source_alpha.shape[0] - 1 - points[:, 1],\n",
    "        color=\"#aaaaaa\",\n",
    "        linewidth=0.4,\n",
    "    )\n",
    "    lines = matplotlib.pyplot.gca().get_lines()\n",
    "    labellines.labelLine(\n",
    "        line=lines[-1],\n",
    "        x=points[len(points) // 2, 0] - 75,\n",
    "        label=f\"RA = {ra:.0f}º\",\n",
    "        backgroundcolor=\"none\",\n",
    "        outline_color=\"#000000\",\n",
    "    )\n",
    "\n",
    "for dec in grid_dec:\n",
    "    points = wcs.all_world2pix([[ra, dec] for ra in grid_ra], 0)\n",
    "    matplotlib.pyplot.plot(\n",
    "        points[:, 0],\n",
    "        source_alpha.shape[0] - 1 - points[:, 1],\n",
    "        color=\"#aaaaaa\",\n",
    "        linewidth=0.4,\n",
    "    )\n",
    "    lines = matplotlib.pyplot.gca().get_lines()\n",
    "    labellines.labelLine(\n",
    "        line=lines[-1],\n",
    "        x=points[len(points) // 2, 0] - 200,\n",
    "        label=f\"Dec = {dec:.0f}º\",\n",
    "        backgroundcolor=\"none\",\n",
    "        outline_color=\"#000000\",\n",
    "    )\n",
    "\n",
    "# plot events centers of mass\n",
    "matplotlib.pyplot.scatter(\n",
    "    x=centers[:, 0],\n",
    "    y=source_alpha.shape[0] - 1 - centers[:, 1],\n",
    "    marker=\"o\",\n",
    "    facecolor=\"none\",\n",
    "    edgecolors=\"#ffffff\",\n",
    "    linewidth=0.4,\n",
    "    s=100,\n",
    ")\n",
    "\n",
    "# plot stars\n",
    "matplotlib.pyplot.scatter(\n",
    "    x=stars[:, 0],\n",
    "    y=source_alpha.shape[0] - 1 - stars[:, 1],\n",
    "    marker=\"o\",\n",
    "    facecolor=\"none\",\n",
    "    edgecolors=[\n",
    "        f\"#de5246{int(magnitude):02x}\"\n",
    "        for magnitude in numpy.floor((scaled_magnitudes**0.5) * 255)\n",
    "    ],\n",
    "    linewidth=1.0,\n",
    "    s=scaled_magnitudes * 200.0,\n",
    ")\n",
    "\n",
    "matplotlib.pyplot.xlim((0.0, source_alpha.shape[1]))\n",
    "matplotlib.pyplot.ylim((source_alpha.shape[0], 0.0))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 ('.venv': venv)",
   "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.17"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "07c4f84319129d1c976bad471c8642b46ea9b7b96fed5777f2c683a8a5a74ac2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}