This page was generated by nbsphinx from docs/notebooks/motion_list/Shadow2DExclusion.ipynb.
Shadow2DExclusion Overview
[1]:
%matplotlib inline
[2]:
import numpy as np
import matplotlib.pyplot as plt
import sys
import xarray as xr
from matplotlib.patches import Polygon
from typing import Tuple
sys.executable;
[3]:
try:
from bapsf_motion.motion_builder.exclusions import (
CircularExclusion,
DividerExclusion,
GovernExclusion,
Shadow2DExclusion,
)
except ModuleNotFoundError:
from pathlib import Path
HERE = Path().cwd()
BAPSF_MOTION = (HERE / ".." / ".." / ".." ).resolve()
sys.path.append(str(BAPSF_MOTION))
from bapsf_motion.motion_builder.exclusions import (
CircularExclusion,
DividerExclusion,
GovernExclusion,
Shadow2DExclusion,
)
[4]:
plt.rcParams.update(
{
# "figure.figsize": [12, 0.56 * 12],
"figure.figsize": [10, 0.8 * 10],
"font.size": 16,
}
)
Usage
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.
Is a GovernExclusion
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.
[5]:
issubclass(Shadow2DExclusion, GovernExclusion)
[5]:
True
Defining a Shadow2DExclusion
Create an initial global mask.
[6]:
size = 111
side = np.linspace(-55, 55, num=size)
ds = xr.Dataset(
{"mask": (("x", "y"), np.ones((size, size), dtype=bool))},
coords={
"x": side,
"y": side,
},
)
Add some initial exlustion layers for Shadow2DExclusion to act on.
[7]:
ex_1 = CircularExclusion(ds, radius=10, center=(0,0), exclude="inside")
ex_2 = CircularExclusion(ds, radius=8, center=(25,30), exclude="inside")
ds.mask.plot(x="x", y="y");
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.
[8]:
ex = Shadow2DExclusion(ds, source_point=[45,-58])
ds.mask.plot(x="x", y="y");
plt.scatter(
ex.source_point[0],
ex.source_point[1],
marker="+",
color="black",
s=8**2,
)
axis = plt.gca()
axis.set_xlim(-60, 60)
axis.set_ylim(-60, 60)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[8], line 1
----> 1 ex = Shadow2DExclusion(ds, source_point=[45,-58])
2
3 ds.mask.plot(x="x", y="y");
4 plt.scatter(
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/bapsf_motion/motion_builder/exclusions/shadow.py:97, in Shadow2DExclusion.__init__(self, ds, source_point, skip_ds_add)
94 self._boundaries = None
95 self._insertion_edge_indices = None
---> 97 super().__init__(
98 ds,
99 source_point=source_point,
100 skip_ds_add=skip_ds_add,
101 )
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/bapsf_motion/motion_builder/exclusions/base.py:70, in BaseExclusion.__init__(self, ds, skip_ds_add, **kwargs)
67 return
69 # store this mask to the Dataset
---> 70 self.regenerate_exclusion()
72 # update the global mask
73 self.update_global_mask()
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/bapsf_motion/motion_builder/exclusions/base.py:212, in BaseExclusion.regenerate_exclusion(self)
204 raise RuntimeError(
205 f"For exclusion {self.name} skip_ds_add={self.skip_ds_add} and thus "
206 f"the exclusion can not be regenerated and updated in the Dataset. "
207 f"To get the exclusion matrix use the 'ex.exclusion' property."
208 )
210 self.composed_exclusions.clear()
--> 212 self._ds[self.name] = self._generate_exclusion()
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/bapsf_motion/motion_builder/exclusions/shadow.py:190, in Shadow2DExclusion._generate_exclusion(self)
186 corner_rays = self._build_corner_rays(edge_pool)
188 # Reduce corner rays
189 # - remove corner arrays that point to points behind nearer edges
--> 190 ray_mask = self._build_corner_ray_mask(corner_rays, edge_pool)
191 corner_rays = corner_rays[ray_mask]
193 # Build an array of fanned rays from all the corner rays
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/bapsf_motion/motion_builder/exclusions/shadow.py:314, in Shadow2DExclusion._build_corner_ray_mask(self, corner_rays, edge_pool)
303 # determine if a corner_ray intersects an edge that is closer
304 # to the insertion point
305 # - solving the eqn:
(...) 311 # closer edge to the insertion point
312 #
313 point_deltas = edge_pool[..., 0, :] - self.source_point
--> 314 denominator = np.cross(corner_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(
315 0, 1
316 )
317 mu_array = np.cross(point_deltas, edge_vectors) / denominator
318 nu_array = (
319 np.cross(point_deltas[:, np.newaxis, ...], corner_rays).swapaxes(0, 1)
320 / denominator
321 )
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/numpy/_core/numeric.py:1679, in cross(a, b, axisa, axisb, axisc, axis)
1677 b = moveaxis(b, axisb, -1)
1678 if a.shape[-1] != 3 or b.shape[-1] != 3:
-> 1679 raise ValueError(
1680 f"Both input arrays must be (arrays of) 3-dimensional vectors, "
1681 f"but they are {a.shape[-1]} and {b.shape[-1]} dimensional instead."
1682 )
1684 # Create the output array
1685 shape = *broadcast(a[..., 0], b[..., 0]).shape, 3
ValueError: Both input arrays must be (arrays of) 3-dimensional vectors, but they are 2 and 2 dimensional instead.
The configuration dictionary for this exclusion looks like…
[9]:
ex.config
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 ex.config
NameError: name 'ex' is not defined
The Underlying Algorithm
Lets first reset the initial mask.
[10]:
size = 111
side = np.linspace(-55, 55, num=size)
ds = xr.Dataset(
{"mask": (("x", "y"), np.ones((size, size), dtype=bool))},
coords={
"x": side,
"y": side,
},
)
ex_1 = CircularExclusion(ds, radius=10, center=(0,0), exclude="inside")
ex_2 = CircularExclusion(ds, radius=8, center=(25,30), exclude="inside")
ds.mask.plot(x="x", y="y");
Lets first define some key parameters
[11]:
# gather some mask parameters
dx, dy = ex_1.mask_resolution
xkey, ykey = ex_1.mspace_dims
x_coord, y_coord = ex_1.mspace_coords[xkey], ex_1.mspace_coords[ykey]
# point source
source_point = np.array([45, -58])
[12]:
# define boundary of the motion space
# - index_0 = the N-th side
# - index_1 = start (0) and stop (1) of the side
# - index_2 = (x, y) coordinate of the associated point
#
boundaries = np.zeros((4, 2, 2))
x_min = x_coord[0] - 0.5 * dx
x_max = x_coord[-1] + 0.5 * dx
y_min =y_coord[0] - 0.5 * dy
y_max = y_coord[-1] + 0.5 * dy
# lower horizontal
boundaries[0, 0, :] = [x_min, y_min]
boundaries[0, 1, :] = [x_max, y_min]
# right vertical
boundaries[1, 0, :] = [x_max, y_min]
boundaries[1, 1, :] = [x_max, y_max]
# upper horizontal
boundaries[2, 0, :] = [x_max, y_max]
boundaries[2, 1, :] = [x_min, y_max]
# left vertical
boundaries[3, 0, :] = [x_min, y_max]
boundaries[3, 1, :] = [x_min, y_min]
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:
Assume there are exclusion layers to shadow. If the motion space mask is all
Trueor allFalsethen the resulting mask is the same.The motion space outside a boundary side is
Trueif the side is a side the point source would shine through (outside-to-inside), otherwise it is consideredFalse.
[13]:
def _determine_insertion_edge_indicies():
x_range = [x_coord[0] - 0.5 * dx, x_coord[-1] + 0.5 * dx]
y_range = [y_coord[0] - 0.5 * dy, y_coord[-1] + 0.5 * dy]
if (
(x_range[0] <= source_point[0] <= x_range[1])
and (y_range[0] <= source_point[1] <= y_range[1])
):
# insertion point is within the motion space
return ()
insertion_edge_indices = []
deltas = boundaries[..., 1, :] - boundaries[..., 0, :]
for _orientation, _index in zip(["horizontal", "vertical"], [1, 0]):
_indices = np.where(np.isclose(deltas[..., _index], 0))[0]
ii_min, ii_max = (
_indices
if (
boundaries[_indices[0], 0, _index]
< boundaries[_indices[1], 0, _index]
)
else (_indices[1], _indices[0])
)
if source_point[_index] > boundaries[ii_max, 0, _index]:
insertion_edge_indices.append(ii_max)
elif source_point[_index] < boundaries[ii_min, 0, _index]:
insertion_edge_indices.append(ii_min)
return tuple(set(insertion_edge_indices))
insertion_edge_indices = _determine_insertion_edge_indicies()
insertion_edge_indices
[13]:
(np.int64(0),)
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.
[14]:
def _add_to_edge_pool(edge, pool=None) -> Tuple[int, np.ndarray]:
# edge.shape == (2, 2)
# index_0 -> start (0) and end (1) point of edge
# index_1 -> edge coordinate (0, 1) = (x, y)
if pool is None:
pool = np.array(edge)[np.newaxis, ...]
else:
pool = np.concatenate(
(pool, np.array(edge)[np.newaxis, ...]),
axis=0,
)
return pool.shape[0] - 1, pool
[15]:
def _build_edge_pool(mask: xr.DataArray) -> np.ndarray:
# Find the (x, y) coordinates for the starting and ending points
# of an edge in the mask array. An edge occurs then neighboring
# cells change values (i.e. switch between True and False)
pool = None
# gather vertical edges
edge_indices = np.where(np.diff(mask, axis=0))
ix_array = np.unique(edge_indices[0])
for ix in ix_array:
iy_array = edge_indices[1][edge_indices[0] == ix]
x = x_coord[ix] + 0.5 * dx
if iy_array.size == 1:
iy = iy_array[0]
edge = np.array(
[
[x, y_coord[iy] - 0.5 * dy],
[x, y_coord[iy] + 0.5 * dy],
]
)
eid, pool = _add_to_edge_pool(edge, pool)
else:
jumps = np.where(np.diff(iy_array) != 1)[0]
starts = np.array([0])
starts = np.concatenate((starts, jumps + 1))
starts = iy_array[starts]
stops = np.concatenate((jumps, [iy_array.size - 1]))
stops = iy_array[stops]
for iy_start, iy_stop in zip(starts, stops):
edge = np.array(
[
[x, y_coord[iy_start] - 0.5 * dy],
[x, y_coord[iy_stop] + 0.5 * dy],
]
)
eid, pool = _add_to_edge_pool(edge, pool)
# gather horizontal edges
edge_indices = np.where(np.diff(mask, axis=1))
iy_array = np.unique(edge_indices[1])
for iy in iy_array:
ix_array = edge_indices[0][edge_indices[1] == iy]
y = y_coord[iy] + 0.5 * dy
if ix_array.size == 1:
ix = ix_array[0]
edge = np.array(
[
[x_coord[ix] - 0.5 * dx, y],
[x_coord[ix] + 0.5 * dx, y],
]
)
eid, pool = _add_to_edge_pool(edge, pool)
else:
jumps = np.where(np.diff(ix_array) != 1)[0]
starts = np.array([0])
starts = np.concatenate((starts, jumps + 1))
starts = ix_array[starts]
stops = np.concatenate((jumps, [ix_array.size - 1]))
stops = ix_array[stops]
for ix_start, ix_stop in zip(starts, stops):
edge = np.array(
[
[x_coord[ix_start] - 0.5 * dx, y],
[x_coord[ix_stop] + 0.5 * dx, y],
]
)
eid, pool = _add_to_edge_pool(edge, pool)
# gather motion space perimeter edges
for ii in range(4):
boundary_side = boundaries[ii, ...]
delta = boundary_side[1, ...] - boundary_side[0, ...]
edge_type = "horizontal" if np.isclose(delta[1], 0) else "vertical"
if edge_type == "horizontal":
edge_vals = mask.sel(**{ykey: boundary_side[0, 1], "method": "nearest"})
else:
edge_vals = mask.sel(**{xkey: boundary_side[0, 0], "method": "nearest"})
compare_val = ii in insertion_edge_indices
_conditional_array = edge_vals if compare_val else np.logical_not(edge_vals)
if np.all(_conditional_array):
# perimeter side is not considered an "edge" (i.e. a boundary
# where True-False state switches
pass
elif np.all(np.logical_not(_conditional_array)):
# whole side is an edge
eid, pool = _add_to_edge_pool(boundary_side, pool)
else:
# array contain edges and non-edges
# False entries are edges
new_edge_indices = np.where(np.diff(_conditional_array))[0] + 1
if not _conditional_array[0]:
# boundary side starts as a new edge ... this is not captured
# by np.diff so manually add the first index
new_edge_indices = np.insert(new_edge_indices, 0, 0)
if not _conditional_array[-1]:
# boundary side ends as a new edge ... this is not captured
# by np.diff so manually add the last index
new_edge_indices = np.append(
new_edge_indices, _conditional_array.size - 1
)
for jj in range(0, new_edge_indices.size, 2):
istart = new_edge_indices[jj]
istop = new_edge_indices[jj + 1] - 1
if edge_type == "horizontal":
new_edge = np.array(
[
[x_coord[istart] - 0.5 * dx, boundary_side[0, 1]],
[x_coord[istop] + 0.5 * dx, boundary_side[0, 1]],
],
)
else:
new_edge = np.array(
[
[boundary_side[0, 0], y_coord[istart] - 0.5 * dy],
[boundary_side[0, 0], y_coord[istop] + 0.5 * dy],
],
)
eid, pool = _add_to_edge_pool(new_edge, pool)
return pool
edge_pool = _build_edge_pool(ds.mask)
edge_pool.shape
[15]:
(99, 2, 2)
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.
[16]:
ds.mask.plot(x="x", y="y", aspect=1.25, size=12)
ax = plt.gca()
ax.set_ylim(-60, 60)
ax.set_xlim(-60, 60)
plt.scatter(
edge_pool[..., 0, 0],
edge_pool[..., 0, 1],
color="blue",
s=4**2,
)
plt.scatter(
edge_pool[..., 1, 0],
edge_pool[..., 1, 1],
facecolors="none",
edgecolors="red",
s=6**2,
)
plt.scatter(
source_point[0],
source_point[1],
marker="+",
color="black",
s=10**2,
);
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:
Calculate all the vectors (
corner_rays) that go between the source point and each of the defining edge points.Calculate the vectors that define an edge (
edge_vectors). This is just the difference between the end and start point of an edge.Determine if a corner ray intersects and closer edge vector.
[17]:
def _build_corner_rays(edge_pool: np.ndarray) -> np.ndarray:
# collect unique edge points (i.e. unique (x,y) coords of edge
# segment start and stop locations)
edge_points = edge_pool.reshape(-1, 2)
edge_points = np.unique(edge_points, axis=0)
corner_rays = edge_points - source_point
# sort corner_rays and edge_points corresponding to the ray angle
delta = edge_points - source_point[np.newaxis, :]
perp_indices = np.where(delta[..., 0] == 0)[0]
if perp_indices.size > 0:
delta[perp_indices, 0] = 1 # dx
delta[perp_indices, 1] = np.inf * (
delta[perp_indices, 1] / np.abs(delta[perp_indices, 1])
) # dy
ray_angles = np.arctan(delta[..., 1] / delta[..., 0])
sort_i = np.argsort(ray_angles)
corner_rays = corner_rays[sort_i]
return corner_rays
corner_rays = _build_corner_rays(edge_pool)
corner_rays.shape
[17]:
(100, 2)
In order to determine if a corner ray crosses and edge lets first two lines defined as:
\[\begin{split}\vec{l}_1 = \vec{p}_s + \mu * \vec{v}_{ray}\\
\vec{l}_2 = \vec{p}_e + \nu * \vec{v}_{edge}\end{split}\]
where
\[\begin{split}\vec{p}_s = \text{vector to the source point}\\
\vec{p}_e = \text{vector to the edge start point}\\
\vec{v}_{ray} = \text{the corner ray vector}\\
\vec{v}_{edge} = \text{the edge vector}\\
\mu, \nu\enspace \text{are real numbers}\end{split}\]
The two lines cross each other when \(\vec{l}_1 = \vec{l}_2\), which allows us to solve for \(\mu\) and \(\nu\).
\[\begin{split}\mu = \frac{(\vec{p}_e - \vec{p}_s)\times\vec{v}_{edge}}{\vec{v}_{ray}\times\vec{v}_{edge}}\\
\nu = \frac{(\vec{p}_e - \vec{p}_s)\times\vec{v}_{ray}}{\vec{v}_{ray}\times\vec{v}_{edge}}\end{split}\]
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.
Now lets filter out all the corner_rays that point to edges behind other edges.
[18]:
def _build_corner_ray_mask(
corner_rays: np.ndarray, edge_pool: np.ndarray
) -> np.ndarray:
edge_vectors = edge_pool[..., 1, :] - edge_pool[..., 0, :]
# determine if a corner_ray intersects an edge that is closer
# to the insertion point
# - solving the eqn:
#
# source_point + mu * corner_ray = edge_pool[..., 0, :] + nu * edge_vector
#
# * mu and nu are scalars
# * if 0 < mu < 1 and 0 < nu < 1, then the corner_ray passes through a
# closer edge to the insertion point
#
point_deltas = edge_pool[..., 0, :] - source_point
denominator = np.cross(corner_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(0, 1)
mu_array = np.cross(point_deltas, edge_vectors) / denominator
nu_array = (
np.cross(point_deltas[:, np.newaxis, ...], corner_rays ).swapaxes(0, 1)
/ denominator
)
mu_condition = np.logical_and(mu_array >= 0, mu_array < 1)
nu_condition = np.logical_and(nu_array >= 0, nu_array <= 1)
_condition = np.logical_and(mu_condition, nu_condition)
_count = np.count_nonzero(_condition, axis=1)
return np.where(_count == 0, True, False)
ray_mask = _build_corner_ray_mask(corner_rays, edge_pool)
corner_rays = corner_rays[ray_mask]
corner_rays.shape
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[18], line 31
27 _count = np.count_nonzero(_condition, axis=1)
28
29 return np.where(_count == 0, True, False)
30
---> 31 ray_mask = _build_corner_ray_mask(corner_rays, edge_pool)
32 corner_rays = corner_rays[ray_mask]
33 corner_rays.shape
Cell In[18], line 18, in _build_corner_ray_mask(corner_rays, edge_pool)
14 # * if 0 < mu < 1 and 0 < nu < 1, then the corner_ray passes through a
15 # closer edge to the insertion point
16 #
17 point_deltas = edge_pool[..., 0, :] - source_point
---> 18 denominator = np.cross(corner_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(0, 1)
19 mu_array = np.cross(point_deltas, edge_vectors) / denominator
20 nu_array = (
21 np.cross(point_deltas[:, np.newaxis, ...], corner_rays ).swapaxes(0, 1)
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/numpy/_core/numeric.py:1679, in cross(a, b, axisa, axisb, axisc, axis)
1677 b = moveaxis(b, axisb, -1)
1678 if a.shape[-1] != 3 or b.shape[-1] != 3:
-> 1679 raise ValueError(
1680 f"Both input arrays must be (arrays of) 3-dimensional vectors, "
1681 f"but they are {a.shape[-1]} and {b.shape[-1]} dimensional instead."
1682 )
1684 # Create the output array
1685 shape = *broadcast(a[..., 0], b[..., 0]).shape, 3
ValueError: Both input arrays must be (arrays of) 3-dimensional vectors, but they are 2 and 2 dimensional instead.
We have now reduced the number of corner rays by 60%.
[19]:
ds.mask.plot(x="x", y="y", aspect=1.25, size=12)
ax = plt.gca()
ax.set_ylim(-60, 60)
ax.set_xlim(-60, 60)
plt.scatter(
corner_rays[..., 0] + source_point[0],
corner_rays[..., 1] + source_point[1],
facecolors="none",
edgecolors="red",
s=6**2,
)
plt.scatter(
source_point[0],
source_point[1],
marker="+",
color="black",
s=10**2,
);
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.
[20]:
def _build_fanned_rays(
edge_pool: np.ndarray, corner_rays: np.ndarray
) -> np.ndarray:
# calculate angles of corner_rays
angles = np.arcsin(corner_rays[..., 1] / np.linalg.norm(corner_rays, axis=1))
corner_ray_angles = np.where(corner_rays[..., 0] >= 0, angles, np.pi - angles)
# generate fanned rays
delta_angle = (
0.01 * np.min([dx, dy])
/ np.linalg.norm(corner_rays, axis=1)
)
fan_plus = np.array(
[
np.cos(corner_ray_angles + delta_angle),
np.sin(corner_ray_angles + delta_angle),
]
).swapaxes(0, 1)
fan_minus = np.array(
[
np.cos(corner_ray_angles - delta_angle),
np.sin(corner_ray_angles - delta_angle),
]
).swapaxes(0, 1)
print(f"fan_minus.shape = {fan_minus.shape}")
fan_rays = np.concatenate((fan_plus, fan_minus), axis=0)
# project fan rays to the nearest edge
edge_vectors = edge_pool[..., 1, :] - edge_pool[..., 0, :]
point_deltas = edge_pool[..., 0, :] - source_point
denominator = np.cross(fan_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(0, 1)
mu_array = np.cross(point_deltas, edge_vectors) / denominator
nu_array = (
np.cross(point_deltas[:, np.newaxis, ...], fan_rays).swapaxes(0, 1)
/ denominator
)
mu_condition = mu_array > 0
nu_condition = np.logical_and(nu_array >= 0, nu_array <= 1)
mask = np.logical_and(mu_condition, nu_condition)
# Note: rays that do not satisfy the mask conditions project to
# infinity. This happens when the insertion point is not
# located in the motion space and a boundary corner is
# fanned.
mu_array[np.logical_not(mask)] = np.inf
adjusted_mu_array = np.nanmin(mu_array, axis=1)
fan_rays = adjusted_mu_array[..., None] * fan_rays
# filter close points before merging
double_corner_rays = np.concatenate((corner_rays, corner_rays), axis=0)
mask = np.logical_and(
np.isclose(fan_rays[..., 0], double_corner_rays[..., 0], atol=.5 * dx),
np.isclose(fan_rays[..., 1], double_corner_rays[..., 1], atol=.5 * dy),
)
fan_rays = fan_rays[np.logical_not(mask)]
# filter out np.inf rays
finite_mask = np.all(np.isfinite(fan_rays), axis=1)
fan_rays = fan_rays[finite_mask]
# sort fan rays
ray_angles = np.arcsin(fan_rays[..., 1] / np.linalg.norm(fan_rays, axis=1))
ray_angles = np.where(fan_rays[..., 0] >= 0, ray_angles, np.pi - ray_angles)
sort_i = np.argsort(ray_angles)
fan_rays = fan_rays[sort_i]
return fan_rays
fan_rays = _build_fanned_rays(edge_pool, corner_rays)
fan_rays.shape
fan_minus.shape = (100, 2)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[20], line 75
71 fan_rays = fan_rays[sort_i]
72
73 return fan_rays
74
---> 75 fan_rays = _build_fanned_rays(edge_pool, corner_rays)
76 fan_rays.shape
Cell In[20], line 36, in _build_fanned_rays(edge_pool, corner_rays)
32 # project fan rays to the nearest edge
33 edge_vectors = edge_pool[..., 1, :] - edge_pool[..., 0, :]
34
35 point_deltas = edge_pool[..., 0, :] - source_point
---> 36 denominator = np.cross(fan_rays, edge_vectors[:, np.newaxis, ...]).swapaxes(0, 1)
37 mu_array = np.cross(point_deltas, edge_vectors) / denominator
38 nu_array = (
39 np.cross(point_deltas[:, np.newaxis, ...], fan_rays).swapaxes(0, 1)
File ~/checkouts/readthedocs.org/user_builds/bapsf-motion/envs/latest/lib/python3.13/site-packages/numpy/_core/numeric.py:1679, in cross(a, b, axisa, axisb, axisc, axis)
1677 b = moveaxis(b, axisb, -1)
1678 if a.shape[-1] != 3 or b.shape[-1] != 3:
-> 1679 raise ValueError(
1680 f"Both input arrays must be (arrays of) 3-dimensional vectors, "
1681 f"but they are {a.shape[-1]} and {b.shape[-1]} dimensional instead."
1682 )
1684 # Create the output array
1685 shape = *broadcast(a[..., 0], b[..., 0]).shape, 3
ValueError: Both input arrays must be (arrays of) 3-dimensional vectors, but they are 2 and 2 dimensional instead.
Merge the corner and fanned rays.
[21]:
rays = np.concatenate((corner_rays, fan_rays), axis=0)
# sort rays by angle
angles = np.arcsin(rays[..., 1] / np.linalg.norm(rays, axis=1))
angles = np.where(rays[..., 0] >= 0, angles, np.pi - angles)
sort_i = np.argsort(angles)
rays = rays[sort_i]
rays.shape
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 1
----> 1 rays = np.concatenate((corner_rays, fan_rays), axis=0)
2
3 # sort rays by angle
4 angles = np.arcsin(rays[..., 1] / np.linalg.norm(rays, axis=1))
NameError: name 'fan_rays' is not defined
[22]:
ds.mask.plot(x="x", y="y", aspect=1.25, size=12)
ax = plt.gca()
ax.set_ylim(-60, 60)
ax.set_xlim(-60, 60)
plt.scatter(
rays[..., 0] + source_point[0],
rays[..., 1] + source_point[1],
facecolors="none",
edgecolors="red",
s=6**2,
)
plt.scatter(
source_point[0],
source_point[1],
marker="+",
color="black",
s=10**2,
)
for ii in range(rays.shape[0]):
plt.plot(
[source_point[0], rays[ii, 0] + source_point[0]],
[source_point[1], rays[ii, 1] + source_point[1]],
color="red",
linewidth=0.5,
);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 8
4 ax.set_ylim(-60, 60)
5 ax.set_xlim(-60, 60)
6
7 plt.scatter(
----> 8 rays[..., 0] + source_point[0],
9 rays[..., 1] + source_point[1],
10 facecolors="none",
11 edgecolors="red",
NameError: name 'rays' is not defined
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?
Using the
source_pointand two neighboring rays inrayswe can iterate throughraysto form and around of triangles.Using Barycentric Coordinates we can say a point
\vec{p}can be written as\[\vec{p} = \lambda_1 \vec{p}_1 + \lambda_2 \vec{p}_2 + \lambda_3 \vec{p}_3\]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.
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.
[23]:
def _paint_mask(mask, rays: np.ndarray) -> xr.DataArray:
rays = np.append(rays, rays[0, ...][None, ...], axis=0)
endpoints = rays + source_point[None, :]
triangles = np.zeros((rays.shape[0] - 1, 3, 2))
triangles[..., 0, :] = source_point
triangles[..., 1, :] = endpoints[:-1, :]
triangles[..., 2, :] = endpoints[1:, :]
grid_points = np.zeros((x_coord.size, y_coord.size, 2))
grid_points[..., 0] = np.repeat(
x_coord.values[..., np.newaxis], y_coord.size, axis=1
)
grid_points[..., 1] = np.repeat(
y_coord.values[np.newaxis, ...], x_coord.size, axis=0
)
# This processes uses Barycentric coordinates to determine if a
# grid point is within the triangle.
#
# https://en.wikipedia.org/wiki/Barycentric_coordinate_system
#
# lambda shape is (x_size, y_size, N_rays)
#
# calculate lambda_3
numerator = np.cross(
grid_points[:, :, None, :] - triangles[None, None, :, 0, :],
(triangles[:, 1, :] - triangles[:, 0, :])[None, None, :, :],
)
denominator = np.cross(
triangles[:, 2, :] - triangles[:, 0, :],
triangles[:, 1, :] - triangles[:, 0, :]
)
zero_mask = denominator == 0
if np.any(zero_mask):
# denominator can be zero if all points on the triangle lie
# on a line
not_zero_mask = np.logical_not(zero_mask)
triangles = triangles[not_zero_mask, ...]
numerator = numerator[..., not_zero_mask]
denominator = denominator[not_zero_mask]
lambda_3 = numerator / denominator[None, None, ...]
# calculate lambda_2
numerator = np.cross(
grid_points[:, :, None, :] - triangles[None, None, :, 0, :],
(triangles[:, 2, :] - triangles[:, 0, :])[None, None, :, :],
)
denominator = -denominator
lambda_2 = numerator / denominator[None, None, ...]
# calculate lambda_1
lambda_1 = 1 - lambda_2 - lambda_3
# generate the conditional for each point in the motion space
# _conditional.shape = (mspace Nx, mspace Ny, N_rays)
#
lambda_1_condition = np.logical_and(lambda_1 >= 0, lambda_1 <= 1)
lambda_2_condition = np.logical_and(lambda_2 >= 0, lambda_2 <= 1)
lambda_3_condition = np.logical_and(lambda_3 >= 0, lambda_3 <= 1)
_condition = np.logical_and(
np.logical_and(lambda_1_condition, lambda_2_condition),
lambda_3_condition,
)
return mask.copy(data=np.any(_condition, axis=2))
shadow_mask = _paint_mask(ds.mask, rays)
shadow_mask
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 70
66 )
67
68 return mask.copy(data=np.any(_condition, axis=2))
69
---> 70 shadow_mask = _paint_mask(ds.mask, rays)
71 shadow_mask
NameError: name 'rays' is not defined
[24]:
shadow_mask.plot(x="x", y="y", aspect=1.25, size=12)
ax = plt.gca()
ax.set_ylim(-60, 60)
ax.set_xlim(-60, 60)
plt.scatter(
rays[..., 0] + source_point[0],
rays[..., 1] + source_point[1],
facecolors="none",
edgecolors="red",
s=6**2,
)
plt.scatter(
source_point[0],
source_point[1],
marker="+",
color="black",
s=10**2,
)
for ii in range(rays.shape[0]):
plt.plot(
[source_point[0], rays[ii, 0] + source_point[0]],
[source_point[1], rays[ii, 1] + source_point[1]],
color="red",
linewidth=0.5,
);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 1
----> 1 shadow_mask.plot(x="x", y="y", aspect=1.25, size=12)
2
3 ax = plt.gca()
4 ax.set_ylim(-60, 60)
NameError: name 'shadow_mask' is not defined
[ ]: