surfacetools module

class clabtoolkit.surfacetools.Surface(surface_file=None, vertices=None, faces=None, color='#f0f0f0', alpha=1.0, hemi=None)

Bases: object

Comprehensive class for loading and visualizing brain surface data.

Provides interface for working with brain surface geometries including loading from files or arrays, managing scalar maps and parcellations, and creating visualizations using PyVista. Supports FreeSurfer and other surface formats.

surf

Path to surface file if loaded from file.

Type:

str or None

mesh

PyVista mesh object containing surface geometry and data.

Type:

pv.PolyData

hemi

Hemisphere designation (‘lh’, ‘rh’, or ‘unknown’).

Type:

str

colortables

Dictionary storing color table information for parcellations.

Type:

dict

Examples

>>> # Load from FreeSurfer surface file
>>> surface = Surface('lh.pial')
>>>
>>> # Create from vertex/face arrays
>>> surface = Surface(vertices=verts, faces=faces, hemi='lh')
>>>
>>> # Load scalar data and parcellations
>>> surface.load_scalar_map('thickness.mgh', 'thickness')
>>> surface.load_annotation('lh.aparc.annot', 'aparc')

__init__(surface_file=None, vertices=None, faces=None, color='#f0f0f0', alpha=1.0, hemi=None)

Initialize Surface object from file, arrays, or create empty instance.

Parameters:

• surface_file (str or Path, optional) – Path to surface file (FreeSurfer .pial, .white, .inflated). Default is None.

• vertices (np.ndarray, optional) – Vertex coordinates array with shape (n_vertices, 3). Default is None.

• faces (np.ndarray, optional) – Face connectivity array with shape (n_faces, 3). Default is None.

• color (str or np.ndarray, optional) – Color for the surface mesh. Can be a hex color string (e.g., ‘#f0f0f0’) or an RGB array in [0, 1] or [0, 255] range. Default is ‘#f0f0f0’.

• alpha (float, optional) – Alpha transparency value in [0, 1] for the surface color. Default is 1.0 (opaque).

• hemi (str, optional) – Hemisphere designation (‘lh’ or ‘rh’). Auto-detected from filename if None. Default is None.

Raises:

• ValueError – If both surface_file and vertices/faces are provided, or if only one of vertices/faces is provided.

• FileNotFoundError – If surface file doesn’t exist.

Examples

>>> # Load from file with auto-detection
>>> surface = Surface('lh.pial')
>>>
>>> # Create from arrays
>>> vertices = np.random.rand(100, 3)
>>> faces = np.array([[0, 1, 2], [1, 2, 3]])
>>> surface = Surface(vertices=vertices, faces=faces, hemi='lh')
>>>
>>> # Create empty instance
>>> surface = Surface()

load_from_file(surface_file, color='#f0f0f0', alpha=1.0, hemi=None)

Load surface geometry from FreeSurfer or compatible surface file.

Parameters:

• surface_file (str or Path) – Path to surface file (e.g., FreeSurfer .pial, .white, .inflated).

• color (str or np.ndarray, optional) – Color for the surface mesh. Can be a hex color string (e.g., ‘#f0f0f0’) or an RGB array in [0, 1] or [0, 255] range. Default is ‘#f0f0f0’.

• alpha (float, optional) – Alpha transparency value in [0, 1] for the surface color. Default is 1.0 (opaque).

• hemi (str, optional) – Hemisphere designation (‘lh’ or ‘rh’). Auto-detected from filename if None. Default is None.

Raises:

• FileNotFoundError – If surface file cannot be found.

• ValueError – If surface file format is unsupported or corrupted.

Notes

• Automatically detects hemisphere from filename if not provided.

• Converts color string to RGB array if provided as hex.

• Adds alpha channel to color for RGBA representation.

• Creates default parcellation data for visualization.

• Uses nibabel to read FreeSurfer geometry files.

• Handles both left (‘lh’) and right (‘rh’) hemisphere surfaces.

• If color is not provided, defaults to light gray (‘#f0f0f0’).

Examples

>>> surface = Surface()
>>> surface.load_from_file('lh.pial')
>>> print(f"Loaded {surface.mesh.n_points} vertices")
>>>
>>> # Explicit hemisphere specification
>>> surface.load_from_file('brain_surface.surf', hemi='rh')

load_from_arrays(vertices, faces, color='#f0f0f0', alpha=1.0, hemi=None, surface_file=None)

Load surface geometry from vertex and face arrays.

Parameters:

• vertices (np.ndarray) – Vertex coordinates with shape (n_vertices, 3).

• faces (np.ndarray) – Face connectivity with shape (n_faces, 3).

• color (str or np.ndarray, optional) – Color for the surface mesh. Can be a hex color string (e.g., ‘#f0f0f0’) or an RGB array in [0, 1] or [0, 255] range. Default is ‘#f0f0f0’.

• alpha (float, optional) – Alpha transparency value in [0, 1] for the surface color. Default is 1.0 (opaque).

• hemi (str, optional) – Hemisphere designation (‘lh’ or ‘rh’). Defaults to ‘lh’.

• surface_file (str, optional) – Associated surface file path for metadata. Default is None.

Raises:

ValueError – If vertices or faces arrays have incorrect shapes.

Examples

>>> # Basic triangle mesh
>>> vertices = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]])
>>> faces = np.array([[0, 1, 2]])
>>> surface = Surface()
>>> surface.load_from_arrays(vertices, faces, hemi='lh')

load_from_mesh(mesh, color='#f0f0f0', alpha=1.0, hemi=None)

Load surface geometry from existing PyVista mesh object.

Parameters:

• mesh (pv.PolyData) – PyVista mesh object containing surface geometry.

• color (str or np.ndarray, optional) – Color for the surface mesh. Can be a hex color string (e.g., ‘#f0f0f0’) or an RGB array in [0, 1] or [0, 255] range. Default is ‘#f0f0f0’.

• alpha (float, optional) – Alpha transparency value in [0, 1] for the surface color. Default is 1.0 (opaque).

• hemi (str, optional) – Hemisphere designation (‘lh’ or ‘rh’). Defaults to ‘lh’.

Notes

Creates a deep copy of the input mesh to avoid modifying the original. Adds default surface colors if not present in the mesh.

Examples

>>> # From existing PyVista mesh
>>> existing_mesh = pv.PolyData(vertices, faces)
>>> surface = Surface()
>>> surface.load_from_mesh(existing_mesh, hemi='rh')
>>>
>>> # From procedural mesh
>>> sphere = pv.Sphere(radius=50)
>>> surface.load_from_mesh(sphere, hemi='lh')

is_loaded()

Check whether surface data has been loaded.

Returns:

True if surface data is loaded, False otherwise.

Return type:

bool

Examples

>>> surface = Surface()
>>> print(surface.is_loaded())  # False
>>> surface.load_from_file('lh.pial')
>>> print(surface.is_loaded())  # True

create_mesh_from_arrays(vertices, faces)

Create PyVista mesh object from vertex and face arrays.

Parameters:

• vertices (np.ndarray) – Vertex coordinates with shape (n_vertices, 3).

• faces (np.ndarray) – Face connectivity with shape (n_faces, 3).

Returns:

PyVista mesh object with vertices, faces, and default surface colors.

Return type:

pv.PolyData

Raises:

ValueError – If arrays have incorrect shapes or face indices are invalid.

Notes

Validates input arrays and creates properly formatted PyVista mesh with default surface colors. Adds normals to point data if provided.

Examples

>>> surface = Surface()
>>> vertices = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]])
>>> faces = np.array([[0, 1, 2]])
>>> mesh = surface.create_mesh_from_arrays(vertices, faces)
>>> print(f"Created mesh with {mesh.n_points} vertices")

