From fcc2cfaa87d58d24e38dd8ac7edaf5181cf117e6 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Mon, 26 Jan 2026 13:51:21 -0800
Subject: [PATCH 01/11] bci application nilearn viz

---
 .../bci_visualization/bci_visualization.py    | 133 ++++----
 .../nilearn_visualization_operator.py         | 312 ++++++++++++++++++
 .../bci_visualization/requirements.txt        |   2 +
 3 files changed, 382 insertions(+), 65 deletions(-)
 create mode 100644 applications/bci_visualization/operators/nilearn_visualization_operator.py

diff --git a/applications/bci_visualization/bci_visualization.py b/applications/bci_visualization/bci_visualization.py
index f1a31e95e8..54fbf0e409 100644
--- a/applications/bci_visualization/bci_visualization.py
+++ b/applications/bci_visualization/bci_visualization.py
@@ -13,6 +13,7 @@
 from holoscan.operators import HolovizOp
 from holoscan.resources import CudaStreamPool, UnboundedAllocator
 from holoscan.schedulers import EventBasedScheduler
+from operators.nilearn_visualization_operator import NilearnVisualizationOperator
 from operators.reconstruction import (
     BuildRHSOperator,
     ConvertToVoxelsOperator,
@@ -109,40 +110,42 @@ def compose(self):
             xyz=pipeline_assets.xyz,
         )
 
-        # ========== Visualization Pipeline Operators ==========
-        # Get volume_renderer kwargs from YAML config to extract density range
-        volume_renderer_kwargs = self.kwargs("volume_renderer")
-        density_min = volume_renderer_kwargs.get("density_min", -100.0)
-        density_max = volume_renderer_kwargs.get("density_max", 100.0)
-
-        voxel_to_volume = VoxelStreamToVolumeOp(
-            self,
-            name="voxel_to_volume",
-            pool=volume_allocator,
-            mask_nifti_path=self._mask_path,
-            density_min=density_min,
-            density_max=density_max,
-            **self.kwargs("voxel_stream_to_volume"),
-        )
-
-        volume_renderer = VolumeRendererOp(
-            self,
-            name="volume_renderer",
-            config_file=self._rendering_config,
-            allocator=volume_allocator,
-            cuda_stream_pool=cuda_stream_pool,
-            **volume_renderer_kwargs,
-        )
-
-        # IMPORTANT changes to avoid deadlocks of volume_renderer and holoviz
-        # when running in multi-threading mode
-        # 1. Set the output port condition to NONE to remove backpressure
-        volume_renderer.spec.outputs["color_buffer_out"].condition(ConditionType.NONE)
-        # 2. Use a passthrough operator to configure queue policy as POP to use latest frame
-        color_buffer_passthrough = ColorBufferPassthroughOp(
-            self,
-            name="color_buffer_passthrough",
-        )
+        nilearn_visualization_operator = NilearnVisualizationOperator(fragment=self)
+
+        # # ========== Visualization Pipeline Operators ==========
+        # # Get volume_renderer kwargs from YAML config to extract density range
+        # volume_renderer_kwargs = self.kwargs("volume_renderer")
+        # density_min = volume_renderer_kwargs.get("density_min", -100.0)
+        # density_max = volume_renderer_kwargs.get("density_max", 100.0)
+
+        # voxel_to_volume = VoxelStreamToVolumeOp(
+        #     self,
+        #     name="voxel_to_volume",
+        #     pool=volume_allocator,
+        #     mask_nifti_path=self._mask_path,
+        #     density_min=density_min,
+        #     density_max=density_max,
+        #     **self.kwargs("voxel_stream_to_volume"),
+        # )
+
+        # volume_renderer = VolumeRendererOp(
+        #     self,
+        #     name="volume_renderer",
+        #     config_file=self._rendering_config,
+        #     allocator=volume_allocator,
+        #     cuda_stream_pool=cuda_stream_pool,
+        #     **volume_renderer_kwargs,
+        # )
+
+        # # IMPORTANT changes to avoid deadlocks of volume_renderer and holoviz
+        # # when running in multi-threading mode
+        # # 1. Set the output port condition to NONE to remove backpressure
+        # volume_renderer.spec.outputs["color_buffer_out"].condition(ConditionType.NONE)
+        # # 2. Use a passthrough operator to configure queue policy as POP to use latest frame
+        # color_buffer_passthrough = ColorBufferPassthroughOp(
+        #     self,
+        #     name="color_buffer_passthrough",
+        # )
 
         holoviz = HolovizOp(
             self,
@@ -183,43 +186,43 @@ def compose(self):
         )
         self.add_flow(
             convert_to_voxels_operator,
-            voxel_to_volume,
+            nilearn_visualization_operator,
             {
                 ("affine_4x4", "affine_4x4"),
                 ("hb_voxel_data", "hb_voxel_data"),
             },
         )
 
-        # ========== Connect Visualization Pipeline ==========
-        # voxel_to_volume → volume_renderer
-        self.add_flow(
-            voxel_to_volume,
-            volume_renderer,
-            {
-                ("volume", "density_volume"),
-                ("spacing", "density_spacing"),
-                ("permute_axis", "density_permute_axis"),
-                ("flip_axes", "density_flip_axes"),
-            },
-        )
-        # Add mask connections to VolumeRendererOp
-        self.add_flow(
-            voxel_to_volume,
-            volume_renderer,
-            {
-                ("mask_volume", "mask_volume"),
-                ("mask_spacing", "mask_spacing"),
-                ("mask_permute_axis", "mask_permute_axis"),
-                ("mask_flip_axes", "mask_flip_axes"),
-            },
-        )
-
-        # volume_renderer ↔ holoviz
-        self.add_flow(
-            volume_renderer, color_buffer_passthrough, {("color_buffer_out", "color_buffer_in")}
-        )
-        self.add_flow(color_buffer_passthrough, holoviz, {("color_buffer_out", "receivers")})
-        self.add_flow(holoviz, volume_renderer, {("camera_pose_output", "camera_pose")})
+        # # ========== Connect Visualization Pipeline ==========
+        # # voxel_to_volume → volume_renderer
+        # self.add_flow(
+        #     voxel_to_volume,
+        #     volume_renderer,
+        #     {
+        #         ("volume", "density_volume"),
+        #         ("spacing", "density_spacing"),
+        #         ("permute_axis", "density_permute_axis"),
+        #         ("flip_axes", "density_flip_axes"),
+        #     },
+        # )
+        # # Add mask connections to VolumeRendererOp
+        # self.add_flow(
+        #     voxel_to_volume,
+        #     volume_renderer,
+        #     {
+        #         ("mask_volume", "mask_volume"),
+        #         ("mask_spacing", "mask_spacing"),
+        #         ("mask_permute_axis", "mask_permute_axis"),
+        #         ("mask_flip_axes", "mask_flip_axes"),
+        #     },
+        # )
+
+        # # volume_renderer ↔ holoviz
+        # self.add_flow(
+        #     volume_renderer, color_buffer_passthrough, {("color_buffer_out", "color_buffer_in")}
+        # )
+        # self.add_flow(color_buffer_passthrough, holoviz, {("color_buffer_out", "receivers")})
+        # self.add_flow(holoviz, volume_renderer, {("camera_pose_output", "camera_pose")})
 
 
 def main():
