diff --git a/.gitignore b/.gitignore
index 74a3feeb..7d8dfea1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -186,3 +186,9 @@ scratch/
 plugins/
 
 .claude/worktrees/sharp-burnell/
+
+# thesis-only generated artifacts (regenerable from validation notebooks)
+examples/head_models/thesis_results_out/
+
+# machine-local scan-data symlink (points to /home/ma7/BA/PG_Subjects)
+examples/head_models/PG_Subjects
diff --git a/README.md b/README.md
index 369cc7e5..f8788f85 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,111 @@
+# Face Anonymization for Photogrammetry Scans
+
+M.Sc. thesis at TU Berlin / IBS Lab: "Development of a Landmark-Aware Face-Removal Algorithm for Photogrammetric Head Scans for Data Protection".
+
+This is a fork of [cedalion](https://github.com/ibs-lab/cedalion) with a geometric face anonymization module for Einstar photogrammetry scans used in fNIRS research. The module lives at:
+
+```
+src/cedalion/geometry/photogrammetry/anonymization/
+```
+
+## Installation
+
+Requires the cedalion conda environment:
+
+```bash
+conda env create -f environment_dev.yml
+conda activate cedalion
+pip install -e .
+```
+
+## Usage
+
+### Interactive single-scan workflow
+
+Open `examples/head_models/51_manual_5pt_anonymization.ipynb`. The notebook:
+1. Loads an Einstar scan (`cedalion.io.read_einstar_obj`)
+2. Picks the five 10-20 landmarks interactively (Nz, Iz, Cz, LPA, RPA)
+3. Calls `anonymize_scan(surface, landmarks)`
+4. Shows a before/after comparison
+5. Saves the anonymized OBJ + sanitized JPG texture bundle via `save_anonymized_scan`
+
+## Module structure
+
+```
+src/cedalion/geometry/photogrammetry/anonymization/
+├── __init__.py        public API, re-exports all functions listed below
+├── pipeline.py        anonymize_scan (entry point)
+├── preprocessing.py   normalize_axes, isolate_head, align_axes_from_landmarks,
+│                      revert_to_einstar_frame
+├── mask.py            detect_cap_boundary, face_mask_from_landmarks,
+│                      delete_masked_vertices, save_anonymized_scan
+└── _utils.py          private helpers shared by preprocessing and mask
+                       (_rebuild_mesh, _copy_visual, _reindex_faces,
+                        _apply_affine, _transform_labeled_points,
+                        _ear_midpoint, _upper_head_centroid,
+                        _resolve_texture_image)
+
+examples/head_models/
+└── 51_manual_5pt_anonymization.ipynb   interactive workflow notebook
+
+tests/
+└── test_anonymization.py               26 unit tests
+```
+
+Pipeline steps inside `anonymize_scan`:
+
+1. `normalize_axes`: rotate so +Y points anterior (handles arbitrary Einstar orientation)
+2. `isolate_head`: remove body, shoulders, and disconnected fragments
+3. `align_axes_from_landmarks`: map to CTF frame (+X anterior, +Y left, +Z up)
+4. `detect_cap_boundary`: locate the front cap-edge height along Z
+5. `face_mask_from_landmarks`: face region union ear spheres, clamped below the cap
+6. Landmark preservation: 8 mm spheres around each landmark + midline nasion strip
+7. `delete_masked_vertices`: drop triangles touching any masked vertex, UVs in sync
+8. `revert_to_einstar_frame`: return to `crs="digitized"` for saving
+
+## Tests
+
+```bash
+pytest tests/test_anonymization.py -v
+```
+
+26 tests covering all eight public functions and the end-to-end pipeline. No real scan data needed: all tests build synthetic geometry with `trimesh.creation.icosphere`, the same approach used elsewhere in cedalion (see `test_geodesics.py`, `test_dataclasses_geometry.py`). Three fixtures are shared:
+
+- `simple_sphere_surface`: unit icosphere for geometry-only checks
+- `head_like_surface`: elongated icosphere (X scaled x1.2), giving the masking step a real face region to remove
+- `axis_normalized_landmarks`: five `LabeledPoints` (Nz, Iz, Cz, LPA, RPA) on the sphere axes in the post-`normalize_axes` frame
+
+## Branch layout
+
+| Branch | Contents |
+|--------|---------|
+| `feature/face-anonymization` | **This branch**: thesis implementation (anonymization module, notebook 51, test suite) |
+| `main` | Upstream cedalion base |
+| `validation/face-anonymization` | Measurement notebooks behind the thesis results tables (see below) |
+| `auxiliary/mediapipe-nasion` | Experimental automatic nasion detection (MediaPipe) |
+
+### Validation branch contents
+
+The notebooks on `validation/face-anonymization` produce the per-subject numbers in the thesis results tables. They run on top of the shipped pipeline's output and do not modify it.
+
+| Notebook | What it checks | Headline result on the 11-subject cohort |
+|---|---|---|
+| `64_batch_validation` | Driver that runs the four validators below on every subject and emits per-subject CSVs | per-subject CSVs feeding the thesis tables |
+| `66_preservation_check` | 10-20 landmark deviation + bit-exact vertex preservation across the surviving surface | 0.000 mm across all 55 landmark measurements (5 landmarks x 11 subjects) |
+| `68_coreg_invariance` | Cedalion `ColoredStickerProcessor` re-run on original vs anonymized mesh (optode-cap subcohort) | 131 matched stickers, 0 mm deviation in sticker centres and scalp-projected optode positions |
+| `69_coreg_visual_check` | Side-by-side render of detected stickers on the original and anonymized mesh; any original-side detection without a match on the anonymized side is drawn in red so the operator can spot it | qualitative; the manual visibility check that backs notebook 68 |
+| `70_auxiliary_nasion` | MediaPipe Face Landmarker auto-nasion vs the manually picked nasion | 13.68 mm cohort mean offset, 30.78 mm worst case |
+| `72_face_detectability_comparison` | MediaPipe Face Detector hit counts under a 21-view sweep: original vs vertex-deletion vs noise-perturbation | 93/231 -> 34/231 (63% reduction); optode-cap subcohort 32/147 -> 4/147 (88%); noise reaches 28/231 (Wilcoxon p=0.375, not distinguishable from deletion) |
+| `73_s8_mediapipe_boxes` | BlazeFace bounding box and six-keypoint inspection on a single subject | qualitative figure for the discussion of why post-deletion residual hits are silhouette artefacts |
+
+Helper modules `_thesis_data.py`, `_thesis_pipeline.py`, `_validator_noise.py`, and `_validator_render.py` hold the cohort walk, pipeline call, noise-perturbation operator, and contact-sheet rendering shared across these notebooks.
+
+---
+
+*Upstream cedalion documentation below.*
+
+---
+
 # Cedalion - fNIRS analysis toolbox
 
 A python-based framework for the data-driven analysis of multimodal fNIRS and DOT in naturalistic environments. Developed by the [Intelligent Biomedical Sensing (IBS) Lab](https://ibs-lab.com/) with and for the community.
diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 1dca972b..18431423 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -1,7 +1,13 @@
 # Changelog
 
 ## Unreleased changes (available on the `dev` branch)
- 
+
+### Additions and Changes
+
+#### Photogrammetry
+
+- Added `cedalion.geometry.photogrammetry.anonymization`, a deletion-based face removal pipeline for photogrammetry scans. From the five anatomical landmarks (Nz, Iz, Cz, LPA, RPA), the pipeline normalizes axes, isolates the head, detects the cap boundary, builds a face mask, deletes the masked vertices, and saves the anonymized mesh together with the landmarks. See `examples/head_models/51_manual_5pt_anonymization.ipynb` for an end-to-end walkthrough.
+
 ## Version 26.04.0 (in preparation)
 
 ### Additions and Changes
diff --git a/examples/head_models/51_manual_5pt_anonymization.ipynb b/examples/head_models/51_manual_5pt_anonymization.ipynb
new file mode 100644
index 00000000..79cef4cc
--- /dev/null
+++ b/examples/head_models/51_manual_5pt_anonymization.ipynb
@@ -0,0 +1,779 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "e505f83d-79aa-467b-a6f2-44d973491f55",
+   "metadata": {},
+   "source": [
+    "# Face Anonymization from 5 Given Landmarks\n",
+    "\n",
+    "Geometric face removal for Einstar photogrammetry scans. Given the 5\n",
+    "landmarks (Nz, Iz, Cz, Lpa, Rpa), `anonymize_scan` deletes the face and\n",
+    "returns the result in the raw Einstar (`crs=digitized`) frame, so the\n",
+    "saved `.obj` lines up with the original at co-registration.\n",
+    "\n",
+    "Flow: load → pick 5 landmarks → `anonymize_scan` → save.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a916040f-02d3-44e3-a1f4-05188c814d1d",
+   "metadata": {},
+   "outputs": [],
+   "source": "import logging\nimport os\nfrom pathlib import Path\n\nimport numpy as np\nimport pyvista as pv\nimport xarray as xr\nfrom PIL import Image\nimport trimesh\n\nimport cedalion\nimport cedalion.dataclasses as cdc\nimport cedalion.io\nimport cedalion.vis.blocks\nfrom cedalion.vtktutils import trimesh_to_vtk_polydata\nfrom cedalion.geometry.photogrammetry.anonymization import (\n    anonymize_scan,\n    save_anonymized_scan,\n)\n\nlogging.getLogger('cedalion').setLevel(logging.WARNING)\npv.set_jupyter_backend('server')\n\n# === Options ===\n# Path to your Einstar .obj scan. The notebook expects a sibling .jpg\n# texture (Einstar exports both).\nSCAN_PATH = './scan.obj'         # edit me\nOUT_PATH = None                  # None -> write '<stem>_anon.obj' next to SCAN_PATH\n\n# True: run the picker. False: load landmarks from a TSV sidecar (set below).\nINTERACTIVE = True\nCACHED_LANDMARKS_TSV = None      # used only when INTERACTIVE = False\n\n\n# === Mask parameters (edit to adjust the anonymized region) ===\nEAR_DELETE_RADIUS_MM = 40.0      # sphere radius around LPA/RPA included in deletion\nLANDMARK_KEEP_RADIUS_MM = 8.0    # sphere radius around each landmark kept intact\nEYEBROW_OFFSET_MM = 10.0         # failsafe cap height above Nz (flush-cap fallback)\nCAP_Z_CEILING_MM = 40.0          # mm above Nz where cap detection is clamped\nHEAD_ISOLATION_RADIUS_MM = 220.0 # sphere radius used to strip body/shoulders\n"
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c7b11266-a722-49ac-9aee-67953820c2ac",
+   "metadata": {},
+   "source": [
+    "## Load the Einstar scan\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "5c99124b-f88c-4aaf-bff2-124cfadd5da1",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Loaded: 943,539 vertices, 1,736,738 faces\n"
+     ]
+    }
+   ],
+   "source": [
+    "surface = cedalion.io.read_einstar_obj(SCAN_PATH)\n",
+    "print(f'Loaded: {surface.nvertices:,} vertices, {surface.nfaces:,} faces')\n",
+    "\n",
+    "if OUT_PATH is None:\n",
+    "    p = Path(SCAN_PATH)\n",
+    "    OUT_PATH = str(p.with_name(p.stem + '_anon' + p.suffix))\n",
+    "\n",
+    "# trimesh 4.6 places the texture image on visual.material.image; if neither\n",
+    "# path has it, attach it from the sibling JPG so save_anonymized_scan can\n",
+    "# sanitize the final texture.\n",
+    "visual = surface.mesh.visual\n",
+    "img = getattr(visual, 'image', None) or getattr(getattr(visual, 'material', None), 'image', None)\n",
+    "if img is None:\n",
+    "    jpg = SCAN_PATH.replace('.obj', '.jpg')\n",
+    "    uv = getattr(visual, 'uv', None)\n",
+    "    assert os.path.exists(jpg) and uv is not None, 'no texture available'\n",
+    "    surface.mesh.visual = trimesh.visual.TextureVisuals(\n",
+    "        uv=uv, image=Image.open(jpg).convert('RGBA'),\n",
+    "    )\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5c282588-d37d-4810-a5c1-104b751d2155",
+   "metadata": {},
+   "source": [
+    "## Pick the 5 landmarks\n",
+    "\n",
+    "Right-click on the mesh to place a sphere; click a sphere to cycle its label\n",
+    "through `Nz -> Iz -> Cz -> Lpa -> Rpa`. Close the window when all 5 are placed.\n",
+    "Skip this cell if `INTERACTIVE = False`.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "4ea3ff12-ca93-4c03-9f50-b84c28c7cfdc",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "8931ae160417497089d51fca6c1cd106",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Widget(value='<iframe src=\"http://localhost:44253/index.html?ui=P_0x7fe28a702550_2&reconnect=auto\" class=\"pyvi…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "if INTERACTIVE:\n",
+    "    pvplt = pv.Plotter()\n",
+    "    get_landmarks = cedalion.vis.blocks.plot_surface(\n",
+    "        pvplt, surface, opacity=1.0, pick_landmarks=True,\n",
+    "    )\n",
+    "    pvplt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f98684b3-69b3-49c8-a212-c726d6223125",
+   "metadata": {},
+   "source": [
+    "## Wrap picked points into a `LabeledPoints`\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "3a09ed60-b599-4191-b394-460b852f4b27",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
+       "<defs>\n",
+       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
+       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
+       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
+       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
+       "</symbol>\n",
+       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
+       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
+       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
+       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
+       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
+       "</symbol>\n",
+       "</defs>\n",
+       "</svg>\n",
+       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
+       " *\n",
+       " */\n",
+       "\n",
+       ":root {\n",
+       "  --xr-font-color0: var(\n",
+       "    --jp-content-font-color0,\n",
+       "    var(--pst-color-text-base rgba(0, 0, 0, 1))\n",
+       "  );\n",
+       "  --xr-font-color2: var(\n",
+       "    --jp-content-font-color2,\n",
+       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.54))\n",
+       "  );\n",
+       "  --xr-font-color3: var(\n",
+       "    --jp-content-font-color3,\n",
+       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.38))\n",
+       "  );\n",
+       "  --xr-border-color: var(\n",
+       "    --jp-border-color2,\n",
+       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 10))\n",
+       "  );\n",
+       "  --xr-disabled-color: var(\n",
+       "    --jp-layout-color3,\n",
+       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 40))\n",
+       "  );\n",
+       "  --xr-background-color: var(\n",
+       "    --jp-layout-color0,\n",
+       "    var(--pst-color-on-background, white)\n",
+       "  );\n",
+       "  --xr-background-color-row-even: var(\n",
+       "    --jp-layout-color1,\n",
+       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 5))\n",
+       "  );\n",
+       "  --xr-background-color-row-odd: var(\n",
+       "    --jp-layout-color2,\n",
+       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 15))\n",
+       "  );\n",
+       "}\n",
+       "\n",
+       "html[theme=\"dark\"],\n",
+       "html[data-theme=\"dark\"],\n",
+       "body[data-theme=\"dark\"],\n",
+       "body.vscode-dark {\n",
+       "  --xr-font-color0: var(\n",
+       "    --jp-content-font-color0,\n",
+       "    var(--pst-color-text-base, rgba(255, 255, 255, 1))\n",
+       "  );\n",
+       "  --xr-font-color2: var(\n",
+       "    --jp-content-font-color2,\n",
+       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.54))\n",
+       "  );\n",
+       "  --xr-font-color3: var(\n",
+       "    --jp-content-font-color3,\n",
+       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.38))\n",
+       "  );\n",
+       "  --xr-border-color: var(\n",
+       "    --jp-border-color2,\n",
+       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 10))\n",
+       "  );\n",
+       "  --xr-disabled-color: var(\n",
+       "    --jp-layout-color3,\n",
+       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 40))\n",
+       "  );\n",
+       "  --xr-background-color: var(\n",
+       "    --jp-layout-color0,\n",
+       "    var(--pst-color-on-background, #111111)\n",
+       "  );\n",
+       "  --xr-background-color-row-even: var(\n",
+       "    --jp-layout-color1,\n",
+       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 5))\n",
+       "  );\n",
+       "  --xr-background-color-row-odd: var(\n",
+       "    --jp-layout-color2,\n",
+       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 15))\n",
+       "  );\n",
+       "}\n",
+       "\n",
+       ".xr-wrap {\n",
+       "  display: block !important;\n",
+       "  min-width: 300px;\n",
+       "  max-width: 700px;\n",
+       "}\n",
+       "\n",
+       ".xr-text-repr-fallback {\n",
+       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-header {\n",
+       "  padding-top: 6px;\n",
+       "  padding-bottom: 6px;\n",
+       "  margin-bottom: 4px;\n",
+       "  border-bottom: solid 1px var(--xr-border-color);\n",
+       "}\n",
+       "\n",
+       ".xr-header > div,\n",
+       ".xr-header > ul {\n",
+       "  display: inline;\n",
+       "  margin-top: 0;\n",
+       "  margin-bottom: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-obj-type,\n",
+       ".xr-array-name {\n",
+       "  margin-left: 2px;\n",
+       "  margin-right: 10px;\n",
+       "}\n",
+       "\n",
+       ".xr-obj-type {\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-sections {\n",
+       "  padding-left: 0 !important;\n",
+       "  display: grid;\n",
+       "  grid-template-columns: 150px auto auto 1fr 0 20px 0 20px;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item {\n",
+       "  display: contents;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input {\n",
+       "  display: inline-block;\n",
+       "  opacity: 0;\n",
+       "  height: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input + label {\n",
+       "  color: var(--xr-disabled-color);\n",
+       "  border: 2px solid transparent !important;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input:enabled + label {\n",
+       "  cursor: pointer;\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input:focus + label {\n",
+       "  border: 2px solid var(--xr-font-color0) !important;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input:enabled + label:hover {\n",
+       "  color: var(--xr-font-color0);\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary {\n",
+       "  grid-column: 1;\n",
+       "  color: var(--xr-font-color2);\n",
+       "  font-weight: 500;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary > span {\n",
+       "  display: inline-block;\n",
+       "  padding-left: 0.5em;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:disabled + label {\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in + label:before {\n",
+       "  display: inline-block;\n",
+       "  content: \"►\";\n",
+       "  font-size: 11px;\n",
+       "  width: 15px;\n",
+       "  text-align: center;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:disabled + label:before {\n",
+       "  color: var(--xr-disabled-color);\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:checked + label:before {\n",
+       "  content: \"▼\";\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:checked + label > span {\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary,\n",
+       ".xr-section-inline-details {\n",
+       "  padding-top: 4px;\n",
+       "  padding-bottom: 4px;\n",
+       "}\n",
+       "\n",
+       ".xr-section-inline-details {\n",
+       "  grid-column: 2 / -1;\n",
+       "}\n",
+       "\n",
+       ".xr-section-details {\n",
+       "  display: none;\n",
+       "  grid-column: 1 / -1;\n",
+       "  margin-bottom: 5px;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
+       "  display: contents;\n",
+       "}\n",
+       "\n",
+       ".xr-array-wrap {\n",
+       "  grid-column: 1 / -1;\n",
+       "  display: grid;\n",
+       "  grid-template-columns: 20px auto;\n",
+       "}\n",
+       "\n",
+       ".xr-array-wrap > label {\n",
+       "  grid-column: 1;\n",
+       "  vertical-align: top;\n",
+       "}\n",
+       "\n",
+       ".xr-preview {\n",
+       "  color: var(--xr-font-color3);\n",
+       "}\n",
+       "\n",
+       ".xr-array-preview,\n",
+       ".xr-array-data {\n",
+       "  padding: 0 5px !important;\n",
+       "  grid-column: 2;\n",
+       "}\n",
+       "\n",
+       ".xr-array-data,\n",
+       ".xr-array-in:checked ~ .xr-array-preview {\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-array-in:checked ~ .xr-array-data,\n",
+       ".xr-array-preview {\n",
+       "  display: inline-block;\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list {\n",
+       "  display: inline-block !important;\n",
+       "  list-style: none;\n",
+       "  padding: 0 !important;\n",
+       "  margin: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list li {\n",
+       "  display: inline-block;\n",
+       "  padding: 0;\n",
+       "  margin: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list:before {\n",
+       "  content: \"(\";\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list:after {\n",
+       "  content: \")\";\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list li:not(:last-child):after {\n",
+       "  content: \",\";\n",
+       "  padding-right: 5px;\n",
+       "}\n",
+       "\n",
+       ".xr-has-index {\n",
+       "  font-weight: bold;\n",
+       "}\n",
+       "\n",
+       ".xr-var-list,\n",
+       ".xr-var-item {\n",
+       "  display: contents;\n",
+       "}\n",
+       "\n",
+       ".xr-var-item > div,\n",
+       ".xr-var-item label,\n",
+       ".xr-var-item > .xr-var-name span {\n",
+       "  background-color: var(--xr-background-color-row-even);\n",
+       "  border-color: var(--xr-background-color-row-odd);\n",
+       "  margin-bottom: 0;\n",
+       "  padding-top: 2px;\n",
+       "}\n",
+       "\n",
+       ".xr-var-item > .xr-var-name:hover span {\n",
+       "  padding-right: 5px;\n",
+       "}\n",
+       "\n",
+       ".xr-var-list > li:nth-child(odd) > div,\n",
+       ".xr-var-list > li:nth-child(odd) > label,\n",
+       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
+       "  background-color: var(--xr-background-color-row-odd);\n",
+       "  border-color: var(--xr-background-color-row-even);\n",
+       "}\n",
+       "\n",
+       ".xr-var-name {\n",
+       "  grid-column: 1;\n",
+       "}\n",
+       "\n",
+       ".xr-var-dims {\n",
+       "  grid-column: 2;\n",
+       "}\n",
+       "\n",
+       ".xr-var-dtype {\n",
+       "  grid-column: 3;\n",
+       "  text-align: right;\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-var-preview {\n",
+       "  grid-column: 4;\n",
+       "}\n",
+       "\n",
+       ".xr-index-preview {\n",
+       "  grid-column: 2 / 5;\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-var-name,\n",
+       ".xr-var-dims,\n",
+       ".xr-var-dtype,\n",
+       ".xr-preview,\n",
+       ".xr-attrs dt {\n",
+       "  white-space: nowrap;\n",
+       "  overflow: hidden;\n",
+       "  text-overflow: ellipsis;\n",
+       "  padding-right: 10px;\n",
+       "}\n",
+       "\n",
+       ".xr-var-name:hover,\n",
+       ".xr-var-dims:hover,\n",
+       ".xr-var-dtype:hover,\n",
+       ".xr-attrs dt:hover {\n",
+       "  overflow: visible;\n",
+       "  width: auto;\n",
+       "  z-index: 1;\n",
+       "}\n",
+       "\n",
+       ".xr-var-attrs,\n",
+       ".xr-var-data,\n",
+       ".xr-index-data {\n",
+       "  display: none;\n",
+       "  border-top: 2px dotted var(--xr-background-color);\n",
+       "  padding-bottom: 20px !important;\n",
+       "  padding-top: 10px !important;\n",
+       "}\n",
+       "\n",
+       ".xr-var-attrs-in + label,\n",
+       ".xr-var-data-in + label,\n",
+       ".xr-index-data-in + label {\n",
+       "  padding: 0 1px;\n",
+       "}\n",
+       "\n",
+       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
+       ".xr-var-data-in:checked ~ .xr-var-data,\n",
+       ".xr-index-data-in:checked ~ .xr-index-data {\n",
+       "  display: block;\n",
+       "}\n",
+       "\n",
+       ".xr-var-data > table {\n",
+       "  float: right;\n",
+       "}\n",
+       "\n",
+       ".xr-var-data > pre,\n",
+       ".xr-index-data > pre,\n",
+       ".xr-var-data > table > tbody > tr {\n",
+       "  background-color: transparent !important;\n",
+       "}\n",
+       "\n",
+       ".xr-var-name span,\n",
+       ".xr-var-data,\n",
+       ".xr-index-name div,\n",
+       ".xr-index-data,\n",
+       ".xr-attrs {\n",
+       "  padding-left: 25px !important;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs,\n",
+       ".xr-var-attrs,\n",
+       ".xr-var-data,\n",
+       ".xr-index-data {\n",
+       "  grid-column: 1 / -1;\n",
+       "}\n",
+       "\n",
+       "dl.xr-attrs {\n",
+       "  padding: 0;\n",
+       "  margin: 0;\n",
+       "  display: grid;\n",
+       "  grid-template-columns: 125px auto;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dt,\n",
+       ".xr-attrs dd {\n",
+       "  padding: 0;\n",
+       "  margin: 0;\n",
+       "  float: left;\n",
+       "  padding-right: 10px;\n",
+       "  width: auto;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dt {\n",
+       "  font-weight: normal;\n",
+       "  grid-column: 1;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dt:hover span {\n",
+       "  display: inline-block;\n",
+       "  background: var(--xr-background-color);\n",
+       "  padding-right: 10px;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dd {\n",
+       "  grid-column: 2;\n",
+       "  white-space: pre-wrap;\n",
+       "  word-break: break-all;\n",
+       "}\n",
+       "\n",
+       ".xr-icon-database,\n",
+       ".xr-icon-file-text2,\n",
+       ".xr-no-icon {\n",
+       "  display: inline-block;\n",
+       "  vertical-align: middle;\n",
+       "  width: 1em;\n",
+       "  height: 1.5em !important;\n",
+       "  stroke-width: 0;\n",
+       "  stroke: currentColor;\n",
+       "  fill: currentColor;\n",
+       "}\n",
+       "\n",
+       ".xr-var-attrs-in:checked + label > .xr-icon-file-text2,\n",
+       ".xr-var-data-in:checked + label > .xr-icon-database,\n",
+       ".xr-index-data-in:checked + label > .xr-icon-database {\n",
+       "  color: var(--xr-font-color0);\n",
+       "  filter: drop-shadow(1px 1px 5px var(--xr-font-color2));\n",
+       "  stroke-width: 0.8px;\n",
+       "}\n",
+       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.DataArray (label: 5, digitized: 3)&gt; Size: 120B\n",
+       "&lt;Quantity([[135.82489764 134.11959626 388.34999789]\n",
+       " [ 83.19674841 -38.36289545 470.02017819]\n",
+       " [226.49436905 -13.33052629 371.30306422]\n",
+       " [144.87231487  60.99825005 490.28171817]\n",
+       " [ 74.20606761  16.62082718 368.93477872]], &#x27;millimeter&#x27;)&gt;\n",
+       "Coordinates:\n",
+       "  * label    (label) &lt;U3 60B &#x27;Nz&#x27; &#x27;Iz&#x27; &#x27;Cz&#x27; &#x27;Lpa&#x27; &#x27;Rpa&#x27;\n",
+       "    type     (label) object 40B PointType.LANDMARK ... PointType.LANDMARK\n",
+       "    group    (label) &lt;U1 20B &#x27;L&#x27; &#x27;L&#x27; &#x27;L&#x27; &#x27;L&#x27; &#x27;L&#x27;\n",
+       "Dimensions without coordinates: digitized</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'></div><ul class='xr-dim-list'><li><span class='xr-has-index'>label</span>: 5</li><li><span>digitized</span>: 3</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-813f730a-31f8-46bf-90d1-411fbb4c4615' class='xr-array-in' type='checkbox' checked><label for='section-813f730a-31f8-46bf-90d1-411fbb4c4615' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>[mm] 135.8 134.1 388.3 83.2 -38.36 ... 61.0 490.3 74.21 16.62 368.9</span></div><div class='xr-array-data'><table><tbody><tr><th>Magnitude</th><td style='text-align:left;'><pre>[[135.82489764260836 134.1195962551057 388.34999788776486] [83.19674841020856 -38.36289544683734 470.02017819060416] [226.49436904915945 -13.33052629089562 371.3030642194328] [144.8723148742921 60.99825004860746 490.2817181736568] [74.20606761144799 16.620827183111732 368.93477871973005]]</pre></td></tr><tr><th>Units</th><td style='text-align:left;'>millimeter</td></tr></tbody></table></div></div></li><li class='xr-section-item'><input id='section-863b5ce2-1e30-44df-a965-c675e3ec206a' class='xr-section-summary-in' type='checkbox'  checked><label for='section-863b5ce2-1e30-44df-a965-c675e3ec206a' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>label</span></div><div class='xr-var-dims'>(label)</div><div class='xr-var-dtype'>&lt;U3</div><div class='xr-var-preview xr-preview'>&#x27;Nz&#x27; &#x27;Iz&#x27; &#x27;Cz&#x27; &#x27;Lpa&#x27; &#x27;Rpa&#x27;</div><input id='attrs-1ae75da9-4f24-4008-a9c8-26fb768aedac' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-1ae75da9-4f24-4008-a9c8-26fb768aedac' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-42f5893b-9bcf-42d6-83af-eb15277be687' class='xr-var-data-in' type='checkbox'><label for='data-42f5893b-9bcf-42d6-83af-eb15277be687' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([&#x27;Nz&#x27;, &#x27;Iz&#x27;, &#x27;Cz&#x27;, &#x27;Lpa&#x27;, &#x27;Rpa&#x27;], dtype=&#x27;&lt;U3&#x27;)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>type</span></div><div class='xr-var-dims'>(label)</div><div class='xr-var-dtype'>object</div><div class='xr-var-preview xr-preview'>PointType.LANDMARK ... PointType...</div><input id='attrs-f6113b8c-f472-4a7c-8e45-0014e5560c6a' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-f6113b8c-f472-4a7c-8e45-0014e5560c6a' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-e6830bc7-aae4-4b53-a423-f215588a0e4a' class='xr-var-data-in' type='checkbox'><label for='data-e6830bc7-aae4-4b53-a423-f215588a0e4a' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([&lt;PointType.LANDMARK: 3&gt;, &lt;PointType.LANDMARK: 3&gt;,\n",
+       "       &lt;PointType.LANDMARK: 3&gt;, &lt;PointType.LANDMARK: 3&gt;,\n",
+       "       &lt;PointType.LANDMARK: 3&gt;], dtype=object)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>group</span></div><div class='xr-var-dims'>(label)</div><div class='xr-var-dtype'>&lt;U1</div><div class='xr-var-preview xr-preview'>&#x27;L&#x27; &#x27;L&#x27; &#x27;L&#x27; &#x27;L&#x27; &#x27;L&#x27;</div><input id='attrs-6c9359df-b959-47dc-a89d-4a4aad3e038f' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-6c9359df-b959-47dc-a89d-4a4aad3e038f' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f0e9a6dc-56be-4e46-b876-bec43da6a675' class='xr-var-data-in' type='checkbox'><label for='data-f0e9a6dc-56be-4e46-b876-bec43da6a675' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([&#x27;L&#x27;, &#x27;L&#x27;, &#x27;L&#x27;, &#x27;L&#x27;, &#x27;L&#x27;], dtype=&#x27;&lt;U1&#x27;)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-6aee2e82-d78d-4eaf-a768-be84733a8bf3' class='xr-section-summary-in' type='checkbox'  ><label for='section-6aee2e82-d78d-4eaf-a768-be84733a8bf3' class='xr-section-summary' >Indexes: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>label</div></div><div class='xr-index-preview'>PandasIndex</div><input type='checkbox' disabled/><label></label><input id='index-0351f477-b136-47da-b290-ed9b4ed249d6' class='xr-index-data-in' type='checkbox'/><label for='index-0351f477-b136-47da-b290-ed9b4ed249d6' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Index([&#x27;Nz&#x27;, &#x27;Iz&#x27;, &#x27;Cz&#x27;, &#x27;Lpa&#x27;, &#x27;Rpa&#x27;], dtype=&#x27;object&#x27;, name=&#x27;label&#x27;))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-e7577a12-af76-4399-ba73-750b2be851de' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-e7577a12-af76-4399-ba73-750b2be851de' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
+      ],
+      "text/plain": [
+       "<xarray.DataArray (label: 5, digitized: 3)> Size: 120B\n",
+       "<Quantity([[135.82489764 134.11959626 388.34999789]\n",
+       " [ 83.19674841 -38.36289545 470.02017819]\n",
+       " [226.49436905 -13.33052629 371.30306422]\n",
+       " [144.87231487  60.99825005 490.28171817]\n",
+       " [ 74.20606761  16.62082718 368.93477872]], 'millimeter')>\n",
+       "Coordinates:\n",
+       "  * label    (label) <U3 60B 'Nz' 'Iz' 'Cz' 'Lpa' 'Rpa'\n",
+       "    type     (label) object 40B PointType.LANDMARK ... PointType.LANDMARK\n",
+       "    group    (label) <U1 20B 'L' 'L' 'L' 'L' 'L'\n",
+       "Dimensions without coordinates: digitized"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "if INTERACTIVE:\n",
+    "    landmarks = get_landmarks()\n",
+    "else:\n",
+    "    assert CACHED_LANDMARKS_TSV, 'set CACHED_LANDMARKS_TSV or INTERACTIVE = True'\n",
+    "    landmarks = cedalion.io.load_tsv(CACHED_LANDMARKS_TSV)\n",
+    "\n",
+    "labels = list(landmarks['label'].values)\n",
+    "assert set(labels) == {'Nz', 'Iz', 'Cz', 'Lpa', 'Rpa'}, f'bad labels: {labels}'\n",
+    "display(landmarks)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d495ae12-8920-4294-aae1-4eb83f205403",
+   "metadata": {},
+   "source": [
+    "## Anonymize\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "acfa6540-281e-4123-8a0d-ee30cd33aae1",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Subject13: 943,539 -> 590,551 verts (-352,988)\n"
+     ]
+    }
+   ],
+   "source": [
+    "surface_anon, landmarks_anon = anonymize_scan(\n",
+    "    surface, landmarks,\n",
+    "    ear_delete_radius_mm=EAR_DELETE_RADIUS_MM,\n",
+    "    landmark_keep_radius_mm=LANDMARK_KEEP_RADIUS_MM,\n",
+    "    eyebrow_offset_mm=EYEBROW_OFFSET_MM,\n",
+    "    cap_z_ceiling_mm=CAP_Z_CEILING_MM,\n",
+    "    head_isolation_radius_mm=HEAD_ISOLATION_RADIUS_MM,\n",
+    ")\n",
+    "stem = Path(SCAN_PATH).stem\n",
+    "print(f'{stem}: {surface.nvertices:,} -> {surface_anon.nvertices:,} verts '\n",
+    "      f'(-{surface.nvertices - surface_anon.nvertices:,})')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "74a83f75-488b-48c9-a660-e5b38eeaf384",
+   "metadata": {},
+   "source": [
+    "## Compare original vs anonymized\n",
+    "\n",
+    "Both meshes live in the `digitized` frame, so the side-by-side view shares one\n",
+    "coordinate system. The same five landmark spheres are overlaid on each.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "309c17ba-a234-438e-bd5a-3372dc328266",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "687724dd208b4edca3d05ca2c1e3db15",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Widget(value='<iframe src=\"http://localhost:44253/index.html?ui=P_0x7fe109fd7890_3&reconnect=auto\" class=\"pyvi…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "stem = Path(SCAN_PATH).stem\n",
+    "lm_colors = {\n",
+    "    'Nz': 'lime', 'Iz': 'magenta', 'Cz': 'cyan',\n",
+    "    'Lpa': 'orange', 'Rpa': 'blue', 'LPA': 'orange', 'RPA': 'blue',\n",
+    "}\n",
+    "n_removed = surface.nvertices - surface_anon.nvertices\n",
+    "\n",
+    "orig_vtk = trimesh_to_vtk_polydata(surface.mesh)\n",
+    "anon_vtk = trimesh_to_vtk_polydata(surface_anon.mesh)\n",
+    "\n",
+    "lm_pos = landmarks_anon.pint.dequantify().values\n",
+    "lm_lbls = landmarks_anon['label'].values\n",
+    "\n",
+    "pvplt = pv.Plotter(shape=(1, 2))\n",
+    "\n",
+    "pvplt.subplot(0, 0)\n",
+    "pvplt.add_mesh(pv.wrap(orig_vtk), rgb=True, smooth_shading=True)\n",
+    "for lbl, pos in zip(lm_lbls, lm_pos):\n",
+    "    pvplt.add_mesh(pv.Sphere(radius=4, center=pos), color=lm_colors.get(lbl, 'yellow'))\n",
+    "pvplt.add_text(f'{stem} Original', position='upper_left', font_size=14)\n",
+    "\n",
+    "pvplt.subplot(0, 1)\n",
+    "pvplt.add_mesh(pv.wrap(anon_vtk), rgb=True, smooth_shading=True)\n",
+    "for lbl, pos in zip(lm_lbls, lm_pos):\n",
+    "    pvplt.add_mesh(pv.Sphere(radius=4, center=pos), color=lm_colors.get(lbl, 'yellow'))\n",
+    "pvplt.add_text(\n",
+    "    f'{stem} Anonymized (-{n_removed:,} verts)',\n",
+    "    position='upper_left', font_size=14,\n",
+    ")\n",
+    "\n",
+    "pvplt.link_views()\n",
+    "pvplt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "75e40eed-bccc-4eb5-a4d5-9fd1c64f0d1b",
+   "metadata": {},
+   "source": [
+    "## Save\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "ac270274-9e08-4e50-90d2-65c4c49aa34c",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "wrote /tmp/Subject13_notebook_test_anon.jpg\n",
+      "wrote /tmp/Subject13_notebook_test_anon.mtl\n",
+      "wrote /tmp/Subject13_notebook_test_anon.obj\n"
+     ]
+    }
+   ],
+   "source": [
+    "written = save_anonymized_scan(surface_anon, OUT_PATH)\n",
+    "for p in written:\n",
+    "    print(f'wrote {p}')\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5d9540ae-a939-4b18-be01-321627b8ad49",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
index 224c6a07..f8300ed2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,7 +61,6 @@ dependencies = [
   "seedir",
 ]
 
-
 [project.urls]
 Home = "https://ibs-lab.com/cedalion/"
 Documentation = "https://doc.ibs.tu-berlin.de/cedalion/doc/dev/"
diff --git a/src/cedalion/geometry/photogrammetry/anonymization/__init__.py b/src/cedalion/geometry/photogrammetry/anonymization/__init__.py
new file mode 100644
index 00000000..2abe8a63
--- /dev/null
+++ b/src/cedalion/geometry/photogrammetry/anonymization/__init__.py
@@ -0,0 +1,65 @@
+"""Face anonymization module for photogrammetry scans.
+
+Anonymizes facial regions in 3D photogrammetry scans while preserving
+optode positions and anatomical landmarks for fNIRS research.
+
+The canonical entry point is :func:`anonymize_scan`. Pass it a raw
+Einstar surface and the 5 landmarks (Nz, Iz, Cz, LPA, RPA) and it returns
+the anonymized surface plus landmarks, ready for
+:func:`save_anonymized_scan`.
+
+Example:
+    >>> from cedalion.geometry.photogrammetry.anonymization import (
+    ...     anonymize_scan, save_anonymized_scan,
+    ... )
+    >>> surface_anon, landmarks_anon = anonymize_scan(surface, landmarks)
+    >>> save_anonymized_scan(surface_anon, "out.obj")
+
+Pipeline steps inside ``anonymize_scan``:
+
+1. :func:`normalize_axes`: rotate around X so Y points anterior.
+2. :func:`isolate_head`: strip body, shoulders, and disconnected fragments.
+3. :func:`align_axes_from_landmarks`: map into the CTF frame.
+4. :func:`detect_cap_boundary`: find the cap front edge along Z.
+5. :func:`face_mask_from_landmarks`: face region + ear spheres, clamped
+   below the cap.
+6. Preserve small spheres around each landmark and a midline nasion strip.
+7. :func:`delete_masked_vertices`: drop triangles touching any masked vertex.
+8. :func:`revert_to_einstar_frame`: return to ``crs="digitized"`` so the
+   output matches :func:`cedalion.io.read_einstar_obj`.
+
+Each of those functions is also exported so callers that need to inspect
+or override an intermediate step can do so without re-implementing the
+whole pipeline.
+"""
+
+from .preprocessing import (
+    normalize_axes,
+    isolate_head,
+    align_axes_from_landmarks,
+    revert_to_einstar_frame,
+)
+from .mask import (
+    detect_cap_boundary,
+    face_mask_from_landmarks,
+    delete_masked_vertices,
+    save_anonymized_scan,
+)
+from .pipeline import anonymize_scan
+
+
+__all__ = [
+    # Top-level orchestrator (canonical entry point)
+    "anonymize_scan",
+    # Preprocessing (axis normalization, head isolation, full alignment,
+    # and the inverse mapping back to the raw Einstar frame)
+    "normalize_axes",
+    "isolate_head",
+    "align_axes_from_landmarks",
+    "revert_to_einstar_frame",
+    # Mask construction and application
+    "detect_cap_boundary",
+    "face_mask_from_landmarks",
+    "delete_masked_vertices",
+    "save_anonymized_scan",
+]
diff --git a/src/cedalion/geometry/photogrammetry/anonymization/_utils.py b/src/cedalion/geometry/photogrammetry/anonymization/_utils.py
new file mode 100644
index 00000000..a8845fb4
--- /dev/null
+++ b/src/cedalion/geometry/photogrammetry/anonymization/_utils.py
@@ -0,0 +1,151 @@
+"""Shared helpers used by ``preprocessing`` and ``mask``.
+
+Consolidates texture-image lookup, mesh rebuilding, face/vertex reindexing,
+affine application, landmark transforms, ear midpoint, and upper-head centroid
+to avoid duplicating these across sibling modules.
+"""
+
+from typing import Callable
+
+import numpy as np
+import trimesh
+
+import cedalion.typing as cdt
+
+
+def _resolve_texture_image(visual):
+    """Return ``visual.image`` if set, else fall back to ``material.image``."""
+    image = getattr(visual, "image", None)
+    if image is not None:
+        return image
+    mat = getattr(visual, "material", None)
+    return getattr(mat, "image", None) if mat is not None else None
+
+
+def _reindex_faces(
+    vertices: np.ndarray,
+    faces: np.ndarray,
+    kept_face_mask: np.ndarray,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Slice a triangle mesh by a face mask and rebuild dense vertex indices.
+
+    Args:
+        vertices: ``(N, 3)`` vertex array.
+        faces: ``(F, 3)`` face index array.
+        kept_face_mask: ``(F,)`` bool array.
+
+    Returns:
+        ``(new_vertices, new_faces, kept_vidx)`` where ``new_faces`` indexes
+        into ``new_vertices`` and ``kept_vidx`` are the original vertex
+        indices that survived (useful for slicing UVs / vertex colors).
+    """
+    kept_faces = faces[kept_face_mask]
+    kept_vidx = np.unique(kept_faces)
+    reindex = -np.ones(len(vertices), dtype=np.int64)
+    reindex[kept_vidx] = np.arange(len(kept_vidx))
+    new_faces = reindex[kept_faces]
+    new_vertices = vertices[kept_vidx]
+    return new_vertices, new_faces, kept_vidx
+
+
+def _apply_affine(points: np.ndarray, M: np.ndarray) -> np.ndarray:
+    """Apply a 4x4 affine to an ``(N, 3)`` point array."""
+    return points @ M[:3, :3].T + M[:3, 3]
+
+
+def _transform_labeled_points(
+    landmarks: cdt.LabeledPoints,
+    transform_xyz: Callable[[np.ndarray], np.ndarray],
+    new_crs: str,
+) -> cdt.LabeledPoints:
+    """Transform a LabeledPoints array and rename its spatial dim.
+
+    Wraps the dequantify -> math -> ``copy(data=...)`` -> rename -> quantify
+    cycle that ``align_axes_from_landmarks`` and ``revert_to_einstar_frame``
+    both need.
+    """
+    dequant = landmarks.pint.dequantify()
+    new_xyz = transform_xyz(dequant.values)
+    old_crs_dim = next(d for d in landmarks.dims if d != "label")
+    return (
+        dequant.copy(data=new_xyz)
+        .rename({old_crs_dim: new_crs})
+        .pint.quantify()
+    )
+
+
+def _ear_midpoint(Lpa: np.ndarray, Rpa: np.ndarray) -> np.ndarray:
+    """Return the LPA/RPA midpoint (CTF origin in the aligned frame)."""
+    return 0.5 * (Lpa + Rpa)
+
+
+def _upper_head_centroid(
+    vertices: np.ndarray,
+) -> tuple[np.ndarray, float]:
+    """Centroid of vertices in the top 40% of the X range.
+
+    Returns ``(centroid_xyz, x_max)``. Callers that only need the YZ
+    components zero out X themselves.
+    """
+    x_max = float(vertices[:, 0].max())
+    x_min = float(vertices[:, 0].min())
+    upper = vertices[vertices[:, 0] > x_min + 0.6 * (x_max - x_min)]
+    centroid = upper.mean(axis=0)
+    return centroid, x_max
+
+
+def _copy_visual(src_mesh, dst_mesh, vertex_index=None) -> None:
+    """Copy a trimesh ``visual`` onto ``dst_mesh``, optionally reindexing.
+
+    Passing ``visual=old.visual`` to a new ``Trimesh(...)`` silently downgrades
+    a ``TextureVisuals`` to ``ColorVisuals`` because the visual still
+    back-references the source mesh's vertex count. Rebuilding it from scratch
+    preserves TextureVisuals.
+
+    Args:
+        src_mesh: Source trimesh.
+        dst_mesh: Destination trimesh (mutated in place).
+        vertex_index: Indices into the source vertex array. If given, UVs /
+            vertex_colors are sliced by this index. Leave None when the
+            vertex count is unchanged.
+    """
+    src_visual = src_mesh.visual
+    uv = getattr(src_visual, "uv", None)
+    n_src = len(src_mesh.vertices)
+
+    if uv is not None and len(uv) == n_src:
+        uv_arr = np.asarray(uv)
+        if vertex_index is not None:
+            uv_arr = uv_arr[vertex_index]
+        image = _resolve_texture_image(src_visual)
+        material = getattr(src_visual, "material", None)
+        dst_mesh.visual = trimesh.visual.TextureVisuals(
+            uv=uv_arr,
+            image=image,
+            material=material,
+        )
+        return
+
+    vcol = getattr(src_visual, "vertex_colors", None)
+    if vcol is not None and len(vcol) == n_src:
+        vcol_arr = np.asarray(vcol)
+        if vertex_index is not None:
+            vcol_arr = vcol_arr[vertex_index]
+        dst_mesh.visual.vertex_colors = vcol_arr
+
+
+def _rebuild_mesh(
+    src_mesh: trimesh.Trimesh,
+    *,
+    vertices: np.ndarray,
+    faces: np.ndarray,
+    vertex_index: np.ndarray | None = None,
+) -> trimesh.Trimesh:
+    """Build a new ``Trimesh`` and re-attach the source visual.
+
+    Used by every transform in the pipeline to keep the
+    ``Trimesh(process=False)`` + ``_copy_visual`` pattern in one place.
+    """
+    new_mesh = trimesh.Trimesh(vertices=vertices, faces=faces, process=False)
+    _copy_visual(src_mesh, new_mesh, vertex_index=vertex_index)
+    return new_mesh
diff --git a/src/cedalion/geometry/photogrammetry/anonymization/mask.py b/src/cedalion/geometry/photogrammetry/anonymization/mask.py
new file mode 100644
index 00000000..f55975ae
--- /dev/null
+++ b/src/cedalion/geometry/photogrammetry/anonymization/mask.py
@@ -0,0 +1,376 @@
+"""Facial-region mask construction and application.
+
+Given the 5 CTF-aligned landmarks, builds a boolean deletion mask and
+applies it to the mesh:
+
+- ``detect_cap_boundary``: finds the Z height where the EEG cap front edge
+  sits so the mask can be clamped below it.
+- ``face_mask_from_landmarks``: unions a forward face region with two ear
+  spheres, clamped below the cap.
+- ``delete_masked_vertices``: drops triangles touching any masked vertex.
+
+All mask math assumes the CTF frame: +X=anterior, +Y=left, +Z=up, origin
+at the LPA-RPA midpoint (see ``align_axes_from_landmarks``).
+"""
+
+import logging
+import os
+
+import numpy as np
+import trimesh
+from PIL import Image, ImageDraw, ImageFilter
+from scipy.signal import savgol_filter
+
+import cedalion.dataclasses as cdc
+
+from ._utils import _ear_midpoint, _rebuild_mesh, _reindex_faces, _resolve_texture_image
+
+logger = logging.getLogger("cedalion")
+
+
+def detect_cap_boundary(
+    verts: np.ndarray,
+    Nz: np.ndarray,
+    Cz: np.ndarray,
+    Lpa: np.ndarray,
+    Rpa: np.ndarray,
+    band_width: float = 15.0,
+    bin_size: float = 1.0,
+    foot_grad_threshold: float = 0.2,
+    cap_z_ceiling_mm: float = 40.0,
+    eyebrow_offset_mm: float = 10.0,
+) -> tuple[float, np.ndarray, np.ndarray, np.ndarray]:
+    """Find the Z height where the EEG cap front edge sits.
+
+    Scans upward from Nz along the midline and records max-X per Z-bin. The
+    cap protrudes anteriorly, so X-max on the midline traces the cap. The
+    cap edge is the foot of the rise leading to the peak: walk back from the
+    peak until the smoothed gradient drops below ``foot_grad_threshold``.
+
+    Failsafe for flush / no-optode caps: when the cap sits flat against the
+    head, the anterior bump vanishes and X-max is roughly monotonic up to
+    ``Cz``, stranding ``cap_z`` near the crown. If detection lands above
+    ``Nz[2] + cap_z_ceiling_mm``, fall back to ``Nz[2] + eyebrow_offset_mm``
+    (just above the supraorbital ridge).
+
+    Expects the CTF frame (+X=anterior, +Y=left, +Z=up).
+
+    Args:
+        verts: Mesh vertices, shape (N, 3).
+        Nz: Nasion position.
+        Cz: Cz position.
+        Lpa: Left preauricular position.
+        Rpa: Right preauricular position.
+        band_width: Y-band half-width for the midline X-profile (mm).
+        bin_size: Z-bin size for the X-profile (mm).
+        foot_grad_threshold: dX/dZ below this value marks the foot of the rise.
+        cap_z_ceiling_mm: Absolute mm above Nz at which the cap-peak
+            detection is considered untrustworthy and the failsafe fires.
+        eyebrow_offset_mm: Failsafe cut height, expressed as mm above Nz.
+            Anatomically just above the supraorbital ridge.
+
+    Returns:
+        Tuple of (cap_z, profile_z, profile_x_raw, profile_x_smooth).
+    """
+    ear_mid = _ear_midpoint(Lpa, Rpa)
+    mid_y = ear_mid[1]
+    in_band = np.abs(verts[:, 1] - mid_y) < band_width
+    above_nz = verts[:, 2] > Nz[2]
+    anterior = verts[:, 0] > ear_mid[0]
+    sel_verts = verts[in_band & above_nz & anterior]
+
+    bins = np.arange(Nz[2], Cz[2], bin_size)
+    bin_centers = bins[:-1] + bin_size / 2
+    max_x = np.full(len(bin_centers), np.nan)
+
+    for i in range(len(bin_centers)):
+        in_bin = (sel_verts[:, 2] >= bins[i]) & (sel_verts[:, 2] < bins[i + 1])
+        if in_bin.any():
+            max_x[i] = sel_verts[in_bin, 0].max()
+
+    valid = ~np.isnan(max_x)
+    if valid.sum() < 7:
+        fallback = 0.5 * (Nz[2] + Cz[2])
+        return fallback, bin_centers[valid], max_x[valid], max_x[valid]
+
+    zv = bin_centers[valid]
+    xv = max_x[valid]
+
+    win = min(11, len(xv) if len(xv) % 2 == 1 else len(xv) - 1)
+    xv_s = savgol_filter(xv, window_length=win, polyorder=2)
+
+    peak_idx = int(np.argmax(xv_s))
+    grad = np.gradient(xv_s, zv)
+    cap_z = zv[0]
+    for i in range(peak_idx - 1, 0, -1):
+        if grad[i] < foot_grad_threshold:
+            cap_z = zv[i]
+            break
+
+    ceiling = Nz[2] + cap_z_ceiling_mm
+    if cap_z > ceiling:
+        fallback = Nz[2] + eyebrow_offset_mm
+        logger.info(
+            f"detect_cap_boundary: cap_z={cap_z:.1f} mm exceeded ceiling "
+            f"Nz+{cap_z_ceiling_mm:.0f}={ceiling:.1f}; assuming flush cap "
+            f"and falling back to Nz+{eyebrow_offset_mm:.0f}={fallback:.1f} mm."
+        )
+        cap_z = fallback
+
+    return cap_z, zv, xv, xv_s
+
+
+def face_mask_from_landmarks(
+    verts: np.ndarray,
+    Nz: np.ndarray,
+    Lpa: np.ndarray,
+    Rpa: np.ndarray,
+    cap_z: float | None = None,
+    ear_delete_radius: float = 40.0,
+) -> tuple[np.ndarray, dict]:
+    """Build face + ear deletion mask from the landmarks (CTF frame).
+
+    Mask is the union of two regions, both clamped below the cap boundary:
+
+    1. Face region: anterior to the ear coronal plane (X > ear_mid_X).
+    2. Ear spheres: ``ear_delete_radius`` mm around Lpa and Rpa.
+
+    Expects the CTF frame (+X=anterior, +Y=left, +Z=up).
+
+    Args:
+        verts: Mesh vertices, shape (N, 3).
+        Nz: Nasion position in the CTF frame.
+        Lpa: Left preauricular position in the CTF frame.
+        Rpa: Right preauricular position in the CTF frame.
+        cap_z: Upper bound Z value (typically from ``detect_cap_boundary``).
+            Defaults to Nz[2] if not provided.
+        ear_delete_radius: Sphere radius around Lpa/Rpa in mm.
+
+    Returns:
+        Tuple of (mask, info). ``mask`` is a boolean array of shape (N,).
+        ``info`` has keys ``upper_bound``, ``ear_mid``, and ``counts``
+        (per-region vertex counts).
+    """
+    ear_mid = _ear_midpoint(Lpa, Rpa)
+    upper_bound = cap_z if cap_z is not None else Nz[2]
+
+    below_cap = verts[:, 2] < upper_bound
+    anterior = verts[:, 0] > ear_mid[0]
+    face_region = below_cap & anterior
+
+    d_lpa = np.linalg.norm(verts - Lpa, axis=1)
+    d_rpa = np.linalg.norm(verts - Rpa, axis=1)
+    ear_region = ((d_lpa < ear_delete_radius) | (d_rpa < ear_delete_radius)) & below_cap
+
+    mask = face_region | ear_region
+
+    info = {
+        "upper_bound": upper_bound,
+        "ear_mid": ear_mid,
+        "counts": {
+            "below_cap": int(below_cap.sum()),
+            "face_region": int(face_region.sum()),
+            "ear_region": int(ear_region.sum()),
+            "all": int(mask.sum()),
+        },
+    }
+    return mask, info
+
+
+@cdc.validate_schemas
+def delete_masked_vertices(
+    surface: cdc.TrimeshSurface,
+    mask: np.ndarray,
+) -> cdc.TrimeshSurface:
+    """Drop triangles touching any masked vertex and strip unreferenced vertices.
+
+    Reindexes ``mesh.visual.uv`` in lockstep with the vertex array.
+    ``Trimesh.remove_unreferenced_vertices`` does not touch UV arrays, so a
+    naive call leaves ``len(uv) != len(vertices)`` and breaks textured export.
+
+    Args:
+        surface: Input TrimeshSurface.
+        mask: Boolean array of shape (n_vertices,). True = vertex to remove.
+
+    Returns:
+        New TrimeshSurface with masked vertices (and the faces touching them)
+        removed. CRS, units, UVs, and texture image are preserved in sync.
+    """
+    old_mesh = surface.mesh
+    old_verts = np.asarray(old_mesh.vertices)
+    old_faces = np.asarray(old_mesh.faces)
+
+    kept_face_mask = ~mask[old_faces].any(axis=1)
+    new_verts, new_faces, kept_vidx = _reindex_faces(
+        old_verts, old_faces, kept_face_mask
+    )
+
+    new_mesh = _rebuild_mesh(
+        old_mesh,
+        vertices=new_verts,
+        faces=new_faces,
+        vertex_index=kept_vidx,
+    )
+
+    return cdc.TrimeshSurface(
+        mesh=new_mesh, crs=surface.crs, units=surface.units,
+    )
+
+
+def _sanitize_texture_from_uv(mesh: trimesh.Trimesh) -> Image.Image | None:
+    """Rebuild the texture image keeping only pixels referenced by mesh UVs.
+
+    Rasterizes every surviving face's UV triangle onto a keep-mask, then
+    blacks out pixels outside the mask so face-region pixels cannot leak
+    through the JPG. The mask is dilated by 1 pixel to absorb anti-aliasing
+    fringes at triangle boundaries.
+
+    OBJ UVs use a bottom-left origin; PIL images use top-left, so V is
+    flipped before rasterizing.
+
+    Args:
+        mesh: A ``trimesh.Trimesh`` with ``mesh.visual`` as ``TextureVisuals``.
+
+    Returns:
+        Sanitized PIL Image in RGB mode, or ``None`` if the mesh has no
+        usable texture (caller should fall back to geometry-only output).
+    """
+    visual = mesh.visual
+    uv = getattr(visual, "uv", None)
+    image = _resolve_texture_image(visual)
+
+    if uv is None or image is None or len(uv) != len(mesh.vertices):
+        return None
+
+    img_rgb = np.asarray(image.convert("RGB"))
+    h, w = img_rgb.shape[:2]
+
+    uv = np.asarray(uv)
+    px = np.mod(uv[:, 0], 1.0) * w
+    py = (1.0 - np.mod(uv[:, 1], 1.0)) * h
+
+    mask_img = Image.new("L", (w, h), 0)
+    draw = ImageDraw.Draw(mask_img)
+    for face in np.asarray(mesh.faces):
+        pts = [(px[i], py[i]) for i in face]
+        draw.polygon(pts, fill=255)
+    mask_img = mask_img.filter(ImageFilter.MaxFilter(3))
+    keep = np.asarray(mask_img) > 0
+
+    out = np.empty_like(img_rgb)
+    out[keep] = img_rgb[keep]
+    out[~keep] = (0, 0, 0)
+    return Image.fromarray(out, mode="RGB")
+
+
+_MTL_TEMPLATE = (
+    "newmtl _texture\n"
+    "Kd 1.00 1.00 1.00\n"
+    "Ka 0.00 0.00 0.00\n"
+    "Tf 1.00 1.00 1.00\n"
+    "Ni 1.00\n"
+    "map_Kd {jpg_name}\n"
+)
+
+
+@cdc.validate_schemas
+def save_anonymized_scan(
+    surface: cdc.TrimeshSurface,
+    out_path: str,
+    strip_texture: bool = False,
+) -> list[str]:
+    """Export an anonymized photogrammetry surface to disk.
+
+    Default path (``strip_texture=False``): writes an ``.obj`` + ``.mtl`` +
+    sanitized ``.jpg`` bundle. The JPG is rebuilt to contain colors *only*
+    for UV regions still referenced by the anonymized mesh; face-region
+    pixels are replaced by the fill color, so opening the JPG alone reveals
+    no face.
+
+    ``strip_texture=True`` writes geometry only (no MTL, no JPG). Same
+    fallback applies when the input mesh has no usable texture, with a
+    warning.
+
+    Args:
+        surface: Anonymized TrimeshSurface (typically output of
+            ``delete_masked_vertices``).
+        out_path: Destination path ending in ``.obj``.
+        strip_texture: If True, skip the MTL + JPG and write geometry only.
+
+    Returns:
+        List of absolute paths written (``.obj`` plus ``.mtl`` + ``.jpg``
+        when a texture was written).
+
+    Raises:
+        ValueError: If ``out_path`` does not end in ``.obj``.
+    """
+    if not out_path.lower().endswith(".obj"):
+        raise ValueError(f"out_path must end in .obj, got: {out_path}")
+
+    out_dir = os.path.dirname(os.path.abspath(out_path)) or "."
+    stem = os.path.splitext(os.path.basename(out_path))[0]
+
+    sanitized = None if strip_texture else _sanitize_texture_from_uv(surface.mesh)
+
+    # Bare rebuild: the visual is rewritten below (with the sanitized
+    # texture) or the mesh is exported geometry-only, so we deliberately
+    # skip ``_rebuild_mesh``'s visual-copy step.
+    mesh_to_write = trimesh.Trimesh(
+        vertices=np.asarray(surface.mesh.vertices),
+        faces=np.asarray(surface.mesh.faces),
+        process=False,
+    )
+
+    written: list[str] = []
+
+    if sanitized is not None:
+        jpg_name = f"{stem}.jpg"
+        mtl_name = f"{stem}.mtl"
+        jpg_path = os.path.join(out_dir, jpg_name)
+        mtl_path = os.path.join(out_dir, mtl_name)
+
+        sanitized.save(jpg_path, format="JPEG", quality=92)
+        with open(mtl_path, "w") as fh:
+            fh.write(_MTL_TEMPLATE.format(jpg_name=jpg_name))
+
+        mesh_to_write.visual = trimesh.visual.TextureVisuals(
+            uv=np.asarray(surface.mesh.visual.uv),
+            image=sanitized,
+        )
+        obj_text = trimesh.exchange.obj.export_obj(
+            mesh_to_write,
+            include_texture=True,
+            write_texture=False,
+            mtl_name=mtl_name,
+        )
+        obj_text = _ensure_mtllib(obj_text, mtl_name)
+        with open(out_path, "w") as fh:
+            fh.write(obj_text)
+
+        written.extend([out_path, mtl_path, jpg_path])
+    else:
+        if not strip_texture:
+            logger.warning(
+                "save_anonymized_scan: input mesh has no usable texture; "
+                "falling back to geometry-only OBJ."
+            )
+        mesh_to_write.export(out_path)
+        written.append(out_path)
+
+    logger.info(
+        f"Saved anonymized scan: "
+        f"{[os.path.basename(p) for p in written]}"
+    )
+    return sorted(written)
+
+
+def _ensure_mtllib(obj_text: str, mtl_name: str) -> str:
+    """Guarantee the OBJ references our MTL via ``mtllib`` and ``usemtl``.
+
+    Trimesh's OBJ exporter emits material lines inconsistently across
+    versions; rewrite deterministically so the MTL pairing is stable.
+    """
+    lines = obj_text.splitlines()
+    lines = [line for line in lines if not line.startswith(("mtllib ", "usemtl "))]
+    header = [f"mtllib {mtl_name}", "usemtl _texture"]
+    return "\n".join(header + lines) + "\n"
diff --git a/src/cedalion/geometry/photogrammetry/anonymization/pipeline.py b/src/cedalion/geometry/photogrammetry/anonymization/pipeline.py
new file mode 100644
index 00000000..66a1ee8b
--- /dev/null
+++ b/src/cedalion/geometry/photogrammetry/anonymization/pipeline.py
@@ -0,0 +1,159 @@
+"""Top-level orchestrator: ``anonymize_scan``.
+
+Single entry point that chains the helper functions from ``preprocessing``
+and ``mask``. The individual functions remain public so downstream code can
+inspect intermediate state (e.g. the cap-detection profile, the head
+isolation mask), but the canonical pipeline lives here and only here.
+"""
+
+from typing import Literal
+
+import numpy as np
+
+import cedalion.dataclasses as cdc
+import cedalion.typing as cdt
+from cedalion.geometry.landmarks import normalize_landmarks_labels
+
+from ._utils import _ear_midpoint, _apply_affine, _transform_labeled_points
+from .mask import (
+    delete_masked_vertices,
+    detect_cap_boundary,
+    face_mask_from_landmarks,
+)
+from .preprocessing import (
+    align_axes_from_landmarks,
+    isolate_head,
+    normalize_axes,
+    revert_to_einstar_frame,
+)
+
+
+_REQUIRED_LABELS = ("Nz", "Iz", "Cz", "LPA", "RPA")
+
+
+@cdc.validate_schemas
+def anonymize_scan(
+    surface: cdc.TrimeshSurface,
+    landmarks: cdt.LabeledPoints,
+    *,
+    head_isolation_radius_mm: float = 220.0,
+    cap_band_width_mm: float = 15.0,
+    cap_bin_size_mm: float = 1.0,
+    cap_foot_grad_threshold: float = 0.2,
+    cap_z_ceiling_mm: float = 40.0,
+    eyebrow_offset_mm: float = 10.0,
+    ear_delete_radius_mm: float = 40.0,
+    landmark_keep_radius_mm: float = 8.0,
+    return_frame: Literal["digitized", "ctf"] = "digitized",
+) -> tuple[cdc.TrimeshSurface, cdt.LabeledPoints]:
+    """Run the full face-anonymization pipeline on a raw Einstar scan.
+
+    Steps (each is also exposed as a standalone function for inspection):
+
+    1. ``normalize_axes``: rotate around X so Y points anterior.
+    2. ``isolate_head``: strip body, shoulders, fragments.
+    3. ``align_axes_from_landmarks``: map into the CTF frame.
+    4. ``detect_cap_boundary``: find the cap front edge along Z.
+    5. ``face_mask_from_landmarks``: union face region + ear spheres,
+       clamped below the cap.
+    6. Preserve ``landmark_keep_radius_mm``-spheres around each landmark
+       and a midline nasion strip from Nz to the cap.
+    7. ``delete_masked_vertices``: drop triangles touching any masked
+       vertex, keeping UVs in sync.
+    8. (default) ``revert_to_einstar_frame``: return to ``crs="digitized"``
+       so the output matches ``read_einstar_obj``'s convention and can be
+       fed to ``save_anonymized_scan`` and downstream co-registration.
+
+    Args:
+        surface: Raw Einstar TrimeshSurface (``crs="digitized"``).
+        landmarks: LabeledPoints with Nz, Iz, Cz, LPA, RPA (mixed-case
+            aliases like Lpa/Rpa accepted; normalized via
+            ``normalize_landmarks_labels``).
+        head_isolation_radius_mm: Sphere radius around the upper-head
+            centroid for ``isolate_head``.
+        cap_band_width_mm: Y-band half-width for the cap X-profile.
+        cap_bin_size_mm: Z-bin size for the cap X-profile.
+        cap_foot_grad_threshold: dX/dZ below which the cap-foot is
+            recognized.
+        cap_z_ceiling_mm: mm above Nz at which cap detection is
+            considered untrustworthy and the failsafe fires.
+        eyebrow_offset_mm: failsafe cap height above Nz.
+        ear_delete_radius_mm: sphere radius around LPA/RPA for the ear
+            region of the deletion mask.
+        landmark_keep_radius_mm: per-landmark preservation sphere radius
+            and half-width of the midline nasion strip.
+        return_frame: ``"digitized"`` (default) reverts back to the raw
+            Einstar frame; ``"ctf"`` keeps the CTF frame.
+
+    Returns:
+        Tuple of (anonymized_surface, anonymized_landmarks). Frame is
+        controlled by ``return_frame``. The surface can be written with
+        ``save_anonymized_scan`` when ``return_frame="digitized"``.
+    """
+    landmarks = normalize_landmarks_labels(landmarks)
+    labels = list(landmarks["label"].values)
+    missing = set(_REQUIRED_LABELS) - set(labels)
+    if missing:
+        raise ValueError(f"Missing landmarks for anonymization: {missing}")
+    idx = {lbl: i for i, lbl in enumerate(labels)}
+    Nz_raw = landmarks.pint.dequantify().values[idx["Nz"]]
+
+    surface_n, _, R_norm = normalize_axes(surface, Nz_raw)
+    R_norm4 = np.eye(4)
+    R_norm4[:3, :3] = R_norm
+    crs_dim = next(d for d in landmarks.dims if d != "label")
+    landmarks_n = _transform_labeled_points(
+        landmarks,
+        lambda p: _apply_affine(p, R_norm4),
+        new_crs=crs_dim,
+    )
+
+    Nz_n = landmarks_n.pint.dequantify().values[idx["Nz"]]
+    surface_n, _ = isolate_head(
+        surface_n, Nz_n, radius=head_isolation_radius_mm
+    )
+
+    surface_h, landmarks_ctf, M_align = align_axes_from_landmarks(
+        surface_n, landmarks_n
+    )
+    lm_ctf = landmarks_ctf.pint.dequantify().values
+    Nz, Iz, Cz, Lpa, Rpa = (
+        lm_ctf[idx[lbl]] for lbl in _REQUIRED_LABELS
+    )
+    ear_mid = _ear_midpoint(Lpa, Rpa)
+
+    verts = np.asarray(surface_h.mesh.vertices)
+    cap_z, *_ = detect_cap_boundary(
+        verts, Nz, Cz, Lpa, Rpa,
+        band_width=cap_band_width_mm,
+        bin_size=cap_bin_size_mm,
+        foot_grad_threshold=cap_foot_grad_threshold,
+        cap_z_ceiling_mm=cap_z_ceiling_mm,
+        eyebrow_offset_mm=eyebrow_offset_mm,
+    )
+
+    mask, _ = face_mask_from_landmarks(
+        verts, Nz, Lpa, Rpa,
+        cap_z=cap_z,
+        ear_delete_radius=ear_delete_radius_mm,
+    )
+
+    for lm in (Nz, Iz, Cz, Lpa, Rpa):
+        near = np.linalg.norm(verts - lm, axis=1) < landmark_keep_radius_mm
+        mask[near] = False
+    nasion_strip = (
+        (verts[:, 2] >= Nz[2])
+        & (verts[:, 2] < cap_z)
+        & (np.abs(verts[:, 1] - Nz[1]) < landmark_keep_radius_mm)
+        & (verts[:, 0] > ear_mid[0])
+    )
+    mask[nasion_strip] = False
+
+    surface_anon = delete_masked_vertices(surface_h, mask)
+
+    if return_frame == "ctf":
+        return surface_anon, landmarks_ctf
+
+    return revert_to_einstar_frame(
+        surface_anon, landmarks_ctf, R_norm, M_align
+    )
diff --git a/src/cedalion/geometry/photogrammetry/anonymization/preprocessing.py b/src/cedalion/geometry/photogrammetry/anonymization/preprocessing.py
new file mode 100644
index 00000000..db2f9d6e
--- /dev/null
+++ b/src/cedalion/geometry/photogrammetry/anonymization/preprocessing.py
@@ -0,0 +1,331 @@
+"""Preprocessing: axis normalization and head isolation.
+
+Gets a raw Einstar photogrammetry scan into a standard reference frame suitable for
+landmark detection and mask building:
+
+- ``normalize_axes`` rotates around X so Y points anterior (preliminary,
+  used before head isolation).
+- ``isolate_head`` strips body/shoulders/chair and disconnected fragments.
+- ``align_axes_from_landmarks`` maps the scan into the CTF frame
+  (+X=anterior, +Y=left, +Z=up, origin at the LPA-RPA midpoint) once all 5
+  landmarks are available.
+"""
+
+import logging
+
+import numpy as np
+import trimesh
+
+import cedalion.dataclasses as cdc
+import cedalion.typing as cdt
+from cedalion.geometry.landmarks import normalize_landmarks_labels
+
+from ._utils import (
+    _apply_affine,
+    _ear_midpoint,
+    _rebuild_mesh,
+    _reindex_faces,
+    _transform_labeled_points,
+    _upper_head_centroid,
+)
+
+logger = logging.getLogger("cedalion")
+
+
+@cdc.validate_schemas
+def normalize_axes(
+    surface: cdc.TrimeshSurface,
+    nasion: np.ndarray,
+) -> tuple[cdc.TrimeshSurface, np.ndarray, np.ndarray]:
+    """Rotate mesh around X-axis so Y=anterior (toward face), Z=left.
+
+    The Einstar scanner fixes X=up (gravity-based) but leaves Y/Z arbitrary
+    per scan. The forward direction is inferred from the nasion against the
+    upper-head centroid.
+
+    Args:
+        surface: TrimeshSurface in raw Einstar coordinates.
+        nasion: Nasion position as numpy array of shape (3,), in mm, in the
+            raw Einstar frame (matching ``surface``).
+
+    Returns:
+        Tuple of (rotated_surface, rotated_nasion, rotation_matrix).
+        ``rotated_nasion`` is in the same units as the input. The 3x3
+        ``rotation_matrix`` can be applied to other points via ``p @ R.T``.
+    """
+    vertices = surface.mesh.vertices
+
+    centroid, _ = _upper_head_centroid(np.asarray(vertices))
+    angle = np.arctan2(nasion[2] - centroid[2], nasion[1] - centroid[1])
+    cos_a = np.cos(-angle)
+    sin_a = np.sin(-angle)
+    R = np.array([
+        [1.0, 0.0, 0.0],
+        [0.0, cos_a, -sin_a],
+        [0.0, sin_a, cos_a],
+    ])
+
+    new_mesh = _rebuild_mesh(
+        surface.mesh,
+        vertices=vertices @ R.T,
+        faces=surface.mesh.faces,
+    )
+    rotated_surface = cdc.TrimeshSurface(new_mesh, crs=surface.crs, units=surface.units)
+
+    rotated_nasion = R @ nasion
+
+    logger.debug(
+        f"Axis normalization: rotated {np.degrees(angle):.1f}deg around X. "
+        f"Y now points anterior."
+    )
+
+    return rotated_surface, rotated_nasion, R
+
+
+def _largest_component_mask(mesh: trimesh.Trimesh) -> np.ndarray:
+    """Boolean vertex mask selecting the largest connected component.
+
+    Einstar scans contain floating fragments (loose triangles, cable shreds,
+    background patches) that drag vertex extrema off into empty space; strip
+    them first.
+
+    Uses ``trimesh.graph.connected_component_labels`` on face adjacency rather
+    than ``mesh.split``: split allocates a ``Trimesh`` per component and can
+    OOM on scans with thousands of fragments.
+
+    Einstar OBJs duplicate vertices along UV seams; ``Trimesh.merge_vertices``
+    preserves seams (to keep textures), so face adjacency on the raw mesh
+    over-fragments the head. ``trimesh.grouping.unique_rows`` does
+    position-only merging for connectivity analysis.
+
+    Args:
+        mesh: A ``trimesh.Trimesh`` instance.
+
+    Returns:
+        Boolean array of shape ``(n_vertices,)``. All-True if the mesh is
+        empty or already a single connected component.
+    """
+    n_verts = len(mesh.vertices)
+    n_faces = len(mesh.faces)
+    if n_faces == 0:
+        return np.ones(n_verts, dtype=bool)
+
+    _, inverse = trimesh.grouping.unique_rows(np.asarray(mesh.vertices))
+    canonical_faces = inverse[mesh.faces]
+    adjacency = trimesh.graph.face_adjacency(faces=canonical_faces)
+    face_labels = trimesh.graph.connected_component_labels(
+        adjacency, node_count=n_faces
+    )
+    counts = np.bincount(face_labels)
+    if len(counts) <= 1:
+        return np.ones(n_verts, dtype=bool)
+    biggest_label = int(np.argmax(counts))
+    face_mask = face_labels == biggest_label
+
+    kept_vidx = np.unique(mesh.faces[face_mask])
+    mask = np.zeros(n_verts, dtype=bool)
+    mask[kept_vidx] = True
+    return mask
+
+
+@cdc.validate_schemas
+def isolate_head(
+    surface: cdc.TrimeshSurface,
+    nasion: np.ndarray,
+    radius: float = 220.0,
+) -> tuple[cdc.TrimeshSurface, np.ndarray]:
+    """Remove shoulders, body, and chair; keep only the head.
+
+    Uses a sphere centered on the upper-head centroid. The radius is
+    chosen to capture the full head (~180mm wide x 230mm tall) while
+    excluding shoulders and body. Scans that are already head-only are
+    returned unchanged (the sphere just contains everything).
+
+    The surface must be axis-normalized first (X=up, Y=anterior, Z=left).
+
+    Args:
+        surface: Axis-normalized TrimeshSurface.
+        nasion: Nasion position as numpy array of shape (3,), in mm,
+            in the axis-normalized frame (matching ``surface``).
+        radius: Sphere radius in mm (default 220). A human head has
+            ~90mm radius; 220mm adds margin for ears and jaw.
+
+    Returns:
+        Tuple of (head_surface, head_mask). head_mask is a boolean
+        array of shape (n_vertices,) indicating which original
+        vertices were kept.
+    """
+    vertices = np.asarray(surface.mesh.vertices)
+    faces = np.asarray(surface.mesh.faces)
+
+    center, x_max = _upper_head_centroid(vertices)
+    midpoint_x = (x_max + nasion[0]) / 2.0
+    center[0] = min(center[0], midpoint_x)
+
+    dist = np.linalg.norm(vertices - center, axis=1)
+    head_mask = dist < radius
+
+    head_mask = head_mask & _largest_component_mask(surface.mesh)
+
+    if head_mask.sum() < 100 or head_mask.mean() > 0.95:
+        logger.debug(
+            f"Head isolation: sphere contains {head_mask.mean()*100:.0f}% "
+            f"of vertices -- scan is already head-only"
+        )
+        return surface, head_mask
+
+    face_mask = head_mask[faces].all(axis=1)
+    new_verts, new_faces, kept_vidx = _reindex_faces(vertices, faces, face_mask)
+
+    new_mesh = _rebuild_mesh(
+        surface.mesh,
+        vertices=new_verts,
+        faces=new_faces,
+        vertex_index=kept_vidx,
+    )
+
+    head_surface = cdc.TrimeshSurface(
+        new_mesh, crs=surface.crs, units=surface.units
+    )
+
+    logger.debug(
+        f"Head isolation: {len(vertices):,} -> {len(new_verts):,} vertices "
+        f"({len(vertices) - len(new_verts):,} removed), "
+        f"center=[{center[0]:.0f},{center[1]:.0f},{center[2]:.0f}], "
+        f"radius={radius:.0f}mm"
+    )
+
+    return head_surface, head_mask
+
+
+@cdc.validate_schemas
+def align_axes_from_landmarks(
+    surface: cdc.TrimeshSurface,
+    landmarks: cdt.LabeledPoints,
+) -> tuple[cdc.TrimeshSurface, cdt.LabeledPoints, np.ndarray]:
+    """Map mesh + landmarks into the CTF anatomical frame.
+
+    CTF convention:
+
+        +X = anterior (toward Nz)
+        +Y = left (toward Lpa)
+        +Z = up (toward Cz)
+        origin = midpoint of Lpa and Rpa (interaural midpoint)
+
+    The returned surface and landmarks carry ``crs="ctf"``.
+
+    Args:
+        surface: Axis-normalized TrimeshSurface (post ``normalize_axes`` and
+            ``isolate_head``).
+        landmarks: LabeledPoints with labels Nz, Iz, Cz, and LPA/RPA (or
+            aliases like Lpa/Rpa, "left ear"/"right ear"; normalized
+            internally via ``normalize_landmarks_labels``)
+            (matching the surface frame).
+
+    Returns:
+        Tuple of (aligned_surface, aligned_landmarks, transform).
+        ``transform`` is a 4x4 homogeneous affine that maps input-frame
+        points into CTF; apply as ``hom = transform @ [x, y, z, 1]``.
+
+    Raises:
+        ValueError: If any of the required landmarks (Nz, Iz, Cz, LPA, RPA)
+            are missing after label normalization.
+    """
+    landmarks = normalize_landmarks_labels(landmarks)
+    lm = landmarks.pint.dequantify().values
+    labels = list(landmarks["label"].values)
+    idx = {lbl: i for i, lbl in enumerate(labels)}
+
+    required = {"Nz", "Iz", "Cz", "LPA", "RPA"}
+    missing = required - set(labels)
+    if missing:
+        raise ValueError(f"Missing landmarks for alignment: {missing}")
+
+    Nz = lm[idx["Nz"]]
+    Cz = lm[idx["Cz"]]
+    Lpa = lm[idx["LPA"]]
+    Rpa = lm[idx["RPA"]]
+    origin = _ear_midpoint(Lpa, Rpa)
+
+    y_ax = Lpa - Rpa
+    y_ax = y_ax / np.linalg.norm(y_ax)
+
+    nz_dir = Nz - origin
+    nz_dir = nz_dir - np.dot(nz_dir, y_ax) * y_ax
+    x_ax = nz_dir / np.linalg.norm(nz_dir)
+
+    z_ax = np.cross(x_ax, y_ax)
+
+    if np.dot(Cz - origin, z_ax) < 0:
+        z_ax = -z_ax
+    y_ax = np.cross(z_ax, x_ax)
+
+    R = np.vstack([x_ax, y_ax, z_ax])
+    M = np.eye(4)
+    M[:3, :3] = R
+    M[:3, 3] = -R @ origin
+
+    aligned_verts = _apply_affine(np.asarray(surface.mesh.vertices), M)
+    new_mesh = _rebuild_mesh(
+        surface.mesh,
+        vertices=aligned_verts,
+        faces=surface.mesh.faces,
+    )
+    aligned_surface = cdc.TrimeshSurface(
+        new_mesh, crs="ctf", units=surface.units,
+    )
+
+    aligned_landmarks = _transform_labeled_points(
+        landmarks, lambda p: _apply_affine(p, M), "ctf"
+    )
+
+    return aligned_surface, aligned_landmarks, M
+
+
+@cdc.validate_schemas
+def revert_to_einstar_frame(
+    surface: cdc.TrimeshSurface,
+    landmarks: cdt.LabeledPoints,
+    R_normalize: np.ndarray,
+    M_align: np.ndarray,
+) -> tuple[cdc.TrimeshSurface, cdt.LabeledPoints]:
+    """Map an aligned surface and landmarks back into the raw Einstar frame.
+
+    Inverse of ``normalize_axes`` composed with ``align_axes_from_landmarks``,
+    so the returned mesh and landmarks carry ``crs="digitized"`` and match
+    the original ``read_einstar_obj`` output.
+
+    Note that ``isolate_head`` is not invertible: the returned mesh is still
+    head-only even though its coordinates are in the digitized frame.
+
+    Args:
+        surface: TrimeshSurface in the CTF frame (post
+            ``align_axes_from_landmarks``, optionally after masking).
+        landmarks: LabeledPoints in the CTF frame.
+        R_normalize: 3x3 rotation returned by ``normalize_axes``.
+        M_align: 4x4 affine returned by ``align_axes_from_landmarks``.
+
+    Returns:
+        Tuple of (surface_digitized, landmarks_digitized). Both carry
+        ``crs="digitized"`` and the mesh preserves UVs / vertex colors.
+    """
+    M_inv = np.linalg.inv(M_align)
+    R_inv4 = np.eye(4)
+    R_inv4[:3, :3] = R_normalize.T
+    M_total = R_inv4 @ M_inv
+
+    raw_verts = _apply_affine(np.asarray(surface.mesh.vertices), M_total)
+    new_mesh = _rebuild_mesh(
+        surface.mesh,
+        vertices=raw_verts,
+        faces=surface.mesh.faces,
+    )
+    raw_surface = cdc.TrimeshSurface(
+        new_mesh, crs="digitized", units=surface.units,
+    )
+
+    raw_landmarks = _transform_labeled_points(
+        landmarks, lambda p: _apply_affine(p, M_total), "digitized"
+    )
+
+    return raw_surface, raw_landmarks
diff --git a/tests/test_anonymization.py b/tests/test_anonymization.py
new file mode 100644
index 00000000..7fbf18e6
--- /dev/null
+++ b/tests/test_anonymization.py
@@ -0,0 +1,347 @@
+"""Unit tests for the face anonymization pipeline."""
+
+import os
+
+import numpy as np
+import pytest
+import trimesh
+from numpy.testing import assert_allclose
+
+import cedalion
+import cedalion.dataclasses as cdc
+
+from cedalion.geometry.photogrammetry.anonymization import (
+    align_axes_from_landmarks,
+    anonymize_scan,
+    delete_masked_vertices,
+    detect_cap_boundary,
+    face_mask_from_landmarks,
+    isolate_head,
+    normalize_axes,
+    revert_to_einstar_frame,
+    save_anonymized_scan,
+)
+
+
+@pytest.fixture
+def simple_sphere_surface():
+    """Unit sphere as a minimal TrimeshSurface for geometry-only tests."""
+    sphere = trimesh.creation.icosphere(subdivisions=3, radius=100)
+    return cdc.TrimeshSurface(
+        mesh=sphere, crs="scanner", units=cedalion.units.millimeter
+    )
+
+
+@pytest.fixture
+def head_like_surface():
+    """Elongated sphere; X=up, Y=anterior, Z=left (axis-normalized frame)."""
+    sphere = trimesh.creation.icosphere(subdivisions=4, radius=100)
+    vertices = sphere.vertices.copy()
+    vertices[:, 0] *= 1.2
+    vertices[:, 1] *= 0.9
+    head_mesh = trimesh.Trimesh(vertices=vertices, faces=sphere.faces)
+    return cdc.TrimeshSurface(
+        mesh=head_mesh, crs="scanner", units=cedalion.units.millimeter
+    )
+
+
+@pytest.fixture
+def axis_normalized_landmarks():
+    """5 anatomical landmarks in the post-normalize_axes frame (X=up, Y=ant, Z=left)."""
+    coords = np.array([
+        [0, 100, 0],     # Nz (anterior)
+        [0, -100, 0],    # Iz (posterior)
+        [100, 0, 0],     # Cz (top, +X)
+        [0, 0, 100],     # LPA (left, +Z)
+        [0, 0, -100],    # RPA (right, -Z)
+    ], dtype=float)
+    return cdc.build_labeled_points(
+        coords,
+        crs="scanner",
+        units="mm",
+        labels=["Nz", "Iz", "Cz", "LPA", "RPA"],
+        types=[cdc.PointType.LANDMARK] * 5,
+    )
+
+
+def test_normalize_rotation_orthogonal(simple_sphere_surface):
+    """Rotation matrix returned by normalize_axes is orthogonal."""
+    nasion = np.array([0, 50, 50])
+    _, _, R = normalize_axes(simple_sphere_surface, nasion)
+    assert_allclose(R @ R.T, np.eye(3), atol=1e-10)
+
+
+def test_normalize_nasion_to_positive_y(simple_sphere_surface):
+    """After normalization the rotated nasion has a positive Y component."""
+    nasion = np.array([0, 50, 50])
+    _, rotated_nasion, _ = normalize_axes(simple_sphere_surface, nasion)
+    assert rotated_nasion[1] > 0
+
+
+def test_normalize_identity_when_aligned(simple_sphere_surface):
+    """No rotation is applied when the nasion already points along +Y."""
+    nasion = np.array([0, 100, 0])
+    _, rotated_nasion, R = normalize_axes(simple_sphere_surface, nasion)
+    assert_allclose(R, np.eye(3), atol=1e-6)
+    assert_allclose(rotated_nasion, nasion, atol=1e-6)
+
+
+def test_isolate_head_only_scan_unchanged(simple_sphere_surface):
+    """A head-only scan (no body) is returned with almost all vertices kept."""
+    nasion = np.array([0, 100, 0])
+    _, mask = isolate_head(simple_sphere_surface, nasion)
+    assert mask.mean() > 0.9
+
+
+def test_isolate_reduces_count_with_body():
+    """Body sphere 400 mm away is removed; only the head sphere survives."""
+    sphere = trimesh.creation.icosphere(subdivisions=3, radius=100)
+    body = trimesh.creation.icosphere(subdivisions=3, radius=80)
+    body.vertices[:, 0] -= 400
+    combined = trimesh.util.concatenate([sphere, body])
+
+    surface = cdc.TrimeshSurface(
+        combined, crs="scanner", units=cedalion.units.millimeter
+    )
+    nasion = np.array([0, 100, 0])
+    head_surface, _ = isolate_head(surface, nasion)
+    assert head_surface.nvertices < surface.nvertices
+
+
+def test_align_origin_at_ear_midpoint(simple_sphere_surface, axis_normalized_landmarks):
+    """CTF origin is placed at the midpoint of LPA and RPA."""
+    _, aligned_lm, _ = align_axes_from_landmarks(
+        simple_sphere_surface, axis_normalized_landmarks
+    )
+    lm = aligned_lm.pint.dequantify()
+    lpa = lm.sel(label="LPA").values
+    rpa = lm.sel(label="RPA").values
+    assert_allclose(0.5 * (lpa + rpa), [0, 0, 0], atol=1e-6)
+
+
+def test_align_axes_orientation(simple_sphere_surface, axis_normalized_landmarks):
+    """Nz points along +X, LPA along +Y, and Cz along +Z in the CTF frame."""
+    _, aligned_lm, _ = align_axes_from_landmarks(
+        simple_sphere_surface, axis_normalized_landmarks
+    )
+    lm = aligned_lm.pint.dequantify()
+    nz = lm.sel(label="Nz").values
+    lpa = lm.sel(label="LPA").values
+    cz = lm.sel(label="Cz").values
+    assert nz[0] > abs(nz[1]) and nz[0] > abs(nz[2])
+    assert lpa[1] > abs(lpa[0]) and lpa[1] > abs(lpa[2])
+    assert cz[2] > abs(cz[0]) and cz[2] > abs(cz[1])
+
+
+def test_align_returns_ctf_crs(simple_sphere_surface, axis_normalized_landmarks):
+    """Aligned surface carries crs='ctf'."""
+    aligned_surface, _, _ = align_axes_from_landmarks(
+        simple_sphere_surface, axis_normalized_landmarks
+    )
+    assert aligned_surface.crs == "ctf"
+
+
+def test_align_raises_on_missing_landmark(
+    simple_sphere_surface, axis_normalized_landmarks
+):
+    """ValueError is raised when a required landmark is absent."""
+    partial = axis_normalized_landmarks.isel(label=slice(0, 4))
+    with pytest.raises(ValueError, match="Missing landmarks"):
+        align_axes_from_landmarks(simple_sphere_surface, partial)
+
+
+def test_cap_boundary_in_sane_range():
+    """cap_z stays inside the documented [Nz[2], Nz[2] + ceiling] window."""
+    sphere = trimesh.creation.icosphere(subdivisions=4, radius=100)
+    verts = np.asarray(sphere.vertices)
+    Nz = np.array([100, 0, 0])
+    Cz = np.array([0, 0, 100])
+    Lpa = np.array([0, 1, 0])
+    Rpa = np.array([0, -1, 0])
+    cap_z, *_ = detect_cap_boundary(verts, Nz, Cz, Lpa, Rpa)
+    assert Nz[2] <= cap_z <= Nz[2] + 40.0
+
+
+def test_cap_boundary_failsafe_clamps_implausible_cap():
+    """Failsafe clamps cap_z to Nz[2] + eyebrow_offset_mm when detection runs hot.
+
+    Builds a midline X-max profile that rises steeply early then creeps gently
+    up to Cz; the gradient walk-back catches the gentle phase at high z, so
+    without the failsafe cap_z would land above the ceiling.
+    """
+    z = np.linspace(0.5, 99.5, 200)
+    x = 50.0 + 50.0 * np.tanh(z / 8.0) + 0.05 * z
+    y = np.zeros_like(z)
+    anterior = np.column_stack([x, y, z])
+    posterior = np.column_stack([-x, y, z])
+    verts = np.vstack([anterior, posterior])
+    Nz = np.array([100.0, 0.0, 0.0])
+    Cz = np.array([0.0, 0.0, 100.0])
+    Lpa = np.array([0.0, 1.0, 0.0])
+    Rpa = np.array([0.0, -1.0, 0.0])
+    cap_z, *_ = detect_cap_boundary(verts, Nz, Cz, Lpa, Rpa)
+    assert cap_z == pytest.approx(Nz[2] + 10.0)
+
+
+def test_face_mask_region_semantics():
+    """Anterior-below-cap vertex is masked; posterior and above-cap are not."""
+    Nz = np.array([100, 0, 10])
+    Lpa = np.array([0, 100, 0])
+    Rpa = np.array([0, -100, 0])
+    verts = np.array([
+        [50, 0, 0],     # anterior, below cap
+        [-50, 0, 0],    # posterior
+        [50, 0, 50],    # anterior, above cap
+    ])
+    mask, _ = face_mask_from_landmarks(verts, Nz, Lpa, Rpa, cap_z=20.0)
+    assert mask.tolist() == [True, False, False]
+
+
+def test_face_mask_ear_sphere():
+    """Vertex near LPA is masked even when posterior to the ear coronal plane."""
+    Nz = np.array([100, 0, 0])
+    Lpa = np.array([0, 100, 0])
+    Rpa = np.array([0, -100, 0])
+    # 10mm posterior to ear midpoint, but 20mm from LPA (within ear_delete_radius=40)
+    verts = np.array([[-10, 80, 0]])
+    mask, _ = face_mask_from_landmarks(
+        verts, Nz, Lpa, Rpa, cap_z=10.0, ear_delete_radius=40.0
+    )
+    assert mask[0]
+
+
+def test_delete_no_op_on_false_mask(simple_sphere_surface):
+    """All-False mask leaves the vertex count unchanged."""
+    mask = np.zeros(simple_sphere_surface.nvertices, dtype=bool)
+    result = delete_masked_vertices(simple_sphere_surface, mask)
+    assert result.nvertices == simple_sphere_surface.nvertices
+
+
+def test_delete_partial_mask_reduces_count(simple_sphere_surface):
+    """Masking one hemisphere reduces the vertex count."""
+    mask = np.asarray(simple_sphere_surface.mesh.vertices)[:, 0] > 0
+    result = delete_masked_vertices(simple_sphere_surface, mask)
+    assert result.nvertices < simple_sphere_surface.nvertices
+
+
+def test_delete_preserves_crs_and_units(simple_sphere_surface):
+    """CRS and units are propagated unchanged after vertex deletion."""
+    mask = np.zeros(simple_sphere_surface.nvertices, dtype=bool)
+    mask[0] = True
+    result = delete_masked_vertices(simple_sphere_surface, mask)
+    assert result.crs == simple_sphere_surface.crs
+    assert result.units == simple_sphere_surface.units
+
+
+def test_revert_round_trip_with_align(simple_sphere_surface, axis_normalized_landmarks):
+    """align then revert recovers the original vertex positions."""
+    aligned_surface, aligned_lm, M = align_axes_from_landmarks(
+        simple_sphere_surface, axis_normalized_landmarks
+    )
+    reverted_surface, _ = revert_to_einstar_frame(
+        aligned_surface, aligned_lm, R_normalize=np.eye(3), M_align=M
+    )
+    assert_allclose(
+        np.asarray(reverted_surface.mesh.vertices),
+        np.asarray(simple_sphere_surface.mesh.vertices),
+        atol=1e-6,
+    )
+
+
+def test_revert_returns_digitized_crs(
+    simple_sphere_surface, axis_normalized_landmarks
+):
+    """Reverted surface carries crs='digitized'."""
+    aligned_surface, aligned_lm, M = align_axes_from_landmarks(
+        simple_sphere_surface, axis_normalized_landmarks
+    )
+    reverted_surface, _ = revert_to_einstar_frame(
+        aligned_surface, aligned_lm, np.eye(3), M
+    )
+    assert reverted_surface.crs == "digitized"
+
+
+def test_save_raises_on_bad_extension(simple_sphere_surface, tmp_path):
+    """ValueError is raised when out_path does not end in .obj."""
+    out = str(tmp_path / "out.txt")
+    with pytest.raises(ValueError, match=".obj"):
+        save_anonymized_scan(simple_sphere_surface, out)
+
+
+def test_save_geometry_only(simple_sphere_surface, tmp_path):
+    """strip_texture=True writes an OBJ file without MTL or JPG."""
+    out = str(tmp_path / "anon.obj")
+    written = save_anonymized_scan(simple_sphere_surface, out, strip_texture=True)
+    assert os.path.exists(out)
+    assert any(p.endswith(".obj") for p in written)
+
+
+def test_anonymize_scan_reduces_vertices(head_like_surface, axis_normalized_landmarks):
+    """Anonymized surface has fewer vertices than the input (face region deleted)."""
+    surface_anon, _ = anonymize_scan(head_like_surface, axis_normalized_landmarks)
+    assert surface_anon.nvertices < head_like_surface.nvertices
+
+
+def test_anonymize_scan_returns_digitized_frame(
+    head_like_surface, axis_normalized_landmarks
+):
+    """Default return frame is 'digitized' for both surface and landmarks."""
+    surface_anon, landmarks_anon = anonymize_scan(
+        head_like_surface, axis_normalized_landmarks
+    )
+    assert surface_anon.crs == "digitized"
+    assert "digitized" in landmarks_anon.dims
+
+
+def test_anonymize_scan_return_frame_ctf(head_like_surface, axis_normalized_landmarks):
+    """return_frame='ctf' keeps the surface in the CTF coordinate frame."""
+    surface_anon, landmarks_anon = anonymize_scan(
+        head_like_surface, axis_normalized_landmarks, return_frame="ctf"
+    )
+    assert surface_anon.crs == "ctf"
+
+
+def test_anonymize_scan_raises_on_missing_landmark(
+    head_like_surface, axis_normalized_landmarks
+):
+    """ValueError is raised when fewer than 5 required landmarks are provided."""
+    partial = axis_normalized_landmarks.isel(label=slice(0, 3))
+    with pytest.raises(ValueError, match="Missing landmarks"):
+        anonymize_scan(head_like_surface, partial)
+
+
+def test_full_anonymization_pipeline(
+    head_like_surface, axis_normalized_landmarks, tmp_path
+):
+    """End-to-end pipeline: normalize, isolate, align, mask, revert, save."""
+    nasion = axis_normalized_landmarks.pint.dequantify().sel(label="Nz").values
+    surface_n, nasion_n, R = normalize_axes(head_like_surface, nasion)
+    lm_arr = axis_normalized_landmarks.pint.dequantify().values
+    landmarks_n = (
+        axis_normalized_landmarks.pint.dequantify()
+        .copy(data=lm_arr @ R.T)
+        .pint.quantify()
+    )
+    surface_n, _ = isolate_head(surface_n, nasion_n)
+
+    surface_h, landmarks_n, M_ctf = align_axes_from_landmarks(surface_n, landmarks_n)
+    verts = np.asarray(surface_h.mesh.vertices)
+    lm_n = landmarks_n.pint.dequantify()
+    Nz = lm_n.sel(label="Nz").values
+    Cz = lm_n.sel(label="Cz").values
+    Lpa = lm_n.sel(label="LPA").values
+    Rpa = lm_n.sel(label="RPA").values
+
+    cap_z, *_ = detect_cap_boundary(verts, Nz, Cz, Lpa, Rpa)
+    mask, _ = face_mask_from_landmarks(verts, Nz, Lpa, Rpa, cap_z=cap_z)
+    surface_anon = delete_masked_vertices(surface_h, mask)
+    surface_anon_dig, _ = revert_to_einstar_frame(
+        surface_anon, landmarks_n, R, M_ctf
+    )
+
+    out = str(tmp_path / "anon.obj")
+    save_anonymized_scan(surface_anon_dig, out, strip_texture=True)
+
+    assert surface_anon.nvertices < surface_h.nvertices
+    assert surface_anon_dig.crs == "digitized"
+    assert os.path.exists(out)