get_vertices()

Get vertex coordinates from the surface mesh.

Returns:

Array of vertex coordinates with shape (n_vertices, 3).

Return type:

np.ndarray

Raises:

RuntimeError – If no surface data has been loaded.

Examples

>>> surface = Surface('lh.pial')
>>> vertices = surface.get_vertices()
>>> print(f"Surface has {len(vertices)} vertices")
>>> print(f"First vertex: {vertices[0]}")

get_faces()

Get face connectivity from the surface mesh.

Returns:

Array of face indices with shape (n_faces, 3). Each row contains three vertex indices forming a triangular face.

Return type:

np.ndarray

Raises:

RuntimeError – If no surface data has been loaded.

Notes

Extracts face connectivity from PyVista’s internal format which stores faces as [n_vertices, vertex_id1, vertex_id2, …]. This method returns only the vertex indices in standard format.

Examples

>>> surface = Surface('lh.pial')
>>> faces = surface.get_faces()
>>> print(f"Surface has {len(faces)} triangular faces")
>>> print(f"First face connects vertices: {faces[0]}")

get_edges(return_counts=False)

Extract unique edges from a triangular mesh using vectorized operations.

This function efficiently extracts all unique edges from a triangular mesh represented as a faces array. Each triangle contributes three edges, and the function automatically removes duplicates that occur when triangles share edges.

Parameters:

return_counts (bool, optional) – If True, also return the count of how many faces each edge belongs to. This is useful for identifying boundary edges (count=1) vs interior edges (count=2). Default is False.

Returns:

• edges (np.ndarray of shape (n_edges, 2)) – Array of unique edges where each row contains two vertex indices [v1, v2] with v1 <= v2. Edges are sorted lexicographically.

• counts (np.ndarray of shape (n_edges,), optional) – Number of faces that contain each edge. Only returned if return_counts=True. Boundary edges have count=1, interior edges have count=2.

Raises:

ValueError – If faces array does not have exactly 3 columns (not triangular). If faces array is empty. If faces array contains negative indices.

Return type:

ndarray

Examples

>>> # Simple triangular mesh with 2 triangles sharing an edge
>>> faces = np.array([[0, 1, 2], [1, 3, 2]])
>>> edges = get_edges(faces)
>>> print(edges)
[[0 1]
[0 2]
[1 2]
[1 3]
[2 3]]

>>> # Get edge counts to identify boundary vs interior edges
>>> edges, counts = get_edges(faces, return_counts=True)
>>> boundary_edges = edges[counts == 1]
>>> interior_edges = edges[counts == 2]
>>> print("Boundary edges:", boundary_edges)
>>> print("Interior edges:", interior_edges)
Boundary edges: [[0 1]
                [0 2]
                [1 3]
                [2 3]]
Interior edges: [[1 2]]

>>> # Cube mesh (8 vertices, 12 triangular faces)
>>> cube_faces = np.array([
...     [0, 1, 2], [0, 2, 3],  # Bottom face
...     [4, 5, 6], [4, 6, 7],  # Top face
...     [0, 1, 5], [0, 5, 4],  # Front face
...     [2, 3, 7], [2, 7, 6],  # Back face
...     [0, 3, 7], [0, 7, 4],  # Left face
...     [1, 2, 6], [1, 6, 5]   # Right face
... ])
>>> edges = get_edges(cube_faces)
>>> print(f"Cube has {len(edges)} unique edges")
Cube has 18 unique edges

Notes

This function uses vectorized NumPy operations for high performance on large meshes. The algorithm:

1. Extracts all three edges from each triangle simultaneously

2. Sorts vertex pairs to canonical form (smaller index first)

3. Uses numpy.unique to efficiently remove duplicates

Time complexity: O(n log n) where n is the number of faces Space complexity: O(n) for intermediate arrays

For non-triangular meshes, use the general extract_edges_from_faces function instead.

The canonical edge representation ensures that edge (i, j) and edge (j, i) are treated as the same edge, with the final representation always having the smaller vertex index first.

See also

extract_edges_from_faces

General version for arbitrary polygon meshes

numpy.unique

Used internally for deduplication

get_boundary_edges()

Extract only the boundary edges from a triangular mesh.

Boundary edges are those that belong to only one triangle, indicating the mesh boundary or holes in the mesh.

Parameters:

faces (np.ndarray of shape (n_faces, 3)) – Triangular mesh faces array.

Returns:

boundary_edges – Array of boundary edges where each edge belongs to only one face.

Return type:

np.ndarray of shape (n_boundary_edges, 2)

Examples

>>> # Mesh with a hole (incomplete sphere)
>>> faces = np.array([[0, 1, 2], [1, 3, 2], [3, 4, 2]])
>>> boundary = get_boundary_edges(faces)
>>> print("Boundary edges:", boundary)

get_manifold_edges()

Extract only the manifold (interior) edges from a triangular mesh.

Manifold edges are those shared by exactly two triangles, indicating proper mesh topology without boundaries or non-manifold geometry.

Parameters:

faces (np.ndarray of shape (n_faces, 3)) – Triangular mesh faces array.

Returns:

manifold_edges – Array of manifold edges where each edge belongs to exactly two faces.