diff --git a/applications/bci_visualization/operators/nilearn_visualization_operator.py b/applications/bci_visualization/operators/nilearn_visualization_operator.py
new file mode 100644
index 0000000000..af48fdf638
--- /dev/null
+++ b/applications/bci_visualization/operators/nilearn_visualization_operator.py
@@ -0,0 +1,312 @@
+"""Visualization operator for HbO voxel data on brain surface."""
+
+from __future__ import annotations
+
+import gzip
+import logging
+import nibabel as nib
+from nilearn.plotting import plot_surf_stat_map
+from nilearn.surface import vol_to_surf
+import struct
+from pathlib import Path
+from time import perf_counter
+from typing import Any, Tuple
+from urllib.request import urlretrieve
+import matplotlib.pyplot as plt
+
+import numpy as np
+from numpy.typing import NDArray
+
+from holoscan.core import (
+    ConditionType,
+    ExecutionContext,
+    InputContext,
+    Operator,
+    OperatorSpec,
+    OutputContext,
+)
+
+logger = logging.getLogger(__name__)
+
+# URL to the MNI152 surface mesh
+MNI152_SURFACE_URL = (
+    "https://raw.githubusercontent.com/neurolabusc/surf-ice/master/sample/mni152_2009.mz3"
+)
+
+# Default cache directory for downloaded assets
+_CACHE_DIR = Path.home() / ".cache" / "myelin"
+
+
+def _read_mz3(filepath: Path) -> Tuple[NDArray[np.int32], NDArray[np.float32]]:
+    """Read an MZ3 mesh file and return faces and vertices.
+
+    Parameters
+    ----------
+    filepath : Path
+        Path to the MZ3 file.
+
+    Returns
+    -------
+    faces : NDArray[np.int32]
+        Face indices with shape (n_faces, 3).
+    vertices : NDArray[np.float32]
+        Vertex coordinates with shape (n_vertices, 3).
+
+    Raises
+    ------
+    ValueError
+        If the file is not a valid MZ3 file.
+    """
+    hdr_bytes = 16
+
+    # First, try to read as raw binary
+    with open(filepath, "rb") as f:
+        data = f.read(hdr_bytes)
+        hdr = struct.unpack_from("<HHIII", data)
+
+    # Check if gzip compressed (magic bytes 0x1f8b) or raw MZ3 (0x4D5A = 23117)
+    is_gzip = hdr[0] != 23117
+
+    if is_gzip:
+        # Re-open with gzip
+        f = gzip.open(filepath, "rb")
+        data = f.read(hdr_bytes)
+        hdr = struct.unpack_from("<HHIII", data)
+        if hdr[0] != 23117:
+            f.close()
+            raise ValueError(f"Not a valid MZ3 file: {filepath}")
+    else:
+        f = open(filepath, "rb")
+        f.read(hdr_bytes)  # Skip header we already read
+
+    try:
+        attr = hdr[1]
+        n_face = hdr[2]
+        n_vert = hdr[3]
+        n_skip = hdr[4]
+
+        is_face = (attr & 1) != 0
+        is_vert = (attr & 2) != 0
+
+        if attr > 127:
+            raise ValueError("Unable to read future version of MZ3 file")
+        if n_vert < 1:
+            raise ValueError("Unable to read MZ3 files without vertices")
+        if n_face < 1 and is_face:
+            raise ValueError("MZ3 files with isFACE must specify NFACE")
+
+        # Read faces
+        f.seek(hdr_bytes + n_skip)
+        faces = np.zeros((0, 3), dtype=np.int32)
+        if is_face:
+            face_data = f.read(12 * n_face)  # 3 * 4 bytes per face
+            faces_flat = np.array(
+                struct.unpack_from(f"<{3 * n_face}I", face_data), dtype=np.int32
+            )
+            faces = faces_flat.reshape((n_face, 3))
+
+        # Read vertices
+        f.seek(hdr_bytes + n_skip + (is_face * n_face * 12))
+        vertices = np.zeros((0, 3), dtype=np.float32)
+        if is_vert:
+            vert_data = f.read(12 * n_vert)  # 3 * 4 bytes per vertex
+            verts_flat = np.array(
+                struct.unpack_from(f"<{3 * n_vert}f", vert_data), dtype=np.float32
+            )
+            vertices = verts_flat.reshape((n_vert, 3))
+    finally:
+        f.close()
+
+    return faces, vertices
+
+
+def _ensure_surface_mesh(cache_dir: Path | None = None) -> Path:
+    """Download the MNI152 surface mesh if not already cached.
+
+    Parameters
+    ----------
+    cache_dir : Path or None
+        Directory to cache the mesh. Defaults to ~/.cache/myelin
+
+    Returns
+    -------
+    Path
+        Path to the cached mesh file.
+    """
+    if cache_dir is None:
+        cache_dir = _CACHE_DIR
+
+    cache_dir.mkdir(parents=True, exist_ok=True)
+    mesh_path = cache_dir / "mni152_2009.mz3"
+
+    if not mesh_path.exists():
+        logger.info("Downloading MNI152 surface mesh to %s", mesh_path)
+        urlretrieve(MNI152_SURFACE_URL, mesh_path)
+        logger.info("Download complete")
+
+    return mesh_path
+
+
+def _load_surface_mesh(
+    cache_dir: Path | None = None,
+) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
+    """Load the MNI152 surface mesh.
+
+    Returns
+    -------
+    coordinates : NDArray[np.float32]
+        Vertex coordinates with shape (n_vertices, 3).
+    faces : NDArray[np.int32]
+        Face indices with shape (n_faces, 3).
+    """
+    mesh_path = _ensure_surface_mesh(cache_dir)
+    faces, vertices = _read_mz3(mesh_path)
+    return vertices, faces
+
+
+class NilearnVisualizationOperator(Operator):
+    """Operator that visualizes HbO voxel data on a brain surface.
+
+    This operator uses nilearn's vol_to_surf to project voxel data onto
+    a brain surface mesh and renders it using matplotlib in real-time.
+    """
+
+    def __init__(
+        self,
+        *,
+        fragment: Any | None = None,
+    ) -> None:
+        """Initialize the visualization operator.
+
+        Parameters
+        ----------
+        fragment : Any or None
+            Parent fragment (Application) for Holoscan wiring.
+        """
+        super().__init__(fragment, name=self.__class__.__name__)
+        self._last_frame_time: float | None = None
+        self._frame_count = 0
+        self._surface_mesh: Tuple[NDArray[np.float32], NDArray[np.int32]] | None = None
+        self._affine: NDArray[np.float32] | None = None
+        # Matplotlib state for real-time updates
+        self._fig: Any = None
+        self._axes: Any = None
+        self._initialized_plot: bool = False
+
+    def setup(self, spec: OperatorSpec) -> None:
+        """Configure operator inputs."""
+        spec.input("hb_voxel_data")
+        spec.input("affine_4x4").condition(ConditionType.NONE)
+
+    def _ensure_mesh_loaded(self) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
+        """Lazy-load the surface mesh."""
+        if self._surface_mesh is None:
+            logger.info("Loading MNI152 surface mesh...")
+            self._surface_mesh = _load_surface_mesh()
+            logger.info(
+                "Surface mesh loaded: %d vertices, %d faces",
+                self._surface_mesh[0].shape[0],
+                self._surface_mesh[1].shape[0],
+            )
+        return self._surface_mesh
+
+    def compute(
+        self,
+        op_input: InputContext,
+        op_output: OutputContext,
+        context: ExecutionContext,
+    ) -> None:
+        """Receive voxel data and render on surface."""
+        del op_output, context
+
+        hb_frame = op_input.receive("hb_voxel_data")
+        affine = op_input.receive("affine_4x4")
+
+        # Store affine on first frame
+        if affine is not None and self._affine is None:
+            self._affine = np.asarray(affine)
+            logger.info("Received affine matrix with shape %s", self._affine.shape)
+
+        now = perf_counter()
+        prev = self._last_frame_time
+        self._last_frame_time = now
+
+        if prev is None:
+            logger.info(
+                "Visualizer received HbO voxel frame with shape %s",
+                getattr(hb_frame, "shape", None),
+            )
+        else:
+            logger.debug(
+                "Visualizer received frame after %.2f ms",
+                (now - prev) * 1000.0,
+            )
+
+        self._frame_count += 1
+        self._render_frame(hb_frame)
+
+    def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
+        """Render the HbO volume on the brain surface.
+
+        Parameters
+        ----------
+        hb_volume : NDArray[np.float32]
+            3D volume of HbO values.
+        """
+        if self._affine is None:
+            logger.warning("No affine matrix received yet, skipping visualization")
+            return
+
+        # Load the surface mesh
+        coords, faces = self._ensure_mesh_loaded()
+
+        # Create a NIfTI image from the volume
+        nifti_img = nib.nifti1.Nifti1Image(hb_volume.astype(np.float32), self._affine)
+
+        # Project volume to surface using vol_to_surf with default args
+        # surf_mesh can be a tuple of (coordinates, faces) or a mesh object
+        surf_mesh = (coords, faces)
+
+        try:
+            surface_data = vol_to_surf(nifti_img, surf_mesh)
+        except Exception as e:
+            logger.warning("vol_to_surf failed: %s. Skipping frame.", e)
+            return
+
+        # Check for valid data
+        if np.all(np.isnan(surface_data)):
+            logger.warning("All surface data is NaN, skipping visualization")
+            return
+
+        # Initialize figure on first frame
+        if not self._initialized_plot:
+            plt.ion()  # Enable interactive mode for real-time updates
+            self._fig, self._axes = plt.subplots(
+                subplot_kw={"projection": "3d"}, figsize=(10, 8)
+            )
+            self._initialized_plot = True
+        else:
+            # Clear axes but preserve camera view
+            elev = self._axes.elev
+            azim = self._axes.azim
+            self._axes.clear()
+            self._axes.view_init(elev=elev, azim=azim)
+
+        # Create the surface plot, reusing existing figure/axes
+        plot_surf_stat_map(
+            surf_mesh=surf_mesh,
+            stat_map=surface_data,
+            hemi="left",  # The MNI152 mesh is a full brain, but we view from left
+            view="lateral",
+            cmap="RdBu_r",
+            colorbar=False,  # Disable to avoid accumulating colorbars
+            title=f"HbO Frame {self._frame_count}",
+            threshold=None,  # Show all values
+            figure=self._fig,
+            axes=self._axes,
+        )
+
+        # Force display update for real-time visualization
+        self._fig.canvas.draw_idle()
+        self._fig.canvas.flush_events()
+        plt.pause(0.001)
diff --git a/applications/bci_visualization/requirements.txt b/applications/bci_visualization/requirements.txt
index ae5827bd3e..9f814da0ee 100644
--- a/applications/bci_visualization/requirements.txt
+++ b/applications/bci_visualization/requirements.txt
@@ -1,4 +1,6 @@
+matplotlib==3.10.8
 nibabel==5.3.2
+nilearn==0.13.0
 numpy==2.3.3
 scipy==1.16.3
 h5py==3.15.1

From 89540b3a8b07f4c9cce8f407eb288a4082c3ce70 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Mon, 26 Jan 2026 14:07:42 -0800
Subject: [PATCH 02/11] make own surf

---
 .../nilearn_visualization_operator.py         | 45 +++++++++++++++----
 1 file changed, 37 insertions(+), 8 deletions(-)

diff --git a/applications/bci_visualization/operators/nilearn_visualization_operator.py b/applications/bci_visualization/operators/nilearn_visualization_operator.py
index af48fdf638..a47987ec00 100644
--- a/applications/bci_visualization/operators/nilearn_visualization_operator.py
+++ b/applications/bci_visualization/operators/nilearn_visualization_operator.py
@@ -5,8 +5,8 @@
 import gzip
 import logging
 import nibabel as nib
+import scipy.ndimage
 from nilearn.plotting import plot_surf_stat_map
-from nilearn.surface import vol_to_surf
 import struct
 from pathlib import Path
 from time import perf_counter
