Refactor spatial activation functions

examples/augmentation/62_synthetic_hrfs_example.ipynb

@@ -277,7 +277,7 @@
    "source": [
     "# create the spatial activation\n",
     "spatial_act = synhrf.build_spatial_activation(\n",
-    "    head_ras,\n",
+    "    head_ras.brain,\n",
     "    c3_seed,\n",
     "    spatial_scale=2 * cedalion.units.cm,\n",
     "    intensity_scale=1 * units.micromolar,\n",
@@ -359,7 +359,7 @@
     "f, ax = plt.subplots(2, 3, figsize=(9, 6))\n",
     "for i, spatial_scale in enumerate([0.5 * units.cm, 2 * units.cm, 3 * units.cm]):\n",
     "    spatial_act = synhrf.build_spatial_activation(\n",
-    "        head_ras,\n",
+    "        head_ras.brain,\n",
     "        c3_seed,\n",
     "        spatial_scale=spatial_scale,\n",
     "        intensity_scale=1 * units.micromolar,\n",
@@ -389,7 +389,7 @@
     "    ]\n",
     "):\n",
     "    spatial_act = synhrf.build_spatial_activation(\n",
-    "        head_ras,\n",
+    "        head_ras.brain,\n",
     "        c3_seed,\n",
     "        spatial_scale=2 * units.cm,\n",
     "        intensity_scale=intensity_scale,\n",
@@ -429,14 +429,14 @@
    "outputs": [],
    "source": [
     "spatial_act_c3 = synhrf.build_spatial_activation(\n",
-    "    head_ras,\n",
+    "    head_ras.brain,\n",
     "    c3_seed,\n",
     "    spatial_scale=2 * cedalion.units.cm,\n",
     "    intensity_scale=1 * units.micromolar,\n",
     "    hbr_scale=-0.4,\n",
     ")\n",
     "spatial_act_c4 = synhrf.build_spatial_activation(\n",
-    "    head_ras,\n",
+    "    head_ras.brain,\n",
     "    c4_seed,\n",
     "    spatial_scale=2 * cedalion.units.cm,\n",
     "    intensity_scale=1 * units.micromolar,\n",

examples/machine_learning/51_multimodal_source_decomposition_on_synthetic_fNIRS-EEG_data.ipynb

@@ -707,8 +707,8 @@
    "outputs": [],
    "source": [
     "# Build fnirs spatial activation\n",
-    "s_spatial_fnirs = build_spatial_activation(head_model=head,\n",
-    "                                           seed_vertex=seed_vertex_fnirs,\n",
+    "s_spatial_fnirs = build_spatial_activation(surface=head.brain,\n",
+    "                                           seed_vertices=seed_vertex_fnirs,\n",
     "                                           spatial_scale= 20 * units.mm,\n",
     "                                           intensity_scale= 1 * units.micromolar,\n",
     "                                           hbr_scale=-0.4)\n",
@@ -2149,7 +2149,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "cedalion_250526",
+   "display_name": "cedalion",
    "language": "python",
    "name": "python3"
   },

src/cedalion/dataclasses/geometry.py

@@ -50,7 +50,7 @@ class Surface(ABC):
     crs: str
     units: pint.Unit
 
-    vertex_coords : dict[str, ArrayLike] = field(default_factory=dict)
+    vertex_coords: dict[str, ArrayLike] = field(default_factory=dict)
 
     @property
     @abstractmethod
@@ -316,7 +316,7 @@ class VTKSurface(Surface):
     def vertices(self) -> cdt.LabeledPointCloud:
         vertices = vtk_to_numpy(self.mesh.GetPoints().GetData())
         coords = {"label": np.arange(len(vertices))}
-        coords.update({k : ("label", v) for k,v in self.vertex_coords.items()})
+        coords.update({k: ("label", v) for k, v in self.vertex_coords.items()})
 
         result = xr.DataArray(
             vertices,
@@ -731,9 +731,11 @@ def geodesic_distance(self, verts, m=1.0, fem=False):
             goodrows = np.nonzero(~np.array(lfac.sum(0) == 0).ravel())[0]
             self._goodrows = goodrows
             self._rlfac_solvers[m] = sparse.linalg.factorized(
-                lfac[goodrows][:, goodrows]
+                lfac.tocsc()[goodrows][:, goodrows].tocsc()
+            )
+            self._nLC_solvers[m] = sparse.linalg.factorized(
+                nLC.tocsc()[goodrows][:, goodrows].tocsc()
             )
-            self._nLC_solvers[m] = sparse.linalg.factorized(nLC[goodrows][:, goodrows])
 
         # I. "Integrate the heat flow ̇u = ∆u for some fixed time t"
         # ---------------------------------------------------------
@@ -898,3 +900,66 @@ def affine_transform_from_numpy(
     units = cedalion.units.Unit(to_units) / cedalion.units.Unit(from_units)
 
     return xr.DataArray(transform, dims=[to_crs, from_crs]).pint.quantify(units)
+
+
+def robust_geodesic_distance(
+    surface: TrimeshSurface, seed_vertices: list[int] | int, m: float = 10.0
+) -> np.ndarray:
+    """Compute geodesic distances from seed vertices robustly across multiple submeshes.
+
+    Args:
+        surface (TrimeshSurface): The surface.
+        seed_vertices (list[int] | int): Seed vertices (single int or list of ints).
+        m (float): Geodesic distance parameter.
+
+    Returns:
+        np.ndarray: Distance array of shape (num_vertices,), inf for vertices
+                    not reachable (not on the same submesh as any seed).
+
+    Initial Contributors:
+        - Thomas Fischer | [email protected] | 2025
+    """
+
+    mesh = surface.mesh
+
+    if isinstance(seed_vertices, int):
+        seed_vertices = [seed_vertices]
+
+    seed_positions = mesh.vertices[seed_vertices]
+    mesh_split = mesh.split(only_watertight=False)
+    distances_from_seed = np.ones(mesh.vertices.shape[0]) * np.inf
+
+    for submesh in mesh_split:
+        submesh_set = set(map(tuple, submesh.vertices))
+        seeds_in_submesh = [pos for pos in seed_positions if tuple(pos) in submesh_set]
+
+        if not seeds_in_submesh:
+            continue
+
+        seed_indices_in_submesh = [
+            np.where((submesh.vertices == pos).all(axis=1))[0][0]
+            for pos in seeds_in_submesh
+        ]
+
+        cortex_surface = PycortexSurface(
+            SimpleMesh(submesh.vertices, submesh.faces),
+            crs=surface.crs,
+            units=mesh.units,
+        )
+
+        distances_on_submesh = cortex_surface.geodesic_distance(
+            seed_indices_in_submesh, m=m
+        )
+
+        submesh_vertex_indices = [
+            i
+            for i, coord in enumerate(map(tuple, mesh.vertices))
+            if coord in submesh_set
+        ]
+
+        distances_from_seed[submesh_vertex_indices] = np.minimum(
+            distances_from_seed[submesh_vertex_indices],
+            distances_on_submesh,
+        )
+
+    return distances_from_seed

src/cedalion/sim/synthetic_hrf.py

@@ -23,84 +23,45 @@
 
 
 def build_spatial_activation(
-    head_model: cfm.TwoSurfaceHeadModel,
-    seed_vertex: int,
+    surface: cdg.TrimeshSurface,
+    seed_vertices: list[int] | int,
     spatial_scale: cdt.QLength = 1 * units.cm,
     intensity_scale: cdt.QConcentration = 1 * units.micromolar,
     hbr_scale: float = None,
     m: float = 10.0,
 ):
-    """Generates a spatial activation at a seed vertex.
+    """Generates a spatial activation at one or multiple seed vertices.
 
     This function generates a blob of activity on the brain surface.
     The blob is centered at the seed vertex.
     Geodesic distances, and therefore also the blob, can be distorded
     due to mesh decimation or unsuitable m value.
 
     Args:
-        head_model (cfm.TwoSurfaceHeadModel): Head model with brain and scalp surfaces.
-        seed_vertex (int): Index of the seed vertex.
+        surface (cdg.TrimeshSurface): Brain or scalp surface.
+        seed_vertices (list[int] | int): Indices of the seed vertices.
         spatial_scale (Quantity): Scale of the spatial size.
         intensity_scale (Quantity): Scaling factor for the intensity of the blob.
-        hbr_scale (float): Scaling factor for HbR relative to HbO. If None, the blob
-            will have no concentration dimension and only represent HbO.
-        m (float): Geodesic distance parameter. Larger values of m will smooth &
-            regularize the distance computation. Smaller values of m will roughen and
-            will usually increase error in the distance computation.
+        hbr_scale (float): Scaling factor for HbR relative to HbO.
+        m (float): Geodesic distance parameter.
 
     Returns:
         xr.DataArray: Spatial image with activation values for each vertex.
 
     Initial Contributors:
         - Thomas Fischer | [email protected] | 2024
-
     """
 
-    spatial_scale_unit = (
-        (spatial_scale / head_model.brain.units).to_base_units().magnitude
-    )
-
-    seed_pos = head_model.brain.mesh.vertices[seed_vertex]
+    if isinstance(seed_vertices, int):
+        seed_vertices = [seed_vertices]
 
-    # if the mesh is not contiguous:
-    # only calculate distances on the submesh on which the seed vertex lies
+    spatial_scale_unit = (spatial_scale / surface.units).to_base_units().magnitude
 
-    # get a list of the distinct submeshes
-    mesh_split = head_model.brain.mesh.split(only_watertight=False)
+    distances_from_seeds = cdg.robust_geodesic_distance(surface, seed_vertices, m=m)
 
-    # check in which submesh the seed vertex is
-    for i, submesh in enumerate(mesh_split):
-        if seed_pos in submesh.vertices:
-            break
-
-    # get index of the seed vertex in the submesh
-    seed_vertex = np.where((submesh.vertices == seed_pos).all(axis=1))[0]
-
-    # create a pycortex surface of the submesh to calculate geodesic distances
-    cortex_surface = cdg.PycortexSurface(
-        cdg.SimpleMesh(submesh.vertices, submesh.faces),
-        crs=head_model.brain.crs,
-        units=head_model.brain.units,
-    )
-    distances_on_submesh = cortex_surface.geodesic_distance([seed_vertex], m=m)
-
-    # find indices of submesh in original mesh
-    # convert meshes into set of tuples for fast lookup
-    submesh_set = set(map(tuple, submesh.vertices))
-    submesh_indices = [
-        i
-        for i, coord in enumerate(map(tuple, head_model.brain.mesh.vertices))
-        if coord in submesh_set
-    ]
-
-    # set distances on vertices outside of the submesh to inf
-    distances_from_seed = np.ones(head_model.brain.mesh.vertices.shape[0]) * np.inf
-    distances_from_seed[submesh_indices] = distances_on_submesh
-
-    # plug the distances in a normal distribution
+    # plug geodesic distances in gaussian pdf
     norm_pdf = stats.norm(scale=spatial_scale_unit).pdf
-
-    blob_img = norm_pdf(distances_from_seed)
+    blob_img = norm_pdf(distances_from_seeds)
     blob_img = blob_img / np.max(blob_img)
     blob_img = xr.DataArray(blob_img, dims=["vertex"])
 
@@ -113,7 +74,6 @@ def build_spatial_activation(
         )
 
     blob_img = blob_img * intensity_scale
-
     blob_img = blob_img.pint.to(units.molar)
 
     return blob_img
@@ -358,6 +318,7 @@ def plot_spatial_activation(
     plt_pv.add_text(title, position="upper_edge", font_size=20)
     plt_pv.show()
 
+
 class RandomGaussianSum(TemporalBasisFunction):
     r"""A single HRF composed of Gaussians with time-dependent random weights."""
 
@@ -405,16 +366,17 @@ def __call__(self, ts: cdt.NDTimeSeries) -> xr.DataArray:
         spread = ((self.t_end - self.t_start) / 2.5).magnitude
         mu_float = mu.to("s").magnitude
         envelope = np.exp(-((mu_float - mid) ** 2) / spread**2)
-        stds = self.weight_std_min + (
-            self.weight_std_max - self.weight_std_min
-            ) * envelope
+        stds = (
+            self.weight_std_min + (self.weight_std_max - self.weight_std_min) * envelope
+        )
 
         weights = self.rng.normal(loc=self.weight_mean, scale=stds)  # <-- changed
 
         # Generate weighted Gaussians
-        gaussians = np.exp(
-            -((t_hrf[:, None] - mu[None, :]) ** 2) / self.t_std**2
-        ) * weights[None, :]
+        gaussians = (
+            np.exp(-((t_hrf[:, None] - mu[None, :]) ** 2) / self.t_std**2)
+            * weights[None, :]
+        )
 
         hrf = np.sum(gaussians, axis=1)
         window = tukey(len(t_hrf), alpha=0.1)

tests/test_synthetic_hrf.py

@@ -24,14 +24,14 @@ def test_build_spatial_activation(head_model):
 
     for seed in seeds:
         blob_small = syn.build_spatial_activation(
-            head_model,
+            head_model.brain,
             seed,
             spatial_scale=scale_small,
             intensity_scale=intensity_scale,
             hbr_scale=-0.4,
         )
         blob_big = syn.build_spatial_activation(
-            head_model,
+            head_model.brain,
             seed,
             spatial_scale=scale_big,
             intensity_scale=intensity_scale,