Return type:

np.ndarray of shape (n_manifold_edges, 2)

Examples

>>> # Closed mesh (tetrahedron)
>>> faces = np.array([[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 2, 3]])
>>> manifold = get_manifold_edges(faces)
>>> print("Manifold edges:", manifold)

compute_normals()

Compute and store vertex normals for the surface mesh.

Calculates unit normal vectors for each vertex and stores them in the mesh point data under the key “Normals”. Normals are automatically normalized to unit length.

Raises:

• RuntimeError – If no surface data has been loaded.

• RuntimeError – If computed normals have zero length and cannot be normalized.

Notes

Uses PyVista’s built-in normal computation which averages face normals at each vertex. The resulting normals are forced to be unit vectors. Overwrites any existing normals in the mesh.

Examples

>>> surface = Surface('lh.pial')
>>> surface.compute_normals()
>>> normals = surface.get_normals()
>>> print(f"Computed {len(normals)} unit normal vectors")
>>>
>>> # Check that normals are unit vectors
>>> norms = np.linalg.norm(normals, axis=1)
>>> print(f"Normal lengths range: {norms.min():.3f} - {norms.max():.3f}")

get_normals()

Get vertex normals from the surface mesh if available.

Returns:

Array of normal vectors with shape (n_vertices, 3) if normals have been computed, None otherwise.

Return type:

np.ndarray or None

Notes

Returns None if normals haven’t been computed yet. Use compute_normals() to calculate normals before calling this method.

Examples

>>> surface = Surface('lh.pial')
>>> normals = surface.get_normals()
>>> if normals is not None:
...     print(f"Found {len(normals)} normal vectors")
... else:
...     print("No normals computed yet")
...     surface.compute_normals()
...     normals = surface.get_normals()

load_annotation(annotation, parc_name=None)

Load parcellation annotation onto surface for visualization.

Loads FreeSurfer annotation files or AnnotParcellation objects, storing labels and color information for region-based visualization.

Parameters:

• annotation (str , Path or AnnotParcellation) – Path to annotation file (.annot) or AnnotParcellation object.

• parc_name (str) – Name for parcellation reference in visualizations.

Raises:

• FileNotFoundError – If annotation file cannot be found.

• ValueError – If invalid input type or vertex count mismatch.

Examples

>>> # Load Desikan-Killiany parcellation
>>> surface.load_annotation('lh.aparc.annot', 'aparc')
>>>
>>> # Load from object
>>> annot = AnnotParcellation('lh.aparc.a2009s.annot')
>>> surface.load_annotation(annot, 'destrieux')

load_scalar_maps(scalar_map, annotation=None, maps_names=None)

Load data from a FreeSurfer vertex-wise map, a numpy array, a CSV file or pandas Dataframe onto surface for visualization.

Handles both vertex-wise data (one value per vertex) and region-wise data (requires annotation for mapping to vertices). It is important that the CSV file has a header row with column names because the first row is used to name the maps.

If it contains region-wise data then the Annotation file is mandatory.

Parameters:

• scalar_map (str, Path, pd.DataFrame) – Path to a FreeSurfer vertex-wise map file, a CSV file, a numpy array or a pandas DataFrame.

• annotation (str, Path or AnnotParcellation, optional) – Annotation file/object for mapping region data to vertices. Required if the Dataframe has region-wise data. Default is None.

• maps_names (str or list, optional) – Names for scalar data. If None, uses column names from CSV. Default is None.

Raises:

• FileNotFoundError – If map file or annotation file cannot be found.

• ValueError – If annot_file required but not provided or invalid type.

• ValueError – If maps_names length does not match number of columns in CSV.

Notes

Automatically detects if the array, the CSV or the array contains vertex-wise or region-wise data based on the number of rows. If the number of rows matches the number of vertices in the mesh, it is treated as vertex-wise data. Otherwise, it is treated as region-wise data and requires an annotation file to map region values to vertices.

The annotation can be provided as a file path, a string or as an AnnotParcellation object. If the annotation file is not found, it will try to load it from the colortables associated with the surface.

If maps_names is not provided, it uses the column names from the CSV file. If the annotation is provided as a string, it will try to load it as an AnnotParcellation object. If it is provided as an AnnotParcellation object, it will use it directly.

Examples

>>> surf_lh = cltsurf.Surface("/opt/freesurfer/subjects/fsaverage/surf/lh.pial")

>>> # Example 1: Reading a region-wise map from a CSV file and selecting a specific column name
>>> print("Example 1: Reading a region-wise map from a CSV file with a specific column name")
>>> surf_lh.load_maps_scalar_maps("/tmp/values.csv",
                                annotation="/opt/freesurfer/subjects/fsaverage/label/lh.aparc.annot",
                                maps_names="region_index")

