Electrode microstructure analysis

The executed script for detector preparation, joint segmentation and physical measurements. The reported vectors use its stacked_three_channel output.

1#!/usr/bin/env python3
2"""Separate three-detector segmentation experiment; original exports are read-only."""
3from __future__ import annotations
4 
5import argparse
6import csv
7import hashlib
8import json
9from pathlib import Path
10import platform
11import shutil
12import subprocess
13import sys
14 
15import matplotlib
16matplotlib.use("Agg")
17import matplotlib.pyplot as plt
18from matplotlib.colors import ListedColormap
19from matplotlib.patches import Patch
20import numpy as np
21from PIL import Image
22import sklearn
23import scipy
24import skimage
25import PIL
26from sklearn.cluster import KMeans
27from skimage.filters import gaussian, threshold_multiotsu
28from skimage.measure import label, regionprops
29from threadpoolctl import threadpool_limits
30 
31PORTAL_ROOT = Path(__file__).resolve().parents[1]
32ROOT = PORTAL_ROOT.parent
33sys.path.insert(0, str(PORTAL_ROOT / "scripts"))
34from extract_physical_features import FEATURE_COLUMNS, load_imagerep
35from compare_physical_batches import FEATURES
36 
37PIXEL_SIZE = .025
38VARIANTS = ["original_bse_otsu", "bse_kmeans_control", "stacked_three_channel"]
39PHASES = ["Pore-labelled", "Graphite-labelled", "Silicon-labelled"]
40COLOURS = ["#202020", "#b7bbc0", "#bd8b43"]
41 
42 
43def digest(path):
44    return hashlib.sha256(Path(path).read_bytes()).hexdigest()
45 
46 
47def inventory(data_dir, batch_names=None):
48    rows = []
49    for batch in batch_names or ["Batch_1", "Batch_2", "Batch_3", "test_1"]:
50        files = sorted((data_dir / batch).glob("img_*_BSE.tif"))
51        if not files:
52            raise ValueError(f"No BSE images in requested directory: {data_dir / batch}")
53        for bse in files:
54            sample = bse.stem.removeprefix("img_").removesuffix("_BSE")
55            inlens = bse.with_name(f"img_{sample}_Inlens.tif")
56            thirds = [bse.with_name(f"img_{sample}_{name}.tif") for name in ["ETD", "SE"]]
57            thirds = [p for p in thirds if p.exists()]
58            if not inlens.exists() or len(thirds) != 1:
59                raise ValueError(f"Expected BSE/Inlens and exactly one ETD or SE file: {sample}")
60            rows.append({"sample_id": sample, "batch": batch, "paths": [bse, inlens, thirds[0]],
61                         "third_channel": thirds[0].stem.rsplit("_", 1)[1]})
62    if not rows:
63        raise ValueError("No samples found")
64    return rows
65 
66 
67def load_channels(paths):
68    images, shapes = [], []
69    for path in paths:
70        with Image.open(path) as im:
71            raw = np.array(im)
72        if raw.ndim == 3:
73            raw = raw[:, :, 0]
74        if raw.ndim != 2:
75            raise ValueError(f"Expected a 2D intensity image: {path}")
76        shapes.append(raw.shape)
77        h = raw.shape[0]
78        images.append(gaussian(raw[int(.1*h):int(.9*h), :], sigma=1.0, preserve_range=True))
79    if len(set(shapes)) != 1:
80        raise ValueError("Channels have different pixel grids; no automatic resize or registration is performed")
81    return images, shapes[0]
82 
83 
84def otsu_mask(bse):
85    low, high = threshold_multiotsu(bse, classes=3)
86    labels = np.zeros(bse.shape, dtype=np.uint8)
87    labels[bse >= low] = 1
88    labels[bse >= high] = 2
89    return labels, [float(low), float(high)]
90 
91 
92def clustered_mask(channels, original, seed=42, max_pixels=100_000):
93    means = np.array([im.mean() for im in channels])
94    scales = np.array([im.std(ddof=0) for im in channels])
95    if not np.isfinite(scales).all() or np.any(scales <= 0):
96        raise ValueError("A detector image is constant or nonfinite")
97    flat = [im.ravel() for im in channels]
98    sampled_indices = np.random.default_rng(seed).choice(original.size, min(max_pixels, original.size), replace=False)
99    sampled = (np.column_stack([im[sampled_indices] for im in flat]) - means) / scales
100    sampled_labels = original.ravel()[sampled_indices]
101    if set(sampled_labels) != {0, 1, 2}:
102        raise ValueError("The sampled pixels do not contain all three initial BSE classes")
103    initial = np.array([sampled[sampled_labels == phase].mean(axis=0) for phase in range(3)])
104    model = KMeans(n_clusters=3, init=initial, n_init=1, random_state=seed,
105                   max_iter=300, tol=1e-4, algorithm="lloyd").fit(sampled)
106    if model.n_iter_ >= 300:
107        raise ValueError("Clustering reached the iteration limit")
108    if len(set(model.labels_)) != 3:
109        raise ValueError("Clustering did not retain three distinct classes")
110    # Map joint clusters to the same intensity-based interpretation as the original.
111    order = np.argsort(model.cluster_centers_[:, 0], kind="stable")
112    mapping = np.empty(3, dtype=np.uint8)
113    mapping[order] = np.arange(3)
114    labels = np.empty(original.size, dtype=np.uint8)
115    for start in range(0, labels.size, 250_000):
116        end = min(start + 250_000, labels.size)
117        pixels = (np.column_stack([im[start:end] for im in flat]) - means) / scales
118        labels[start:end] = mapping[model.predict(pixels)]
119    return labels.reshape(original.shape), {
120        "channel_means": means.tolist(), "channel_population_sd": scales.tolist(),
121        "initial_normalised_centres_in_original_phase_order": initial.tolist(),
122        "normalised_cluster_centres_in_phase_order": model.cluster_centers_[order].tolist(),
123        "raw_cluster_centres_in_phase_order": (model.cluster_centers_[order] * scales + means).tolist(),
124        "sampled_pixels": len(sampled_indices), "seed": seed, "iterations": int(model.n_iter_),
125        "normalisation": "Per-image mean and population SD; each included channel has equal weight",
126        "initialisation": "Means of sampled pixels in the original three BSE Multi-Otsu masks",
127        "phase_mapping": "Increasing BSE cluster-centre intensity: pore, graphite, silicon",
128    }
129 
130 
131def measure_masks(labels, inlens, core):
132    """Original downstream calculations, with the segmentation masks as inputs."""
133    pore, matrix, silicon = labels == 0, labels == 1, labels == 2
134    if not 0 < pore.mean() < 1:
135        raise ValueError("Pore mask cannot be empty or fill the image")
136    median = np.median(inlens[pore])
137    cbd = pore & (inlens > median)
138    uncertainty = core.make_error_prediction(pore.astype(int), confidence=.95, target_error=.05)
139    chords = []
140    for axis in ["y", "x"]:
141        lines = (pore[:, col] for col in range(0, pore.shape[1], 40)) if axis == "y" else (pore[row, :] for row in range(0, pore.shape[0], 40))
142        lengths = []
143        for line in lines:
144            diffs = np.diff(np.concatenate(([0], line.astype(int), [0])))
145            starts, ends = np.where(diffs == 1)[0], np.where(diffs == -1)[0]
146            lengths.extend((ends - starts) * PIXEL_SIZE)
147        chords.append(float(np.mean(lengths)) if lengths else 0.0)
148    ly, lx = chords
149    props = [p for p in regionprops(label(silicon)) if p.area >= 10]
150    diameters = [p.equivalent_diameter_area * PIXEL_SIZE for p in props]
151    ratios = [p.axis_major_length / max(p.axis_minor_length, 1) for p in props if p.axis_minor_length > 0]
152    d10, d50, d90 = np.percentile(diameters, [10, 50, 90]) if diameters else [0., 0., 0.]
153    tiles = [100*np.mean(tile) for band in np.array_split(silicon, 4, axis=0) for tile in np.array_split(band, 4, axis=1)]
154    values = [100*float(pore.mean()), 100*float(uncertainty["abs_err"]),
155              100*float(matrix.mean()), 100*float(cbd.mean()), 100*float(silicon.mean()),
156              PIXEL_SIZE*float(uncertainty["integral_range"]), ly, lx, lx/ly if ly > 0 else 1.,
157              float(np.mean(ratios)) if ratios else 1., float(d10), float(d50), float(d90), float(np.std(tiles, ddof=0))]
158    if not np.isfinite(values).all():
159        raise ValueError("A physical calculation returned a nonfinite value")
160    return dict(zip(FEATURE_COLUMNS, values))
161 
162 
163def saved_features():
164    result = {}
165    aliases = {f["key"]: f["raw_key"] for f in FEATURES}
166    aliases["porosity_ci95_pct"] = "ci95_abs_pct"
167    for filename, canonical in [("qc_dataset_features.csv", False), ("test_1_physical_features.csv", True)]:
168        with (PORTAL_ROOT / "public" / filename).open(newline="") as handle:
169            for row in csv.DictReader(handle):
170                result[row["sample_id"]] = {key: float(row[key if canonical else aliases[key]]) for key in FEATURE_COLUMNS}
171    return result
172 
173 
174def mask_comparison(original, candidate):
175    output = {"changed_pixel_fraction": float(np.mean(original != candidate))}
176    for k, name in enumerate(["pore", "graphite", "silicon"]):
177        a, b = original == k, candidate == k
178        output[name + "_iou"] = float(np.sum(a & b) / np.sum(a | b))
179        output[name + "_dice"] = float(2 * np.sum(a & b) / (np.sum(a) + np.sum(b)))
180        output[name + "_changed_pixel_fraction"] = float(np.mean(a != b))
181    return output
182 
183 
184def draw_masks(row, channels, masks, measurements, output):
185    # Display a fixed centre window; segmentation and measurements use the full crop.
186    h, w = channels[0].shape
187    size = min(900, h, w)
188    y, x = (h-size)//2, (w-size)//2
189    view = (slice(y, y+size), slice(x, x+size))
190    fig, axes = plt.subplots(2, 3, figsize=(12, 9.8))
191    for ax, im, name in zip(axes[0], channels, ["BSE", "Inlens", row["third_channel"]]):
192        ax.imshow(im[view], cmap="gray", vmin=float(im.min()), vmax=float(im.max()))
193        ax.set_title(name, fontweight="bold")
194    for ax, mask, variant, title in zip(axes[1], masks, VARIANTS,
195                                      ["Original BSE Multi-Otsu", "BSE-only clustering control", "Three-channel clustering"]):
196        value = measurements[variant]
197        ax.imshow(mask[view], cmap=ListedColormap(COLOURS), vmin=0, vmax=2, interpolation="nearest")
198        ax.set_title(f"{title}\nPore {value['porosity_pct']:.2f}%; silicon {value['inclusion_pct']:.2f}%", fontsize=11)
199    for ax in axes.flat:
200        ax.set_axis_off()
201    fig.suptitle(f"{row['batch']} / {row['sample_id']}: matched detector views and masks", fontsize=14, fontweight="bold")
202    fig.legend(handles=[Patch(facecolor=c, label=p) for c, p in zip(COLOURS, PHASES)],
203               loc="lower center", bbox_to_anchor=(.5, .055), ncol=3, frameon=False)
204    fig.text(.5, .022, f"Same central {size} × {size} pixel display window. Percentages use the entire {h} × {w} pixel crop.\nIntensity clusters are candidate phase masks; CBD median allocation is not shown.", ha="center", fontsize=9)
205    fig.subplots_adjust(left=.025, right=.975, top=.94, bottom=.13, wspace=.04, hspace=.12)
206    fig.savefig(output, dpi=150, facecolor="white")
207    plt.close(fig)
208 
209 
210def main():
211    parser = argparse.ArgumentParser(description=__doc__)
212    parser.add_argument("--data-dir", type=Path, default=ROOT / "data")
213    parser.add_argument("--run-dir", type=Path, required=True)
214    parser.add_argument("--imagerep-path", type=Path, default=ROOT / "ImageRep")
215    parser.add_argument("--batch-name", action="append", help="Folder under --data-dir; repeat to select several folders")
216    parser.add_argument("--sample-id", action="append", default=[])
217    args = parser.parse_args()
218    if args.run_dir.exists():
219        parser.error("Run directory already exists; choose a new directory")
220    rows = inventory(args.data_dir, args.batch_name)
221    if args.sample_id:
222        rows = [r for r in rows if r["sample_id"] in args.sample_id]
223        if set(args.sample_id) != {r["sample_id"] for r in rows}:
224            parser.error("Requested sample not found")
225    args.run_dir.mkdir(parents=True)
226    (args.run_dir / "figures").mkdir()
227    (args.run_dir / "masks").mkdir()
228    original_values = saved_features()
229    core = load_imagerep(args.imagerep_path)
230    source_files = [Path(__file__), PORTAL_ROOT / "scripts/extract_physical_features.py",
231                    PORTAL_ROOT / "scripts/compare_physical_batches.py", Path(core.__file__)]
232    source_dir = args.run_dir / "source_snapshot"
233    source_dir.mkdir()
234    source_hashes = {str(p): digest(p) for p in source_files}
235    for path in source_files:
236        shutil.copy2(path, source_dir / path.name)
237    imagerep_revision = subprocess.check_output(["git", "-C", str(args.imagerep_path), "rev-parse", "HEAD"], text=True).strip()
238    records, features = [], []
239    for index, row in enumerate(rows):
240        print(f"[{index+1}/{len(rows)}] {row['batch']}/{row['sample_id']} ({row['third_channel']})", flush=True)
241        channels, native_shape = load_channels(row["paths"])
242        original, thresholds = otsu_mask(channels[0])
243        # Validate the unchanged downstream calculation against saved full precision.
244        baseline = measure_masks(original, channels[1], core)
245        saved = original_values.get(row["sample_id"])
246        errors = {key: abs(baseline[key] - saved[key]) for key in FEATURE_COLUMNS} if saved else None
247        if errors is not None and max(errors.values()) > 1e-10:
248            raise ValueError(f"Original-method reproduction failed for {row['sample_id']}: {errors}")
249        with threadpool_limits(limits=1):
250            bse, bse_fit = clustered_mask(channels[:1], original)
251            stacked, stack_fit = clustered_mask(channels, original)
252        masks = [original, bse, stacked]
253        measured = {VARIANTS[0]: baseline, VARIANTS[1]: measure_masks(bse, channels[1], core),
254                    VARIANTS[2]: measure_masks(stacked, channels[1], core)}
255        for variant, mask in zip(VARIANTS, masks):
256            features.append({"batch": row["batch"], "sample_id": row["sample_id"],
257                             "third_channel": row["third_channel"], "variant": variant, **measured[variant]})
258            Image.fromarray(mask).save(args.run_dir / "masks" / f"{row['sample_id']}_{variant}.png")
259        draw_masks(row, channels, masks, measured, args.run_dir / "figures" / f"{row['sample_id']}.png")
260        records.append({"batch": row["batch"], "sample_id": row["sample_id"], "third_channel": row["third_channel"],
261                        "input_files": [{"path": str(p), "sha256": digest(p)} for p in row["paths"]],
262                        "native_shape": list(native_shape), "crop_shape": list(original.shape),
263                        "original_thresholds": thresholds, "original_max_absolute_feature_error": max(errors.values()) if errors is not None else None,
264                        "original_verification": "Compared with saved values" if saved else "New sample; no previous export to compare",
265                        "fits": {VARIANTS[1]: bse_fit, VARIANTS[2]: stack_fit}, "measurements": measured,
266                        "agreement_with_original": {VARIANTS[1]: mask_comparison(original, bse), VARIANTS[2]: mask_comparison(original, stacked)},
267                        "added_channel_effect": mask_comparison(bse, stacked)})
268        (args.run_dir / "progress.json").write_text(json.dumps(records, allow_nan=False) + "\n")
269        print("  pore fractions: " + ", ".join(f"{v}={measured[v]['porosity_pct']:.2f}%" for v in VARIANTS), flush=True)
270    with (args.run_dir / "features.csv").open("w", newline="") as handle:
271        writer = csv.DictWriter(handle, fieldnames=["batch", "sample_id", "third_channel", "variant", *FEATURE_COLUMNS])
272        writer.writeheader()
273        writer.writerows(features)
274    report = {"experiment": "Three-channel intensity segmentation v1", "status": "Exploratory side experiment; original results unchanged",
275              "method": {"pixel_size_um": PIXEL_SIZE, "crop": "central 80% rows, all columns", "gaussian_sigma_px": 1,
276                         "third_channel": "ETD or SE, treated as the third detector at the user's direction",
277                         "registration": "No registration or resampling; matched pixel grids assumed, separate alignment audit",
278                         "clustering": "Three KMeans intensity classes; per-image channel z-scores; 100000 sampled pixels; BSE Multi-Otsu initial centres; seed42; n_init1",
279                         "downstream": "Same 14 saved physical fields, including the separate heuristic Inlens median CBD allocation",
280                         "silicon": "Bright phase identified as silicon by the user; intensity-derived boundaries still require validation",
281                         "agreement": "Mask IoU, Dice and changed-pixel fractions measure agreement, not accuracy against ground truth",
282                         "test_labels": "Not used in segmentation, measurement or settings"},
283              "script_sha256": digest(__file__), "original_extractor_sha256": digest(PORTAL_ROOT / "scripts/extract_physical_features.py"),
284              "source_hashes": source_hashes, "imagerep_revision": imagerep_revision,
285              "versions": {"python": platform.python_version(), "numpy": np.__version__, "sklearn": sklearn.__version__,
286                           "scipy": scipy.__version__, "skimage": skimage.__version__, "pillow": PIL.__version__,
287                           "matplotlib": matplotlib.__version__},
288              "samples": records}
289    (args.run_dir / "report.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
290    print(f"Wrote {len(rows)} × 3 feature records to {args.run_dir}")
291 
292 
293if __name__ == "__main__":
294    main()

The executed assignment script. The reported confidence bands use the empirical_precision_v2 option specified in the protocol.

1#!/usr/bin/env python3
2"""Assign multichannel physical vectors with a fixed, transparent centroid rule.
3 
4First run without --test-csv to freeze the method and evaluate known crops.
5Then supply --test-csv with the same output directory. No test labels are used.
6Confidence bands describe heuristic assignment support, never probabilities or
7manufacturing acceptance. All evaluation remains provisional at crop level.
8"""
9from __future__ import annotations
10 
11import argparse
12import csv
13from datetime import datetime, timezone
14import hashlib
15import json
16from pathlib import Path
17import platform
18 
19import numpy as np
20 
21 
22CLASSES = ("Batch_1", "Batch_2", "Batch_3")
23VARIANT = "stacked_three_channel"
24FEATURES = (
25    ("porosity_pct", "Pore-labelled area fraction", "%"),
26    ("porosity_ci95_pct", "Porosity uncertainty half-width", "percentage points"),
27    ("active_material_pct", "Graphite-labelled area fraction", "%"),
28    ("cbd_pct", "Carbon-binder allocation", "%"),
29    ("inclusion_pct", "Silicon-labelled area fraction", "%"),
30    ("char_length_scale_cls_um", "Pore correlation length", "µm"),
31    ("throat_through_plane_ly_um", "Vertical pore chord", "µm"),
32    ("chord_in_plane_lx_um", "Horizontal pore chord", "µm"),
33    ("pore_anisotropy_ratio", "Horizontal/vertical pore chord ratio", "ratio"),
34    ("particle_aspect_ratio", "Silicon-labelled domain aspect ratio", "ratio"),
35    ("particle_d50_um", "Silicon-labelled domain diameter D50", "µm"),
36    ("particle_d90_um", "Silicon-labelled domain diameter D90", "µm"),
37    ("slurry_dispersion_index", "Silicon-labelled area spatial variation", "percentage points"),
38)
39KEYS = tuple(item[0] for item in FEATURES)
40# Fixed before inspecting the full test vectors. These are design choices,
41# not empirically calibrated probability cutoffs or optimized hyperparameters.
42THRESHOLDS = {
43    "High": {"oof_precision_min": 0.85, "oof_prediction_count_min": 5,
44             "relative_margin_min": 0.25, "leave_one_crop_out_agreement_min": 0.90},
45    "Medium": {"oof_precision_min": 0.60, "oof_prediction_count_min": 5,
46               "relative_margin_min": 0.10, "leave_one_crop_out_agreement_min": 0.75},
47}
48SUPPORT_QUANTILE = 0.95
49POLICIES = ("conservative_v1", "empirical_precision_v2")
50CONFIDENCE_POLICY_V2 = {
51    "requested_by": "User", "threshold_source": "Explicit user instruction; no outcome-based threshold search",
52    "revision_context": "Reporting bands revised after the v1 assignment summary. The fitted classifier, features and assignments are unchanged.",
53    "score": "100 * correctly assigned known crops / all known crops predicted as the assigned class in leave-one-crop-out evaluation",
54    "score_scope": "Class-level observed validation precision; not an individual calibrated probability",
55    "thresholds": {"High": "strictly greater than70%", "Medium": "50% through70%, inclusive", "Low": "below50%"},
56    "other_checks": "Distance margin, omission agreement and empirical support envelope remain diagnostics and do not change this user-defined band",
57}
58LIMITATIONS = [
59    "The ratings are heuristic assignment support, not calibrated probabilities.",
60    "Known crops may share original images; parent-image groups are unavailable. Leave-one-crop-out evaluation can therefore be optimistic.",
61    "Per-class precision estimates are based on few crop predictions and unknown test class proportions.",
62    "Thirteen equally standardized exported fields retain correlated and derived variables, including porosity uncertainty and the carbon-binder allocation. They are not thirteen independent physical mechanisms.",
63    "The support envelope is an empirical distance summary, not a conformal prediction region or a 95% probability statement.",
64    "Leave-one-crop-out agreement measures sensitivity to individual known crops, not correctness or uncertainty in phase boundaries.",
65    "Silicon-labelled mask boundaries remain unvalidated and can include rims or lines; pore chords do not establish three-dimensional transport.",
66    "Three earlier test identities have already been revealed. This is not an untouched external validation exercise; test labels are nevertheless excluded from fitting and confidence calculations.",
67    "Batch identity does not establish manufacturing quality. No accept or reject decision is produced.",
68]
69 
70 
71def sha256(path):
72    return hashlib.sha256(Path(path).read_bytes()).hexdigest()
73 
74 
75def canonical_hash(value):
76    return hashlib.sha256(json.dumps(value, sort_keys=True, separators=(",", ":"),
77                                     allow_nan=False).encode()).hexdigest()
78 
79 
80def write_json(path, value):
81    Path(path).write_text(json.dumps(value, indent=2, allow_nan=False) + "\n")
82 
83 
84def read_vectors(path, *, known):
85    """Read only the requested method and split; never read test-label columns."""
86    records = []
87    with Path(path).open(newline="", encoding="utf-8-sig") as handle:
88        reader = csv.DictReader(handle)
89        header = reader.fieldnames or []
90        required = {"batch", "sample_id", "variant", *KEYS}
91        if len(header) != len(set(header)) or not required.issubset(header):
92            raise ValueError("Unique CSV columns including batch, sample_id, variant and all13 fixed fields are required")
93        seen = set()
94        for row in reader:
95            if None in row or any(value is None for value in row.values()):
96                raise ValueError("Malformed CSV row")
97            if row["variant"] != VARIANT:
98                continue
99            batch = row["batch"].strip()
100            if known and batch not in CLASSES:
101                continue
102            if not known and batch in CLASSES:
103                raise ValueError("Test input includes known-batch rows; provide the separate test extraction")
104            sample = row["sample_id"].strip()
105            if not sample or sample in seen:
106                raise ValueError("Sample IDs must be nonempty and unique within the selected variant")
107            seen.add(sample)
108            values = [float(row[key]) for key in KEYS]
109            if not np.isfinite(values).all():
110                raise ValueError(f"Nonfinite vector: {sample}")
111            records.append({"sample_id": sample, "batch": batch,
112                            "third_channel": row.get("third_channel", ""), "values": values})
113    if not records:
114        raise ValueError("No vectors for the selected method and split")
115    if known and (len(records) != 31 or {label: sum(r["batch"] == label for r in records)
116                                          for label in CLASSES} != {"Batch_1": 7, "Batch_2": 7, "Batch_3": 17}):
117        raise ValueError("This frozen study requires the original31 known crops (7,7,17)")
118    return records
119 
120 
121def fit_rule(x, y):
122    mean = x.mean(axis=0)
123    scale = x.std(axis=0, ddof=1)
124    if not np.isfinite(scale).all() or np.any(scale <= 0):
125        raise ValueError("All13 frozen fields must have positive training sample SD")
126    centers = np.stack([x[y == label].mean(axis=0) for label in CLASSES])
127    return {"mean": mean, "scale": scale, "centroids": centers}
128 
129 
130def distances(rule, x):
131    x = np.atleast_2d(x)
132    squared = ((x[:, None, :] - rule["centroids"][None, :, :]) / rule["scale"]) ** 2
133    return np.sqrt(squared.mean(axis=2)), squared
134 
135 
136def relative_margin(values):
137    order = np.argsort(values, kind="stable")
138    best, runner = order[:2]
139    # Equal zero distances are an exact tie, not maximal alignment.
140    margin = float((values[runner] - values[best]) / values[runner]) if values[runner] > 0 else 0.0
141    return int(best), int(runner), margin
142 
143 
144def wilson_interval(correct, total):
145    """Descriptive binomial interval only; shared parent crops violate its model."""
146    if total == 0:
147        return [None, None]
148    z = 1.959963984540054
149    p = correct / total
150    denominator = 1 + z * z / total
151    center = (p + z * z / (2 * total)) / denominator
152    half = z * np.sqrt(p * (1 - p) / total + z * z / (4 * total * total)) / denominator
153    return [float(center - half), float(center + half)]
154 
155 
156def evaluate_known(records):
157    x = np.array([r["values"] for r in records])
158    y = np.array([r["batch"] for r in records])
159    rows, omission_rules = [], []
160    matrix = np.zeros((3, 3), dtype=int)
161    for i, record in enumerate(records):
162        keep = np.arange(len(records)) != i
163        rule = fit_rule(x[keep], y[keep])
164        omission_rules.append(rule)
165        ds, _ = distances(rule, x[i])
166        best, runner, margin = relative_margin(ds[0])
167        actual = CLASSES.index(record["batch"])
168        matrix[actual, best] += 1
169        rows.append({"sample_id": record["sample_id"], "actual_batch": record["batch"],
170                     "assigned_batch": CLASSES[best], "runner_up": CLASSES[runner],
171                     "correct": best == actual, "relative_margin": margin,
172                     "true_class_distance": float(ds[0, actual]),
173                     "distances": dict(zip(CLASSES, ds[0].tolist()))})
174    per_class = {}
175    for j, label in enumerate(CLASSES):
176        n_predicted, n_actual, correct = int(matrix[:, j].sum()), int(matrix[j].sum()), int(matrix[j, j])
177        true_distances = [r["true_class_distance"] for r in rows if r["actual_batch"] == label]
178        per_class[label] = {"actual_crop_count": n_actual, "oof_prediction_count": n_predicted,
179                            "oof_correct_predictions": correct,
180                            "oof_precision": correct / n_predicted if n_predicted else None,
181                            "oof_recall": correct / n_actual,
182                            "descriptive_binomial_precision_interval95": wilson_interval(correct, n_predicted),
183                            "true_class_loo_distance_quantile95": float(np.quantile(true_distances, SUPPORT_QUANTILE, method="linear")),
184                            "true_class_loo_distances": true_distances}
185    evaluation = {"class_order": list(CLASSES), "confusion_matrix_actual_rows_predicted_columns": matrix.tolist(),
186                  "correct": int(matrix.trace()), "crop_count": len(records),
187                  "accuracy": float(matrix.trace() / len(records)),
188                  "balanced_accuracy": float(np.mean(matrix.diagonal() / matrix.sum(axis=1))),
189                  "per_class": per_class, "predictions": rows,
190                  "precision_interval_note": "Wilson intervals assume independent Bernoulli outcomes; shown descriptively only because crop dependence is unresolved.",
191                  "evaluation_scope": "Fixed-rule leave-one-known-crop-out; scaling and centroids refitted using the other30 crops in each fold."}
192    return evaluation, fit_rule(x, y), omission_rules
193 
194 
195def confidence_band(precision, count, margin, stability, within_support):
196    checks = {}
197    for name, thresholds in THRESHOLDS.items():
198        checks[name] = {
199            "class_oof_precision": precision is not None and precision >= thresholds["oof_precision_min"],
200            "class_oof_prediction_count": count >= thresholds["oof_prediction_count_min"],
201            "relative_margin": margin >= thresholds["relative_margin_min"],
202            "leave_one_crop_out_agreement": stability >= thresholds["leave_one_crop_out_agreement_min"],
203            "within_assigned_class_support": bool(within_support),
204        }
205    band = next((name for name in THRESHOLDS if all(checks[name].values())), "Low")
206    return band, checks
207 
208 
209def precision_rating(precision):
210    if precision is None:
211        return "Low"
212    return "High" if precision > 0.70 else "Medium" if precision >= 0.50 else "Low"
213 
214 
215def assign_tests(test, known, evaluation, rule, omission_rules, policy="conservative_v1"):
216    x = np.array([r["values"] for r in test])
217    known_x = np.array([r["values"] for r in known])
218    known_y = np.array([r["batch"] for r in known])
219    ds, squared = distances(rule, x)
220    omission_predictions = np.array([np.argmin(distances(r, x)[0], axis=1) for r in omission_rules])
221    assignments = []
222    for i, record in enumerate(test):
223        best, runner, margin = relative_margin(ds[i])
224        label = CLASSES[best]
225        validation = evaluation["per_class"][label]
226        stability = float(np.mean(omission_predictions[:, i] == best))
227        envelope = validation["true_class_loo_distance_quantile95"]
228        inside = bool(ds[i, best] <= envelope)
229        band, checks = confidence_band(validation["oof_precision"], validation["oof_prediction_count"], margin, stability, inside)
230        # The summed contributions equal runner-up squared RMS distance minus
231        # assigned squared RMS distance. Positive values favour the assignment.
232        support = (squared[i, runner] - squared[i, best]) / len(KEYS)
233        features = []
234        for j, (key, name, unit) in enumerate(FEATURES):
235            summaries = {}
236            for batch in CLASSES:
237                values = known_x[known_y == batch, j]
238                summaries[batch] = {"mean": float(values.mean()), "sample_sd": float(values.std(ddof=1)),
239                                    "min": float(values.min()), "max": float(values.max())}
240            features.append({"key": key, "label": name, "unit": unit, "value": record["values"][j],
241                             "known_batches": summaries,
242                             "assigned_minus_runner_squared_distance_support": float(support[j]),
243                             "standardized_difference_from_assigned_mean": float((x[i, j] - rule["centroids"][best, j]) / rule["scale"][j]),
244                             "difference_from_batch3_mean": float(x[i, j] - rule["centroids"][2, j])})
245        ranked = sorted(features, key=lambda f: f["assigned_minus_runner_squared_distance_support"], reverse=True)
246        positive = [f for f in ranked if f["assigned_minus_runner_squared_distance_support"] > 0][:5]
247        counter = [f for f in reversed(ranked) if f["assigned_minus_runner_squared_distance_support"] < 0][:3]
248        record_out = {"sample_id": record["sample_id"], "third_channel": record["third_channel"],
249                      "assigned_batch": label, "runner_up": CLASSES[runner],
250                      "confidence": band, "confidence_meaning": "Heuristic assignment support, not probability or QC release",
251                      "distances": dict(zip(CLASSES, ds[i].tolist())),
252                      "relative_margin": margin,
253                      "squared_distance_margin": float(ds[i, runner] ** 2 - ds[i, best] ** 2),
254                      "leave_one_crop_out_agreement": stability,
255                      "leave_one_crop_out_votes": {label_: int(np.sum(omission_predictions[:, i] == j)) for j, label_ in enumerate(CLASSES)},
256                      "assigned_class_oof_precision": validation["oof_precision"],
257                      "assigned_class_oof_prediction_count": validation["oof_prediction_count"],
258                      "assigned_class_support_distance95": envelope,
259                      "within_assigned_class_support": inside,
260                      "confidence_checks": checks,
261                      "unmet_medium_conditions": [key for key, passed in checks["Medium"].items() if not passed],
262                      "top_supporting_features": positive,
263                      "top_counter_features": counter,
264                      "features": features,
265                      "manufacturing_decision": "Not evaluated"}
266        if policy == "empirical_precision_v2":
267            precision = validation["oof_precision"]
268            correct = validation["oof_correct_predictions"]
269            total = validation["oof_prediction_count"]
270            votes = record_out["leave_one_crop_out_votes"][label]
271            explanation = (f"{correct} of {total} held-out known crops assigned to {label.replace('_', ' ')} were correctly identified "
272                           f"({100 * precision:.1f}%). This is a class-level validation rate, not an individual probability. "
273                           f"{votes} of 31 crop-omission refits retained this assignment.")
274            if margin < 0.10:
275                explanation += f" The nearest alternatives are close (relative distance margin {100 * margin:.1f}%)."
276            if not inside:
277                explanation += " This image is outside the assigned class's empirical distance envelope."
278            record_out.update({"diagnostic_confidence_v1": band, "confidence": precision_rating(precision),
279                               "confidence_meaning": "User-defined band of class-level observed validation precision; not individual probability",
280                               "confidence_pct": 100 * precision,
281                               "validation_rate_pct": 100 * precision,
282                               "validation_n_correct": correct, "validation_n_predictions": total,
283                               "explanation": explanation})
284        if not np.isclose(sum(f["assigned_minus_runner_squared_distance_support"] for f in features),
285                          record_out["squared_distance_margin"], rtol=1e-10, atol=1e-12):
286            raise AssertionError("Feature contributions do not sum to the squared-distance margin")
287        assignments.append(record_out)
288    return assignments
289 
290 
291def export_assignment_csv(path, assignments, policy="conservative_v1"):
292    if policy == "empirical_precision_v2":
293        fields = ["sample_id", "assigned_batch", "confidence", "validation_rate_pct", "explanation"]
294        with Path(path).open("w", newline="") as handle:
295            writer = csv.DictWriter(handle, fieldnames=fields)
296            writer.writeheader()
297            writer.writerows({key: record[key] for key in fields} for record in assignments)
298        return
299    fields = ["sample_id", "third_channel", "assigned_batch", "runner_up", "confidence",
300              "distance_Batch_1", "distance_Batch_2", "distance_Batch_3", "relative_margin",
301              "leave_one_crop_out_agreement", "assigned_class_oof_precision", "assigned_class_oof_prediction_count",
302              "assigned_class_support_distance95", "within_assigned_class_support", "unmet_medium_conditions",
303              "top_supporting_features", *KEYS]
304    with Path(path).open("w", newline="") as handle:
305        writer = csv.DictWriter(handle, fieldnames=fields)
306        writer.writeheader()
307        for assignment in assignments:
308            row = {key: assignment[key] for key in fields if key in assignment}
309            row.update({f"distance_{label}": assignment["distances"][label] for label in CLASSES})
310            row.update({f["key"]: f["value"] for f in assignment["features"]})
311            row["unmet_medium_conditions"] = "; ".join(assignment["unmet_medium_conditions"])
312            row["top_supporting_features"] = "; ".join(f'{f["label"]}: {f["value"]:.6g} {f["unit"]}' for f in assignment["top_supporting_features"])
313            writer.writerow(row)
314 
315 
316def export_feature_csv(path, assignments):
317    fields = ["sample_id", "assigned_batch", "runner_up", "feature", "label", "unit", "value",
318              "assigned_minus_runner_squared_distance_support", "difference_from_batch3_mean"]
319    fields += [f"{batch}_{measure}" for batch in CLASSES for measure in ("mean", "sample_sd", "min", "max")]
320    with Path(path).open("w", newline="") as handle:
321        writer = csv.DictWriter(handle, fieldnames=fields)
322        writer.writeheader()
323        for record in assignments:
324            for feature in record["features"]:
325                row = {key: record[key] for key in ("sample_id", "assigned_batch", "runner_up")}
326                row.update({key: feature[key] for key in ("label", "unit", "value", "assigned_minus_runner_squared_distance_support", "difference_from_batch3_mean")})
327                row["feature"] = feature["key"]
328                for batch, summary in feature["known_batches"].items():
329                    row.update({f"{batch}_{key}": value for key, value in summary.items()})
330                writer.writerow(row)
331 
332 
333def main():
334    parser = argparse.ArgumentParser(description=__doc__)
335    parser.add_argument("--known-csv", type=Path, required=True)
336    parser.add_argument("--test-csv", type=Path)
337    parser.add_argument("--output-dir", type=Path, required=True)
338    parser.add_argument("--confidence-policy", choices=POLICIES, default="conservative_v1")
339    args = parser.parse_args()
340    args.output_dir.mkdir(parents=True, exist_ok=True)
341    specification = {
342        "method": "Fixed nearest batch centroid using equally standardized RMS distance",
343        "segmentation_variant": VARIANT, "class_order": list(CLASSES),
344        "features": [{"key": key, "label": label, "unit": unit} for key, label, unit in FEATURES],
345        "excluded_field": "particle_d10_um, omitted in both methods to retain the previously compared common13 fields",
346        "scaling": "Sample standard deviation (ddof=1) across known training crops only; refitted inside each leave-one-out fold",
347        "tie_rule": "First class in Batch_1, Batch_2, Batch_3 order; exact ties have zero relative margin",
348        "relative_margin_formula": "(runner_up_distance - nearest_distance) / runner_up_distance; zero for an all-zero tie",
349        "confidence_thresholds": THRESHOLDS,
350        "confidence_combination": "Every condition must pass; High checked before Medium; otherwise Low",
351        "support_envelope": "Assigned-class distance <=95th percentile (linear interpolation) of known true-class leave-one-crop-out distances",
352        "stability": "Fraction of31 leave-one-known-crop-out refits retaining the full-fit assignment",
353        "feature_contributions": "(runner-up squared standardized deviation minus assigned squared standardized deviation)/13; sum equals runner-up squared RMS distance minus assigned squared RMS distance",
354        "limitations": LIMITATIONS,
355    }
356    specification["confidence_policy"] = args.confidence_policy
357    if args.confidence_policy == "empirical_precision_v2":
358        specification["confidence_policy_v2"] = CONFIDENCE_POLICY_V2
359        specification["diagnostic_thresholds_v1"] = specification.pop("confidence_thresholds")
360        specification["confidence_thresholds"] = CONFIDENCE_POLICY_V2["thresholds"]
361        specification["confidence_combination"] = "Band the class-level observed validation precision using the user's cutoffs; geometry diagnostics do not alter the band"
362    freeze_path = args.output_dir / "frozen_method.json"
363    signature = {"specification_sha256": canonical_hash(specification), "known_csv_sha256": sha256(args.known_csv),
364                 "script_sha256": sha256(__file__)}
365    if freeze_path.exists():
366        frozen = json.loads(freeze_path.read_text())
367        if frozen["signature"] != signature:
368            raise ValueError("Frozen method/source/known input changed; use a separate versioned output directory, never silently retune")
369    else:
370        if args.test_csv:
371            raise ValueError("Freeze and evaluate using known data first: run once without --test-csv")
372        frozen = {"frozen_utc": datetime.now(timezone.utc).isoformat(), "signature": signature,
373                  "specification": specification, "runtime": {"python": platform.python_version(), "numpy": np.__version__}}
374        write_json(freeze_path, frozen)
375    known = read_vectors(args.known_csv, known=True)
376    evaluation, rule, omission_rules = evaluate_known(known)
377    write_json(args.output_dir / "known_evaluation.json", evaluation)
378    with (args.output_dir / "known_loo.csv").open("w", newline="") as handle:
379        fields = ["sample_id", "actual_batch", "assigned_batch", "runner_up", "correct", "relative_margin", "true_class_distance"] + [f"distance_{label}" for label in CLASSES]
380        writer = csv.DictWriter(handle, fieldnames=fields)
381        writer.writeheader()
382        for record in evaluation["predictions"]:
383            row = {key: record[key] for key in fields if key in record}
384            row.update({f"distance_{label}": value for label, value in record["distances"].items()})
385            writer.writerow(row)
386    fitted = {key: value.tolist() for key, value in rule.items()}
387    fitted.update({"feature_order": list(KEYS), "class_order": list(CLASSES), "known_sample_ids": [r["sample_id"] for r in known]})
388    write_json(args.output_dir / "fitted_centroids.json", fitted)
389    print(json.dumps({"known_correct": evaluation["correct"], "known_total": evaluation["crop_count"],
390                      "accuracy": evaluation["accuracy"], "balanced_accuracy": evaluation["balanced_accuracy"],
391                      "per_class": {key: {field: value for field, value in values.items() if field != "true_class_loo_distances"}
392                                    for key, values in evaluation["per_class"].items()}}, indent=2))
393    if not args.test_csv:
394        return
395    test = read_vectors(args.test_csv, known=False)
396    if set(r["sample_id"] for r in known) & set(r["sample_id"] for r in test):
397        raise ValueError("Known and test sample IDs overlap")
398    assignments = assign_tests(test, known, evaluation, rule, omission_rules, args.confidence_policy)
399    result = {"frozen_method": frozen, "test_csv_sha256": sha256(args.test_csv),
400              "test_sample_count": len(test), "known_evaluation": evaluation,
401              "assignment_counts": {label: sum(r["assigned_batch"] == label for r in assignments) for label in CLASSES},
402              "confidence_counts": {level: sum(r["confidence"] == level for r in assignments) for level in ("High", "Medium", "Low")},
403              "assignments": assignments}
404    if args.confidence_policy == "empirical_precision_v2":
405        result["confidence_policy_v2"] = CONFIDENCE_POLICY_V2
406    write_json(args.output_dir / "assignments.json", result)
407    export_assignment_csv(args.output_dir / "assignments.csv", assignments, args.confidence_policy)
408    export_feature_csv(args.output_dir / "feature_evidence.csv", assignments)
409    print(json.dumps({"test_samples": len(test), "assignment_counts": result["assignment_counts"],
410                      "confidence_counts": result["confidence_counts"]}, indent=2))
411 
412 
413if __name__ == "__main__":
414    main()

Current three-detector analysis: 31 known images, nine test images and 14 saved measurement fields per image.

• Test assignments and confidence (CSV)

  All nine images; percentages are class-level validation rates.

• Test measurements (CSV)

  All nine images, with 13 physical descriptors and porosity uncertainty at full precision.

• Known-batch measurements (CSV)

  31 known images processed by the same three-detector segmentation.

• Measurement summaries (JSON)

  Known-batch means, sample standard deviations and all nine individual test vectors.

• Assignment measurements and validation (JSON)

  Per-image measurements, validation rates and the assignment specification.

• Assignment feature evidence (CSV)

  Every assignment input compared with the known batch means.

• Protocol (Markdown)

  Current measurement and assignment methods, with executable instructions.

• Feature vector (LaTeX)

  The combined-method example for Batch 3, sample 0grcilhi; porosity uncertainty is separate.

Project repository