Add avoidance and intrusion number computation

notebooks/user_guide.ipynb

@@ -4495,6 +4495,93 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Intrusion and Avoidance numbers\n",
+    "\n",
+    "The dimensionless Intrusion number $\\mathcal{I}n$ and Avoidance number $\\mathcal{A}v$ characterise a crowd by how much pedestrians encroach on each other's personal space and how imminent collisions are ([Cordes et al., *PNAS Nexus*, 2024](https://doi.org/10.1093/pnasnexus/pgae120)). Both are computed per pedestrian and frame.\n",
+    "\n",
+    "The **Intrusion number** adds up the encroachments of all neighbours $j$ within $3\\,r_{\\text{soc}}$ of pedestrian $i$ (Eq. 1 of the paper):\n",
+    "\n",
+    "$$\n",
+    "\\mathcal{I}n_i = \\sum_{j \\in \\mathcal{N}_i}\n",
+    "\\left(\\frac{r_{\\text{soc}} - \\ell_{\\text{min}}}{r_{ij} - \\ell_{\\text{min}}}\\right)^2\n",
+    "$$\n",
+    "\n",
+    "with the social radius $r_{\\text{soc}} = 0.8\\,$m, the pedestrian diameter $\\ell_{\\text{min}} = 0.2\\,$m and the centre-to-centre distance $r_{ij}$. `IntrusionMethod.SUM` (the default) is this definition: encroachments add up, as mechanical forces do. `IntrusionMethod.MAX` is the variant of Eq. S6 in the paper's SI Appendix B: the intrusion is dominated by the closest neighbour and does not grow with the number of neighbours. The authors report that the delineation of regimes is robust under this change.\n",
+    "\n",
+    "The **Avoidance number** is based on the time-to-collision $\\tau_{ij}$, the time until $i$ and $j$ first touch if both keep their current velocity, with $\\tau_{ij} = \\infty$ if they never touch (Eq. 2 of the paper):\n",
+    "\n",
+    "$$\n",
+    "\\mathcal{A}v_i = \\sum_{j \\in \\mathcal{N}'_i} \\frac{\\tau_0}{\\tau_{ij}}\n",
+    "$$\n",
+    "\n",
+    "with the reference time $\\tau_0 = 3\\,$s. The equation is written as a sum, but the paper restricts the neighbour set $\\mathcal{N}'_i$ to the single neighbour with the shortest $\\tau_{ij}$, the most imminent risk. The sum therefore has one term, and because $\\tau_0/\\tau_{ij}$ decreases with $\\tau_{ij}$ it equals the maximum of $\\tau_0/\\tau_{ij}$ over all neighbours in the frame. That maximum is what `compute_avoidance` returns; a pedestrian on no collision course gets $\\mathcal{A}v_i = 0$.\n",
+    "\n",
+    "Recorded head trajectories produce apparent overlaps of the discs. Following the paper, both functions cap each per-pair term, at 400 for $\\mathcal{I}n_{ij}$ and 60 for $\\mathcal{A}v_{ij}$ (`cap=None` disables it). The paper also smooths the trajectories with a 4th-order Butterworth filter at 0.5 Hz and computes velocities over about 1 s; the functions do neither, so we choose `frame_step` for a 1 s window below.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pedpy import IntrusionMethod, compute_avoidance, compute_intrusion\n",
+    "\n",
+    "intrusion = compute_intrusion(traj_data=traj)\n",
+    "intrusion.head()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "frame_step = int(traj.frame_rate / 2)  # velocity window 2 * frame_step / frame_rate = 1 s, as in the paper\n",
+    "avoidance = compute_avoidance(traj_data=traj, frame_step=frame_step, tau_0=3.0)\n",
+    "avoidance.head()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Temporal evolution\n",
+    "\n",
+    "We can visualize how these numbers evolve over time by averaging over the pedestrians in each frame. We follow the paper's rules: the intrusion number is averaged over all pedestrians, isolated ones counting as 0, while the avoidance number is averaged only over pedestrians on a collision course, $\\mathcal{A}v_i > 0$. The paper samples these averages every 0.5 s over the stationary state; here we show every frame.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "mean_in = intrusion.groupby(\"frame\")[\"intrusion\"].mean()\n",
+    "mean_av = avoidance[avoidance[\"avoidance\"] > 0].groupby(\"frame\")[\"avoidance\"].mean()\n",
+    "\n",
+    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "\n",
+    "ax1.plot(mean_in.index / traj.frame_rate, mean_in.values)\n",
+    "ax1.set_xlabel(\"Time / s\")\n",
+    "ax1.set_ylabel(r\"$\\langle \\mathcal{I}n \\rangle$\")\n",
+    "ax1.set_title(\"Mean Intrusion Number\")\n",
+    "\n",
+    "ax2.plot(mean_av.index / traj.frame_rate, mean_av.values)\n",
+    "ax2.set_xlabel(\"Time / s\")\n",
+    "ax2.set_ylabel(r\"$\\langle \\mathcal{A}v \\rangle$\")\n",
+    "ax2.set_title(\"Mean Avoidance Number\")\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {

pedpy/__init__.py

@@ -12,6 +12,7 @@
 
 from .column_identifier import (
     ACC_COL,
+    AVOIDANCE_COL,
     A_X_COL,
     A_Y_COL,
     COUNT_COL,
@@ -25,6 +26,7 @@
     FRAME_COL,
     ID_COL,
     INTERSECTION_COL,
+    INTRUSION_COL,
     LAST_FRAME_COL,
     MEAN_SPEED_COL,
     MID_FRAME_COL,
@@ -96,6 +98,7 @@
     Cutoff,
     SpeedCalculation,
     compute_frame_range_in_area,
+    compute_individual_distances,
     compute_individual_voronoi_polygons,
     compute_intersecting_polygons,
     compute_neighbor_distance,
@@ -117,7 +120,12 @@
     compute_speed_profile,
     get_grid_cells,
 )
-from .methods.spatial_analysis import compute_pair_distribution_function
+from .methods.spatial_analysis import (
+    IntrusionMethod,
+    compute_avoidance,
+    compute_intrusion,
+    compute_pair_distribution_function,
+)
 from .methods.speed_calculator import (
     compute_individual_speed,
     compute_line_speed,
@@ -187,6 +195,7 @@
     "compute_frame_range_in_area",
     "compute_individual_voronoi_polygons",
     "compute_intersecting_polygons",
+    "compute_individual_distances",
     "compute_neighbors",
     "compute_neighbor_distance",
     "compute_time_distance_line",
@@ -216,6 +225,9 @@
     "compute_mean_acceleration_per_frame",
     "compute_voronoi_acceleration",
     "correct_invalid_trajectories",
+    "compute_intrusion",
+    "IntrusionMethod",
+    "compute_avoidance",
     "detect_anomalies_in_trajectories",
     "PEDPY_BLUE",
     "PEDPY_GREEN",
@@ -271,6 +283,8 @@
     "MID_POSITION_COL",
     "END_POSITION_COL",
     "WINDOW_SIZE_COL",
+    "INTRUSION_COL",
+    "AVOIDANCE_COL",
     "__version__",
     "AccelerationError",
     "GeometryError",

pedpy/column_identifier.py

@@ -44,3 +44,6 @@
 SPEED_SP2_COL: Final = "s_sp-1"
 FLOW_SP1_COL: Final = "j_sp+1"
 FLOW_SP2_COL: Final = "j_sp-1"
+
+INTRUSION_COL: Final = "intrusion"
+AVOIDANCE_COL: Final = "avoidance"

pedpy/methods/method_utils.py

@@ -1309,3 +1309,35 @@ def is_individual_speed_valid(
         return DataValidationStatus.ENTRY_MISSING
 
     return DataValidationStatus.DATA_CORRECT
+
+
+def compute_individual_distances(*, traj_data: TrajectoryData) -> pd.DataFrame:
+    """Compute pairwise distances between all pedestrians per frame.
+
+    Args:
+        traj_data: trajectory data
+
+    Returns:
+        DataFrame with columns 'id', 'frame', 'neighbor_id', and 'distance'
+    """
+    neighbor_point_col = f"{POINT_COL}_{NEIGHBOR_ID_COL}"
+    points = traj_data.data[[ID_COL, FRAME_COL, POINT_COL]]
+    neighbor_points = points.rename(
+        columns={
+            ID_COL: NEIGHBOR_ID_COL,
+            POINT_COL: neighbor_point_col,
+        }
+    )
+    matrix = points.merge(
+        neighbor_points,
+        how="inner",
+        on=FRAME_COL,
+    )
+    matrix = matrix[matrix[ID_COL] != matrix[NEIGHBOR_ID_COL]]
+    matrix[DISTANCE_COL] = np.linalg.norm(
+        shapely.get_coordinates(matrix[POINT_COL]) - shapely.get_coordinates(matrix[neighbor_point_col]),
+        axis=1,
+    )
+
+    matrix = matrix.drop(columns=[POINT_COL, neighbor_point_col])
+    return matrix.reset_index(drop=True)