@@ -188,6 +188,7 @@ def __init__(
         self._frame_count = 0
         self._surface_mesh: Tuple[NDArray[np.float32], NDArray[np.int32]] | None = None
         self._affine: NDArray[np.float32] | None = None
+        self._voxel_coords: NDArray[np.float32] | None = None
         # Matplotlib state for real-time updates
         self._fig: Any = None
         self._axes: Any = None
@@ -257,20 +258,48 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
             logger.warning("No affine matrix received yet, skipping visualization")
             return
 
+        # Ensure we have a numpy array (handle CuPy arrays from GPU-based operators)
+        if hasattr(hb_volume, "get"):
+            hb_volume = hb_volume.get()
+
         # Load the surface mesh
         coords, faces = self._ensure_mesh_loaded()
+        surf_mesh = (coords, faces)
 
-        # Create a NIfTI image from the volume
-        nifti_img = nib.nifti1.Nifti1Image(hb_volume.astype(np.float32), self._affine)
+        # Precompute voxel coordinates for the mesh vertices if not done yet
+        if self._voxel_coords is None:
+            try:
+                inv_affine = np.linalg.inv(self._affine)
 
-        # Project volume to surface using vol_to_surf with default args
-        # surf_mesh can be a tuple of (coordinates, faces) or a mesh object
-        surf_mesh = (coords, faces)
+                # Transform vertices to voxel space: [x, y, z, 1] . inv_affine.T
+                # coords is (N, 3)
+                ones = np.ones((coords.shape[0], 1), dtype=coords.dtype)
+                vertices_homo = np.hstack((coords, ones))  # (N, 4)
+
+                # Result is (N, 4)
+                voxel_coords_homo = vertices_homo @ inv_affine.T
+
+                # Store (3, N) for map_coordinates
+                self._voxel_coords = voxel_coords_homo[:, :3].T
+
+                logger.info("Precomputed voxel mapping for %d vertices", coords.shape[0])
+
+            except Exception as e:
+                logger.error("Failed to compute voxel coordinates: %s", e)
+                return
 
         try:
-            surface_data = vol_to_surf(nifti_img, surf_mesh)
+            # Map volume data to surface using trilinear interpolation
+            # map_coordinates expects (ndim, n_points) coords
+            surface_data = scipy.ndimage.map_coordinates(
+                hb_volume.astype(np.float32),
+                self._voxel_coords,
+                order=1,
+                mode="constant",
+                cval=0.0,
+            )
         except Exception as e:
-            logger.warning("vol_to_surf failed: %s. Skipping frame.", e)
+            logger.warning("Visualization mapping failed: %s. Skipping frame.", e)
             return
 
         # Check for valid data

From c5d03a940ae019f78637bcdfde163c3d14e9ef38 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Mon, 26 Jan 2026 14:46:11 -0800
Subject: [PATCH 03/11] switch to mpl toolkit and faces

---
 .../bci_visualization/bci_visualization.py    |   6 +-
 ...rator.py => mpl_visualization_operator.py} | 140 +++++++++++++-----
 .../bci_visualization/requirements.txt        |   1 -
 3 files changed, 106 insertions(+), 41 deletions(-)
 rename applications/bci_visualization/operators/{nilearn_visualization_operator.py => mpl_visualization_operator.py} (70%)

diff --git a/applications/bci_visualization/bci_visualization.py b/applications/bci_visualization/bci_visualization.py
index 54fbf0e409..c7c91afbc6 100644
--- a/applications/bci_visualization/bci_visualization.py
+++ b/applications/bci_visualization/bci_visualization.py
@@ -13,7 +13,7 @@
 from holoscan.operators import HolovizOp
 from holoscan.resources import CudaStreamPool, UnboundedAllocator
 from holoscan.schedulers import EventBasedScheduler
-from operators.nilearn_visualization_operator import NilearnVisualizationOperator
+from operators.mpl_visualization_operator import MplVisualizationOperator
 from operators.reconstruction import (
     BuildRHSOperator,
     ConvertToVoxelsOperator,
@@ -110,7 +110,7 @@ def compose(self):
             xyz=pipeline_assets.xyz,
         )
 
-        nilearn_visualization_operator = NilearnVisualizationOperator(fragment=self)
+        mpl_visualization_operator = MplVisualizationOperator(fragment=self)
 
         # # ========== Visualization Pipeline Operators ==========
         # # Get volume_renderer kwargs from YAML config to extract density range
