{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "15505fcf-22cc-4d8e-8b83-4a9210e0f2ed",
   "metadata": {},
   "source": [
    "# `Shadow2DExclusion` Overview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e0a519a-1b16-4242-99f1-2cfd036802fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46416d10-4fd6-4454-a1c1-f3b9f1826ad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import xarray as xr\n",
    "\n",
    "from matplotlib.patches import Polygon\n",
    "from typing import Tuple\n",
    "\n",
    "sys.executable;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c757daf-b507-4119-acbe-135a6a91d197",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    from bapsf_motion.motion_builder.exclusions import (\n",
    "        CircularExclusion,\n",
    "        DividerExclusion,\n",
    "        GovernExclusion,\n",
    "        Shadow2DExclusion,\n",
    "    )\n",
    "except ModuleNotFoundError:\n",
    "    from pathlib import Path\n",
    "\n",
    "    HERE = Path().cwd()\n",
    "    BAPSF_MOTION = (HERE / \"..\" / \"..\" / \"..\" ).resolve()\n",
    "    sys.path.append(str(BAPSF_MOTION))\n",
    "    \n",
    "    from bapsf_motion.motion_builder.exclusions import (\n",
    "        CircularExclusion,\n",
    "        DividerExclusion,\n",
    "        GovernExclusion,\n",
    "        Shadow2DExclusion,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ec82b7e-8d2c-4d1c-a4b6-b68d0c359109",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update(\n",
    "    {\n",
    "        # \"figure.figsize\": [12, 0.56 * 12],\n",
    "        \"figure.figsize\": [10, 0.8 * 10],\n",
    "        \"font.size\": 16,\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a682e855-dceb-4d11-97b7-72584acda59f",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Usage\n",
    "\n",
    "Direct usage should never be needed, since the `MotionBuilder` will handle this given the correct configuration is given to `MotionBuilder`.  The appropriate TOML or dictionary like configurations can be found in the documentation for `Shadow2DExclusion`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dc35d2b-f2f3-4f46-a08e-06af38014c98",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Is a `GovernExclusion`\n",
    "\n",
    "`Shadow2DExclusion` is a subclass of `GovernExclusion`.  This means the mask generated by `Shadow2DExclusion` will examine the existing global mask to generate its own mask, and that generated mask will replace the global mask.  As a result, there should only be one `GovernExclusion` used wheneve constructing a motion space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d125417f-116a-434c-ba53-2a58312092f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "issubclass(Shadow2DExclusion, GovernExclusion)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a75ae2d-ff4c-4d3b-bd97-ef57886c2cba",
   "metadata": {},
   "source": [
    "### Defining a `Shadow2DExclusion`\n",
    "\n",
    "Create an initial global mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e2a92bc-0dd7-4957-85d6-59c54635b7c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 111\n",
    "side = np.linspace(-55, 55, num=size)\n",
    "ds = xr.Dataset(\n",
    "    {\"mask\": ((\"x\", \"y\"), np.ones((size, size), dtype=bool))},\n",
    "    coords={\n",
    "        \"x\": side,\n",
    "        \"y\": side,\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e2fe97b-9127-4024-9239-4241b4e02169",
   "metadata": {},
   "source": [
    "Add some initial exlustion layers for `Shadow2DExclusion` to act on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e333ca09-ac3d-444c-a0e3-752bd1a94adf",
   "metadata": {},
   "outputs": [],
   "source": [
    "ex_1 = CircularExclusion(ds, radius=10, center=(0,0), exclude=\"inside\")\n",
    "ex_2 = CircularExclusion(ds, radius=8, center=(25,30), exclude=\"inside\")\n",
    "\n",
    "ds.mask.plot(x=\"x\", y=\"y\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a619a35b-1c20-430f-8956-d7155a7b2484",
   "metadata": {},
   "source": [
    "Now add the `Shadow2DExclusion` Layer, which has only one configuration parameter `source_point`.  `source_point` can be view as the \"light\" point source for the shadowing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77c8e1dc-5c10-4302-afb0-6a13356f0fc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ex = Shadow2DExclusion(ds, source_point=[45,-58])\n",
    "\n",
    "ds.mask.plot(x=\"x\", y=\"y\");\n",
    "plt.scatter(\n",
    "    ex.source_point[0],\n",
    "    ex.source_point[1],\n",
    "    marker=\"+\",\n",
    "    color=\"black\",\n",
    "    s=8**2,\n",
    ")\n",
    "axis = plt.gca()\n",
    "axis.set_xlim(-60, 60)\n",
    "axis.set_ylim(-60, 60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d03fd68-ca68-4570-9362-d5c8f4cf2e7b",
   "metadata": {},
   "source": [
    "The configuration dictionary for this exclusion looks like..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86ef6963-c6bf-4860-a8f3-fb3f91f61a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ex.config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd0451bc-2f59-4bc5-8fbc-5a2ce3026834",
   "metadata": {},
   "source": [
    "## The Underlying Algorithm\n",
    "\n",
    "Lets first reset the initial mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9fd4a86-4d8c-4b1c-8455-e8e6121a4756",
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 111\n",
    "side = np.linspace(-55, 55, num=size)\n",
    "ds = xr.Dataset(\n",
    "    {\"mask\": ((\"x\", \"y\"), np.ones((size, size), dtype=bool))},\n",
    "    coords={\n",
    "        \"x\": side,\n",
    "        \"y\": side,\n",
    "    },\n",
    ")\n",
    "\n",
    "ex_1 = CircularExclusion(ds, radius=10, center=(0,0), exclude=\"inside\")\n",
    "ex_2 = CircularExclusion(ds, radius=8, center=(25,30), exclude=\"inside\")\n",
    "\n",
    "ds.mask.plot(x=\"x\", y=\"y\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1cdc8c7-f101-46c9-b850-ff4e039f2c53",
   "metadata": {},
   "source": [
    "Lets first define some key parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d08f9f0d-b826-43b8-89be-941078e191dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gather some mask parameters\n",
    "dx, dy = ex_1.mask_resolution\n",
    "xkey, ykey = ex_1.mspace_dims\n",
    "x_coord, y_coord = ex_1.mspace_coords[xkey], ex_1.mspace_coords[ykey]\n",
    "\n",
    "# point source\n",
    "source_point = np.array([45, -58])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a83aab6-74e9-4e05-82f2-cda8853cfd9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define boundary of the motion space\n",
    "# - index_0 = the N-th side\n",
    "# - index_1 = start (0) and stop (1) of the side\n",
    "# - index_2 = (x, y) coordinate of the associated point\n",
    "#\n",
    "boundaries = np.zeros((4, 2, 2))\n",
    "\n",
    "x_min = x_coord[0] - 0.5 * dx\n",
    "x_max = x_coord[-1] + 0.5 * dx\n",
    "y_min =y_coord[0] - 0.5 * dy\n",
    "y_max = y_coord[-1] + 0.5 * dy\n",
    "\n",
    "# lower horizontal\n",
    "boundaries[0, 0, :] = [x_min, y_min]\n",
    "boundaries[0, 1, :] = [x_max, y_min]\n",
    "\n",
    "# right vertical\n",
    "boundaries[1, 0, :] = [x_max, y_min]\n",
    "boundaries[1, 1, :] = [x_max, y_max]\n",
    "\n",
    "# upper horizontal\n",
    "boundaries[2, 0, :] = [x_max, y_max]\n",
    "boundaries[2, 1, :] = [x_min, y_max]\n",
    "\n",
    "# left vertical\n",
    "boundaries[3, 0, :] = [x_min, y_max]\n",
    "boundaries[3, 1, :] = [x_min, y_min]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd3ef8e4-be8f-41f2-a314-243711e9dc4f",
   "metadata": {},
   "source": [
    "Now we need to determine if the `source_point` is inside the defined motion space or outside, and which boundary sides the rays would pass through if the source point is outside.  A few assumptions:\n",
    "\n",
    "1. Assume there are exclusion layers to shadow.  If the motion space mask is all `True` or all `False` then the resulting mask is the same.\n",
    "2. The motion space outside a boundary side is `True` if the side is a side the point source would shine through (outside-to-inside), otherwise it is considered `False`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "259d46f4-1767-4d71-b31b-16912b326a9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _determine_insertion_edge_indicies():\n",
    "    x_range = [x_coord[0] - 0.5 * dx, x_coord[-1] + 0.5 * dx]\n",
    "    y_range = [y_coord[0] - 0.5 * dy, y_coord[-1] + 0.5 * dy]\n",
    "\n",
    "    if (\n",
    "        (x_range[0] <= source_point[0] <= x_range[1])\n",
    "        and (y_range[0] <= source_point[1] <= y_range[1])\n",
    "    ):\n",
    "        # insertion point is within the motion space\n",
    "        return ()\n",
    "\n",
    "    insertion_edge_indices = []\n",
    "\n",
    "    deltas = boundaries[..., 1, :] - boundaries[..., 0, :]\n",
    "\n",
    "    for _orientation, _index in zip([\"horizontal\", \"vertical\"], [1, 0]):\n",
    "        _indices = np.where(np.isclose(deltas[..., _index], 0))[0]\n",
    "        ii_min, ii_max = (\n",
    "            _indices\n",
    "            if (\n",
    "                boundaries[_indices[0], 0, _index]\n",
    "                < boundaries[_indices[1], 0, _index]\n",
    "            )\n",
    "            else (_indices[1], _indices[0])\n",
    "        )\n",
    "        if source_point[_index] > boundaries[ii_max, 0, _index]:\n",
    "            insertion_edge_indices.append(ii_max)\n",
    "        elif source_point[_index] < boundaries[ii_min, 0, _index]:\n",
    "            insertion_edge_indices.append(ii_min)\n",
    "\n",
    "    return tuple(set(insertion_edge_indices))\n",
    "\n",
    "insertion_edge_indices = _determine_insertion_edge_indicies()\n",
    "insertion_edge_indices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6da6c7d-c8e8-49c1-9fb5-52e649eb55d7",
   "metadata": {},
   "source": [
    "Now we need to determine where every edge is located.  An edge is a boundary where the global mask switches between `True` and `False`.  To gather these edges we define a pool (array) that contains the (x, y) locations of every starting and ending point of an edge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aeb94e5-6ad8-4b68-bf55-fba1909b6019",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _add_to_edge_pool(edge, pool=None) -> Tuple[int, np.ndarray]:\n",
    "    # edge.shape == (2, 2)\n",
    "    # index_0 -> start (0) and end (1) point of edge\n",
    "    # index_1 -> edge coordinate (0, 1) = (x, y)\n",
    "    if pool is None:\n",
    "        pool = np.array(edge)[np.newaxis, ...]\n",
    "    else:\n",
    "        pool = np.concatenate(\n",
    "            (pool, np.array(edge)[np.newaxis, ...]),\n",
    "            axis=0,\n",
    "        )\n",
    "\n",
    "    return pool.shape[0] - 1, pool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c97de5f-bff5-423e-8f27-5bad0500ffff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _build_edge_pool(mask: xr.DataArray) -> np.ndarray:\n",
    "    # Find the (x, y) coordinates for the starting and ending points\n",
    "    # of an edge in the mask array.  An edge occurs then neighboring\n",
    "    # cells change values (i.e. switch between True and False)\n",
    "    pool = None\n",
    "\n",
    "    # gather vertical edges\n",
    "    edge_indices = np.where(np.diff(mask, axis=0))\n",
    "    ix_array = np.unique(edge_indices[0])\n",
    "\n",
    "    for ix in ix_array:\n",
    "        iy_array = edge_indices[1][edge_indices[0] == ix]\n",
    "\n",
    "        x = x_coord[ix] + 0.5 * dx\n",
    "\n",
    "        if iy_array.size == 1:\n",
    "            iy = iy_array[0]\n",
    "\n",
    "            edge = np.array(\n",
    "                [\n",
    "                    [x, y_coord[iy] - 0.5 * dy],\n",
    "                    [x, y_coord[iy] + 0.5 * dy],\n",
    "                ]\n",
    "            )\n",
    "            eid, pool = _add_to_edge_pool(edge, pool)\n",
    "        else:\n",
    "            jumps = np.where(np.diff(iy_array) != 1)[0]\n",
    "\n",
    "            starts = np.array([0])\n",
    "            starts = np.concatenate((starts, jumps + 1))\n",
    "            starts = iy_array[starts]\n",
    "\n",
    "            stops = np.concatenate((jumps, [iy_array.size - 1]))\n",
    "            stops = iy_array[stops]\n",
    "\n",
    "            for iy_start, iy_stop in zip(starts, stops):\n",
    "                edge = np.array(\n",
    "                    [\n",
    "                        [x, y_coord[iy_start] - 0.5 * dy],\n",
    "                        [x, y_coord[iy_stop] + 0.5 * dy],\n",
    "                    ]\n",
    "                )\n",
    "                eid, pool = _add_to_edge_pool(edge, pool)\n",
    "\n",
    "    # gather horizontal edges\n",
    "    edge_indices = np.where(np.diff(mask, axis=1))\n",
    "    iy_array = np.unique(edge_indices[1])\n",
    "\n",
    "    for iy in iy_array:\n",
    "        ix_array = edge_indices[0][edge_indices[1] == iy]\n",
    "\n",
    "        y = y_coord[iy] + 0.5 * dy\n",
    "\n",
    "        if ix_array.size == 1:\n",
    "            ix = ix_array[0]\n",
    "\n",
    "            edge = np.array(\n",
    "                [\n",
    "                    [x_coord[ix] - 0.5 * dx, y],\n",
    "                    [x_coord[ix] + 0.5 * dx, y],\n",
    "                ]\n",
    "            )\n",
    "            eid, pool = _add_to_edge_pool(edge, pool)\n",
    "        else:\n",
    "            jumps = np.where(np.diff(ix_array) != 1)[0]\n",
    "\n",
    "            starts = np.array([0])\n",
    "            starts = np.concatenate((starts, jumps + 1))\n",
    "            starts = ix_array[starts]\n",
    "\n",
    "            stops = np.concatenate((jumps, [ix_array.size - 1]))\n",
    "            stops = ix_array[stops]\n",
    "\n",
    "            for ix_start, ix_stop in zip(starts, stops):\n",
    "                edge = np.array(\n",
    "                    [\n",
    "                        [x_coord[ix_start] - 0.5 * dx, y],\n",
    "                        [x_coord[ix_stop] + 0.5 * dx, y],\n",
    "                    ]\n",
    "                )\n",
    "                eid, pool = _add_to_edge_pool(edge, pool)\n",
    "\n",
    "    # gather motion space perimeter edges\n",
    "    for ii in range(4):\n",
    "        boundary_side = boundaries[ii, ...]\n",
    "        delta = boundary_side[1, ...] - boundary_side[0, ...]\n",
    "        edge_type = \"horizontal\" if np.isclose(delta[1], 0) else \"vertical\"\n",
    "\n",
    "        if edge_type == \"horizontal\":\n",
    "            edge_vals = mask.sel(**{ykey: boundary_side[0, 1], \"method\": \"nearest\"})\n",
    "        else:\n",
    "            edge_vals = mask.sel(**{xkey: boundary_side[0, 0], \"method\": \"nearest\"})\n",
    "\n",
    "        compare_val = ii in insertion_edge_indices\n",
    "        _conditional_array = edge_vals if compare_val else np.logical_not(edge_vals)\n",
    "        if np.all(_conditional_array):\n",
    "            # perimeter side is not considered an \"edge\" (i.e. a boundary\n",
    "            # where True-False state switches\n",
    "            pass\n",
    "        elif np.all(np.logical_not(_conditional_array)):\n",
    "            # whole side is an edge\n",
    "            eid, pool = _add_to_edge_pool(boundary_side, pool)\n",
    "        else:\n",
    "            # array contain edges and non-edges\n",
    "            # False entries are edges\n",
    "            new_edge_indices = np.where(np.diff(_conditional_array))[0] + 1\n",
    "            if not _conditional_array[0]:\n",
    "                # boundary side starts as a new edge ... this is not captured\n",
    "                # by np.diff so manually add the first index\n",
    "                new_edge_indices = np.insert(new_edge_indices, 0, 0)\n",
    "            if not _conditional_array[-1]:\n",
    "                # boundary side ends as a new edge ... this is not captured\n",
    "                # by np.diff so manually add the last index\n",
    "                new_edge_indices = np.append(\n",
    "                    new_edge_indices, _conditional_array.size - 1\n",
    "                )\n",
    "\n",
    "            for jj in range(0, new_edge_indices.size, 2):\n",
    "                istart = new_edge_indices[jj]\n",
    "                istop = new_edge_indices[jj + 1] - 1\n",
    "\n",
    "                if edge_type == \"horizontal\":\n",
    "                    new_edge = np.array(\n",
    "                        [\n",
    "                            [x_coord[istart] - 0.5 * dx, boundary_side[0, 1]],\n",
    "                            [x_coord[istop] + 0.5 * dx, boundary_side[0, 1]],\n",
    "                        ],\n",
    "                    )\n",
    "                else:\n",
    "                    new_edge = np.array(\n",
    "                        [\n",
    "                            [boundary_side[0, 0], y_coord[istart] - 0.5 * dy],\n",
    "                            [boundary_side[0, 0], y_coord[istop] + 0.5 * dy],\n",
    "                        ],\n",
    "                    )\n",
    "\n",
    "                eid, pool = _add_to_edge_pool(new_edge, pool)\n",
    "\n",
    "    return pool\n",
    "\n",
    "edge_pool = _build_edge_pool(ds.mask)\n",
    "edge_pool.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ad7c7d7-cf09-42d0-8d73-0aa77bb559f4",
   "metadata": {},
   "source": [
    "The edge pool has a total of 99 identified edges.  This includes all edges of the currently defined exclusion layers, as well as the boundary edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b02deee5-481e-4934-9457-e8cfd96fd062",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.mask.plot(x=\"x\", y=\"y\", aspect=1.25, size=12)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_ylim(-60, 60)\n",
    "ax.set_xlim(-60, 60)\n",
    "\n",
    "plt.scatter(\n",
    "    edge_pool[..., 0, 0],\n",
    "    edge_pool[..., 0, 1],\n",
    "    color=\"blue\",\n",
    "    s=4**2,\n",
    ")\n",
    "plt.scatter(\n",
    "    edge_pool[..., 1, 0],\n",
    "    edge_pool[..., 1, 1],\n",
    "    facecolors=\"none\",\n",
    "    edgecolors=\"red\",\n",
    "    s=6**2,\n",
    ")\n",
    "plt.scatter(\n",
    "    source_point[0], \n",
    "    source_point[1],\n",
    "    marker=\"+\",\n",
    "    color=\"black\",\n",
    "    s=10**2,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8b9ad3f-d357-4f21-95c8-3c0b9e66891f",
   "metadata": {},
   "source": [
    "A lot of these edges are redundant when it comes calculating the shadow.  That is, an edge behind another edge does nothing for determinging the overall shadow.  Thus we need to filter out all the edge points that are behind other edges.  To do this we:\n",
    "\n",
    "1. Calculate all the vectors (`corner_rays`) that go between the source point and each of the defining edge points.\n",
    "2. Calculate the vectors that define an edge (`edge_vectors`).  This is just the difference between the end and start point of an edge.\n",
    "3. Determine if a corner ray intersects and closer edge vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd508450-794e-44ba-97b2-651151dcd76e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _build_corner_rays(edge_pool: np.ndarray) -> np.ndarray:\n",
    "    # collect unique edge points (i.e. unique (x,y) coords of edge\n",
    "    # segment start and stop locations)\n",
    "    edge_points = edge_pool.reshape(-1, 2)\n",
    "    edge_points = np.unique(edge_points, axis=0)\n",
    "\n",
    "    corner_rays = edge_points - source_point\n",
    "\n",
    "    # sort corner_rays and edge_points corresponding to the ray angle\n",
    "    delta = edge_points - source_point[np.newaxis, :]\n",
    "    perp_indices = np.where(delta[..., 0] == 0)[0]\n",
    "    if perp_indices.size > 0:\n",
    "        delta[perp_indices, 0] = 1  # dx\n",
    "        delta[perp_indices, 1] = np.inf * (\n",
    "                delta[perp_indices, 1] / np.abs(delta[perp_indices, 1])\n",
    "        )  # dy\n",
    "    ray_angles = np.arctan(delta[..., 1] / delta[..., 0])\n",
    "    sort_i = np.argsort(ray_angles)\n",
    "    corner_rays = corner_rays[sort_i]\n",
    "\n",
    "    return corner_rays\n",
    "\n",
    "corner_rays = _build_corner_rays(edge_pool)\n",
    "corner_rays.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c257e19-2aaa-4c4a-aec1-0eeb90effa2b",
   "metadata": {},
   "source": [
    "In order to determine if a corner ray crosses and edge lets first two lines defined as:\n",
    "\n",
    "$$\n",
    "\\vec{l}_1 = \\vec{p}_s + \\mu * \\vec{v}_{ray}\\\\\n",
    "\\vec{l}_2 = \\vec{p}_e + \\nu * \\vec{v}_{edge}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\vec{p}_s = \\text{vector to the source point}\\\\\n",
    "\\vec{p}_e = \\text{vector to the edge start point}\\\\\n",
    "\\vec{v}_{ray} = \\text{the corner ray vector}\\\\\n",
    "\\vec{v}_{edge} = \\text{the edge vector}\\\\\n",
    "\\mu, \\nu\\enspace \\text{are real numbers}\n",
    "$$\n",
    "\n",
    "The two lines cross each other when $\\vec{l}_1 = \\vec{l}_2$, which allows us to solve for $\\mu$ and $\\nu$.\n",
    "\n",
    "$$\n",
    "\\mu = \\frac{(\\vec{p}_e - \\vec{p}_s)\\times\\vec{v}_{edge}}{\\vec{v}_{ray}\\times\\vec{v}_{edge}}\\\\\n",
    "\\nu = \\frac{(\\vec{p}_e - \\vec{p}_s)\\times\\vec{v}_{ray}}{\\vec{v}_{ray}\\times\\vec{v}_{edge}}\n",
    "$$\n",
    "\n",
    "If $0 \\le \\nu \\le 1$ then the corner ray crosses the edge.  If $0 \\le \\mu < 1$, then that crossed edge is closer than the edge that defined the ray.\n",
    "\n",
    "Now lets filter out all the `corner_rays` that point to edges behind other edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07e254b3-305b-4dcc-a633-1b291b0c336c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _build_corner_ray_mask(\n",
    "    corner_rays: np.ndarray, edge_pool: np.ndarray\n",
    ") -> np.ndarray:\n",
    "\n",
    "    edge_vectors = edge_pool[..., 1, :] - edge_pool[..., 0, :]\n",
    "\n",
    "    # determine if a corner_ray intersects an edge that is closer\n",
    "    # to the insertion point\n",
    "    # - solving the eqn:\n",
    "    #\n",
    "    #   source_point + mu * corner_ray = edge_pool[..., 0, :] + nu * edge_vector\n",
    "    #\n",
    "    #   * mu and nu are scalars\n",
    "    #   * if 0 < mu < 1 and 0 < nu < 1, then the corner_ray passes through a\n",
    "    #     closer edge to the insertion point\n",
    "    #\n",
    "    point_deltas = edge_pool[..., 0, :] - source_point\n",
    "    denominator = np.cross(corner_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(0, 1)\n",
    "    mu_array = np.cross(point_deltas, edge_vectors) / denominator\n",
    "    nu_array = (\n",
    "        np.cross(point_deltas[:, np.newaxis, ...], corner_rays ).swapaxes(0, 1)\n",
    "        / denominator\n",
    "    )\n",
    "    mu_condition = np.logical_and(mu_array >= 0, mu_array < 1)\n",
    "    nu_condition = np.logical_and(nu_array >= 0, nu_array <= 1)\n",
    "    _condition = np.logical_and(mu_condition, nu_condition)\n",
    "    _count = np.count_nonzero(_condition, axis=1)\n",
    "\n",
    "    return np.where(_count == 0, True, False)\n",
    "\n",
    "ray_mask = _build_corner_ray_mask(corner_rays, edge_pool)\n",
    "corner_rays = corner_rays[ray_mask]\n",
    "corner_rays.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5438c505-c051-413b-b518-6c9d721162f2",
   "metadata": {},
   "source": [
    "We have now reduced the number of corner rays by 60%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8e66621-065b-4e5d-8a4d-13eecaa207cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.mask.plot(x=\"x\", y=\"y\", aspect=1.25, size=12)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_ylim(-60, 60)\n",
    "ax.set_xlim(-60, 60)\n",
    "\n",
    "plt.scatter(\n",
    "    corner_rays[..., 0] + source_point[0],\n",
    "    corner_rays[..., 1] + source_point[1],\n",
    "    facecolors=\"none\",\n",
    "    edgecolors=\"red\",\n",
    "    s=6**2,\n",
    ")\n",
    "plt.scatter(\n",
    "    source_point[0], \n",
    "    source_point[1],\n",
    "    marker=\"+\",\n",
    "    color=\"black\",\n",
    "    s=10**2,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "577b2215-1140-425f-9d01-1997b5df8b46",
   "metadata": {},
   "source": [
    "Now we need to fan each corner ray to see if one of those fanned rays project to a boundary or exclusion layer in the \"background\".  If it does, then that fanned ray needs to be added to the collection of corner rays.  This collection of rays will eventually allows us to paint the motion space `True` in the areas \"lit\" by the point source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bd7ada2-6a6b-4717-bedc-89c4c32c308c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _build_fanned_rays(\n",
    "    edge_pool: np.ndarray, corner_rays: np.ndarray\n",
    ") -> np.ndarray:\n",
    "\n",
    "    # calculate angles of corner_rays\n",
    "    angles = np.arcsin(corner_rays[..., 1] / np.linalg.norm(corner_rays, axis=1))\n",
    "    corner_ray_angles = np.where(corner_rays[..., 0] >= 0, angles, np.pi - angles)\n",
    "\n",
    "    # generate fanned rays\n",
    "    delta_angle = (\n",
    "        0.01 * np.min([dx, dy])\n",
    "        / np.linalg.norm(corner_rays, axis=1)\n",
    "    )\n",
    "\n",
    "    fan_plus = np.array(\n",
    "        [\n",
    "            np.cos(corner_ray_angles + delta_angle),\n",
    "            np.sin(corner_ray_angles + delta_angle),\n",
    "        ]\n",
    "    ).swapaxes(0, 1)\n",
    "\n",
    "    fan_minus = np.array(\n",
    "        [\n",
    "            np.cos(corner_ray_angles - delta_angle),\n",
    "            np.sin(corner_ray_angles - delta_angle),\n",
    "        ]\n",
    "    ).swapaxes(0, 1)\n",
    "    print(f\"fan_minus.shape = {fan_minus.shape}\")\n",
    "\n",
    "    fan_rays = np.concatenate((fan_plus, fan_minus), axis=0)\n",
    "\n",
    "    # project fan rays to the nearest edge\n",
    "    edge_vectors = edge_pool[..., 1, :] - edge_pool[..., 0, :]\n",
    "\n",
    "    point_deltas = edge_pool[..., 0, :] - source_point\n",
    "    denominator = np.cross(fan_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(0, 1)\n",
    "    mu_array = np.cross(point_deltas, edge_vectors) / denominator\n",
    "    nu_array = (\n",
    "        np.cross(point_deltas[:, np.newaxis, ...], fan_rays).swapaxes(0, 1)\n",
    "        / denominator\n",
    "    )\n",
    "\n",
    "    mu_condition = mu_array > 0\n",
    "    nu_condition = np.logical_and(nu_array >= 0, nu_array <= 1)\n",
    "    mask = np.logical_and(mu_condition, nu_condition)\n",
    "\n",
    "    # Note: rays that do not satisfy the mask conditions project to\n",
    "    #       infinity.  This happens when the insertion point is not\n",
    "    #       located in the motion space and a boundary corner is\n",
    "    #       fanned.\n",
    "    mu_array[np.logical_not(mask)] = np.inf\n",
    "    adjusted_mu_array = np.nanmin(mu_array, axis=1)\n",
    "    fan_rays = adjusted_mu_array[..., None] * fan_rays\n",
    "\n",
    "    # filter close points before merging\n",
    "    double_corner_rays = np.concatenate((corner_rays, corner_rays), axis=0)\n",
    "    mask = np.logical_and(\n",
    "        np.isclose(fan_rays[..., 0], double_corner_rays[..., 0], atol=.5 * dx),\n",
    "        np.isclose(fan_rays[..., 1], double_corner_rays[..., 1], atol=.5 * dy),\n",
    "    )\n",
    "    fan_rays = fan_rays[np.logical_not(mask)]\n",
    "    \n",
    "    # filter out np.inf rays\n",
    "    finite_mask = np.all(np.isfinite(fan_rays), axis=1)\n",
    "    fan_rays = fan_rays[finite_mask]\n",
    "    \n",
    "    # sort fan rays\n",
    "    ray_angles = np.arcsin(fan_rays[..., 1] / np.linalg.norm(fan_rays, axis=1))\n",
    "    ray_angles = np.where(fan_rays[..., 0] >= 0, ray_angles, np.pi - ray_angles)\n",
    "    sort_i = np.argsort(ray_angles)\n",
    "    fan_rays = fan_rays[sort_i]\n",
    "\n",
    "    return fan_rays\n",
    "\n",
    "fan_rays = _build_fanned_rays(edge_pool, corner_rays)\n",
    "fan_rays.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae5e09b8-6ebd-4dac-bb86-3d4b49e1a99a",
   "metadata": {},
   "source": [
    "Merge the corner and fanned rays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55f40392-578a-48f9-81d4-f2e722471791",
   "metadata": {},
   "outputs": [],
   "source": [
    "rays = np.concatenate((corner_rays, fan_rays), axis=0)\n",
    "\n",
    "# sort rays by angle\n",
    "angles = np.arcsin(rays[..., 1] / np.linalg.norm(rays, axis=1))\n",
    "angles = np.where(rays[..., 0] >= 0, angles, np.pi - angles)\n",
    "sort_i = np.argsort(angles)\n",
    "rays = rays[sort_i]\n",
    "rays.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5cc817d-b4fb-44b5-a9d7-f6ac486a9208",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.mask.plot(x=\"x\", y=\"y\", aspect=1.25, size=12)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_ylim(-60, 60)\n",
    "ax.set_xlim(-60, 60)\n",
    "\n",
    "plt.scatter(\n",
    "    rays[..., 0] + source_point[0],\n",
    "    rays[..., 1] + source_point[1],\n",
    "    facecolors=\"none\",\n",
    "    edgecolors=\"red\",\n",
    "    s=6**2,\n",
    ")\n",
    "plt.scatter(\n",
    "    source_point[0], \n",
    "    source_point[1],\n",
    "    marker=\"+\",\n",
    "    color=\"black\",\n",
    "    s=10**2,\n",
    ")\n",
    "for ii in range(rays.shape[0]):\n",
    "    plt.plot(\n",
    "        [source_point[0], rays[ii, 0] + source_point[0]],\n",
    "        [source_point[1], rays[ii, 1] + source_point[1]],\n",
    "        color=\"red\",\n",
    "        linewidth=0.5,\n",
    "    );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0417fb7a-3554-43f9-8f56-4f05249113d2",
   "metadata": {},
   "source": [
    "[Barycentric Coordinates]: https://en.wikipedia.org/wiki/Barycentric_coordinate_system\n",
    "\n",
    "The last step is to paint the motion space mask true using all the triangles formed by the `rays` and `source_point`.  How is this done?\n",
    "\n",
    "1. Using the `source_point` and two neighboring rays in `rays` we can iterate through `rays` to form and around of triangles.\n",
    "2. Using [Barycentric Coordinates] we can say a point `\\vec{p}` can be written as\n",
    "\n",
    "   $$\n",
    "   \\vec{p} = \\lambda_1 \\vec{p}_1 + \\lambda_2 \\vec{p}_2 + \\lambda_3 \\vec{p}_3\n",
    "   $$\n",
    "\n",
    "    where $(\\vec{p}_1, \\vec{p}_2, \\vec{p}_3)$ are the three points that define a triangle and $(\\lambda_1, \\lambda_2, \\lambda_3)$ are real scalars.\n",
    "3. If all $0 \\le \\lambda_i \\le 1$ and $1=\\lambda_1 + \\lambda_2 +\\lambda_3$, then we can say point $\\vec{p}$ is inside the triangle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a73b959b-eb9e-40c1-b859-e8dd24e6f75e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _paint_mask(mask, rays: np.ndarray) -> xr.DataArray:\n",
    "    rays = np.append(rays, rays[0, ...][None, ...], axis=0)\n",
    "    endpoints = rays + source_point[None, :]\n",
    "\n",
    "    triangles = np.zeros((rays.shape[0] - 1, 3, 2))\n",
    "    triangles[..., 0, :] = source_point\n",
    "    triangles[..., 1, :] = endpoints[:-1, :]\n",
    "    triangles[..., 2, :] = endpoints[1:, :]\n",
    "\n",
    "    grid_points = np.zeros((x_coord.size, y_coord.size, 2))\n",
    "    grid_points[..., 0] = np.repeat(\n",
    "        x_coord.values[..., np.newaxis], y_coord.size, axis=1\n",
    "    )\n",
    "    grid_points[..., 1] = np.repeat(\n",
    "        y_coord.values[np.newaxis, ...], x_coord.size, axis=0\n",
    "    )\n",
    "\n",
    "    # This processes uses Barycentric coordinates to determine if a\n",
    "    # grid point is within the triangle.\n",
    "    #\n",
    "    # https://en.wikipedia.org/wiki/Barycentric_coordinate_system\n",
    "    #\n",
    "    # lambda shape is (x_size, y_size, N_rays)\n",
    "    #\n",
    "    # calculate lambda_3\n",
    "    numerator = np.cross(\n",
    "        grid_points[:, :, None, :] - triangles[None, None, :, 0, :],\n",
    "        (triangles[:, 1, :] - triangles[:, 0, :])[None, None, :, :],\n",
    "    )\n",
    "    denominator = np.cross(\n",
    "        triangles[:, 2, :] - triangles[:, 0, :],\n",
    "        triangles[:, 1, :] - triangles[:, 0, :]\n",
    "    )\n",
    "    \n",
    "    zero_mask = denominator == 0\n",
    "    if np.any(zero_mask):\n",
    "        # denominator can be zero if all points on the triangle lie\n",
    "        # on a line\n",
    "        not_zero_mask = np.logical_not(zero_mask)\n",
    "        triangles = triangles[not_zero_mask, ...]\n",
    "        numerator = numerator[..., not_zero_mask]\n",
    "        denominator = denominator[not_zero_mask]\n",
    "\n",
    "    lambda_3 = numerator / denominator[None, None, ...]\n",
    "\n",
    "    # calculate lambda_2\n",
    "    numerator = np.cross(\n",
    "        grid_points[:, :, None, :] - triangles[None, None, :, 0, :],\n",
    "        (triangles[:, 2, :] - triangles[:, 0, :])[None, None, :, :],\n",
    "    )\n",
    "    denominator = -denominator\n",
    "    lambda_2 = numerator / denominator[None, None, ...]\n",
    "\n",
    "    # calculate lambda_1\n",
    "    lambda_1 = 1 - lambda_2 - lambda_3\n",
    "\n",
    "    # generate the conditional for each point in the motion space\n",
    "    # _conditional.shape = (mspace Nx, mspace Ny, N_rays)\n",
    "    #\n",
    "    lambda_1_condition = np.logical_and(lambda_1 >= 0, lambda_1 <= 1)\n",
    "    lambda_2_condition = np.logical_and(lambda_2 >= 0, lambda_2 <= 1)\n",
    "    lambda_3_condition = np.logical_and(lambda_3 >= 0, lambda_3 <= 1)\n",
    "    _condition = np.logical_and(\n",
    "        np.logical_and(lambda_1_condition, lambda_2_condition),\n",
    "        lambda_3_condition,\n",
    "    )\n",
    "\n",
    "    return mask.copy(data=np.any(_condition, axis=2))\n",
    "\n",
    "shadow_mask = _paint_mask(ds.mask, rays)\n",
    "shadow_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70cb91cf-92a1-4cda-91a8-00b970fea169",
   "metadata": {},
   "outputs": [],
   "source": [
    "shadow_mask.plot(x=\"x\", y=\"y\", aspect=1.25, size=12)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_ylim(-60, 60)\n",
    "ax.set_xlim(-60, 60)\n",
    "\n",
    "plt.scatter(\n",
    "    rays[..., 0] + source_point[0],\n",
    "    rays[..., 1] + source_point[1],\n",
    "    facecolors=\"none\",\n",
    "    edgecolors=\"red\",\n",
    "    s=6**2,\n",
    ")\n",
    "plt.scatter(\n",
    "    source_point[0], \n",
    "    source_point[1],\n",
    "    marker=\"+\",\n",
    "    color=\"black\",\n",
    "    s=10**2,\n",
    ")\n",
    "for ii in range(rays.shape[0]):\n",
    "    plt.plot(\n",
    "        [source_point[0], rays[ii, 0] + source_point[0]],\n",
    "        [source_point[1], rays[ii, 1] + source_point[1]],\n",
    "        color=\"red\",\n",
    "        linewidth=0.5,\n",
    "    );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe8bd710-a43f-4cc2-af50-4503dec413f6",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}