diff --git a/README.md b/README.md index 5de3cc3..c09df79 100644 --- a/README.md +++ b/README.md @@ -93,6 +93,40 @@ ggseg.plot_jhu(data, background='k', edgecolor='w', cmap='Spectral', The comprehensive list of applicable regions can be found in this [folder](https://github.com/ggseg/python-ggseg/tree/main/ggseg/data/jhu). + + +### R ggseg polygon atlas interoperability + +`python-ggseg` can also render 2D polygon atlas tables exported from R `ggseg`. +This is useful when strict R/Python visual parity is needed: both backends can +render the same atlas coordinates, and an optional `fill_hex` column can be used +as identity colors for exact cross-backend color matching. + +```python +import ggseg + +atlas = ggseg.read_polygon_atlas("dk_polygons.csv") + +# Join region-level values by the R ggseg label column. +plot_data = ggseg.join_polygon_values( + atlas, + {"bankssts_left": 1.1, "bankssts_right": 1.4}, + match_column="label", +) + +# Render by numeric value. +ggseg.plot_polygons(plot_data, value_column="value", show=False) + +# Or render exact precomputed colors, if the table contains fill_hex. +ggseg.plot_polygons(plot_data, fill_column="fill_hex", show=False) +``` + +Expected polygon atlas columns are: + +```text +label, region, hemi, view, x, y, .group, subgroup, .feature_id +``` + ## Tests The current development version of `python-ggseg` has a coverage rate close to 100%. diff --git a/doc/r_ggseg_polygon_interop_proposal.md b/doc/r_ggseg_polygon_interop_proposal.md new file mode 100644 index 0000000..3599ffe --- /dev/null +++ b/doc/r_ggseg_polygon_interop_proposal.md @@ -0,0 +1,264 @@ +# Proposal: interoperable R ggseg polygon atlas support for python-ggseg + +## Executive summary + +This document describes and accompanies a small interoperability layer for `python-ggseg`: the ability to load, validate, join, and render cortical atlas polygon tables exported from R `ggseg`. The goal is not to replace the native Python atlas workflow, but to provide a reproducible bridge between the larger R `ggseg` atlas ecosystem and Python-based neuroimaging pipelines. + +The key idea is simple: use a shared polygon atlas file as the source of truth for geometry, and optionally use a precomputed `fill_hex` column as the source of truth for colors. This makes Python and R outputs comparable at the data and figure level. + +```text +R ggseg atlas object + -> canonical polygon CSV / JSON + -> python-ggseg render + -> R ggseg / ggplot render +``` + + + +## Implemented minimal API in this PR + +This PR includes a first minimal implementation of the proposed interoperability path: + +- `read_polygon_atlas(path)` reads a canonical R `ggseg` polygon CSV and converts `x`/`y` to numeric coordinates. +- `validate_polygon_atlas(atlas)` checks the required schema. +- `join_polygon_values(atlas, data, match_column="label", value_column="value")` joins region-level values by an explicit key. +- `plot_polygons(...)` renders the shared polygon table in Matplotlib. +- `fill_column="fill_hex"` uses identity colors and avoids backend-specific colormap remapping. + +The implementation is intentionally additive and does not change existing `plot_dk`, `plot_aseg`, or `plot_jhu` behavior. + +## Background + +R `ggseg` and `ggseg3d` provide a mature ecosystem for visualizing region-level brain statistics on predefined atlas segmentations. Many neuroimaging workflows, however, are Python-first. At present, Python and R ggseg-style plots may differ because each backend can use different atlas objects, geometry preparation steps, plotting defaults, or color-mapping logic. + +For mixed Python/R projects, this creates a reproducibility gap: even when the same values and nominal atlas are used, the final figures may not be guaranteed to share the same geometry and color assignments. + +## Problem statement + +A Python user should be able to render an atlas exported from R `ggseg` without manually rewriting polygon handling code or accepting silent geometry differences. A project that supports both Python and R should be able to produce matched figures from the same atlas geometry and, when required, the same assigned colors. + +## Goals + +- Support R `ggseg`-exported polygon atlas files as a first-class input path in `python-ggseg`. +- Provide a documented canonical schema for 2D polygon atlas interchange. +- Allow users to join region-level data by an explicit key such as `label` or `region`. +- Allow exact cross-backend color parity when `fill_hex` is supplied. +- Preserve raw atlas topology by default so adjacent parcel boundaries continue to fit. +- Add validation and tests for schema, joins, missing regions, and deterministic rendering. + +## Non-goals + +- Do not replace native `python-ggseg` atlas APIs. +- Do not require Python users to install R for normal plotting once an atlas CSV/JSON has been exported. +- Do not make independent per-parcel smoothing a default rendering step. +- Do not attempt to solve 3D mesh interoperability; this proposal is for 2D polygon atlases. + +## Proposed canonical polygon schema + +A minimal 2D atlas polygon table should include: + +| Column | Type | Required | Description | +|---|---:|---:|---| +| `label` | string | yes | Stable parcel label; preferred join key. | +| `region` | string | yes | Human-readable region name. | +| `hemi` | string | yes | Hemisphere, e.g. `left`, `right`. | +| `view` | string | yes | Display view, e.g. `lateral`, `medial`. | +| `x` | numeric | yes | Positioned 2D x coordinate. | +| `y` | numeric | yes | Positioned 2D y coordinate. | +| `.group` | string/int | yes | Polygon group identifier from the ggseg-prepared object. | +| `subgroup` | string/int | yes | Ring/subpolygon identifier for multipart features. | +| `.feature_id` | string/int | yes | Unique displayed polygon feature identifier. | + +Optional columns: + +| Column | Type | Description | +|---|---:|---| +| `value` | numeric | Region-level statistic after joining user data. | +| `fill_hex` | string | Precomputed color, e.g. `#3B4CC0`; if present, plotting should not remap colors unless explicitly requested. | +| `atlas` | string | Atlas name, e.g. `dk`. | +| `source_package` | string | Source R package, e.g. `ggseg`. | +| `source_version` | string | Source package version. | +| `export_date` | string | Export date or timestamp. | + +## Proposed Python API + +A minimal API could look like this: + +```python +import ggseg + +atlas = ggseg.read_polygon_atlas("dk_polygons.csv") + +fig = ggseg.plot_polygons( + atlas, + data=values_df, + match_column="label", + value_column="value", + palette="diverging", + vmin=-1, + midpoint=0, + vmax=1, + edgecolor="#2E2E2E", + linewidth=0.55, +) +fig.savefig("dk_python.svg") +``` + +For exact color parity with R or another backend: + +```python +plot_data = ggseg.join_values( + atlas, + values_df, + match_column="label", + value_column="value", +) +plot_data = ggseg.assign_fill_hex( + plot_data, + palette="diverging", + vmin=-1, + midpoint=0, + vmax=1, +) + +ggseg.plot_polygons(plot_data, fill_column="fill_hex") +``` + +If a table already contains `fill_hex`, the renderer should treat it as an identity color column by default. + +## Suggested R export recipe + +This could be implemented in documentation, a helper script, or an R-side convenience function: + +```r +export_ggseg_polygon_atlas <- function(atlas, file) { + flat <- ggseg:::prepare_polygon_atlas( + atlas = atlas, + hemi = c("left", "right"), + view = c("lateral", "medial"), + position = ggseg::position_brain(hemi ~ view), + context = TRUE, + focus = NULL + ) + + keep <- c("label", "region", "hemi", "view", "x", "y", ".group", "subgroup", ".feature_id") + utils::write.csv(flat[, keep], file, row.names = FALSE) +} + +export_ggseg_polygon_atlas(ggseg::dk(), "dk_polygons.csv") +``` + +The exact implementation may need to follow `ggseg` maintainers' preferred public API boundaries. If relying on an internal function is not desirable, this proposal can instead document a supported export path or request a small R-side public export helper. + +## Boundary and smoothing policy + +Raw ggseg cortex polygons can look angular. In local experiments, a one-pass Chaikin display asset produced smoother cortical outlines and better visual fit than a simplify + Catmull-Rom experiment. However, smoothing is a display-quality decision and should not be silently mixed with raw atlas geometry. + +This proposal therefore recommends two safe options: + +1. preserve the raw exported polygon atlas by default; or +2. if smoothing is desired, generate a separate shared display atlas such as `dk_polygons_chaikin.csv` and render that same smoothed geometry in both Python and R. + +In other words, smoothing should be explicit and reproducible. A renderer should not apply backend-specific smoothing by default if cross-language parity is the goal. If future smoothing is implemented inside python-ggseg, it should ideally be topology-aware or accompanied by tests for gaps, overlaps, and invalid polygons. + +## Validation and tests + +A minimal test set for an upstream PR could include: + +1. Load an R-exported DK atlas polygon CSV. +2. Assert required columns are present. +3. Assert expected row, label, and feature counts for the fixture. +4. Join a small values table by `label`. +5. Report unmatched input regions and missing atlas regions. +6. Render SVG/PDF/PNG without errors. +7. Verify that `fill_hex` is treated as an identity color column. +8. Confirm that no smoothing is applied unless explicitly requested. + +Optional visual regression testing can compare a rendered SVG/PNG against a stored reference, but schema-level and join-level tests should be the required base. + +## Backward compatibility + +This can be added as an optional input path without changing existing `python-ggseg` behavior. Native Python atlas rendering remains supported. Users who do not need R interop do not need to change their workflow. + +## Documentation additions + +Recommended docs/vignette topics: + +- "Rendering an R ggseg atlas export in Python" +- "Matching Python and R ggseg-style figures" +- "Using `fill_hex` for exact cross-backend colors" +- "Why per-parcel smoothing is not topology-safe by default" + +## Suggested contribution strategy + +The GitHub operation for submitting this upstream is a **Pull Request (PR)**. A careful sequence would be: + +1. Open an upstream **Issue** or **Discussion/RFC** first, because this is an interoperability/API proposal. +2. Fork the upstream repository. +3. Create a feature branch, for example `feature/r-ggseg-polygon-interop`. +4. Implement the smallest useful slice: reader + validator + identity `fill_hex` plotting + tests. +5. Push the branch to the fork. +6. Open a **Pull Request** from the fork branch to the upstream repository. + +## Draft issue / discussion title + +```text +Proposal: support R ggseg-exported polygon atlas files for cross-language reproducibility +``` + +## Draft issue / discussion body + +```markdown +Hi, thank you for maintaining python-ggseg. We are building a Python/R dual-backend cortical visualization workflow and found that exact cross-language parity is easiest when both backends render the same polygon geometry. + +Would you be open to supporting R ggseg-exported 2D polygon atlas tables as an optional input path in python-ggseg? + +The proposed minimum schema is: + +label, region, hemi, view, x, y, .group, subgroup, .feature_id + +Optional columns such as value and fill_hex would allow joined data and exact identity-color rendering. This would not replace native python-ggseg atlases; it would add an interoperability bridge for users who want to render R ggseg atlas exports in Python. + +I would be happy to prepare a small PR with: + +- a polygon atlas CSV reader and schema validator; +- explicit join-by-label / join-by-region behavior; +- identity fill_hex rendering; +- a small R-exported DK fixture; +- tests for schema, joins, unmatched labels, and rendering. + +One design point: I would not propose per-parcel smoothing as a default, because independent smoothing can break shared boundaries between adjacent parcels. Raw polygons plus vector/high-DPI export are safer by default. +``` + +## Draft PR title + +```text +Add support for R ggseg-exported polygon atlas tables +``` + +## Draft PR summary + +```markdown +This PR adds an optional interoperability path for rendering polygon atlas tables exported from R ggseg. It includes: + +- a reader for canonical polygon atlas CSV files; +- schema validation for required ggseg-style columns; +- explicit data joining by label or region; +- identity color rendering via fill_hex; +- tests covering schema validation, joins, unmatched labels, and rendering. + +The feature is additive and does not change existing native python-ggseg atlas behavior. +``` + +## Local evidence from this repository + +In the local cortex-visualization skill, a shared DK polygon export was rendered by both Python and R. The shared plot-data table contained: + +```text +rows = 3233 +labels = 70 +features = 204 +fill_hex colors = 63 +``` + +The practical lesson is that cross-language parity is reliable when geometry and colors are shared explicitly rather than reconstructed independently in each backend. diff --git a/ggseg/__init__.py b/ggseg/__init__.py index 54549e9..5257f43 100644 --- a/ggseg/__init__.py +++ b/ggseg/__init__.py @@ -304,3 +304,247 @@ def plot_aseg(data, cmap='Spectral', background='k', edgecolor='w', ylabel='', _add_colorbar_(ax, cmap, norm, edgecolor, fontsize*0.75, ylabel) plt.show() + +# ----------------------------------------------------------------------------- +# Polygon atlas interoperability API +# ----------------------------------------------------------------------------- + +POLYGON_ATLAS_REQUIRED_COLUMNS = ( + 'label', 'region', 'hemi', 'view', 'x', 'y', '.group', 'subgroup', + '.feature_id' +) + + +def read_polygon_atlas(path): + """Read an R ggseg-style 2D polygon atlas table. + + Parameters + ---------- + path : str or path-like + CSV file containing at least the canonical polygon atlas columns: + label, region, hemi, view, x, y, .group, subgroup, .feature_id. + + Returns + ------- + list of dict + One row per polygon coordinate. Numeric x/y coordinates are converted + to float; all other columns are preserved as strings. + """ + import csv + + rows = [] + with open(path, newline='', encoding='utf-8-sig') as f: + reader = csv.DictReader(f) + for row in reader: + if row.get('x') not in (None, ''): + row['x'] = float(row['x']) + if row.get('y') not in (None, ''): + row['y'] = float(row['y']) + rows.append(row) + validate_polygon_atlas(rows) + return rows + + +def validate_polygon_atlas(atlas): + """Validate the canonical R ggseg polygon atlas schema. + + Raises + ------ + ValueError + If the atlas is empty or missing required columns. + """ + if not atlas: + raise ValueError('Polygon atlas is empty') + cols = set(atlas[0].keys()) + missing = [c for c in POLYGON_ATLAS_REQUIRED_COLUMNS if c not in cols] + if missing: + raise ValueError('Polygon atlas missing required columns: %s' % + ', '.join(missing)) + return True + + +def join_polygon_values(atlas, data, match_column='label', value_column='value', + strict=False): + """Join region-level values to a polygon atlas table. + + Parameters + ---------- + atlas : list of dict + Polygon atlas rows, typically returned by ``read_polygon_atlas``. + data : dict or iterable of dict + Region-level values. A dict is interpreted as {match_key: value}. If an + iterable of dict rows is supplied, each row must contain match_column + and value_column. + match_column : str + Atlas/data key used for joining, commonly 'label' or 'region'. + value_column : str + Name of the value column written to returned rows. + strict : bool + If True, raise when input keys are not present in the atlas. + + Returns + ------- + list of dict + Copy of atlas rows with a joined value column. + """ + validate_polygon_atlas(atlas) + if match_column not in atlas[0]: + raise ValueError('Atlas missing match column: %s' % match_column) + + if isinstance(data, dict): + values = data + else: + values = {} + for row in data: + values[row[match_column]] = row[value_column] + + atlas_keys = set(row[match_column] for row in atlas) + unmatched = set(values.keys()).difference(atlas_keys) + if unmatched and strict: + raise ValueError('Input keys not found in atlas: %s' % + ', '.join(sorted(unmatched))) + + joined = [] + for row in atlas: + out = dict(row) + key = row[match_column] + out[value_column] = values.get(key, None) + joined.append(out) + return joined + + +def _polygon_value_colors_(rows, value_column, cmap, vminmax, na_color): + import math + import matplotlib + + values = [] + for row in rows: + value = row.get(value_column) + if value not in (None, ''): + try: + value = float(value) + except (TypeError, ValueError): + continue + if not math.isnan(value): + values.append(value) + + if not values: + return None, None, [na_color for _ in rows] + + cmap_obj = matplotlib.cm.get_cmap(cmap) + if vminmax == [] or vminmax is None: + vmin, vmax = min(values), max(values) + else: + vmin, vmax = vminmax + norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax) + + colors = [] + for row in rows: + value = row.get(value_column) + if value in (None, ''): + colors.append(na_color) + else: + try: + colors.append(cmap_obj(norm(float(value)))) + except (TypeError, ValueError): + colors.append(na_color) + return cmap_obj, norm, colors + + +def _group_polygon_rows_(rows): + from collections import OrderedDict + + features = OrderedDict() + xs = [] + ys = [] + for row in rows: + fid = row['.feature_id'] + subgroup = row.get('subgroup') or '1' + features.setdefault(fid, OrderedDict()) + features[fid].setdefault(subgroup, []) + x = float(row['x']) + y = float(row['y']) + features[fid][subgroup].append((x, y, row)) + xs.append(x) + ys.append(y) + return features, xs, ys + + +def plot_polygons(atlas, data=None, match_column='label', value_column='value', + fill_column=None, cmap='Spectral', vminmax=[], + background='w', edgecolor='k', na_color='gray', + linewidth=1, figsize=(12, 8), title='', fontsize=15, + add_colorbar=True, ylabel='', ax=None, show=True): + """Plot a shared 2D polygon atlas table. + + This function supports R ggseg-exported polygon atlas tables and is intended + for strict R/Python visual parity workflows. If ``fill_column`` is provided + (for example ``fill_hex``), colors are used as identity colors and no + colormap remapping is performed. + + Parameters + ---------- + atlas : list of dict + Polygon atlas rows. + data : dict or iterable of dict, optional + Region-level values to join before plotting. + match_column : str + Join key, commonly 'label' or 'region'. + value_column : str + Numeric value column to color by when fill_column is not supplied. + fill_column : str, optional + Identity color column such as 'fill_hex'. + show : bool + If True, call matplotlib.pyplot.show(). Set False in tests/scripts. + + Returns + ------- + matplotlib.axes.Axes + Axis containing the polygon plot. + """ + import matplotlib.pyplot as plt + import matplotlib.patches as patches + + validate_polygon_atlas(atlas) + rows = join_polygon_values(atlas, data, match_column, value_column) if data is not None else [dict(r) for r in atlas] + + features, xs, ys = _group_polygon_rows_(rows) + if ax is None: + fig, ax = plt.subplots(figsize=figsize, facecolor=background) + else: + fig = ax.figure + fig.set_facecolor(background) + ax.set_facecolor(background) + + if fill_column is not None: + colors = [row.get(fill_column) or na_color for row in rows] + cmap_obj, norm = None, None + else: + cmap_obj, norm, colors = _polygon_value_colors_(rows, value_column, cmap, vminmax, na_color) + + row_color = {id(row): color for row, color in zip(rows, colors)} + for rings in features.values(): + for points_rows in rings.values(): + points = [(x, y) for x, y, _row in points_rows] + if len(points) < 3: + continue + color = row_color[id(points_rows[0][2])] + ax.add_patch(patches.Polygon(points, closed=True, facecolor=color, + edgecolor=edgecolor, linewidth=linewidth, + joinstyle='round')) + + padx = (max(xs) - min(xs)) * 0.04 if xs else 1 + pady = (max(ys) - min(ys)) * 0.04 if ys else 1 + ax.set_xlim(min(xs) - padx, max(xs) + padx) + ax.set_ylim(min(ys) - pady, max(ys) + pady) + ax.set_aspect('equal') + ax.axis('off') + ax.set_title(title, fontsize=fontsize, color=edgecolor) + + if add_colorbar and cmap_obj is not None and norm is not None: + _add_colorbar_(ax, cmap_obj, norm, edgecolor, fontsize * 0.75, ylabel) + + if show: + plt.show() + return ax + diff --git a/ggseg/tests/test_polygon_interop.py b/ggseg/tests/test_polygon_interop.py new file mode 100644 index 0000000..1ea2ec3 --- /dev/null +++ b/ggseg/tests/test_polygon_interop.py @@ -0,0 +1,103 @@ +import csv +import tempfile + +import matplotlib +matplotlib.use("Agg") + +import ggseg + + +ATLAS_ROWS = [ + { + "label": "region_a_left", "region": "region_a", "hemi": "left", + "view": "lateral", "x": "0", "y": "0", ".group": "1", + "subgroup": "1", ".feature_id": "a" + }, + { + "label": "region_a_left", "region": "region_a", "hemi": "left", + "view": "lateral", "x": "1", "y": "0", ".group": "1", + "subgroup": "1", ".feature_id": "a" + }, + { + "label": "region_a_left", "region": "region_a", "hemi": "left", + "view": "lateral", "x": "0", "y": "1", ".group": "1", + "subgroup": "1", ".feature_id": "a" + }, + { + "label": "region_b_left", "region": "region_b", "hemi": "left", + "view": "lateral", "x": "2", "y": "0", ".group": "2", + "subgroup": "1", ".feature_id": "b" + }, + { + "label": "region_b_left", "region": "region_b", "hemi": "left", + "view": "lateral", "x": "3", "y": "0", ".group": "2", + "subgroup": "1", ".feature_id": "b" + }, + { + "label": "region_b_left", "region": "region_b", "hemi": "left", + "view": "lateral", "x": "2", "y": "1", ".group": "2", + "subgroup": "1", ".feature_id": "b" + }, +] + + +def _write_atlas_csv(rows): + f = tempfile.NamedTemporaryFile("w", newline="", suffix=".csv", delete=False) + with f: + writer = csv.DictWriter(f, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + return f.name + + +def test_read_polygon_atlas_and_validate_schema(): + path = _write_atlas_csv(ATLAS_ROWS) + atlas = ggseg.read_polygon_atlas(path) + assert len(atlas) == 6 + assert atlas[0]["x"] == 0.0 + assert atlas[2]["y"] == 1.0 + assert ggseg.validate_polygon_atlas(atlas) is True + + +def test_join_polygon_values_by_label(): + atlas = [dict(row) for row in ATLAS_ROWS] + joined = ggseg.join_polygon_values( + atlas, + {"region_a_left": 1.5, "region_b_left": 2.5}, + match_column="label", + ) + values = set(row["value"] for row in joined) + assert values == {1.5, 2.5} + + +def test_join_polygon_values_strict_reports_unmatched(): + atlas = [dict(row) for row in ATLAS_ROWS] + try: + ggseg.join_polygon_values( + atlas, + {"missing_region": 1.5}, + match_column="label", + strict=True, + ) + except ValueError as exc: + assert "missing_region" in str(exc) + else: + raise AssertionError("strict join should fail for unmatched labels") + + +def test_plot_polygons_with_fill_hex_identity_colors(): + atlas = [] + for row in ATLAS_ROWS: + out = dict(row) + out["fill_hex"] = "#FF0000" if out["label"] == "region_a_left" else "#0000FF" + atlas.append(out) + + ax = ggseg.plot_polygons( + atlas, + fill_column="fill_hex", + background="w", + edgecolor="k", + add_colorbar=False, + show=False, + ) + assert len(ax.patches) == 2