@@ -186,7 +186,7 @@ def compose(self):
         )
         self.add_flow(
             convert_to_voxels_operator,
-            nilearn_visualization_operator,
+            mpl_visualization_operator,
             {
                 ("affine_4x4", "affine_4x4"),
                 ("hb_voxel_data", "hb_voxel_data"),
diff --git a/applications/bci_visualization/operators/nilearn_visualization_operator.py b/applications/bci_visualization/operators/mpl_visualization_operator.py
similarity index 70%
rename from applications/bci_visualization/operators/nilearn_visualization_operator.py
rename to applications/bci_visualization/operators/mpl_visualization_operator.py
index a47987ec00..b7f17c27d4 100644
--- a/applications/bci_visualization/operators/nilearn_visualization_operator.py
+++ b/applications/bci_visualization/operators/mpl_visualization_operator.py
@@ -4,10 +4,9 @@
 
 import gzip
 import logging
-import nibabel as nib
 import scipy.ndimage
-from nilearn.plotting import plot_surf_stat_map
 import struct
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 from pathlib import Path
 from time import perf_counter
 from typing import Any, Tuple
@@ -164,13 +163,59 @@ def _load_surface_mesh(
     return vertices, faces
 
 
-class NilearnVisualizationOperator(Operator):
+def _decimate_mesh(
+    vertices: NDArray[np.float32],
+    faces: NDArray[np.int32],
+    target_faces: int = 10000,
+) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
+    """Decimate a mesh to reduce face count for faster rendering.
+
+    Uses uniform random sampling of faces. This is a simple approach
+    that works well for visualization purposes.
+
+    Parameters
+    ----------
+    vertices : NDArray[np.float32]
+        Vertex coordinates with shape (n_vertices, 3).
+    faces : NDArray[np.int32]
+        Face indices with shape (n_faces, 3).
+    target_faces : int
+        Target number of faces after decimation.
+
+    Returns
+    -------
+    vertices : NDArray[np.float32]
+        Vertex coordinates (may be reduced).
+    faces : NDArray[np.int32]
+        Decimated face indices.
+    """
+    n_faces = faces.shape[0]
+    if n_faces <= target_faces:
+        return vertices, faces
+
+    # Randomly sample faces
+    rng = np.random.default_rng(42)  # Fixed seed for reproducibility
+    selected_idx = rng.choice(n_faces, size=target_faces, replace=False)
+    selected_faces = faces[selected_idx]
+
+    # Find unique vertices used by selected faces
+    unique_verts, inverse = np.unique(selected_faces.ravel(), return_inverse=True)
+    new_vertices = vertices[unique_verts]
+    new_faces = inverse.reshape(-1, 3).astype(np.int32)
+
+    return new_vertices, new_faces
+
+
+class MplVisualizationOperator(Operator):
     """Operator that visualizes HbO voxel data on a brain surface.
 
-    This operator uses nilearn's vol_to_surf to project voxel data onto
-    a brain surface mesh and renders it using matplotlib in real-time.
+    This operator projects voxel data onto a decimated brain surface mesh
+    and renders it using matplotlib's Poly3DCollection for real-time updates.
     """
 
+    # Target face count for decimated mesh (controls performance vs quality)
+    TARGET_FACES = 12000
+
     def __init__(
         self,
         *,
@@ -192,6 +237,8 @@ def __init__(
         # Matplotlib state for real-time updates
         self._fig: Any = None
         self._axes: Any = None
+        self._poly_collection: Poly3DCollection | None = None
+        self._colormap = plt.cm.RdBu_r
         self._initialized_plot: bool = False
 
     def setup(self, spec: OperatorSpec) -> None:
@@ -200,14 +247,22 @@ def setup(self, spec: OperatorSpec) -> None:
         spec.input("affine_4x4").condition(ConditionType.NONE)
 
     def _ensure_mesh_loaded(self) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
-        """Lazy-load the surface mesh."""
+        """Lazy-load and decimate the surface mesh."""
         if self._surface_mesh is None:
             logger.info("Loading MNI152 surface mesh...")
-            self._surface_mesh = _load_surface_mesh()
+            vertices, faces = _load_surface_mesh()
+            logger.info(
+                "Original mesh: %d vertices, %d faces",
+                vertices.shape[0],
+                faces.shape[0],
+            )
+            # Decimate for real-time performance
+            vertices, faces = _decimate_mesh(vertices, faces, self.TARGET_FACES)
+            self._surface_mesh = (vertices, faces)
             logger.info(
-                "Surface mesh loaded: %d vertices, %d faces",
-                self._surface_mesh[0].shape[0],
-                self._surface_mesh[1].shape[0],
+                "Decimated mesh: %d vertices, %d faces",
+                vertices.shape[0],
+                faces.shape[0],
             )
         return self._surface_mesh
 
@@ -302,40 +357,51 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
             logger.warning("Visualization mapping failed: %s. Skipping frame.", e)
             return
 
-        # Check for valid data
-        if np.all(np.isnan(surface_data)):
-            logger.warning("All surface data is NaN, skipping visualization")
-            return
+        # Compute per-face colors by averaging vertex values
+        face_values = surface_data[faces].mean(axis=1)
 
-        # Initialize figure on first frame
+        # Normalize to [0, 1] for colormap
+        vmin, vmax = np.nanmin(face_values), np.nanmax(face_values)
+        if vmax - vmin < 1e-10:
+            vmax = vmin + 1e-10
+        normalized = (face_values - vmin) / (vmax - vmin)
+        face_colors = self._colormap(normalized)
+
+        # Initialize figure and mesh on first frame
         if not self._initialized_plot:
-            plt.ion()  # Enable interactive mode for real-time updates
+            plt.ion()
             self._fig, self._axes = plt.subplots(
                 subplot_kw={"projection": "3d"}, figsize=(10, 8)
             )
+
+            # Build polygon vertices for each face: (n_faces, 3, 3)
+            triangles = coords[faces]
+
+            self._poly_collection = Poly3DCollection(
+                triangles,
+                facecolors=face_colors,
+                edgecolors="none",
+                linewidths=0,
+            )
+            self._axes.add_collection3d(self._poly_collection)
+
+            # Set axis limits based on mesh bounds
+            self._axes.set_xlim(coords[:, 0].min(), coords[:, 0].max())
+            self._axes.set_ylim(coords[:, 1].min(), coords[:, 1].max())
+            self._axes.set_zlim(coords[:, 2].min(), coords[:, 2].max())
+            self._axes.set_box_aspect([1, 1, 1])
+            self._axes.axis("off")
+            self._axes.view_init(elev=0, azim=-90)  # Lateral view
+
             self._initialized_plot = True
+            logger.info("Initialized real-time visualization")
         else:
-            # Clear axes but preserve camera view
-            elev = self._axes.elev
-            azim = self._axes.azim
-            self._axes.clear()
-            self._axes.view_init(elev=elev, azim=azim)
-
-        # Create the surface plot, reusing existing figure/axes
-        plot_surf_stat_map(
-            surf_mesh=surf_mesh,
-            stat_map=surface_data,
-            hemi="left",  # The MNI152 mesh is a full brain, but we view from left
-            view="lateral",
-            cmap="RdBu_r",
-            colorbar=False,  # Disable to avoid accumulating colorbars
-            title=f"HbO Frame {self._frame_count}",
-            threshold=None,  # Show all values
-            figure=self._fig,
-            axes=self._axes,
-        )
-
-        # Force display update for real-time visualization
+            # Update only face colors (fast path)
+            self._poly_collection.set_facecolors(face_colors)
+
+        self._axes.set_title(f"HbO Frame {self._frame_count}")
+
+        # Force display update
         self._fig.canvas.draw_idle()
         self._fig.canvas.flush_events()
         plt.pause(0.001)
diff --git a/applications/bci_visualization/requirements.txt b/applications/bci_visualization/requirements.txt
index 9f814da0ee..42f6e4cfe3 100644
--- a/applications/bci_visualization/requirements.txt
+++ b/applications/bci_visualization/requirements.txt
@@ -1,6 +1,5 @@
 matplotlib==3.10.8
 nibabel==5.3.2
-nilearn==0.13.0
 numpy==2.3.3
 scipy==1.16.3
 h5py==3.15.1

From ab0c1beacd662d83d40fab456ced65917adef427 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Mon, 26 Jan 2026 15:00:19 -0800
Subject: [PATCH 04/11] only flush

---
 .../bci_visualization/operators/mpl_visualization_operator.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/applications/bci_visualization/operators/mpl_visualization_operator.py b/applications/bci_visualization/operators/mpl_visualization_operator.py
index b7f17c27d4..7d6c96141a 100644
--- a/applications/bci_visualization/operators/mpl_visualization_operator.py
+++ b/applications/bci_visualization/operators/mpl_visualization_operator.py
@@ -401,7 +401,5 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
 
         self._axes.set_title(f"HbO Frame {self._frame_count}")
 
-        # Force display update
-        self._fig.canvas.draw_idle()
+        # Process GUI events to update display
         self._fig.canvas.flush_events()
-        plt.pause(0.001)

From eaa7c33f6ea40e1468c0538081711d207d79c923 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Tue, 27 Jan 2026 11:42:54 -0800
Subject: [PATCH 05/11] simpler mesh

---
 applications/bci_visualization/Dockerfile     |   3 +
 .../operators/mpl_visualization_operator.py   | 328 ++++++------------
 .../bci_visualization/requirements.txt        |   1 +
 3 files changed, 103 insertions(+), 229 deletions(-)

diff --git a/applications/bci_visualization/Dockerfile b/applications/bci_visualization/Dockerfile
index 8b744b701d..bafae50df4 100644
--- a/applications/bci_visualization/Dockerfile
+++ b/applications/bci_visualization/Dockerfile
@@ -24,6 +24,9 @@ ARG DEBIAN_FRONTEND=noninteractive
 # Install curl for downloading data files
 RUN apt-get update && apt-get install -y curl && rm -rf /var/lib/apt/lists/*
 
+# Install python3-tk
+RUN apt-get update && apt-get install -y python3-tk && rm -rf /var/lib/apt/lists/*
+
 ENV HOLOSCAN_INPUT_PATH=/workspace/holohub/data/bci_visualization
 
 # Install Python dependencies
diff --git a/applications/bci_visualization/operators/mpl_visualization_operator.py b/applications/bci_visualization/operators/mpl_visualization_operator.py
index 7d6c96141a..59b893590d 100644
--- a/applications/bci_visualization/operators/mpl_visualization_operator.py
+++ b/applications/bci_visualization/operators/mpl_visualization_operator.py
@@ -2,19 +2,17 @@
 
 from __future__ import annotations
 
-import gzip
 import logging
 import scipy.ndimage
-import struct
 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
-from pathlib import Path
-from time import perf_counter
+from time import perf_counter, sleep
 from typing import Any, Tuple
-from urllib.request import urlretrieve
 import matplotlib.pyplot as plt
+import multiprocessing as mp
 
 import numpy as np
 from numpy.typing import NDArray
+from nilearn import datasets, surface
 
 from holoscan.core import (
     ConditionType,
@@ -27,129 +25,13 @@
 
 logger = logging.getLogger(__name__)
 
-# URL to the MNI152 surface mesh
-MNI152_SURFACE_URL = (
-    "https://raw.githubusercontent.com/neurolabusc/surf-ice/master/sample/mni152_2009.mz3"
-)
-
-# Default cache directory for downloaded assets
-_CACHE_DIR = Path.home() / ".cache" / "myelin"
+HB_MIN = -0.0001
+HB_MAX =  0.0001
+HB_RANGE = HB_MAX - HB_MIN
 
 
-def _read_mz3(filepath: Path) -> Tuple[NDArray[np.int32], NDArray[np.float32]]:
-    """Read an MZ3 mesh file and return faces and vertices.
-
-    Parameters
-    ----------
-    filepath : Path
-        Path to the MZ3 file.
-
-    Returns
-    -------
-    faces : NDArray[np.int32]
-        Face indices with shape (n_faces, 3).
-    vertices : NDArray[np.float32]
-        Vertex coordinates with shape (n_vertices, 3).
-
-    Raises
-    ------
-    ValueError
-        If the file is not a valid MZ3 file.
-    """
-    hdr_bytes = 16
-
-    # First, try to read as raw binary
-    with open(filepath, "rb") as f:
-        data = f.read(hdr_bytes)
-        hdr = struct.unpack_from("<HHIII", data)
-
-    # Check if gzip compressed (magic bytes 0x1f8b) or raw MZ3 (0x4D5A = 23117)
-    is_gzip = hdr[0] != 23117
-
-    if is_gzip:
-        # Re-open with gzip
-        f = gzip.open(filepath, "rb")
-        data = f.read(hdr_bytes)
-        hdr = struct.unpack_from("<HHIII", data)
-        if hdr[0] != 23117:
-            f.close()
-            raise ValueError(f"Not a valid MZ3 file: {filepath}")
-    else:
-        f = open(filepath, "rb")
-        f.read(hdr_bytes)  # Skip header we already read
-
-    try:
-        attr = hdr[1]
-        n_face = hdr[2]
-        n_vert = hdr[3]
-        n_skip = hdr[4]
-
-        is_face = (attr & 1) != 0
-        is_vert = (attr & 2) != 0
-
-        if attr > 127:
-            raise ValueError("Unable to read future version of MZ3 file")
-        if n_vert < 1:
-            raise ValueError("Unable to read MZ3 files without vertices")
-        if n_face < 1 and is_face:
-            raise ValueError("MZ3 files with isFACE must specify NFACE")
-
-        # Read faces
-        f.seek(hdr_bytes + n_skip)
-        faces = np.zeros((0, 3), dtype=np.int32)
-        if is_face:
-            face_data = f.read(12 * n_face)  # 3 * 4 bytes per face
-            faces_flat = np.array(
-                struct.unpack_from(f"<{3 * n_face}I", face_data), dtype=np.int32
-            )
-            faces = faces_flat.reshape((n_face, 3))
-
-        # Read vertices
-        f.seek(hdr_bytes + n_skip + (is_face * n_face * 12))
-        vertices = np.zeros((0, 3), dtype=np.float32)
-        if is_vert:
-            vert_data = f.read(12 * n_vert)  # 3 * 4 bytes per vertex
-            verts_flat = np.array(
-                struct.unpack_from(f"<{3 * n_vert}f", vert_data), dtype=np.float32
-            )
-            vertices = verts_flat.reshape((n_vert, 3))
-    finally:
-        f.close()
-
-    return faces, vertices
-
-
-def _ensure_surface_mesh(cache_dir: Path | None = None) -> Path:
-    """Download the MNI152 surface mesh if not already cached.
-
-    Parameters
-    ----------
-    cache_dir : Path or None
-        Directory to cache the mesh. Defaults to ~/.cache/myelin
-
-    Returns
-    -------
-    Path
-        Path to the cached mesh file.
-    """
-    if cache_dir is None:
-        cache_dir = _CACHE_DIR
-
-    cache_dir.mkdir(parents=True, exist_ok=True)
-    mesh_path = cache_dir / "mni152_2009.mz3"
-
-    if not mesh_path.exists():
-        logger.info("Downloading MNI152 surface mesh to %s", mesh_path)
-        urlretrieve(MNI152_SURFACE_URL, mesh_path)
-        logger.info("Download complete")
-
-    return mesh_path
-
-
-def _load_surface_mesh(
-    cache_dir: Path | None = None,
-) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
-    """Load the MNI152 surface mesh.
+def _load_surface_mesh() -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
+    """Load the fsaverage surface mesh using nilearn.
 
     Returns
     -------
@@ -158,64 +40,72 @@ def _load_surface_mesh(
     faces : NDArray[np.int32]
         Face indices with shape (n_faces, 3).
     """
-    mesh_path = _ensure_surface_mesh(cache_dir)
-    faces, vertices = _read_mz3(mesh_path)
-    return vertices, faces
-
-
-def _decimate_mesh(
-    vertices: NDArray[np.float32],
-    faces: NDArray[np.int32],
-    target_faces: int = 10000,
-) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
-    """Decimate a mesh to reduce face count for faster rendering.
-
-    Uses uniform random sampling of faces. This is a simple approach
-    that works well for visualization purposes.
-
-    Parameters
-    ----------
-    vertices : NDArray[np.float32]
-        Vertex coordinates with shape (n_vertices, 3).
-    faces : NDArray[np.int32]
-        Face indices with shape (n_faces, 3).
-    target_faces : int
-        Target number of faces after decimation.
-
-    Returns
-    -------
-    vertices : NDArray[np.float32]
-        Vertex coordinates (may be reduced).
-    faces : NDArray[np.int32]
-        Decimated face indices.
-    """
-    n_faces = faces.shape[0]
-    if n_faces <= target_faces:
-        return vertices, faces
-
-    # Randomly sample faces
-    rng = np.random.default_rng(42)  # Fixed seed for reproducibility
-    selected_idx = rng.choice(n_faces, size=target_faces, replace=False)
-    selected_faces = faces[selected_idx]
-
-    # Find unique vertices used by selected faces
-    unique_verts, inverse = np.unique(selected_faces.ravel(), return_inverse=True)
-    new_vertices = vertices[unique_verts]
-    new_faces = inverse.reshape(-1, 3).astype(np.int32)
-
-    return new_vertices, new_faces
+    fsaverage = datasets.fetch_surf_fsaverage('fsaverage3')
+
+    # Load left and right pial surfaces
+    coords_left, faces_left = surface.load_surf_mesh(fsaverage['pial_left'])
+    coords_right, faces_right = surface.load_surf_mesh(fsaverage['pial_right'])
+
+    # Combine into single mesh
+    combined_coords = np.vstack([coords_left, coords_right])
+    combined_faces = np.vstack([faces_left, faces_right + coords_left.shape[0]])
+
+    return combined_coords, combined_faces
+
+
+class PlotServer(mp.Process):
+    def __init__(self, vertices, faces, queue):
+        super().__init__(daemon=True)
+        self.vertices = vertices
+        self.faces = faces
+        self.queue = queue
+
+    def run(self):
+        # This runs in a fresh process whose main thread owns the GUI loop.
+        import matplotlib.pyplot as plt
+        from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+        plt.ion()
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(10, 8))
+        triangles = self.vertices[self.faces]
+        default_colors = np.full((self.faces.shape[0], 4), [0.5, 0.5, 0.5, 1.0])
+        poly = Poly3DCollection(triangles, facecolors=default_colors, edgecolors=[0.4, 0.4, 0.4])
+        ax.add_collection3d(poly)
+        ax.set_xlim(self.vertices[:,0].min(), self.vertices[:,0].max())
+        ax.set_ylim(self.vertices[:,1].min(), self.vertices[:,1].max())
+        ax.set_zlim(self.vertices[:,2].min(), self.vertices[:,2].max())
+        ax.set_box_aspect([1,1,1])
+        ax.axis("off")
+        ax.view_init(elev=0, azim=-90)
+
+        # event loop: listen for arrays of shape (n_faces, 4) (RGBA) or None to exit
+        try:
+            while True:
+                try:
+                    item = self.queue.get(timeout=0.1)  # small timeout so GUI stays responsive
+                except Exception:
+                    # no data, still allow GUI to process events
+                    fig.canvas.flush_events()
+                    sleep(0.01)
+                    continue
+
+                if item is None:
+                    break  # shutdown signal
+
+                # item expected to be face_colors (n_faces, 4)
+                poly.set_facecolors(item)
+                ax.set_title(f"HbO Frame")
+                fig.canvas.flush_events()
+        finally:
+            plt.close(fig)
 
 
 class MplVisualizationOperator(Operator):
     """Operator that visualizes HbO voxel data on a brain surface.
 
-    This operator projects voxel data onto a decimated brain surface mesh
+    This operator projects voxel data onto a fsaverage surface mesh
     and renders it using matplotlib's Poly3DCollection for real-time updates.
     """
 
-    # Target face count for decimated mesh (controls performance vs quality)
-    TARGET_FACES = 12000
-
     def __init__(
         self,
         *,
@@ -239,31 +129,41 @@ def __init__(
         self._axes: Any = None
         self._poly_collection: Poly3DCollection | None = None
         self._colormap = plt.cm.RdBu_r
-        self._initialized_plot: bool = False
 
     def setup(self, spec: OperatorSpec) -> None:
         """Configure operator inputs."""
         spec.input("hb_voxel_data")
         spec.input("affine_4x4").condition(ConditionType.NONE)
 
+    def start(self) -> None:
+        """Initialize the visualization plot."""
+        coords, faces = self._ensure_mesh_loaded()
+
+        # create queue and start plot process
+        self._plot_queue = mp.Queue(maxsize=4)
+        # send numpy arrays (these are pickled/transferred once)
+        self._plot_proc = PlotServer(coords, faces, self._plot_queue)
+        self._plot_proc.start()
+
+    def stop(self):
+        try:
+            self._plot_queue.put(None)  # sentinel to ask process to exit
+        except Exception:
+            pass
+        if hasattr(self, "_plot_proc"):
+            self._plot_proc.join(timeout=1.0)
+
     def _ensure_mesh_loaded(self) -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
-        """Lazy-load and decimate the surface mesh."""
+        """Lazy-load the surface mesh."""
         if self._surface_mesh is None:
-            logger.info("Loading MNI152 surface mesh...")
+            logger.info("Loading  surface mesh...")
             vertices, faces = _load_surface_mesh()
             logger.info(
-                "Original mesh: %d vertices, %d faces",
+                "Loaded mesh: %d vertices, %d faces",
                 vertices.shape[0],
                 faces.shape[0],
             )
-            # Decimate for real-time performance
-            vertices, faces = _decimate_mesh(vertices, faces, self.TARGET_FACES)
             self._surface_mesh = (vertices, faces)
-            logger.info(
-                "Decimated mesh: %d vertices, %d faces",
-                vertices.shape[0],
-                faces.shape[0],
-            )
         return self._surface_mesh
 
     def compute(
@@ -293,7 +193,7 @@ def compute(
                 getattr(hb_frame, "shape", None),
             )
         else:
-            logger.debug(
+            logger.info(
                 "Visualizer received frame after %.2f ms",
                 (now - prev) * 1000.0,
             )
@@ -317,9 +217,10 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
         if hasattr(hb_volume, "get"):
             hb_volume = hb_volume.get()
 
+        print("hb_volume", hb_volume.min(), hb_volume.max())
+
         # Load the surface mesh
         coords, faces = self._ensure_mesh_loaded()
-        surf_mesh = (coords, faces)
 
         # Precompute voxel coordinates for the mesh vertices if not done yet
         if self._voxel_coords is None:
@@ -361,45 +262,14 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
         face_values = surface_data[faces].mean(axis=1)
 
         # Normalize to [0, 1] for colormap
-        vmin, vmax = np.nanmin(face_values), np.nanmax(face_values)
-        if vmax - vmin < 1e-10:
-            vmax = vmin + 1e-10
-        normalized = (face_values - vmin) / (vmax - vmin)
+        print("face_values", face_values.min(), face_values.max())
+        clamped = np.clip(face_values, HB_MIN, HB_MAX)
+        normalized = (clamped - HB_MIN) / HB_RANGE
         face_colors = self._colormap(normalized)
 
-        # Initialize figure and mesh on first frame
-        if not self._initialized_plot:
-            plt.ion()
-            self._fig, self._axes = plt.subplots(
-                subplot_kw={"projection": "3d"}, figsize=(10, 8)
-            )
-
-            # Build polygon vertices for each face: (n_faces, 3, 3)
-            triangles = coords[faces]
-
-            self._poly_collection = Poly3DCollection(
-                triangles,
-                facecolors=face_colors,
-                edgecolors="none",
-                linewidths=0,
-            )
-            self._axes.add_collection3d(self._poly_collection)
-
-            # Set axis limits based on mesh bounds
-            self._axes.set_xlim(coords[:, 0].min(), coords[:, 0].max())
-            self._axes.set_ylim(coords[:, 1].min(), coords[:, 1].max())
-            self._axes.set_zlim(coords[:, 2].min(), coords[:, 2].max())
-            self._axes.set_box_aspect([1, 1, 1])
-            self._axes.axis("off")
-            self._axes.view_init(elev=0, azim=-90)  # Lateral view
-
-            self._initialized_plot = True
-            logger.info("Initialized real-time visualization")
-        else:
-            # Update only face colors (fast path)
-            self._poly_collection.set_facecolors(face_colors)
-
-        self._axes.set_title(f"HbO Frame {self._frame_count}")
-
-        # Process GUI events to update display
-        self._fig.canvas.flush_events()
+        try:
+            # non-blocking put (drop frame if queue is full)
+            self._plot_queue.put_nowait(face_colors)
+        except Exception:
+            # queue full: drop frame to avoid blocking operator
+            logger.debug("Plot queue full — dropping frame")
diff --git a/applications/bci_visualization/requirements.txt b/applications/bci_visualization/requirements.txt
index 42f6e4cfe3..5a4bfebbc0 100644
--- a/applications/bci_visualization/requirements.txt
+++ b/applications/bci_visualization/requirements.txt
@@ -4,3 +4,4 @@ numpy==2.3.3
 scipy==1.16.3
 h5py==3.15.1
 kernel-sdk==6.6.0
+nilearn==0.13.0
\ No newline at end of file

From be3aa7b4a43fc81301de3aa010d2498b8c618753 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Tue, 27 Jan 2026 13:16:58 -0800
Subject: [PATCH 06/11] adjust by min idx

---
 .../convert_to_voxels_operator.py             | 25 ++++++++++++++-----
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/applications/bci_visualization/operators/reconstruction/convert_to_voxels_operator.py b/applications/bci_visualization/operators/reconstruction/convert_to_voxels_operator.py
index 5d74c6c3ec..ee817d71e6 100644
--- a/applications/bci_visualization/operators/reconstruction/convert_to_voxels_operator.py
+++ b/applications/bci_visualization/operators/reconstruction/convert_to_voxels_operator.py
@@ -139,10 +139,23 @@ def compute(
             self._cum_hbo = data_hbo if self._cum_hbo is None else self._cum_hbo + data_hbo
             self._cum_hbr = data_hbr if self._cum_hbr is None else self._cum_hbr + data_hbr
 
-            layout = self._compute_voxel_layout(result.voxel_metadata)
-            hb_volume = self._voxelize_hbo(self._cum_hbo, layout)
+            scatter_coords, normalized_shape, min_idx = self._compute_voxel_layout(result.voxel_metadata)
+
+            # Adjust the affine ONLY ONCE to account for the shift
+            if not self._affine_sent:
+                # Mathematical correction: A_new = A_raw * Translation(min_idx)
+                # This maps the 0-indexed volume correctly into world space
+                R = self._affine[:3, :3]
+                t = self._affine[:3, 3]
+                
+                # The new translation is the world-space position of the local (0,0,0)
+                corrected_translation = R @ min_idx + t
+                self._affine[:3, 3] = corrected_translation
+                
+                self._emit_affine_once(op_output)
+
+            hb_volume = self._voxelize_hbo(self._cum_hbo, scatter_coords, normalized_shape)
 
-        self._emit_affine_once(op_output)
         op_output.emit(hb_volume, "hb_voxel_data")
 
     def _emit_affine_once(self, op_output: OutputContext) -> None:
@@ -155,9 +168,9 @@ def _emit_affine_once(self, op_output: OutputContext) -> None:
     def _voxelize_hbo(
         self,
         data_hbo: NDArray[np.float32],
-        layout: Tuple[NDArray[np.int_], Tuple[int, int, int], NDArray[np.int_]],
+        scatter_coords: NDArray[np.int_],
+        normalized_shape: Tuple[int, int, int],
     ) -> NDArray[np.float32]:
-        scatter_coords, normalized_shape, _ijk_int = layout
         scatter_coords = scatter_coords.astype(np.int32, copy=False)  # for indexing
 
         num_voxels = data_hbo.shape[0]
@@ -190,4 +203,4 @@ def _compute_voxel_layout(
         shape = tuple(int(axis_max) + 1 for axis_max in normalized.max(axis=0))
         if not all(dim > 0 for dim in shape):
             raise ValueError(f"Shape must be positive, got shape {shape}")
-        return normalized, cast(Tuple[int, int, int], shape), ijk_int
+        return normalized, cast(Tuple[int, int, int], shape), min_idx.get()

From 9642449d068ff60dd4c15b5d109e17e71697fbde Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Wed, 28 Jan 2026 13:43:28 -0800
Subject: [PATCH 07/11] add exponential moving average high pass

---
 .../operators/mpl_visualization_operator.py   |  2 +-
 .../reconstruction/build_rhs_operator.py      |  6 +++++
 .../utils/reconstruction/ema_high_pass.py     | 26 +++++++++++++++++++
 3 files changed, 33 insertions(+), 1 deletion(-)
 create mode 100644 applications/bci_visualization/utils/reconstruction/ema_high_pass.py

diff --git a/applications/bci_visualization/operators/mpl_visualization_operator.py b/applications/bci_visualization/operators/mpl_visualization_operator.py
index 59b893590d..38662e7a26 100644
--- a/applications/bci_visualization/operators/mpl_visualization_operator.py
+++ b/applications/bci_visualization/operators/mpl_visualization_operator.py
@@ -128,7 +128,7 @@ def __init__(
         self._fig: Any = None
         self._axes: Any = None
         self._poly_collection: Poly3DCollection | None = None
-        self._colormap = plt.cm.RdBu_r
+        self._colormap = plt.cm.coolwarm
 
     def setup(self, spec: OperatorSpec) -> None:
         """Configure operator inputs."""
diff --git a/applications/bci_visualization/operators/reconstruction/build_rhs_operator.py b/applications/bci_visualization/operators/reconstruction/build_rhs_operator.py
index 7bad1629b5..6dbe258b8d 100644
--- a/applications/bci_visualization/operators/reconstruction/build_rhs_operator.py
+++ b/applications/bci_visualization/operators/reconstruction/build_rhs_operator.py
@@ -12,6 +12,7 @@
 import numpy as np
 from holoscan.core import ExecutionContext, InputContext, Operator, OperatorSpec, OutputContext
 from utils.reconstruction.assets import Assets
+from utils.reconstruction.ema_high_pass import EmaHighPass
 
 from ..stream import SampleOutput
 from .types import BuildRHSOutput, VoxelMetadata
@@ -42,6 +43,8 @@ def __init__(
         )
         self._wavelengths = assets.wavelengths
 
+        self._ema_high_pass = EmaHighPass()
+
         # GPU caches (CuPy arrays on the propagated CUDA stream)
         self._mega_jacobians_gpu = None
         self._model_optical_properties_gpu = None
@@ -110,6 +113,9 @@ def compute(
             # shape is (moments, channels, wavelengths)
             cp.log(realtime_moments[0], out=realtime_moments[0])
 
+            # apply EMA high-pass filter
+            realtime_moments = self._ema_high_pass.step(realtime_moments)
+
             realtime_moments = self._apply_baseline(realtime_moments)
             if realtime_moments is None:
                 logger.info("Skipping RHS build for first frame (baseline capture)")
diff --git a/applications/bci_visualization/utils/reconstruction/ema_high_pass.py b/applications/bci_visualization/utils/reconstruction/ema_high_pass.py
new file mode 100644
index 0000000000..1ca709e85b
--- /dev/null
+++ b/applications/bci_visualization/utils/reconstruction/ema_high_pass.py
@@ -0,0 +1,26 @@
+"""
+SPDX-FileCopyrightText: Copyright (c) 2026 Kernel.
+SPDX-License-Identifier: Apache-2.0
+"""
+from math import pi
+
+FS = 4.76 # Streaming frequency of Kernel Flow in Hz
+
+class EmaHighPass:
+    """
+    Exponential Moving Average (EMA) High-Pass Filter.
+    """
+
+    def __init__(self, fs=FS, fc=0.01):
+        self.alpha = (2 * pi * fc) / (fs + 2 * pi * fc)
+        self.ema = None
+
+    def step(self, x):
+        """
+        Process one sample (or vector of channels).
+        """
+        if self.ema is None:
+            self.ema = x
+        else:
+            self.ema = self.alpha * x + (1 - self.alpha) * self.ema
+        return x - self.ema
\ No newline at end of file

From bf688188ae5242a998d4cc816cd72ceecfaf7307 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Wed, 28 Jan 2026 14:01:23 -0800
Subject: [PATCH 08/11] different viz modes

---
 applications/bci_visualization/Dockerfile     |  10 +-
 .../bci_visualization/bci_visualization.py    | 164 ++++++++++--------
 applications/bci_visualization/metadata.json  |  23 ++-
 .../requirements-matplotlib.txt               |   2 +
 .../bci_visualization/requirements.txt        |   4 +-
 5 files changed, 120 insertions(+), 83 deletions(-)
 create mode 100644 applications/bci_visualization/requirements-matplotlib.txt

diff --git a/applications/bci_visualization/Dockerfile b/applications/bci_visualization/Dockerfile
index bafae50df4..3578eda1ab 100644
--- a/applications/bci_visualization/Dockerfile
+++ b/applications/bci_visualization/Dockerfile
@@ -24,12 +24,16 @@ ARG DEBIAN_FRONTEND=noninteractive
 # Install curl for downloading data files
 RUN apt-get update && apt-get install -y curl && rm -rf /var/lib/apt/lists/*
 
-# Install python3-tk
-RUN apt-get update && apt-get install -y python3-tk && rm -rf /var/lib/apt/lists/*
+# Install python3-tk (needed for matplotlib GUI support)
+ARG WITH_MATPLOTLIB=no
+RUN if [ "$WITH_MATPLOTLIB" = "yes" ]; then apt-get update && apt-get install -y python3-tk && rm -rf /var/lib/apt/lists/*; fi
 
 ENV HOLOSCAN_INPUT_PATH=/workspace/holohub/data/bci_visualization
 
 # Install Python dependencies
 COPY applications/bci_visualization/requirements.txt /tmp/requirements.txt
 RUN pip install -r /tmp/requirements.txt --no-cache-dir
-    
\ No newline at end of file
+
+# Install optional matplotlib dependencies
+COPY applications/bci_visualization/requirements-matplotlib.txt /tmp/requirements-matplotlib.txt
+RUN if [ "$WITH_MATPLOTLIB" = "yes" ]; then pip install -r /tmp/requirements-matplotlib.txt --no-cache-dir; fi
\ No newline at end of file
diff --git a/applications/bci_visualization/bci_visualization.py b/applications/bci_visualization/bci_visualization.py
index c7c91afbc6..54c68fbc02 100644
--- a/applications/bci_visualization/bci_visualization.py
+++ b/applications/bci_visualization/bci_visualization.py
@@ -50,10 +50,12 @@ def __init__(
         coefficients_path: Path | str,
         mask_path=None,
         reg: float = RegularizedSolverOperator.REG_DEFAULT,
+        viz: str = "clara_viz",
         **kwargs,
     ):
         self._rendering_config = render_config_file
         self._mask_path = mask_path
+        self._viz = viz
 
         # Reconstruction pipeline parameters
         self._stream = stream
@@ -110,42 +112,43 @@ def compose(self):
             xyz=pipeline_assets.xyz,
         )
 
-        mpl_visualization_operator = MplVisualizationOperator(fragment=self)
-
-        # # ========== Visualization Pipeline Operators ==========
-        # # Get volume_renderer kwargs from YAML config to extract density range
-        # volume_renderer_kwargs = self.kwargs("volume_renderer")
-        # density_min = volume_renderer_kwargs.get("density_min", -100.0)
-        # density_max = volume_renderer_kwargs.get("density_max", 100.0)
-
-        # voxel_to_volume = VoxelStreamToVolumeOp(
-        #     self,
-        #     name="voxel_to_volume",
-        #     pool=volume_allocator,
-        #     mask_nifti_path=self._mask_path,
-        #     density_min=density_min,
-        #     density_max=density_max,
-        #     **self.kwargs("voxel_stream_to_volume"),
-        # )
-
-        # volume_renderer = VolumeRendererOp(
-        #     self,
-        #     name="volume_renderer",
-        #     config_file=self._rendering_config,
-        #     allocator=volume_allocator,
-        #     cuda_stream_pool=cuda_stream_pool,
-        #     **volume_renderer_kwargs,
-        # )
-
-        # # IMPORTANT changes to avoid deadlocks of volume_renderer and holoviz
-        # # when running in multi-threading mode
-        # # 1. Set the output port condition to NONE to remove backpressure
-        # volume_renderer.spec.outputs["color_buffer_out"].condition(ConditionType.NONE)
-        # # 2. Use a passthrough operator to configure queue policy as POP to use latest frame
-        # color_buffer_passthrough = ColorBufferPassthroughOp(
-        #     self,
-        #     name="color_buffer_passthrough",
-        # )
+        if self._viz == "matplotlib":
+            mpl_visualization_operator = MplVisualizationOperator(fragment=self)
+        else:
+            # ========== Visualization Pipeline Operators ==========
+            # Get volume_renderer kwargs from YAML config to extract density range
+            volume_renderer_kwargs = self.kwargs("volume_renderer")
+            density_min = volume_renderer_kwargs.get("density_min", -100.0)
+            density_max = volume_renderer_kwargs.get("density_max", 100.0)
+
+            voxel_to_volume = VoxelStreamToVolumeOp(
+                self,
+                name="voxel_to_volume",
+                pool=volume_allocator,
+                mask_nifti_path=self._mask_path,
+                density_min=density_min,
+                density_max=density_max,
+                **self.kwargs("voxel_stream_to_volume"),
+            )
+
+            volume_renderer = VolumeRendererOp(
+                self,
+                name="volume_renderer",
+                config_file=self._rendering_config,
+                allocator=volume_allocator,
+                cuda_stream_pool=cuda_stream_pool,
+                **volume_renderer_kwargs,
+            )
+
+            # IMPORTANT changes to avoid deadlocks of volume_renderer and holoviz
+            # when running in multi-threading mode
+            # 1. Set the output port condition to NONE to remove backpressure
+            volume_renderer.spec.outputs["color_buffer_out"].condition(ConditionType.NONE)
+            # 2. Use a passthrough operator to configure queue policy as POP to use latest frame
+            color_buffer_passthrough = ColorBufferPassthroughOp(
+                self,
+                name="color_buffer_passthrough",
+            )
 
         holoviz = HolovizOp(
             self,
@@ -184,45 +187,47 @@ def compose(self):
                 ("result", "result"),
             },
         )
-        self.add_flow(
-            convert_to_voxels_operator,
-            mpl_visualization_operator,
-            {
-                ("affine_4x4", "affine_4x4"),
-                ("hb_voxel_data", "hb_voxel_data"),
-            },
-        )
 
-        # # ========== Connect Visualization Pipeline ==========
-        # # voxel_to_volume → volume_renderer
-        # self.add_flow(
-        #     voxel_to_volume,
-        #     volume_renderer,
-        #     {
-        #         ("volume", "density_volume"),
-        #         ("spacing", "density_spacing"),
-        #         ("permute_axis", "density_permute_axis"),
-        #         ("flip_axes", "density_flip_axes"),
-        #     },
-        # )
-        # # Add mask connections to VolumeRendererOp
-        # self.add_flow(
-        #     voxel_to_volume,
-        #     volume_renderer,
-        #     {
-        #         ("mask_volume", "mask_volume"),
-        #         ("mask_spacing", "mask_spacing"),
-        #         ("mask_permute_axis", "mask_permute_axis"),
-        #         ("mask_flip_axes", "mask_flip_axes"),
-        #     },
-        # )
-
-        # # volume_renderer ↔ holoviz
-        # self.add_flow(
-        #     volume_renderer, color_buffer_passthrough, {("color_buffer_out", "color_buffer_in")}
-        # )
-        # self.add_flow(color_buffer_passthrough, holoviz, {("color_buffer_out", "receivers")})
-        # self.add_flow(holoviz, volume_renderer, {("camera_pose_output", "camera_pose")})
+        if self._viz == "matplotlib":
+            self.add_flow(
+                convert_to_voxels_operator,
+                mpl_visualization_operator,
+                {
+                    ("affine_4x4", "affine_4x4"),
+                    ("hb_voxel_data", "hb_voxel_data"),
+                },
+            )
+        else:
+            # ========== Connect Visualization Pipeline ==========
+            # voxel_to_volume → volume_renderer
+            self.add_flow(
+                voxel_to_volume,
+                volume_renderer,
+                {
+                    ("volume", "density_volume"),
+                    ("spacing", "density_spacing"),
+                    ("permute_axis", "density_permute_axis"),
+                    ("flip_axes", "density_flip_axes"),
+                },
+            )
+            # Add mask connections to VolumeRendererOp
+            self.add_flow(
+                voxel_to_volume,
+                volume_renderer,
+                {
+                    ("mask_volume", "mask_volume"),
+                    ("mask_spacing", "mask_spacing"),
+                    ("mask_permute_axis", "mask_permute_axis"),
+                    ("mask_flip_axes", "mask_flip_axes"),
+                },
+            )
+
+            # volume_renderer ↔ holoviz
+            self.add_flow(
+                volume_renderer, color_buffer_passthrough, {("color_buffer_out", "color_buffer_in")}
+            )
+            self.add_flow(color_buffer_passthrough, holoviz, {("color_buffer_out", "receivers")})
+            self.add_flow(holoviz, volume_renderer, {("camera_pose_output", "camera_pose")})
 
 
 def main():
@@ -247,6 +252,16 @@ def main():
         help="Path to the mask NIfTI file containing 3D integer labels (I, J, K). Optional.",
     )
 
+    parser.add_argument(
+        "-v",
+        "--viz",
+        action="store",
+        type=str,
+        dest="viz",
+        default="clara_viz",
+        help="Visualization backend to use (e.g., 'clara_viz', 'matplotlib'). Optional.",
+    )
+
     parser.add_argument(
         "-h", "--help", action="help", default=argparse.SUPPRESS, help="Help message"
     )
@@ -272,6 +287,7 @@ def main():
         coefficients_path=kernel_data / "extinction_coefficients_mua.csv",
         mask_path=args.mask_path,
         reg=RegularizedSolverOperator.REG_DEFAULT,
+        viz=args.viz,
     )
 
     # Load YAML configuration
diff --git a/applications/bci_visualization/metadata.json b/applications/bci_visualization/metadata.json
index 7b1553dc2e..22b526890f 100644
--- a/applications/bci_visualization/metadata.json
+++ b/applications/bci_visualization/metadata.json
@@ -69,9 +69,26 @@
 				}
 			]
 		},
-		"run": {
-			"command": "python3 bci_visualization.py --renderer_config <holohub_app_source>/config.json --mask_path <holohub_data_dir>/bci_visualization/anatomy_labels_high_res.nii.gz",
-			"workdir": "holohub_app_source"
+		"default_mode": "clara_viz",
+		"modes": {
+			"clara_viz": {
+				"description": "Run the application with ClaraViz for volume rendering.",
+				"run": {
+					"command": "python3 bci_visualization.py --renderer_config <holohub_app_source>/config.json --mask_path <holohub_data_dir>/bci_visualization/anatomy_labels_high_res.nii.gz",
+					"workdir": "holohub_app_source"
+				}
+			},
+			"matplotlib": {
+				"description": "Run the application with matplotlib for volume rendering.",
+				"build": {
+					"docker_build_args": ["--build-arg", "WITH_MATPLOTLIB=yes"]
+				},
+				"run": {
+					"command": "python3 bci_visualization.py --renderer_config <holohub_app_source>/config.json --mask_path <holohub_data_dir>/bci_visualization/anatomy_labels_high_res.nii.gz --viz matplotlib",
+					"workdir": "holohub_app_source"
+				}
+			}
 		}
+		
 	}
 }
diff --git a/applications/bci_visualization/requirements-matplotlib.txt b/applications/bci_visualization/requirements-matplotlib.txt
new file mode 100644
index 0000000000..40e898820d
--- /dev/null
+++ b/applications/bci_visualization/requirements-matplotlib.txt
@@ -0,0 +1,2 @@
+matplotlib==3.10.8
+nilearn==0.13.0
\ No newline at end of file
diff --git a/applications/bci_visualization/requirements.txt b/applications/bci_visualization/requirements.txt
index 5a4bfebbc0..8e2290fac7 100644
--- a/applications/bci_visualization/requirements.txt
+++ b/applications/bci_visualization/requirements.txt
@@ -1,7 +1,5 @@
-matplotlib==3.10.8
 nibabel==5.3.2
 numpy==2.3.3
 scipy==1.16.3
 h5py==3.15.1
-kernel-sdk==6.6.0
-nilearn==0.13.0
\ No newline at end of file
+kernel-sdk==6.6.0
\ No newline at end of file

From a462d0ec17ea0b16580754bab1d6d3eaf570d7c2 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Wed, 28 Jan 2026 14:05:29 -0800
Subject: [PATCH 09/11] conditional import of operator

---
 applications/bci_visualization/bci_visualization.py          | 2 +-
 .../operators/mpl_visualization_operator.py                  | 5 +----
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/applications/bci_visualization/bci_visualization.py b/applications/bci_visualization/bci_visualization.py
index 54c68fbc02..4ba38dbc24 100644
--- a/applications/bci_visualization/bci_visualization.py
+++ b/applications/bci_visualization/bci_visualization.py
@@ -13,7 +13,6 @@
 from holoscan.operators import HolovizOp
 from holoscan.resources import CudaStreamPool, UnboundedAllocator
 from holoscan.schedulers import EventBasedScheduler
-from operators.mpl_visualization_operator import MplVisualizationOperator
 from operators.reconstruction import (
     BuildRHSOperator,
     ConvertToVoxelsOperator,
@@ -113,6 +112,7 @@ def compose(self):
         )
 
         if self._viz == "matplotlib":
+            from operators.mpl_visualization_operator import MplVisualizationOperator
             mpl_visualization_operator = MplVisualizationOperator(fragment=self)
         else:
             # ========== Visualization Pipeline Operators ==========
diff --git a/applications/bci_visualization/operators/mpl_visualization_operator.py b/applications/bci_visualization/operators/mpl_visualization_operator.py
index 38662e7a26..5518b1b2d0 100644
--- a/applications/bci_visualization/operators/mpl_visualization_operator.py
+++ b/applications/bci_visualization/operators/mpl_visualization_operator.py
@@ -193,7 +193,7 @@ def compute(
                 getattr(hb_frame, "shape", None),
             )
         else:
-            logger.info(
+            logger.debug(
                 "Visualizer received frame after %.2f ms",
                 (now - prev) * 1000.0,
             )
@@ -217,8 +217,6 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
         if hasattr(hb_volume, "get"):
             hb_volume = hb_volume.get()
 
-        print("hb_volume", hb_volume.min(), hb_volume.max())
-
         # Load the surface mesh
         coords, faces = self._ensure_mesh_loaded()
 
@@ -262,7 +260,6 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
         face_values = surface_data[faces].mean(axis=1)
 
         # Normalize to [0, 1] for colormap
-        print("face_values", face_values.min(), face_values.max())
         clamped = np.clip(face_values, HB_MIN, HB_MAX)
         normalized = (clamped - HB_MIN) / HB_RANGE
         face_colors = self._colormap(normalized)

From 46036f6411bd4204c819b797b850f474ae98ce78 Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Wed, 28 Jan 2026 14:07:32 -0800
Subject: [PATCH 10/11] missing flow link

---
 applications/bci_visualization/bci_visualization.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/applications/bci_visualization/bci_visualization.py b/applications/bci_visualization/bci_visualization.py
index 4ba38dbc24..adda36d532 100644
--- a/applications/bci_visualization/bci_visualization.py
+++ b/applications/bci_visualization/bci_visualization.py
@@ -198,6 +198,15 @@ def compose(self):
                 },
             )
         else:
+            self.add_flow(
+                convert_to_voxels_operator,
+                voxel_to_volume,
+                {
+                    ("affine_4x4", "affine_4x4"),
+                    ("hb_voxel_data", "hb_voxel_data"),
+                },
+            )
+
             # ========== Connect Visualization Pipeline ==========
             # voxel_to_volume → volume_renderer
             self.add_flow(

From 7d8169010e3509ea020625c5d7d218dece5c176c Mon Sep 17 00:00:00 2001
From: Gabe Lerner <52436903+gabelerner-kernel@users.noreply.github.com>
Date: Wed, 4 Feb 2026 12:16:09 -0800
Subject: [PATCH 11/11] nicer viz

---
 .../operators/mpl_visualization_operator.py   | 35 ++++++++++++++-----
 1 file changed, 26 insertions(+), 9 deletions(-)

diff --git a/applications/bci_visualization/operators/mpl_visualization_operator.py b/applications/bci_visualization/operators/mpl_visualization_operator.py
index 5518b1b2d0..b34ddf3114 100644
--- a/applications/bci_visualization/operators/mpl_visualization_operator.py
+++ b/applications/bci_visualization/operators/mpl_visualization_operator.py
@@ -25,9 +25,7 @@
 
 logger = logging.getLogger(__name__)
 
-HB_MIN = -0.0001
-HB_MAX =  0.0001
-HB_RANGE = HB_MAX - HB_MIN
+HB_AVERAGE_WINDOW = 100  # number of frames for running average of min/max
 
 
 def _load_surface_mesh() -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
@@ -40,7 +38,7 @@ def _load_surface_mesh() -> Tuple[NDArray[np.float32], NDArray[np.int32]]:
     faces : NDArray[np.int32]
         Face indices with shape (n_faces, 3).
     """
-    fsaverage = datasets.fetch_surf_fsaverage('fsaverage3')
+    fsaverage = datasets.fetch_surf_fsaverage('fsaverage4')
 
     # Load left and right pial surfaces
     coords_left, faces_left = surface.load_surf_mesh(fsaverage['pial_left'])
@@ -66,6 +64,7 @@ def run(self):
         from mpl_toolkits.mplot3d.art3d import Poly3DCollection
         plt.ion()
         fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(10, 8))
+        fig.canvas.manager.set_window_title("Brain Visualization")
         triangles = self.vertices[self.faces]
         default_colors = np.full((self.faces.shape[0], 4), [0.5, 0.5, 0.5, 1.0])
         poly = Poly3DCollection(triangles, facecolors=default_colors, edgecolors=[0.4, 0.4, 0.4])
@@ -75,7 +74,7 @@ def run(self):
         ax.set_zlim(self.vertices[:,2].min(), self.vertices[:,2].max())
         ax.set_box_aspect([1,1,1])
         ax.axis("off")
-        ax.view_init(elev=0, azim=-90)
+        ax.view_init(elev=-1, azim=31, roll=25)
 
         # event loop: listen for arrays of shape (n_faces, 4) (RGBA) or None to exit
         try:
@@ -93,7 +92,6 @@ def run(self):
 
                 # item expected to be face_colors (n_faces, 4)
                 poly.set_facecolors(item)
-                ax.set_title(f"HbO Frame")
                 fig.canvas.flush_events()
         finally:
             plt.close(fig)
@@ -129,6 +127,8 @@ def __init__(
         self._axes: Any = None
         self._poly_collection: Poly3DCollection | None = None
         self._colormap = plt.cm.coolwarm
+        self._hb_min: list[float] = []
+        self._hb_max: list[float] = []
 
     def setup(self, spec: OperatorSpec) -> None:
         """Configure operator inputs."""
@@ -259,9 +259,26 @@ def _render_frame(self, hb_volume: NDArray[np.float32]) -> None:
         # Compute per-face colors by averaging vertex values
         face_values = surface_data[faces].mean(axis=1)
 
-        # Normalize to [0, 1] for colormap
-        clamped = np.clip(face_values, HB_MIN, HB_MAX)
-        normalized = (clamped - HB_MIN) / HB_RANGE
+        hb_min = np.percentile(face_values, 2)
+        self._hb_min.append(hb_min)
+        if len(self._hb_min) > HB_AVERAGE_WINDOW:
+            self._hb_min.pop(0)
+        hb_min = np.mean(self._hb_min)
+
+        hb_max = np.percentile(face_values, 98)
+        self._hb_max.append(hb_max)
+        if len(self._hb_max) > HB_AVERAGE_WINDOW:
+            self._hb_max.pop(0)
+        hb_max = np.mean(self._hb_max)
+
+        # Normalize to [-1, 1] for colormap, centered at 0
+        hb_abs_max = max(abs(hb_min), abs(hb_max))
+        clipped = np.clip(face_values, -hb_abs_max, hb_abs_max)
+        normalized = clipped / hb_abs_max
+
+        # Shift to [0, 1] for colormap
+        normalized = (normalized + 1) / 2
+
         face_colors = self._colormap(normalized)
 
         try: