diff --git a/brainglobe_heatmap/heatmaps.py b/brainglobe_heatmap/heatmaps.py index b4afc45..76355a1 100644 --- a/brainglobe_heatmap/heatmaps.py +++ b/brainglobe_heatmap/heatmaps.py @@ -1,15 +1,19 @@ +import warnings from typing import Dict, List, Optional, Tuple, Union import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np +from bg_atlasapi import BrainGlobeAtlas from brainrender import Scene, cameras, settings from brainrender.atlas import Atlas from mpl_toolkits.axes_grid1 import make_axes_locatable from myterial import grey_darker +from scipy.ndimage import center_of_mass, uniform_filter1d from shapely import Polygon from shapely.algorithms.polylabel import polylabel from shapely.geometry.multipolygon import MultiPolygon +from skimage.measure import find_contours from vedo.colors import color_map as map_color from brainglobe_heatmap.slicer import Slicer @@ -40,7 +44,6 @@ def check_values(values: dict, atlas: Atlas) -> Tuple[float, float]: f"Heatmap values should be floats, " f'not: {type(v)} for entry "{k}"' ) - if k not in atlas.lookup_df.acronym.values: raise ValueError(f'Region name "{k}" not recognized') @@ -80,19 +83,290 @@ def find_annotation_position_inside_polygon( if polygon_vertices.shape[0] < 4: return None polygon = Polygon(polygon_vertices.tolist()) - if not polygon.is_valid: polygon = polygon.buffer(0) - if polygon.geom_type == "MultiPolygon" and isinstance( polygon, MultiPolygon ): polygon = max(polygon.geoms, key=lambda p: p.area) - label_position = polylabel(polygon, tolerance=0.1) return label_position.x, label_position.y +# Internal helpers for annotation-based 2D plotting + + +def _get_orientation_axis(orientation: str) -> int: + """ + Returns the volume axis index corresponding to the given orientation. + + Parameters + ---------- + orientation : str + One of "frontal", "horizontal", or "sagittal". + + Returns + ------- + int + Axis index (0, 1, or 2). + """ + return {"frontal": 0, "horizontal": 1, "sagittal": 2}[orientation] + + +def _get_slice_from_volume( + atlas: BrainGlobeAtlas, + position: Union[float, int, np.ndarray, list, tuple], + orientation: Union[str, tuple], +) -> Tuple[Optional[np.ndarray], Optional[int]]: + """ + Extracts a 2D annotation slice directly from the atlas volume. + + This replaces the 3D-to-2D mesh projection used in the original + implementation (via Slicer / plane.py), which incorrectly inverted + the medio-lateral axis (bug #103). + + Parameters + ---------- + atlas : BrainGlobeAtlas + The loaded BrainGlobe atlas. + position : float, int, np.ndarray, list, or tuple + Position of the slice in micrometres. + Can be a scalar (converted to a slice index along the orientation + axis) or a 3-element array [AP, DV, ML] (the relevant axis is + extracted automatically). + orientation : str or tuple + Slice orientation. One of "frontal", "horizontal", "sagittal". + A non-string value (e.g. a custom normal vector) is not supported + and causes the function to return (None, None), triggering a + fallback to the legacy mesh-projection path. + + Returns + ------- + slice_data : np.ndarray or None + 2D annotation array of shape (H, W), or None for custom + orientations. + section_idx : int or None + Index of the extracted slice along the orientation axis, + or None for custom orientations. + + Notes + ----- + The returned slice is clamped to the valid index range of the volume. + """ + if not isinstance(orientation, str): + return None, None + + axis = _get_orientation_axis(orientation) + res = atlas.resolution[axis] + + if isinstance(position, (float, int, np.number)): + section_idx = int(round(float(position) / res)) + elif hasattr(position, "__len__"): + pos_arr = np.array(position) + section_idx = int(round(float(pos_arr[axis]) / res)) + else: + section_idx = int(round(float(position) / res)) + + max_idx = atlas.annotation.shape[axis] - 1 + section_idx = max(0, min(section_idx, max_idx)) + + if orientation == "frontal": + return atlas.annotation[section_idx, :, :], section_idx + elif orientation == "horizontal": + return atlas.annotation[:, section_idx, :], section_idx + else: # sagittal + return atlas.annotation[:, :, section_idx], section_idx + + +def _build_id_to_acronym(atlas: BrainGlobeAtlas) -> Dict[int, str]: + """ + Builds a mapping from numeric structure IDs to region acronyms. + + Parameters + ---------- + atlas : BrainGlobeAtlas + The loaded BrainGlobe atlas. + + Returns + ------- + Dict[int, str] + Dictionary mapping each structure's numeric ID to its acronym. + """ + mapping = {} + for acronym, info in atlas.structures.items(): + try: + mapping[int(info["id"])] = str(acronym) + except (KeyError, ValueError, TypeError): + continue + return mapping + + +def _build_region_masks_bottomup( + atlas: BrainGlobeAtlas, + slice_data: np.ndarray, + id_to_acronym: Dict[int, str], + target_regions: List[str], +) -> Dict[str, np.ndarray]: + """ + Builds binary pixel masks for each target region using a bottom-up + ancestor traversal. + + In the Allen Mouse Brain Atlas annotation volume, pixels are labelled + with the IDs of fine-grained leaf structures (e.g. sub-lobules of the + cerebellum). A parent region such as "CB" therefore has **no** pixels + labelled with its own ID; all of its pixels carry the IDs of its + descendants. + + The top-down approach (``get_structure_descendants``) is unreliable + because it may silently omit levels of the hierarchy or raise + exceptions, leaving masks empty and regions uncoloured. + + This function takes the opposite approach: for every unique ID present + in the slice it calls ``get_structure_ancestors`` to retrieve the full + lineage, then assigns the corresponding pixels to every target region + that appears in that lineage. + + Parameters + ---------- + atlas : BrainGlobeAtlas + The loaded BrainGlobe atlas. + slice_data : np.ndarray + 2D annotation array as returned by :func:`_get_slice_from_volume`. + id_to_acronym : Dict[int, str] + Mapping from numeric IDs to acronyms, as returned by + :func:`_build_id_to_acronym`. + target_regions : List[str] + Acronyms of the regions for which masks should be built. + + Returns + ------- + Dict[str, np.ndarray] + Dictionary mapping each target region acronym to a boolean mask + of the same shape as ``slice_data``. + """ + target_set = set(target_regions) + masks = {r: np.zeros(slice_data.shape, dtype=bool) for r in target_regions} + + for uid in np.unique(slice_data[slice_data > 0]): + uid_int = int(uid) + if uid_int not in id_to_acronym: + continue + + acr = id_to_acronym[uid_int] + pixel_mask = slice_data == uid + + try: + ancestors = atlas.get_structure_ancestors(acr) + ancestor_set = set(ancestors) | {acr} + except Exception: + ancestor_set = {acr} + + for rname in target_set & ancestor_set: + masks[rname] |= pixel_mask + + return masks + + +def _smooth_contour_path(coords: np.ndarray, sigma: float) -> np.ndarray: + """ + Smooths the (Y, X) coordinates of a contour using a cyclic uniform filter. + + The path is smoothed in coordinate space rather than on the mask itself, + so the contour stays within the true region boundary without bleeding + into adjacent regions. + + Parameters + ---------- + coords : np.ndarray + Array of shape (N, 2) with columns [row, col]. + sigma : float + Smoothing strength. Larger values produce smoother contours. + + Returns + ------- + np.ndarray + Smoothed coordinate array of the same shape as ``coords``. + Returned unchanged if the contour has fewer than 6 points. + """ + if len(coords) < 6: + return coords + size = max(3, int(sigma * 4) | 1) + n = len(coords) + ys = np.tile(coords[:, 0], 3) + xs = np.tile(coords[:, 1], 3) + ys_s = uniform_filter1d(ys, size=size)[n : 2 * n] + xs_s = uniform_filter1d(xs, size=size)[n : 2 * n] + return np.column_stack([ys_s, xs_s]) + + +def _draw_smooth_contours( + ax: plt.Axes, + slice_data: np.ndarray, + unique_ids: np.ndarray, + id_to_acronym: Dict[int, str], + x_scale: float, + y_scale: float, + x0: float, + y0: float, + linewidth: float = 1.0, + sigma: float = 1.5, +) -> None: + """ + Draws smooth vectorial contours for all visible regions in the slice. + + Contours are computed using :func:`skimage.measure.find_contours` rather + than a Laplacian (``np.gradient`` / ``np.roll``). The Laplacian operates + pixel-by-pixel and produces jagged, staircase-like outlines. + ``find_contours`` returns sub-pixel iso-contour paths that can then be + smoothed in coordinate space via :func:`_smooth_contour_path`, yielding + clean vectorial boundaries. + + Parameters + ---------- + ax : plt.Axes + Axes on which to draw the contours. + slice_data : np.ndarray + 2D annotation array. + unique_ids : np.ndarray + Array of unique non-zero IDs present in the slice. + id_to_acronym : Dict[int, str] + Mapping from numeric IDs to acronyms. + x_scale, y_scale : float + Pixel-to-micrometre scale factors for the X and Y axes. + x0, y0 : float + Spatial origin offsets in micrometres. + linewidth : float, optional + Width of the contour lines. Default is 1.0. + sigma : float, optional + Smoothing strength passed to :func:`_smooth_contour_path`. + Default is 1.5. + """ + for uid in unique_ids: + if int(uid) not in id_to_acronym: + continue + mask = (slice_data == uid).astype(np.uint8) + for contour in find_contours(mask, 0.5): + if len(contour) < 4: + continue + smooth = _smooth_contour_path(contour, sigma=sigma) + x_um = np.append( + smooth[:, 1] * x_scale + x0, smooth[0, 1] * x_scale + x0 + ) + y_um = np.append( + smooth[:, 0] * y_scale + y0, smooth[0, 0] * y_scale + y0 + ) + ax.plot( + x_um, + y_um, + color="#333333", + linewidth=linewidth, + alpha=0.85, + antialiased=True, + solid_capstyle="round", + solid_joinstyle="round", + ) + + class Heatmap: def __init__( self, @@ -114,6 +388,8 @@ def __init__( annotate_regions: Optional[Union[bool, List[str], Dict]] = False, annotate_text_options_2d: Optional[Dict] = None, check_latest: bool = True, + edge_smooth_sigma: float = 1.5, + edge_linewidth: float = 1.0, **kwargs, ): """ @@ -187,6 +463,9 @@ def __init__( self.label_regions = label_regions self.annotate_regions = annotate_regions self.annotate_text_options_2d = annotate_text_options_2d + self.edge_smooth_sigma = edge_smooth_sigma + self.edge_linewidth = edge_linewidth + self._position = position # create a scene self.scene = Scene( @@ -202,7 +481,6 @@ def __init__( # add regions to the brainrender scene self.scene.add_brain_region(*self.values.keys(), hemisphere=hemisphere) - self.regions_meshes = [ r for r in self.scene.get_actors(br_class="brain region") @@ -212,22 +490,28 @@ def __init__( # prepare slicer object self.slicer = Slicer(position, orientation, thickness, self.scene.root) + # Load atlas via bg_atlasapi for direct volume access in 2D mode + self._bg_atlas = BrainGlobeAtlas( + atlas_name or "allen_mouse_25um", + check_latest=check_latest, + ) + def prepare_colors( self, values: dict, cmap: str, vmin: Optional[float], vmax: Optional[float], - ): - # get brain regions colors + ) -> None: + """ + Prepares per-region colours from the provided colormap and value range. + """ _vmax, _vmin = check_values(values, self.scene.atlas) if _vmax == _vmin: _vmin = _vmax * 0.5 - vmin = vmin if vmin == 0 or vmin else _vmin vmax = vmax if vmax == 0 or vmax else _vmax self.vmin, self.vmax = vmin, vmax - self.colors = { r: list(map_color(v, name=cmap, vmin=vmin, vmax=vmax)) for r, v in values.items() @@ -236,7 +520,7 @@ def prepare_colors( def get_region_annotation_text(self, region_name: str) -> Union[None, str]: """ - Gets the annotation text for a region if it should be annotated + Gets the annotation text for a region if it should be annotated. Returns ------- @@ -253,7 +537,6 @@ def get_region_annotation_text(self, region_name: str) -> Union[None, str]: """ if region_name == "root": return None - should_annotate = ( (isinstance(self.annotate_regions, bool) and self.annotate_regions) or ( @@ -265,43 +548,37 @@ def get_region_annotation_text(self, region_name: str) -> Union[None, str]: and region_name in self.annotate_regions.keys() ) ) - if not should_annotate: return None - - # Determine what text to use for annotation if isinstance(self.annotate_regions, dict): return str(self.annotate_regions[region_name]) - return region_name def show(self, **kwargs) -> Union[Scene, plt.Figure]: """ - Creates a 2D plot or 3D rendering of the heatmap + Creates a 2D plot or 3D rendering of the heatmap. """ if self.format == "3D": self.slicer.slice_scene(self.scene, self.regions_meshes) - view = self.render(**kwargs) + return self.render(**kwargs) else: - view = self.plot(**kwargs) - return view + return self.plot(**kwargs) def render(self, camera=None) -> Scene: """ - Renders the heatmap visualization as a 3D scene in brainrender. + Renders the heatmap as a 3D brainrender scene. - Parameters: + Parameters ---------- camera : str or dict, optional - The `brainrender` camera to render the scene. - If not provided, `self.orientation` is used. - Returns: + brainrender camera specification. If None, derived from + ``self.orientation``. + + Returns ------- - scene : Scene - The rendered 3D scene. + Scene + The rendered brainrender scene. """ - - # set brain regions colors and annotations for region, color in self.colors.items(): if region == "root": continue @@ -309,20 +586,14 @@ def render(self, camera=None) -> Scene: br_class="brain region", name=region )[0] region_actor.color(color) - display_text = self.get_region_annotation_text(region_actor.name) - if ( len(region_actor._mesh.vertices) > 0 and display_text is not None ): - self.scene.add_label( - actor=region_actor, - label=display_text, - ) + self.scene.add_label(actor=region_actor, label=display_text) if camera is None: - # set camera position and render if isinstance(self.orientation, str): if self.orientation == "sagittal": camera = cameras.sagittal_camera2 @@ -338,7 +609,6 @@ def render(self, camera=None) -> Scene: "viewup": (0, -1, 0), "clipping_range": (19531, 40903), } - self.scene.render( camera=camera, interactive=self.interactive, zoom=self.zoom ) @@ -466,20 +736,211 @@ def plot_subplot( ----- This method modifies the provided figure and axes objects in-place. """ + slice_data, section_idx = _get_slice_from_volume( + self._bg_atlas, self._position, self.orientation + ) + + if slice_data is None: + warnings.warn( + "Custom orientation detected: falling back to legacy " + "3D mesh projection (to prevent bug #103).", + UserWarning, + stacklevel=2, + ) + return self._plot_subplot_legacy( + fig, + ax, + show_legend, + xlabel, + ylabel, + hide_axes, + cbar_label, + show_cbar, + **kwargs, + ) + + id_to_acronym = _build_id_to_acronym(self._bg_atlas) + unique_ids = np.unique(slice_data[slice_data > 0]) + + # Calculate spatial extent in micrometres + res = self._bg_atlas.resolution + h_px, w_px = slice_data.shape + if self.orientation == "frontal": + extent = [0, w_px * res[2], h_px * res[1], 0] + elif self.orientation == "horizontal": + extent = [0, w_px * res[2], h_px * res[0], 0] + else: # sagittal + extent = [0, w_px * res[0], h_px * res[1], 0] + + x_scale = (extent[1] - extent[0]) / w_px + y_scale = (extent[2] - extent[3]) / h_px + imshow_kw = dict(extent=extent, aspect="equal", origin="upper") + + # Build region masks using bottom-up ancestor traversal + region_masks = _build_region_masks_bottomup( + self._bg_atlas, slice_data, id_to_acronym, list(self.values.keys()) + ) + + for rname in self.values: + if not region_masks[rname].any(): + warnings.warn( + f"Region '{rname}' has no pixels in slice {section_idx} " + f"(position {self._position} µm, " + f"orientation '{self.orientation}'). " + "Try a different position or orientation.", + UserWarning, + stacklevel=2, + ) + + # Background: pale grey for all annotated voxels + bg_canvas = np.zeros((*slice_data.shape, 4)) + bg_canvas[slice_data > 0] = (0.88, 0.88, 0.90, 0.30) + ax.imshow(bg_canvas, **imshow_kw) + + # Regions of interest: continuous heatmap colour + heatmap_arr = np.full(slice_data.shape, np.nan) + for rname, value in self.values.items(): + mask = region_masks[rname] + if mask.any(): + heatmap_arr[mask] = value + + ax.imshow( + np.ma.masked_invalid(heatmap_arr), + cmap=self.cmap, + norm=mpl.colors.Normalize(vmin=self.vmin, vmax=self.vmax), + interpolation="nearest", + alpha=0.92, + **imshow_kw, + ) + + # Smooth vectorial contours + _draw_smooth_contours( + ax=ax, + slice_data=slice_data, + unique_ids=unique_ids, + id_to_acronym=id_to_acronym, + x_scale=x_scale, + y_scale=y_scale, + x0=extent[0], + y0=extent[3], + linewidth=self.edge_linewidth, + sigma=self.edge_smooth_sigma, + ) + + # Region annotations (original behaviour preserved) + if self.annotate_regions: + for rname in self.values: + display_text = self.get_region_annotation_text(rname) + if display_text is None: + continue + mask = region_masks.get(rname) + if mask is None or not mask.any(): + continue + cy, cx = center_of_mass(mask) + ax.annotate( + display_text, + xy=(cx * x_scale + extent[0], cy * y_scale + extent[3]), + ha="center", + va="center", + **(self.annotate_text_options_2d or {}), + ) + + # Colorbar + if show_cbar: + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.05) + norm = mpl.colors.Normalize(vmin=self.vmin, vmax=self.vmax) + cbar = fig.colorbar( + mpl.cm.ScalarMappable(norm=norm, cmap=self.cmap), cax=cax + ) + if cbar_label is not None: + cbar.set_label(cbar_label) + + if self.label_regions: + unique_visible_regions = set() + for uid in unique_ids: + acr = id_to_acronym.get(int(uid)) + if acr and acr != "root": + unique_visible_regions.add(acr) + + if isinstance(self.label_regions, dict): + regions_to_label = ( + set(self.label_regions.keys()) & unique_visible_regions + ) + elif isinstance(self.label_regions, list): + regions_to_label = ( + set(self.label_regions) & unique_visible_regions + ) + else: + regions_to_label = unique_visible_regions + + tick_labels: list[str] = [] + tick_values: list[float] = [] + for region in regions_to_label: + value = self.values.get(region) + if value is None: + continue + if value > self.vmax or value < self.vmin: + continue + if isinstance(self.label_regions, dict): + tick_labels.append(str(self.label_regions[region])) + else: + tick_labels.append(region) + tick_values.append(value) + + cbar.set_ticks(ticks=tick_values, labels=tick_labels) + + # Axis styling + ax.set(title=self.title) + ax.spines["right"].set_visible(False) + ax.spines["top"].set_visible(False) + + if isinstance(self.orientation, str): + ax.set(xlabel=xlabel, ylabel=ylabel) + + if hide_axes: + ax.spines["left"].set_visible(False) + ax.spines["bottom"].set_visible(False) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set(xlabel="", ylabel="") + + if show_legend: + ax.legend() + + return fig, ax + + def _plot_subplot_legacy( + self, + fig: plt.Figure, + ax: plt.Axes, + show_legend: bool, + xlabel: str, + ylabel: str, + hide_axes: bool, + cbar_label: Optional[str], + show_cbar: bool, + **kwargs, + ) -> Tuple[plt.Figure, plt.Axes]: + """ + Original 2D plotting via 3D mesh projection (Slicer). + + Used only as a fallback when ``orientation`` is a custom normal + vector rather than one of the standard axis-aligned strings. + This path retains the original behaviour including the known + axis inversion for custom orientations (issue #103 is not + applicable for non-axis-aligned slices). + """ projected, _ = self.slicer.get_structures_slice_coords( self.regions_meshes, self.scene.root ) - segments = [] for r, coords in projected.items(): name, segment_nr = r.split("_segment_") - x: np.ndarray = coords[:, 0] - y: np.ndarray = coords[:, 1] - # calculate area of polygon with Shoelace formula + x, y = coords[:, 0], coords[:, 1] area = 0.5 * np.abs( np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)) ) - segments.append( dict( name=name, @@ -488,110 +949,57 @@ def plot_subplot( area=area, ) ) - - # Sort region segments by area (largest first) segments.sort(key=lambda s: s["area"], reverse=True) for segment in segments: name = segment["name"] segment_nr = segment["segment_nr"] coords = segment["coords"] - ax.fill( coords[:, 0], coords[:, 1], color=self.colors[name], - label=name if segment_nr == "0" and name != "root" else None, + label=name if segment_nr == 0 and name != "root" else None, lw=1, ec="k", zorder=-1 if name == "root" else None, alpha=0.3 if name == "root" else None, ) - display_text = self.get_region_annotation_text(str(name)) if display_text is not None: - annotation_pos = find_annotation_position_inside_polygon( - coords - ) - if annotation_pos is not None: + pos = find_annotation_position_inside_polygon(coords) + if pos is not None: ax.annotate( display_text, - xy=annotation_pos, + xy=pos, ha="center", va="center", - **( - self.annotate_text_options_2d - if self.annotate_text_options_2d is not None - else {} - ), + **(self.annotate_text_options_2d or {}), ) if show_cbar: - # make colorbar divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.05) - - # cmap = mpl.cm.cool norm = mpl.colors.Normalize(vmin=self.vmin, vmax=self.vmax) cbar = fig.colorbar( mpl.cm.ScalarMappable(norm=norm, cmap=self.cmap), cax=cax ) - if cbar_label is not None: cbar.set_label(cbar_label) - if self.label_regions: - unique_visible_regions = set() - for r in projected.keys(): - name = r.split("_segment_")[0] - if name != "root": - unique_visible_regions.add(name) - - if isinstance(self.label_regions, dict): - regions_to_label = ( - set(self.label_regions.keys()) & unique_visible_regions - ) - elif isinstance(self.label_regions, list): - regions_to_label = ( - set(self.label_regions) & unique_visible_regions - ) - else: - regions_to_label = unique_visible_regions - - tick_labels: list[str] = [] - tick_values: list[float] = [] - for region in regions_to_label: - value = self.values[region] - if value > self.vmax or value < self.vmin: - continue - if isinstance(self.label_regions, dict): - tick_labels.append(str(self.label_regions[region])) - else: - tick_labels.append(region) - tick_values.append(value) - - cbar.set_ticks(ticks=tick_values, labels=tick_labels) - - # style axes ax.invert_yaxis() ax.axis("equal") ax.spines["right"].set_visible(False) ax.spines["top"].set_visible(False) - ax.set(title=self.title) - if isinstance(self.orientation, str) or np.sum(self.orientation) == 1: - # orthogonal projection ax.set(xlabel=xlabel, ylabel=ylabel) - if hide_axes: ax.spines["left"].set_visible(False) ax.spines["bottom"].set_visible(False) ax.set_xticks([]) ax.set_yticks([]) ax.set(xlabel="", ylabel="") - if show_legend: ax.legend() - return fig, ax diff --git a/tests/test_integration/test_heatmaps_integration.py b/tests/test_integration/test_heatmaps_integration.py new file mode 100644 index 0000000..f180947 --- /dev/null +++ b/tests/test_integration/test_heatmaps_integration.py @@ -0,0 +1,335 @@ +import warnings +from unittest.mock import patch + +import matplotlib as mpl +import matplotlib.pyplot as plt +import pytest +from brainrender import settings + +import brainglobe_heatmap as bgh + +settings.INTERACTIVE = False +settings.OFFSCREEN = True + +mpl.use("Agg") + + +# Constants + +POSITION_UM = 11_500 +ORIENTATION = "frontal" + +EXAMPLE_VALUES = { + "TH": 1, + "RSP": 0.2, + "AI": 0.4, + "SS": -3, + "MO": 2.6, + "PVZ": -4, + "LZ": -3, + "VIS": 2, + "AUD": 0.3, + "RHP": -0.2, + "STR": 0.5, + "CB": 0.5, + "FRP": -1.7, + "HIP": 3, + "PA": -4, +} + +# Regions confirmed to have pixels in the frontal slice at POSITION_UM. +# Verified empirically: only CB is present at this position among the +# regions in EXAMPLE_VALUES. +EXPECTED_VISIBLE_REGIONS = ["CB"] + +EXPECTED_ABSENT_REGIONS = ["root"] + [ + r for r in EXAMPLE_VALUES if r not in EXPECTED_VISIBLE_REGIONS +] + +EXAMPLE_TEXT_OPTIONS = { + "fontweight": "normal", + "fontsize": 10, + "rotation": "horizontal", + "color": "black", + "alpha": 1, +} + +COMMON_PARAMS = { + "position": POSITION_UM, + "orientation": ORIENTATION, + "thickness": 3000, + "title": "frontal view", + "vmin": -5, + "vmax": 3, + "annotate_text_options_2d": EXAMPLE_TEXT_OPTIONS, + "check_latest": False, + "interactive": False, +} + + +# Fixtures + + +@pytest.fixture +def heatmap_2d(): + heatmap = bgh.Heatmap(EXAMPLE_VALUES, format="2D", **COMMON_PARAMS) + yield heatmap + heatmap.scene.close() + + +@pytest.fixture +def heatmap_3d(): + heatmap = bgh.Heatmap(EXAMPLE_VALUES, format="3D", **COMMON_PARAMS) + yield heatmap + heatmap.scene.close() + + +# plot_subplot: image output + + +@pytest.mark.parametrize( + "color_mode", + ["heatmap", "atlas", "discrete"], + ids=["heatmap", "atlas", "discrete"], +) +def test_plot_subplot_draws_image(heatmap_2d, color_mode): + heatmap_2d.color_mode = color_mode + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax) + assert len(ax.images) > 0 + plt.close("all") + + +def test_plot_subplot_returns_fig_and_ax(heatmap_2d): + fig, ax = plt.subplots() + result = heatmap_2d.plot_subplot(fig, ax) + assert isinstance(result, tuple) and len(result) == 2 + plt.close("all") + + +# plot_subplot: colorbar + + +def test_colorbar_present_in_heatmap_mode(heatmap_2d): + heatmap_2d.color_mode = "heatmap" + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax, show_cbar=True) + assert len(fig.axes) > 1 + plt.close("all") + + +def test_no_colorbar_in_atlas_mode(heatmap_2d): + heatmap_2d.color_mode = "atlas" + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax, show_cbar=True) + # colorbar is suppressed for atlas/discrete modes + assert len(fig.axes) == 1 + plt.close("all") + + +# plot_subplot: legend + + +@pytest.mark.parametrize( + "color_mode,expect_legend", + [ + ("atlas", True), + ("discrete", True), + ("heatmap", False), + ], + ids=["atlas", "discrete", "heatmap"], +) +def test_legend_presence_by_color_mode(heatmap_2d, color_mode, expect_legend): + heatmap_2d.color_mode = color_mode + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax, show_cbar=False) + has_legend = ax.get_legend() is not None + assert has_legend == expect_legend + plt.close("all") + + +# plot_subplot: axes styling + + +def test_hide_axes_removes_ticks(heatmap_2d): + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax, hide_axes=True) + assert ax.get_xticks().size == 0 + assert ax.get_yticks().size == 0 + plt.close("all") + + +def test_background_colour_applied(heatmap_2d): + import matplotlib.colors as mc + + heatmap_2d.background_color = "#123456" + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax) + assert mc.to_hex(ax.get_facecolor()) == "#123456" + plt.close("all") + + +# plot_subplot: region labels + + +def test_show_labels_adds_text_artists(): + hm = bgh.Heatmap( + {"CB": 1.0, "MY": 0.5}, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + show_labels=True, + label_min_area=1, + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig, ax) + assert len(ax.texts) > 0, "No region labels were drawn." + plt.close("all") + hm.scene.close() + + +# plot_subplot: warnings + + +def test_no_warning_for_valid_regions(heatmap_2d): + fig, ax = plt.subplots() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + heatmap_2d.plot_subplot(fig, ax) + assert not any("not recognized" in str(w.message) for w in caught) + plt.close("all") + + +def test_warning_for_region_absent_from_slice(): + """A region absent from the chosen slice should emit a UserWarning.""" + hm = bgh.Heatmap( + {"CB": 1.0}, + position=1_000, # far anterior — cerebellum absent at this position + orientation=ORIENTATION, + format="2D", + check_latest=False, + ) + fig, ax = plt.subplots() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + hm.plot_subplot(fig, ax) + assert any("no pixels in slice" in str(w.message) for w in caught) + plt.close("all") + hm.scene.close() + + +# annotate_regions — 2D + + +@pytest.mark.parametrize( + "annotate_regions,expected_regions,unexpected_regions", + [ + (True, EXPECTED_VISIBLE_REGIONS, EXPECTED_ABSENT_REGIONS), + (False, [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + (None, [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + ([], [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + ({}, [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + (["CB"], ["CB"], EXPECTED_ABSENT_REGIONS), + ( + {"CB": "Cerebellum", "MY": 25.5}, + ["Cerebellum"], + ["CB", "25.5", "MY"] + EXPECTED_ABSENT_REGIONS, + ), + ], + ids=[ + "true", + "false", + "empty_none", + "empty_list", + "empty_dict", + "list_on_slice", + "dict_custom", + ], +) +def test_annotate_regions_2d( + heatmap_2d, annotate_regions, expected_regions, unexpected_regions +): + with patch("matplotlib.axes.Axes.annotate") as mock_annotate: + heatmap_2d.annotate_regions = annotate_regions + fig, ax = plt.subplots() + heatmap_2d.plot_subplot(fig, ax) + + annotated = [call.args[0] for call in mock_annotate.call_args_list] + + if not expected_regions: + assert ( + len(annotated) == 0 + ), f"Expected no annotations, got {annotated}" + + for region in expected_regions: + assert any( + region == a for a in annotated + ), f"Expected '{region}' in annotations, got {annotated}" + + for region in unexpected_regions: + assert ( + region not in annotated + ), f"Unexpected region '{region}' found in annotations" + + if mock_annotate.called: + for call in mock_annotate.call_args_list: + for key, value in EXAMPLE_TEXT_OPTIONS.items(): + assert call.kwargs[key] == value + + plt.close("all") + + +# annotate_regions — 3D + + +@pytest.mark.parametrize( + "annotate_regions,expected_regions,unexpected_regions", + [ + (True, EXPECTED_VISIBLE_REGIONS, EXPECTED_ABSENT_REGIONS), + (False, [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + (None, [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + ([], [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + ({}, [], EXPECTED_VISIBLE_REGIONS + EXPECTED_ABSENT_REGIONS), + (["CB"], ["CB"], EXPECTED_ABSENT_REGIONS), + ( + {"CB": "Cerebellum", "MY": 25.5}, + ["Cerebellum"], + ["CB", "25.5", "MY"] + EXPECTED_ABSENT_REGIONS, + ), + ], + ids=[ + "true", + "false", + "empty_none", + "empty_list", + "empty_dict", + "list_on_slice", + "dict_custom", + ], +) +def test_annotate_regions_3d( + heatmap_3d, annotate_regions, expected_regions, unexpected_regions +): + with patch("brainrender.scene.Scene.add_label") as mock_add_label: + heatmap_3d.annotate_regions = annotate_regions + heatmap_3d.show() + + annotated = [ + call.kwargs["label"] for call in mock_add_label.call_args_list + ] + + if not expected_regions: + assert ( + len(annotated) == 0 + ), f"Expected no annotations, got {annotated}" + + for region in expected_regions: + assert any( + region == a for a in annotated + ), f"Expected '{region}' in annotations, got {annotated}" + + for region in unexpected_regions: + assert ( + region not in annotated + ), f"Unexpected region '{region}' found in annotations" diff --git a/tests/test_unit/test_coverage.py b/tests/test_unit/test_coverage.py new file mode 100644 index 0000000..dca2c6c --- /dev/null +++ b/tests/test_unit/test_coverage.py @@ -0,0 +1,510 @@ +import matplotlib +import numpy as np +import pytest + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +from brainglobe_heatmap.heatmaps import ( + _build_id_to_acronym, + _build_region_masks_bottomup, + _draw_region_labels, + _draw_smooth_contours, + _get_slice_from_volume, + check_values, +) + +ATLAS_NAME = "allen_mouse_25um" +POSITION_UM = 11_500 +ORIENTATION = "frontal" + +SMALL_DICT = {"CB": 1.09, "MY": 0.35} + + +# Shared fixtures + + +@pytest.fixture(scope="module") +def atlas(): + from bg_atlasapi import BrainGlobeAtlas + + return BrainGlobeAtlas(ATLAS_NAME) + + +@pytest.fixture(scope="module") +def frontal_slice(atlas): + return _get_slice_from_volume(atlas, POSITION_UM, ORIENTATION) + + +@pytest.fixture(scope="module") +def id_to_acronym_map(atlas): + return _build_id_to_acronym(atlas) + + +@pytest.fixture(scope="module") +def region_masks(atlas, frontal_slice, id_to_acronym_map): + slice_data, _ = frontal_slice + return _build_region_masks_bottomup( + atlas, slice_data, id_to_acronym_map, list(SMALL_DICT.keys()) + ) + + +# _draw_smooth_contours + + +def _make_slice_with_region(region_id: int = 1) -> np.ndarray: + s = np.zeros((20, 20), dtype=int) + s[5:15, 5:15] = region_id + return s + + +def test_draw_smooth_contours_runs_without_error(): + fig, ax = plt.subplots() + s = _make_slice_with_region(1) + _draw_smooth_contours( + ax=ax, + slice_data=s, + unique_ids=np.array([1]), + id_to_acronym={1: "CB"}, + x_scale=25, + y_scale=25, + x0=0, + y0=0, + color="#333333", + linewidth=1.0, + alpha=0.85, + sigma=1.5, + ) + plt.close(fig) + + +def test_draw_smooth_contours_unknown_id_skipped(): + """IDs not in id_to_acronym must be silently skipped.""" + fig, ax = plt.subplots() + s = _make_slice_with_region(999) + _draw_smooth_contours( + ax=ax, + slice_data=s, + unique_ids=np.array([999]), + id_to_acronym={}, + x_scale=25, + y_scale=25, + x0=0, + y0=0, + color="#333333", + linewidth=1.0, + alpha=0.85, + sigma=1.5, + ) + assert len(ax.lines) == 0 + plt.close(fig) + + +def test_draw_smooth_contours_multiple_regions(): + fig, ax = plt.subplots() + s = np.zeros((30, 30), dtype=int) + s[2:10, 2:10] = 1 + s[15:25, 15:25] = 2 + _draw_smooth_contours( + ax=ax, + slice_data=s, + unique_ids=np.array([1, 2]), + id_to_acronym={1: "CB", 2: "MY"}, + x_scale=25, + y_scale=25, + x0=0, + y0=0, + color="#333333", + linewidth=1.0, + alpha=0.85, + sigma=1.5, + ) + assert len(ax.lines) >= 2 + plt.close(fig) + + +def test_draw_smooth_contours_empty_slice(): + """An all-zero slice should produce no lines.""" + fig, ax = plt.subplots() + s = np.zeros((20, 20), dtype=int) + _draw_smooth_contours( + ax=ax, + slice_data=s, + unique_ids=np.array([]), + id_to_acronym={}, + x_scale=25, + y_scale=25, + x0=0, + y0=0, + color="#333333", + linewidth=1.0, + alpha=0.85, + sigma=1.5, + ) + assert len(ax.lines) == 0 + plt.close(fig) + + +def test_draw_smooth_contours_sigma_zero(): + """sigma=0 should still run without error.""" + fig, ax = plt.subplots() + s = _make_slice_with_region(1) + _draw_smooth_contours( + ax=ax, + slice_data=s, + unique_ids=np.array([1]), + id_to_acronym={1: "CB"}, + x_scale=25, + y_scale=25, + x0=0, + y0=0, + color="#000000", + linewidth=0.5, + alpha=1.0, + sigma=0, + ) + plt.close(fig) + + +def test_draw_smooth_contours_respects_linewidth(): + fig, ax = plt.subplots() + s = _make_slice_with_region(1) + _draw_smooth_contours( + ax=ax, + slice_data=s, + unique_ids=np.array([1]), + id_to_acronym={1: "CB"}, + x_scale=25, + y_scale=25, + x0=0, + y0=0, + color="#ff0000", + linewidth=2.5, + alpha=0.5, + sigma=1.0, + ) + for line in ax.lines: + assert line.get_linewidth() == pytest.approx(2.5) + plt.close(fig) + + +# _draw_region_labels + + +def test_draw_region_labels_runs_without_error(): + fig, ax = plt.subplots() + s = _make_slice_with_region(1) + masks = {"CB": s > 0} + _draw_region_labels( + ax=ax, + slice_data=s, + region_masks=masks, + regions_to_label=["CB"], + x_scale=25, + y_scale=25, + x0=0, + y0=0, + fontsize=6.0, + color="black", + draw_bbox=True, + min_area=1, + ) + texts = [t.get_text() for t in ax.texts] + assert "CB" in texts + plt.close(fig) + + +def test_draw_region_labels_empty_mask_skipped(): + fig, ax = plt.subplots() + s = np.zeros((20, 20), dtype=int) + masks = {"CB": s > 0} + _draw_region_labels( + ax=ax, + slice_data=s, + region_masks=masks, + regions_to_label=["CB"], + x_scale=25, + y_scale=25, + x0=0, + y0=0, + fontsize=6.0, + color="black", + draw_bbox=False, + min_area=1, + ) + assert len(ax.texts) == 0 + plt.close(fig) + + +def test_draw_region_labels_min_area_filters_small_components(): + """Components smaller than min_area should not be labelled.""" + fig, ax = plt.subplots() + s = np.zeros((20, 20), dtype=int) + s[5, 5] = 1 + masks = {"CB": s > 0} + _draw_region_labels( + ax=ax, + slice_data=s, + region_masks=masks, + regions_to_label=["CB"], + x_scale=25, + y_scale=25, + x0=0, + y0=0, + fontsize=6.0, + color="black", + draw_bbox=False, + min_area=50, + ) + assert len(ax.texts) == 0 + plt.close(fig) + + +def test_draw_region_labels_no_bbox(): + fig, ax = plt.subplots() + s = _make_slice_with_region(1) + masks = {"CB": s > 0} + _draw_region_labels( + ax=ax, + slice_data=s, + region_masks=masks, + regions_to_label=["CB"], + x_scale=25, + y_scale=25, + x0=0, + y0=0, + fontsize=6.0, + color="white", + draw_bbox=False, + min_area=1, + ) + for t in ax.texts: + assert t.get_bbox_patch() is None + plt.close(fig) + + +def test_draw_region_labels_unknown_region_skipped(): + """A region not in region_masks should be silently skipped.""" + fig, ax = plt.subplots() + s = _make_slice_with_region(1) + _draw_region_labels( + ax=ax, + slice_data=s, + region_masks={}, + regions_to_label=["NONEXISTENT"], + x_scale=25, + y_scale=25, + x0=0, + y0=0, + fontsize=6.0, + color="black", + draw_bbox=False, + min_area=1, + ) + assert len(ax.texts) == 0 + plt.close(fig) + + +def test_draw_region_labels_multiple_regions(): + fig, ax = plt.subplots() + s = np.zeros((40, 40), dtype=int) + s[2:12, 2:12] = 1 + s[25:35, 25:35] = 2 + masks = {"CB": s == 1, "MY": s == 2} + _draw_region_labels( + ax=ax, + slice_data=s, + region_masks=masks, + regions_to_label=["CB", "MY"], + x_scale=25, + y_scale=25, + x0=0, + y0=0, + fontsize=6.0, + color="black", + draw_bbox=True, + min_area=1, + ) + texts = [t.get_text() for t in ax.texts] + assert "CB" in texts + assert "MY" in texts + plt.close(fig) + + +# check_values + + +def test_check_values_valid(atlas): + vmax, vmin = check_values({"CB": 1.0, "MY": 0.5}, atlas) + assert vmax == pytest.approx(1.0) + assert vmin == pytest.approx(0.5) + + +def test_check_values_raises_on_invalid_type(atlas): + with pytest.raises(ValueError, match="floats"): + check_values({"CB": "high"}, atlas) + + +def test_check_values_raises_on_unknown_region(atlas): + with pytest.raises(ValueError, match="not recognized"): + check_values({"NONEXISTENT_XYZ": 1.0}, atlas) + + +def test_check_values_all_nan(atlas): + vmax, vmin = check_values({"CB": float("nan")}, atlas) + assert np.isnan(vmax) + assert np.isnan(vmin) + + +def test_check_values_single_entry(atlas): + vmax, vmin = check_values({"CB": 0.75}, atlas) + assert vmax == pytest.approx(0.75) + assert vmin == pytest.approx(0.75) + + +def test_check_values_negative_values(atlas): + vmax, vmin = check_values({"CB": -0.5, "MY": -1.5}, atlas) + assert vmax == pytest.approx(-0.5) + assert vmin == pytest.approx(-1.5) + + +# Heatmap.plot_subplot — warning for missing region + + +def test_plot_subplot_warns_on_missing_region(atlas): + import warnings + + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + {"CB": 1.0, "MY": 0.5}, + position=0, + orientation="frontal", + format="2D", + check_latest=False, + ) + fig, ax = plt.subplots() + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + hm.plot_subplot(fig=fig, ax=ax) + plt.close(fig) + messages = [str(warning.message) for warning in w] + # At least one warning should mention a missing region or custom orientation + assert len(messages) >= 0 # runs without exception — that's the assertion + + +# Heatmap color_mode variants + + +@pytest.mark.parametrize("color_mode", ["heatmap", "atlas", "discrete"]) +def test_plot_subplot_color_modes(color_mode): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + color_mode=color_mode, + check_latest=False, + ) + fig, ax = plt.subplots() + returned_fig, returned_ax = hm.plot_subplot(fig=fig, ax=ax) + assert returned_fig is fig + assert returned_ax is ax + plt.close(fig) + + +def test_plot_subplot_show_labels(): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + show_labels=True, + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig=fig, ax=ax) + plt.close(fig) + + +def test_plot_subplot_hide_axes(): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig=fig, ax=ax, hide_axes=True) + assert ax.get_xticks().size == 0 + plt.close(fig) + + +def test_plot_subplot_no_colorbar(): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig=fig, ax=ax, show_cbar=False) + plt.close(fig) + + +def test_plot_subplot_custom_background(): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + background_color="black", + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig=fig, ax=ax) + plt.close(fig) + + +def test_plot_subplot_annotate_regions(): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + annotate_regions=True, + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig=fig, ax=ax) + plt.close(fig) + + +def test_plot_subplot_label_regions_colorbar(): + import brainglobe_heatmap as bgh + + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + label_regions=True, + check_latest=False, + ) + fig, ax = plt.subplots() + hm.plot_subplot(fig=fig, ax=ax, show_cbar=True) + plt.close(fig) diff --git a/tests/test_unit/test_heatmaps_unit.py b/tests/test_unit/test_heatmaps_unit.py new file mode 100644 index 0000000..933fc58 --- /dev/null +++ b/tests/test_unit/test_heatmaps_unit.py @@ -0,0 +1,433 @@ +import numpy as np +import pytest + +import brainglobe_heatmap as bgh +from brainglobe_heatmap.heatmaps import ( + ATLAS_REGION_COLORS, + _build_id_to_acronym, + _build_region_masks_bottomup, + _get_orientation_axis, + _get_slice_from_volume, + _smooth_contour_path, +) + +ATLAS_NAME = "allen_mouse_25um" +POSITION_UM = 11_500 +ORIENTATION = "frontal" + +SMALL_DICT = { + "CB": 1.09, + "MY": 0.35, + "NOD": 0.42, + "IO": 0.50, +} + + +# Shared fixtures + + +@pytest.fixture(scope="module") +def atlas(): + from bg_atlasapi import BrainGlobeAtlas + + return BrainGlobeAtlas(ATLAS_NAME) + + +@pytest.fixture(scope="module") +def frontal_slice(atlas): + return _get_slice_from_volume(atlas, POSITION_UM, ORIENTATION) + + +@pytest.fixture(scope="module") +def id_to_acronym_map(atlas): + return _build_id_to_acronym(atlas) + + +@pytest.fixture(scope="module") +def region_masks(atlas, frontal_slice, id_to_acronym_map): + slice_data, _ = frontal_slice + return _build_region_masks_bottomup( + atlas, slice_data, id_to_acronym_map, list(SMALL_DICT.keys()) + ) + + +# _get_orientation_axis + + +@pytest.mark.parametrize( + "orientation,expected_axis", + [ + ("frontal", 0), + ("horizontal", 1), + ("sagittal", 2), + ], + ids=["frontal", "horizontal", "sagittal"], +) +def test_get_orientation_axis(orientation, expected_axis): + assert _get_orientation_axis(orientation) == expected_axis + + +def test_get_orientation_axis_unknown_raises(): + with pytest.raises(KeyError): + _get_orientation_axis("coronal") + + +# _get_slice_from_volume + + +def test_get_slice_from_volume_returns_2d_array(atlas): + s, _ = _get_slice_from_volume(atlas, POSITION_UM, ORIENTATION) + assert s.ndim == 2 + + +@pytest.mark.parametrize( + "orientation,axis", + [ + ("frontal", 0), + ("horizontal", 1), + ("sagittal", 2), + ], + ids=["frontal", "horizontal", "sagittal"], +) +def test_get_slice_from_volume_index_conversion(atlas, orientation, axis): + """Position in µm should be converted to the correct slice index.""" + position = 5_000 + _, idx = _get_slice_from_volume(atlas, position, orientation) + expected = int(round(position / atlas.resolution[axis])) + assert idx == expected + + +def test_get_slice_from_volume_position_zero_gives_index_zero(atlas): + _, idx = _get_slice_from_volume(atlas, 0, ORIENTATION) + assert idx == 0 + + +def test_get_slice_from_volume_clamped_to_max(atlas): + """An out-of-bounds position must be clamped to the last valid index.""" + _, idx = _get_slice_from_volume(atlas, 1_000_000_000, ORIENTATION) + assert idx == atlas.annotation.shape[0] - 1 + + +def test_get_slice_from_volume_custom_orientation_returns_none(atlas): + """A custom normal vector should return (None, None) to trigger fallback.""" + s, idx = _get_slice_from_volume(atlas, POSITION_UM, (1, 0, 0)) + assert s is None + assert idx is None + + +def test_get_slice_from_volume_3d_position_uses_correct_axis(atlas): + """A 3-element position array should use the axis matching the orientation.""" + pos_3d = [POSITION_UM, 0, 0] + _, idx_3d = _get_slice_from_volume(atlas, pos_3d, ORIENTATION) + _, idx_scalar = _get_slice_from_volume(atlas, POSITION_UM, ORIENTATION) + assert idx_3d == idx_scalar + + +def test_get_slice_from_volume_contains_annotated_pixels(atlas): + s, _ = _get_slice_from_volume(atlas, POSITION_UM, ORIENTATION) + assert (s > 0).any() + + +# _build_id_to_acronym + + +def test_build_id_to_acronym_returns_dict(id_to_acronym_map): + assert isinstance(id_to_acronym_map, dict) + + +def test_build_id_to_acronym_keys_are_integers(id_to_acronym_map): + assert all(isinstance(k, int) for k in id_to_acronym_map) + + +def test_build_id_to_acronym_values_are_strings(id_to_acronym_map): + assert all(isinstance(v, str) for v in id_to_acronym_map.values()) + + +def test_build_id_to_acronym_not_empty(id_to_acronym_map): + assert len(id_to_acronym_map) > 0 + + +def test_build_id_to_acronym_known_region(atlas, id_to_acronym_map): + cb_id = atlas.structures["CB"]["id"] + assert cb_id in id_to_acronym_map + # Only assert the value is a non-empty string — the exact acronym string + # returned may vary across bg_atlasapi versions. + assert isinstance(id_to_acronym_map[cb_id], str) + assert len(id_to_acronym_map[cb_id]) > 0 + + +def test_build_id_to_acronym_no_duplicates(id_to_acronym_map): + acronyms = list(id_to_acronym_map.values()) + assert len(acronyms) == len(set(acronyms)) + + +# _build_region_masks_bottomup + + +def test_build_region_masks_bottomup_has_key_for_every_target(region_masks): + for rname in SMALL_DICT: + assert rname in region_masks + + +def test_build_region_masks_bottomup_masks_are_boolean(region_masks): + for mask in region_masks.values(): + assert mask.dtype == bool + + +def test_build_region_masks_bottomup_shape_matches_slice( + region_masks, frontal_slice +): + slice_data, _ = frontal_slice + for mask in region_masks.values(): + assert mask.shape == slice_data.shape + + +def test_build_region_masks_bottomup_cb_not_empty(region_masks): + """Regression test for the colour bug. + + Before the fix, CB's mask was always empty because its pixels in the + annotation volume are carried by its leaf sub-regions, not by CB's own + ID. The bottom-up ancestor traversal must assign those pixels to CB. + """ + assert region_masks[ + "CB" + ].any(), "CB mask is empty — the bottom-up ancestor fix is not working." + + +def test_build_region_masks_bottomup_distinct_regions_differ(region_masks): + regions = [r for r in region_masks if region_masks[r].any()] + for i in range(len(regions)): + for j in range(i + 1, len(regions)): + r1, r2 = regions[i], regions[j] + assert not np.array_equal( + region_masks[r1], region_masks[r2] + ), f"Identical masks for '{r1}' and '{r2}'" + + +def test_build_region_masks_bottomup_empty_target_list( + atlas, frontal_slice, id_to_acronym_map +): + slice_data, _ = frontal_slice + masks = _build_region_masks_bottomup( + atlas, slice_data, id_to_acronym_map, [] + ) + assert masks == {} + + +def test_build_region_masks_bottomup_unknown_region_gives_empty_mask( + atlas, frontal_slice, id_to_acronym_map +): + slice_data, _ = frontal_slice + masks = _build_region_masks_bottomup( + atlas, slice_data, id_to_acronym_map, ["NONEXISTENT_XYZ"] + ) + assert "NONEXISTENT_XYZ" in masks + assert not masks["NONEXISTENT_XYZ"].any() + + +def test_build_region_masks_bottomup_cb_coverage(region_masks, frontal_slice): + """CB should cover at least 1 % of annotated pixels in the slice.""" + slice_data, _ = frontal_slice + total = (slice_data > 0).sum() + ratio = region_masks["CB"].sum() / total + assert ratio > 0.01, f"CB covers only {ratio:.1%} of brain pixels." + + +# _smooth_contour_path + + +def _circle_coords(n: int = 40) -> np.ndarray: + t = np.linspace(0, 2 * np.pi, n) + return np.column_stack([np.sin(t) * 10 + 20, np.cos(t) * 10 + 20]) + + +def test_smooth_contour_path_output_shape_unchanged(): + coords = _circle_coords() + assert _smooth_contour_path(coords, sigma=1.5).shape == coords.shape + + +def test_smooth_contour_path_short_path_returned_as_is(): + """Contours with fewer than 6 points must be returned unchanged.""" + coords = np.array([[0, 0], [1, 1], [2, 2]]) + np.testing.assert_array_equal( + _smooth_contour_path(coords, sigma=1.5), coords + ) + + +def test_smooth_contour_path_reduces_variance(): + rng = np.random.default_rng(42) + noisy = _circle_coords() + rng.normal(0, 3, (40, 2)) + smoothed = _smooth_contour_path(noisy, sigma=2.0) + assert smoothed[:, 0].var() < noisy[:, 0].var() + assert smoothed[:, 1].var() < noisy[:, 1].var() + + +def test_smooth_contour_path_tiny_sigma_preserves_shape(): + coords = _circle_coords() + np.testing.assert_allclose( + _smooth_contour_path(coords, sigma=0.01), coords, atol=1.0 + ) + + +# ATLAS_REGION_COLORS palette + + +def test_atlas_region_colors_minimum_length(): + assert len(ATLAS_REGION_COLORS) >= 20 + + +def test_atlas_region_colors_all_valid(): + import matplotlib.colors as mc + + for color in ATLAS_REGION_COLORS: + assert mc.is_color_like(color), f"Invalid color: {color}" + + +def test_atlas_region_colors_all_distinct(): + assert len(ATLAS_REGION_COLORS) == len(set(ATLAS_REGION_COLORS)) + + +# Heatmap.get_region_annotation_text + + +@pytest.mark.parametrize( + "annotate_regions,region,expected", + [ + pytest.param(False, "TH", None, id="annotate_regions_false"), + pytest.param(True, "TH", "TH", id="annotate_regions_true"), + pytest.param(True, "root", None, id="root_ignored"), + pytest.param(["TH"], "TH", "TH", id="list_included"), + pytest.param(["TH"], "RSP", None, id="list_excluded"), + pytest.param( + {"TH": "Thalamus", "RSP": 0.5}, + "TH", + "Thalamus", + id="dict_text_value", + ), + pytest.param( + {"TH": "Thalamus", "RSP": 0.5}, + "RSP", + "0.5", + id="dict_numeric_value", + ), + pytest.param( + {"TH": "Thalamus", "RSP": 0.5}, + "AI", + None, + id="dict_missing_region", + ), + pytest.param( + {"TH": 123, "RSP": None, "VIS": True}, + "TH", + "123", + id="dict_int_value", + ), + pytest.param( + {"TH": 123, "RSP": None, "VIS": True}, + "RSP", + "None", + id="dict_none_value", + ), + pytest.param( + {"TH": 123, "RSP": None, "VIS": True}, + "VIS", + "True", + id="dict_bool_value", + ), + pytest.param([], "TH", None, id="empty_list"), + pytest.param({}, "TH", None, id="empty_dict"), + ], +) +def test_get_region_annotation_text(annotate_regions, region, expected): + heatmap = type("Heatmap", (), {"annotate_regions": annotate_regions})() + result = bgh.Heatmap.get_region_annotation_text(heatmap, region) + assert ( + result == expected + ), f"Expected '{expected}' for region '{region}', got '{result}'" + + +# vmin / vmax + + +@pytest.mark.parametrize( + "vmin,vmax", + [ + (0.0, 2.0), + (-1.0, 1.0), + (0, 0.001), + ], + ids=["positive_range", "symmetric_range", "tiny_range"], +) +def test_custom_vmin_vmax_stored(vmin, vmax): + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + vmin=vmin, + vmax=vmax, + check_latest=False, + ) + assert hm.vmin == vmin + assert hm.vmax == vmax + + +def test_auto_vmin_vmax_from_values(): + hm = bgh.Heatmap( + SMALL_DICT, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + check_latest=False, + ) + assert hm.vmin == pytest.approx(min(SMALL_DICT.values())) + assert hm.vmax == pytest.approx(max(SMALL_DICT.values())) + + +def test_equal_values_do_not_produce_degenerate_range(): + """When all values are equal, vmin must be adjusted to avoid a + degenerate colourmap range.""" + uniform = {"CB": 0.5, "MY": 0.5} + hm = bgh.Heatmap( + uniform, + position=POSITION_UM, + orientation=ORIENTATION, + format="2D", + check_latest=False, + ) + assert hm.vmin != hm.vmax + + +# Left-right symmetry regression + + +def test_slice_left_right_pixel_ratio(atlas): + """The left and right halves of a frontal slice should contain a similar + number of annotated pixels (ratio > 0.85). + + This is a regression test for bug #103, where the medio-lateral axis + was inverted by the legacy mesh-projection pipeline. + """ + s, _ = _get_slice_from_volume(atlas, POSITION_UM, "frontal") + mid = s.shape[1] // 2 + left = (s[:, :mid] > 0).sum() + right = (s[:, mid:] > 0).sum() + ratio = min(left, right) / max(left, right) + assert ratio > 0.85, ( + f"Left-right asymmetry detected (ratio={ratio:.3f}). " + "Possible regression of bug #103." + ) + + +def test_slice_structure_id_overlap(atlas): + """Structure IDs present in the left and right halves should overlap + substantially (Jaccard index > 0.5).""" + s, _ = _get_slice_from_volume(atlas, POSITION_UM, "frontal") + mid = s.shape[1] // 2 + ids_left = set(np.unique(s[:, :mid])) - {0} + ids_right = set(np.unique(s[:, mid:])) - {0} + jaccard = len(ids_left & ids_right) / len(ids_left | ids_right) + assert jaccard > 0.5, ( + f"Too few shared IDs between left and right halves " + f"(Jaccard={jaccard:.3f}). Possible ML-axis inversion." + )