>>> print("Loaded maps from CSV file for an specific column name:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> # Example 2: Reading a region-wise map from a dataframe and selecting a column with specified name
>>> import pandas as pd
>>> values_df = pd.read_csv("/tmp/values.csv")
>>> print("Example 2: Reading a region-wise map from a DataFrame with specified names")
>>> surf_lh.load_maps_scalar_maps(values_df,
                                annotation="/opt/freesurfer/subjects/fsaverage/label/lh.aparc.annot",
                                maps_names=["value"])

>>> print(" Loaded maps from DataFrame with specified names:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> # Example 3: Reading a region-wise map from a numpy array without specifiying names
>>> import pandas as pd
>>> print("Example 3: Reading a region-wise map from a numpy array without specifying names")
>>> values_df = pd.read_csv("/tmp/values.csv")
>>> surf_lh.load_maps_scalar_maps(values_df.to_numpy(),
                                annotation="/opt/freesurfer/subjects/fsaverage/label/lh.aparc.annot")

>>> print("Loaded maps from numpy array without specifying names:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> ######### Creating a csv with values as the number of vertices
>>> import pandas as pd
>>> n_points = surf_lh.mesh.n_points
>>> values_df = pd.DataFrame({'vertex_index': np.arange(n_points), 'vertex_value': np.random.rand(n_points)})

>>> values_df.to_csv("/tmp/values-vertexwise.csv", index=False)

>>> # Example 4: Loading vertex-wise maps from a CSV file
>>> print("Example 4: Loading vertex-wise maps from a CSV file")
>>> surf_lh.load_maps_scalar_maps("/tmp/values-vertexwise.csv"
                                )
>>> print("Loaded vertex-wise maps from CSV file:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> # Example 5: Reading a region-wise map from a numpy array and an specified name
>>> print("Example 5: Reading a region-wise map from a numpy array with specified names")
>>> import numpy as np
>>> values_array = np.random.rand(n_points)
>>> surf_lh.load_maps_scalar_maps(values_array,
                                maps_names=["ex5_vertex_value_array"])
>>> print("Loaded vertex-wise maps from numpy array:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> # Example 6: Creating a numpy array without specifying names
>>> values_array = np.random.rand(n_points)
>>> print("Example 6: Creating a numpy array with values as the number of vertices without specifying names")
>>> surf_lh.load_maps_scalar_maps(values_array)
>>> print("Loaded vertex-wise maps from numpy array without specifying names:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> # Example 7: Reading a FreeSurfer map file
>>> print("Example 7: Reading a FreeSurfer map file")
>>> surf_lh.load_maps_scalar_maps("/opt/freesurfer/subjects/fsaverage/surf/lh.thickness",
                                maps_names=["cthickness"])
>>> print("Loaded vertex-wise maps from FreeSurfer map file:")
>>> print(surf_lh.list_overlays())
>>> print("")

>>> # Example 8: Reading multiple FreeSurfer map files
>>> print("Example 8: Reading multiple FreeSurfer map files specifying names")
>>> list_of_maps = [
    "/opt/freesurfer/subjects/fsaverage/surf/lh.thickness",
    "/opt/freesurfer/subjects/fsaverage/surf/lh.curv",
    "/opt/freesurfer/subjects/fsaverage/surf/lh.sulc"]

>>> maps_names = ["thickness", "curvature", "sulc"]

>>> for i, map_file in enumerate(list_of_maps):
    surf_lh.load_maps_scalar_maps(map_file, maps_names=maps_names[i])
>>> print("Loaded vertex-wise maps from FreeSurfer map files:")
>>> print(surf_lh.list_overlays())
>>> print("")

separate_mesh_components(component_labels='components', labels_to_extract=None, clean_mesh=True, preserve_order=False)

Separate a mesh into independent submeshes based on connected component labels.

This method extracts disconnected components from a mesh, creating separate PolyData objects for each component. Each submesh contains only the vertices, faces, and point data belonging to that specific component.

Parameters:

• mesh (pv.PolyData) – Input PyVista mesh containing multiple disconnected components. Must be a triangulated surface mesh (PolyData).

• component_labels (str, default="components") – Name of the point data array containing component labels for each vertex. This array should contain integer labels identifying which component each vertex belongs to.

• labels_to_extract (List[int], optional) – Specific component labels to extract. If None, extracts all unique labels found in the component_labels array. Use this to extract only specific components of interest.

• clean_mesh (bool, default=True) – If True, removes unused vertices and ensures faces use consecutive vertex indices in each submesh. If False, preserves original vertex indices which may result in sparse vertex arrays.

• preserve_order (bool, default=False) – If True, submeshes are returned in the same order as labels_to_extract or sorted label order. If False, may return in arbitrary order for better performance.

Returns:

List of PyVista PolyData objects, one for each connected component. Each submesh contains: - Vertices belonging to that component - Faces using only those vertices (with reindexed vertex references) - All original point data arrays, filtered for the component’s vertices - Original mesh properties and metadata where applicable

Return type:

List[pv.PolyData]

Raises:

• TypeError – If mesh is not a PyVista PolyData object.

• ValueError – If component_labels field doesn’t exist, contains invalid data, or if the mesh has no faces.

• KeyError – If the specified component_labels field is not found in point data.

Examples

>>> import pyvista as pv
>>> import numpy as np
>>>
>>> # Create a mesh with two disconnected triangles
>>> points1 = np.array([[0, 0, 0], [1, 0, 0], [0.5, 1, 0]])
>>> points2 = np.array([[2, 0, 0], [3, 0, 0], [2.5, 1, 0]])
>>> points = np.vstack([points1, points2])
>>> faces1 = np.array([[3, 0, 1, 2]])  # First triangle
>>> faces2 = np.array([[3, 3, 4, 5]])  # Second triangle
>>> faces = np.vstack([faces1, faces2])
>>>
>>> mesh = pv.PolyData(points, faces)
>>>
>>> # Add component labels (would normally come from connected_components)
>>> mesh.point_data["components"] = np.array([1, 1, 1, 2, 2, 2])
>>> mesh.point_data["temperature"] = np.random.rand(6)  # Additional data
>>>
>>> # Separate into independent meshes
>>> submeshes = separate_mesh_components(mesh)
>>> print(f"Original mesh: {mesh.n_points} points, {mesh.n_faces} faces")
>>> for i, submesh in enumerate(submeshes):
...     comp_label = submesh.point_data["components"][0]
...     print(f"Component {comp_label}: {submesh.n_points} points, {submesh.n_faces} faces")
Original mesh: 6 points, 2 faces
Component 1: 3 points, 1 faces
Component 2: 3 points, 1 faces

>>> # Extract only specific components
>>> component_1_only = separate_mesh_components(mesh, labels_to_extract=[1])
>>> print(f"Component 1 only: {len(component_1_only)} submesh(es)")
Component 1 only: 1 submesh(es)

>>> # Check that point data is preserved
>>> original_temp = mesh.point_data["temperature"]
>>> submesh_temps = [sm.point_data["temperature"] for sm in submeshes]
>>> print("Temperature data preserved in submeshes:",
...       len(submesh_temps[0]), len(submesh_temps[1]))
Temperature data preserved in submeshes: 3 3

Notes

• Only triangular faces are currently supported

• Point data arrays are automatically filtered and copied to submeshes

• Cell data is not preserved as faces are restructured

• Vertex indices in faces are automatically remapped to be consecutive

• Empty components (no faces) are excluded from results

list_overlays()

List all available surface overlays and their data types.

Categorizes loaded data based on array dimensions and properties to identify scalar maps, color data, normals, and other overlay types.

Returns:

Dictionary mapping overlay names to their types: - ‘scalar’: 1D arrays of scalar values per vertex - ‘color’: 2D arrays with RGB color values (shape: n_vertices, 3) - ‘normals’: 2D arrays with unit normal vectors (shape: n_vertices, 3) - ‘unknown’: Arrays with other dimensions or unrecognized format

Return type:

dict

Notes

Automatically detects data type based on: - 1D arrays: Classified as scalar data - 2D arrays with 3 columns: Checked for unit vectors (normals) vs colors - Other dimensions: Classified as unknown

Normal vectors are identified by having unit length (norm ≈ 1) and containing negative values.

Examples

>>> # Load various data types
>>> surface.load_scalar_map('thickness.mgh', 'thickness')
>>> surface.load_annotation('aparc.annot', 'aparc')
>>> surface.compute_normals()
>>>
>>> # List all overlays
>>> overlays = surface.list_overlays()
>>> print(overlays)
{'surface': 'color', 'thickness': 'scalar', 'aparc': 'scalar', 'Normals': 'normals'}
>>>
>>> # Filter for scalar maps only
>>> scalar_maps = {k: v for k, v in overlays.items() if v == 'scalar'}
>>> print(f"Available scalar maps: {list(scalar_maps.keys())}")

set_active_overlay(overlay_name)

Set the active overlay for visualization.

Parameters:

overlay_name (str) – Name of the overlay to set as active

Return type:

None

Raises:

ValueError – If the specified overlay is not found in mesh point data

Examples

>>> surface.set_active_overlay("thickness")
>>> surface.set_active_overlay("aparc")

remove_overlay(overlay_name)

Set the active overlay for visualization.

Designates which data array should be used as the primary scalar field for coloring and visualization in PyVista plots. This affects how the surface is colored when rendered.

Parameters:

overlay_name (str) – Name of the overlay to set as active. Must exist in mesh point data.

Raises:

ValueError – If the specified overlay is not found in mesh point data.

Notes

The active overlay determines which data is used for: - Surface coloring in visualizations - Colormap application - Scalar value display in interactive plots

PyVista uses the active scalars for automatic coloring unless explicitly overridden in visualization methods.

Examples

>>> # Set thickness as active for visualization
>>> surface.set_active_overlay('thickness')
>>>
>>> # Switch to parcellation display
>>> surface.set_active_overlay('aparc')
>>>
>>> # Check available overlays first
>>> overlays = surface.list_overlays()
>>> if 'curvature' in overlays:
...     surface.set_active_overlay('curvature')

get_overlay_info(overlay_name)

Get information about a specific surface overlay.

Parameters:

overlay_name (str) – Name of the overlay to query.

Returns:

Dictionary containing overlay metadata with keys:

• ’name’str

  Name of the overlay.

• ’data_shape’tuple

  Shape of the overlay data array.

• ’data_type’str

  NumPy data type of the overlay values.

• ’has_colortable’bool

  Whether the overlay has an associated color table.

• ’num_regions’int, optional

  Number of regions (if parcellation overlay).

• ’region_names’list of str, optional

  Names of regions (if parcellation overlay).

• ’has_annot_object’bool, optional

  Whether annotation object is available (if parcellation overlay).

Return type:

Dict

Raises:

ValueError – If the overlay is not found.

Examples

>>> surface = Surface()
>>> info = surface.get_overlay_info("aparc")
>>> print(f"Overlay has {info['num_regions']} regions")
>>> print(f"Data type: {info['data_type']}")

get_region_vertices(parc_name, region_name)

Get vertices indices for a specific region in a parcellation.

Parameters:

• parc_name (str) – Name of the parcellation

• region_name (str) – Name of the region

Returns:

Array of vertices indices belonging to the region

Return type:

np.ndarray

Raises:

• ValueError – If the parcellation is not found

• ValueError – If the region is not found in the parcellation

Examples

>>> # Get vertices in the precentral gyrus
>>> vertices = surface.get_region_vertices("aparc", "precentral")
>>> print(f"Precentral region has {len(vertices)} vertices")
>>>
>>> # Get all vertices in superior frontal region
>>> vertices = surface.get_region_vertices("aparc", "superiorfrontal")

get_vertexwise_colors(overlay_name='default', colormap='viridis', vmin=None, vmax=None, range_min=None, range_max=None, range_color=(128, 128, 128, 255))

Compute vertices colors for visualization based on the specified overlay.

This method processes the overlay data and creates appropiate vertices colors for visualization, handling both scalar data (with colormaps) and categorical data (with discrete color tables).

Parameters:

• overlay_name (str, optional) – Name of the overlay to visualize. If None, the first available overlay is used.

• colormap (str, optional) – Colormap to use for scalar overlays. If None, uses parcellation color table for categorical data or ‘viridis’ for scalar data.

• vmin (np.float64, optional) – Minimum value for scaling the colormap. If None, uses the minimum value of the overlay

• vmax (np.float64, optional) – Maximum value for scaling the colormap. If None, uses the maximum value of the overlay

• None (If both vmin and vmax are)

• values. (the colormap will be applied to the full range of the overlay)

• provided (If both are)

• colormap. (they will be used to scale the)

• range_min (np.float64, optional) – Minimum value for defining a special color range. Values below this will be colored with range_color

• range_max (np.float64, optional) – Maximum value for defining a special color range. Values above this will be colored with range_color

• range_color (Tuple, optional) – RGBA color to use for values outside the defined range (range_min, range_max)

Returns:

vertices_colors – Array of RGBA colors for each vertex in the mesh.

Return type:

np.ndarray

Raises:

• ValueError – If the specified overlay is not found in the mesh point data

• ValueError – If no overlays are available

Notes

This method sets the vertices colors based on the specified overlay.

Examples

>>> # Prepare colors for a parcellation (uses discrete colors)
>>> surface.get_vertexwise_colors(overlay_name="aparc")
>>>
>>> # Prepare colors for scalar data with custom colormap
>>> surface.get_vertexwise_colors(overlay_name="thickness", colormap="hot")
>>>
>>> # Prepare colors for the surface overlay
>>> surface.get_vertexwise_colors()

prepare_colors(overlay_name=None, cmap='viridis', vmin=None, vmax=None, range_min=None, range_max=None, range_color=(128, 128, 128, 255))

Prepare vertices colors for visualization based on the specified overlay.

This method processes the overlay data and creates appropiate vertices colors for visualization, handling both scalar data (with colormaps) and categorical data (with discrete color tables).

Parameters:

• overlay_name (str, optional) – Name of the overlay to visualize. If None, the first available overlay is used.

• cmap (str, optional) – Colormap to use for scalar overlays. If None, uses parcellation color table for categorical data or ‘viridis’ for scalar data.

• vmin (np.float64, optional) – Minimum value for scaling the colormap. If None, uses the minimum value of the overlay

• vmax (np.float64, optional) – Maximum value for scaling the colormap. If None, uses the maximum value of the overlay

• None (If both vmin and vmax are)

• values. (the colormap will be applied to the full range of the overlay)

• provided (If both are)

• colormap. (they will be used to scale the)

• range_min (np.float64, optional) – Minimum value for defining a special color range. Values below this will be colored with range_color

• range_max (np.float64, optional) – Maximum value for defining a special color range. Values above this will be colored with range_color

• range_color (Tuple, optional) – RGBA color to use for values outside the defined range (range_min, range_max)

Return type:

None

Raises:

• ValueError – If the specified overlay is not found in the mesh point data

• ValueError – If no overlays are available

Notes

This method sets the vertices colors in the mesh based on the specified overlay. The colors are stored in the mesh’s point_data under the key “RGB” and set as the active scalars for visualization.

Examples

>>> # Prepare colors for a parcellation (uses discrete colors)
>>> surface.prepare_colors(overlay_name="aparc")
>>>
>>> # Prepare colors for scalar data with custom colormap
>>> surface.prepare_colors(overlay_name="thickness", cmap="hot")
>>>
>>> # Prepare colors for first available overlay
>>> surface.prepare_colors()

add_surface(surf2add)

Merge this surface with others into a single surface.

This method merges multiple Surface objects by combining their geometries and point data. Only point_data fields that are present in ALL surfaces are retained in the merged result.

Parameters:

surf2add (str, Path or Surface) – Surface to add. It can be a file path (str or Path) to a surface file or another Surface object. If a file path is provided, it will be loaded.

Returns:

New merged Surface object with hemisphere set to “unknown”

Return type:

Surface

Raises:

• TypeError – If surfaces is not a list or contains non-Surface objects

• ValueError – If the surfaces list is empty

Examples

>>> # Merge left and right hemisphere surfaces
>>> lh_surf = Surface("lh.pial")
>>> rh_surf = Surface("rh.pial")
>>> merged = lh_surf.add_surface([rh_surf])
>>> print(f"Merged surface has {merged.mesh.n_points} vertices")
>>>
>>> # Merge multiple surfaces from file paths
>>> surf1 = Surface("surface1.pial")
>>> merged = surf1.add_surface("surface2.pial")

map_volume_to_surface(image, method='nilearn', interp_method='linear', overlay_name=None)

Map volumetric neuroimaging data onto a surface mesh.

This function projects 3D or 4D volumetric data (e.g., structural or functional MRI) onto vertices of a surface mesh using spatial interpolation. It supports multiple projection methods and can handle both file paths and in-memory data arrays.

Parameters:

• image (str, np.ndarray, or nibabel.Nifti1Image) – Input volumetric data to project. Can be: - String: File path to NIfTI image (.nii, .nii.gz, .mgz) - Path: Path object to NIfTI image - numpy.ndarray: 3D or 4D array of volumetric data - nibabel.Nifti1Image: Loaded NIfTI image object

• method (str, default "nilearn") – Projection method to use: - “nilearn”: Uses nilearn.surface.vol_to_surf (recommended for neuroimaging) - “clabtoolkit”: Uses clabtoolkit interpolation functions

• interp_method (str, default "linear") – Interpolation method: - “linear”: Trilinear interpolation (smooth, good for continuous data) - “nearest”: Nearest-neighbor interpolation (preserves discrete values)

• overlay_name (str, optional) – Name to store projected data in surf_obj.mesh.point_data dictionary. Only used with method=”clabtoolkit”. If None, data is not stored.

Returns:

Projected surface data: - For 3D input: 1D array of shape (n_vertices,) - For 4D input: 2D array of shape (n_vertices, n_timepoints) Values correspond to interpolated intensity at each surface vertex.

Return type:

np.ndarray

Raises:

• FileNotFoundError – If image file path does not exist.

• ValueError – If image format is unsupported or has incompatible dimensions.

• ImportError – If required dependencies (nilearn) are not installed.

• RuntimeError – If projection fails due to coordinate system mismatch.

Notes

Coordinate Systems: - nilearn method expects vertices in world coordinates (mm) and handles coordinate transformations internally using the image affine matrix - clabtoolkit method expects vertices in voxel coordinates and requires manual coordinate conversion for NIfTI images

Performance: - nilearn is optimized for neuroimaging data and handles 4D efficiently - clabtoolkit may be slower for 4D data as it processes timepoints sequentially

4D Data Support: - Both methods support 4D fMRI data (x, y, z, time) - Returns (vertices, timepoints) array for temporal analysis

Examples

Project a structural image onto cortical surface:

>>> surface_data = map_volume_to_surface(
...     surf_obj,
...     "T1w.nii.gz",
...     method="nilearn"
... )
>>> print(f"Projected {len(surface_data)} vertex values")

Project 4D fMRI data for time series analysis:

>>> fmri_surface = map_volume_to_surface(
...     surf_obj,
...     "task_fmri.nii.gz",
...     method="nilearn",
...     interp_method="linear"
... )
>>> print(f"Shape: {fmri_surface.shape}")  # (vertices, timepoints)
>>> vertex_timeseries = fmri_surface[1000, :]  # Time series for vertex 1000

Use numpy array input with clabtoolkit method:

>>> import numpy as np
>>> volume_data = np.random.rand(64, 64, 30)  # Synthetic 3D data
>>> surface_data = map_volume_to_surface(
...     surf_obj,
...     volume_data,
...     method="clabtoolkit",
...     overlay_name="random_data"
... )
>>> # Data is now stored in surf_obj.mesh.point_data["random_data"]

Compare interpolation methods:

>>> linear_data = map_volume_to_surface(surf_obj, "mask.nii.gz",
...                                     interp_method="linear")
>>> nearest_data = map_volume_to_surface(surf_obj, "mask.nii.gz",
...                                      interp_method="nearest")

See also

nilearn.surface.vol_to_surf

Underlying nilearn projection function

clabtoolkit.imagetools.interpolate

Underlying clabtoolkit interpolation

save_surface(filename, format='freesurfer', save_annotation=None, map_name=None, overwrite=False)

Save the surface mesh to a file in the specified format.

Exports the surface geometry (vertices and faces) and optionally associated data to various file formats including FreeSurfer, VTK, PLY, STL, and OBJ.

Parameters:

• filename (str) – Output filename with or without extension. Extension will be added automatically if missing for some formats.

• format (str, default "freesurfer") – Output format: ‘freesurfer’, ‘vtk’, ‘ply’, ‘stl’, or ‘obj’.

• save_annotation (str, optional) – Path to save annotation file (for parcellation data). Only applicable for FreeSurfer.

• map_name (str, optional) – Name of overlay/parcellation to include with the surface data.

• overwrite (bool, default False) – Whether to overwrite existing files.

Raises:

• ValueError – If filename is invalid, format is unsupported, or file exists and overwrite is False.

• FileNotFoundError – If the output directory does not exist.

Examples

>>> surface.save_surface("lh.pial.vtk", format="vtk")
>>> surface.save_surface("cortex.ply", format="ply", overwrite=True)

export_to_obj(filename, save_annotation=None, map_name=None, overwrite=False)

Export the surface mesh to an OBJ file format.

Writes surface geometry as a Wavefront OBJ file, which stores vertices and triangular faces in a simple text format widely supported by 3D software and visualization tools.

Parameters:

• filename (str) – Output filename, should end with .obj extension.

• save_annotation (str, optional) – Path to save associated annotation file in FreeSurfer format.

• map_name (str, optional) – Name of parcellation/overlay to export alongside the geometry.

• overwrite (bool, default False) – Whether to overwrite existing files.

Raises:

• ValueError – If filename is invalid or file exists and overwrite is False.

• FileNotFoundError – If the output directory does not exist.

Notes

OBJ format uses 1-based indexing for face connectivity. The exported file includes vertex coordinates and triangular face definitions.

Examples

>>> surface.export_to_obj("brain_surface.obj")
>>> surface.export_to_obj("lh.pial.obj", save_annotation="lh.aparc.annot", map_name="aparc")

export_to_pyvista(filename, save_annotation=None, map_name=None, overwrite=False)

Export surface to VTK, STL, or PLY format using PyVista.

Saves the surface mesh in formats supported by PyVista, preserving geometry and optionally scalar data or colors. The file format is determined by the filename extension.

Parameters:

• filename (str) – Output filename with extension (.vtk, .ply, or .stl).

• save_annotation (str, optional) – Path to save annotation file in FreeSurfer format.

• map_name (str, optional) – Name of overlay to include as scalar data or vertex colors.

• overwrite (bool, default False) – Whether to overwrite existing files.

Raises:

• ValueError – If filename is invalid, map_name not found, or file exists and overwrite is False.

• FileNotFoundError – If the output directory does not exist.

Notes

VTK format can store additional scalar data and colors. PLY and STL formats primarily store geometry. When map_name is specified, the overlay data is prepared as vertex colors using associated colortables.

Examples

>>> surface.export_to_pyvista("brain.vtk")
>>> surface.export_to_pyvista("surface.ply", map_name="thickness")

export_to_freesurfer(filename, save_annotation=None, map_name=None, overwrite=False)

Export surface to FreeSurfer binary format.

Saves the surface mesh in FreeSurfer’s native binary geometry format, which efficiently stores vertex coordinates and triangular face connectivity for neuroimaging applications.

Parameters:

• filename (str) – Output filename, typically without extension (e.g., ‘lh.pial’).

• save_annotation (str, optional) – Path to save annotation file containing parcellation data.

• map_name (str, optional) – Name of parcellation to export with the annotation file.

• overwrite (bool, default False) – Whether to overwrite existing files.

Raises:

• ValueError – If filename is invalid or file exists and overwrite is False.

• FileNotFoundError – If the output directory does not exist.

Notes

FreeSurfer format is a compact binary representation optimized for neuroimaging workflows. The format stores only geometry data; additional data like parcellations are saved separately as .annot files.

Examples

>>> surface.export_to_freesurfer("lh.pial")
>>> surface.export_to_freesurfer("rh.white", save_annotation="rh.aparc.annot", map_name="aparc")

export_annotation(filename, parc_name, overwrite=False)

Export parcellation data to a FreeSurfer annotation file.

Saves vertex-wise parcellation labels, associated color lookup table, and region names in FreeSurfer’s .annot format for use with FreeSurfer tools and visualization software.

Parameters:

• filename (str) – Output filename for annotation file (should end with .annot).

• parc_name (str) – Name of parcellation overlay to export from the surface data.

• overwrite (bool, default False) – Whether to overwrite existing files.

Raises:

• ValueError – If filename or parc_name is invalid, parcellation not found, or file exists and overwrite is False.

• FileNotFoundError – If the output directory does not exist.

Notes

Requires the parcellation to have an associated colortable with region names and colors. The annotation format preserves the mapping between vertex labels, region names, and visualization colors.

Examples

>>> surface.export_annotation("lh.aparc.annot", "aparc")
>>> surface.export_annotation("rh.destrieux.annot", "destrieux", overwrite=True)

plot(overlay_name=None, cmap='viridis', vmin=None, vmax=None, range_min=None, range_max=None, range_color=(128, 128, 128, 255), use_opacity=False, views=['lateral'], views_orientation='grid', hemi='lh', notebook=False, show_colorbar=False, colorbar_title=None, colorbar_position='bottom', save_path=None)

Plot the surface with specified overlay and visualization parameters.

Renders the surface mesh with optional overlays using PyVista, supporting multiple camera views, custom colormaps, and interactive or static output. Handles both categorical parcellation data and continuous scalar overlays.

Parameters:

• overlay_name (str, default "default") – Name of the overlay to visualize from the surface’s point data.

• cmap (str, optional) – Colormap for scalar data. If None, uses parcellation colors for categorical data or ‘viridis’ for scalar data.

• vmin (float, optional) – Minimum value for colormap scaling. If None, uses data minimum.

• vmax (float, optional) – Maximum value for colormap scaling. If None, uses data maximum.

• range_min (float, optional) – Minimum data value to include in the visualization. Values below this threshold will be displayed in range_color.

• range_max (float, optional) – Maximum data value to include in the visualization. Values above this threshold will be displayed in range_color.

• range_color (Tuple, default (128, 128, 128, 255)) – RGBA color for values outside the specified data range.

• use_opacity (bool, default False) – Whether to apply opacity based on data values.

• views (str or List[str], default ["lateral"]) – Camera view(s): ‘lateral’, ‘medial’, ‘dorsal’, ‘ventral’, ‘anterior’, ‘posterior’, or multiple views like [‘lateral’, ‘medial’]. Also supports preset layouts: ‘4_views’, ‘6_views’, ‘8_views’ with optional orientation.

• hemi (str, default "lh") – Hemisphere to visualize: ‘lh’ (left) or ‘rh’ (right).

• notebook (bool, default False) – Whether to display in Jupyter notebook. If False, opens interactive window.

• show_colorbar (bool, default False) – Whether to display colorbar. Automatically determined if None.

• colorbar_title (str, optional) – Title for the colorbar. Uses overlay name if None.

• colorbar_position (str, default "bottom") – Colorbar position: ‘bottom’, ‘top’, ‘left’, or ‘right’.

• save_path (str, optional) – Path to save plot as image. If None, displays interactively.

Returns:

PyVista plotter object for further customization.

Return type:

Plotter

Raises:

ValueError – If overlay not found or invalid view parameter.

Examples

>>> surface.plot(overlay_name="aparc")
>>> surface.plot(overlay_name="thickness", cmap="hot", views="medial", show_colorbar=True)

clabtoolkit.surfacetools.merge_surfaces(surfaces, color_table=None, map_name='surf_id')

Merge multiple surface meshes into a single surface with distinct region IDs. Combines multiple surface meshes into one, assigning unique IDs to each surface region for parcellation and visualization.

Parameters:

• surfaces (List[Union[str, Path, Surface]]) – List of surface meshes to merge. Each item can be: - Surface object - File path (str or Path) to a surface file

• color_table (dict, optional) –

  Colortable dictionary defining names and colors for each surface. If None, a default colortable with distinguishable colors is created.

  The dictionary should have the following structure: {

  ”names”: List[str], # List of surface names “color_table”: np.ndarray, # Nx5 array of RGBA colors and IDs “lookup_table”: None # Placeholder for future use

  }

• map_name (str, default "surf_id") – Name of the overlay map to store surface IDs in the merged surface.

Returns:

Merged Surface object with distinct region IDs, or None if merging fails.

Return type:

Union[Surface, None]

Raises:

• TypeError – If surfaces is not a list or contains invalid items.

• ValueError – If color_table is invalid or does not match number of surfaces.

Examples

Merge multiple surface files with a custom colortable: >>> surfaces = [“lh.pial”, “rh.pial”, “lh.white”, “rh.white”] >>> color_table = { … “names”: [“Left Pial”, “Right Pial”, “Left White”, “Right White”], … “color_table”: np.array([[255, 0, 0, 255, 1], … [0, 255, 0, 255, 2], … [0, 0, 255, 255, 3], … [255, 255, 0, 255, 4]]) … } >>> merged_surface = merge_surfaces(surfaces, color_table=color_table, map_name=”region_id”) >>> print(merged_surface) Merged Surface with 4 regions and distinct IDs.

clabtoolkit.surfacetools.create_surface_colortable(colors, struct_names=None, alpha=1.0)

Create a colortable dictionary used for surface parcellation and visualization. Generates a colortable mapping region names and colors for surface parcellations.

Parameters:

• colors (str, List[str], or np.ndarray) – Colors for each region. Can be: - String: Single color name or hex code (applied to all regions) - List of strings: Color names or hex codes for each region - numpy.ndarray: Nx3 or Nx4 array of RGB(A) values (0-255)

• struct_names (List[str], optional) – List of region names corresponding to colors. If None, generic names like ‘region_1’, ‘region_2’, etc. are assigned. Length must match number of colors if provided.

• alpha (float, default 1.0) – Alpha transparency for colors (0.0 to 1.0). Applied uniformly if colors do not include alpha channel.

Returns:

Dictionary with keys: - ‘struct_names’: List of region names - ‘color_table’: Nx4 numpy array of RGBA colors (0-255) - ‘lookup_table’: None (placeholder for future use)

raises ValueError:

raises If alpha is out of range or struct_names length does not match colors.:

Examples

Create a colortable with specified colors and names: >>> colors = array([[255, 0, 0, 204],

[ 0, 255, 0, 204], [ 0, 0, 255, 204]])

>>> names = ['Region A', 'Region B', 'Region C']
>>> colortable = create_surface_colortable(colors, struct_names=names, alpha=0.8)
>>> print(colortable)
{'struct_names': ['Region A', 'Region B', 'Region C'],
'color_table': array([[255,   0,   0, 204],
                    [  0, 255,   0, 204],
                    [  0,   0, 255, 204]]),
'lookup_table': None}

Return type:

dict

The surfacetools module provides advanced brain surface mesh processing and visualization capabilities using PyVista for 3D rendering and FreeSurfer surface format support.

Key Features

• Load FreeSurfer surface files (.pial, .white, .inflated, .sphere)

• Scalar data overlay and visualization

• Parcellation integration and display

• Interactive 3D plotting with PyVista

• Surface-based analysis tools

• Publication-quality visualization

Main Classes

Surface

The primary class for surface mesh processing and visualization.

Key Methods: - load_from_file(): Load FreeSurfer surface files - load_scalar_data(): Load and attach scalar data to surface - load_parcellation(): Load parcellation overlays - plot(): Interactive 3D visualization with customization options - get_mesh_info(): Extract surface geometry information - map_volume_to_surface(): Project 3D/4D volumetric data onto surface mesh vertices

Common Usage Examples

Basic surface visualization:

from clabtoolkit.surfacetools import Surface

# Load surface
surface = Surface("/path/to/lh.pial")

# Simple surface plot
surface.plot()

# Load and visualize scalar data
surface.load_scalar_data("/path/to/lh.thickness.mgh")
surface.plot(scalar_map=True, colormap='viridis')

Advanced visualization with parcellation:

# Load surface with parcellation overlay
surface = Surface("/path/to/lh.pial")
surface.load_parcellation("/path/to/lh.aparc.annot")

# Plot with parcellation boundaries
surface.plot(
    show_parcellation=True,
    background_color='black',
    lighting=True
)

Multi-view visualization:

# Create multiple views of the same surface
surface.plot(
    views=['lateral', 'medial', 'dorsal'],
    scalar_map=True,
    save_figure="/path/to/output.png"
)

Volume-to-surface mapping:

# Project volumetric data onto surface
from pathlib import Path

# Load surface
surface = Surface("/path/to/lh.pial")

# Map 3D structural image to surface
surface_values = surface.map_volume_to_surface(
    image="/path/to/structural.nii.gz",
    method="nilearn",
    interp_method="linear"
)

# Map 4D functional data to surface
functional_values = surface.map_volume_to_surface(
    image="/path/to/functional_4d.nii.gz",
    method="custom",
    interp_method="linear",
    overlay_name="activation"
)