"""Sweep ChangeStar probability and component-area thresholds without rerunning inference.""" from __future__ import annotations import argparse import csv import json from pathlib import Path import cv2 import numpy as np import rasterio from PIL import Image, ImageDraw, ImageFont import importlib.util def load_demo(): path = Path(__file__).with_name("run_change_detection.py") spec = importlib.util.spec_from_file_location("change_detection_demo", path) if spec is None or spec.loader is None: raise RuntimeError("Unable to load change-detection helpers.") module = importlib.util.module_from_spec(spec) spec.loader.exec_module(module) return module def clean_mask(raw: np.ndarray, minimum_area: int) -> tuple[np.ndarray, int]: cleaned = raw.copy() count, labels, stats, _ = cv2.connectedComponentsWithStats((cleaned > 0).astype(np.uint8), connectivity=8) for component in range(1, count): if int(stats[component, cv2.CC_STAT_AREA]) < minimum_area: cleaned[labels == component] = 0 return cleaned, count - 1 def preview(image: np.ndarray, mask: np.ndarray, generic_mask: np.ndarray | None = None, max_dimension: int = 720) -> Image.Image: height, width = image.shape[:2] scale = min(1.0, max_dimension / max(height, width)) size = (max(1, int(round(width * scale))), max(1, int(round(height * scale)))) base = Image.fromarray(image).resize(size, Image.Resampling.LANCZOS) small_mask = Image.fromarray((mask > 0).astype(np.uint8) * 255).resize(size, Image.Resampling.NEAREST) base_array = np.asarray(base).copy() mask_array = np.asarray(small_mask) > 0 if generic_mask is not None: small_generic = Image.fromarray((generic_mask > 0).astype(np.uint8) * 255).resize(size, Image.Resampling.NEAREST) generic_array = np.asarray(small_generic) > 0 cyan = np.zeros_like(base_array) cyan[..., 1] = 220 cyan[..., 2] = 255 base_array[generic_array] = (base_array[generic_array].astype(np.float32) * 0.45 + cyan[generic_array].astype(np.float32) * 0.55).astype(np.uint8) red = np.zeros_like(base_array) red[..., 0] = 255 base_array[mask_array] = (base_array[mask_array].astype(np.float32) * 0.45 + red[mask_array].astype(np.float32) * 0.55).astype(np.uint8) return Image.fromarray(base_array) def run_scan(run_dir: Path, output_dir: Path, thresholds: list[float], areas: list[int]) -> dict: demo = load_demo() probability_path = run_dir / "change_probability.tif" input_path = run_dir.parents[3] / "data" / "processed" / "00-change-detection" / run_dir.name / "after_registered.tif" if not probability_path.is_file() or not input_path.is_file(): raise FileNotFoundError("The source run must contain change_probability.tif and processed after_registered.tif.") with rasterio.open(probability_path) as probability_dataset: probability = probability_dataset.read(1).astype(np.float32) transform = probability_dataset.transform crs = probability_dataset.crs.to_string() if probability_dataset.crs else None with rasterio.open(input_path) as input_dataset: image = np.transpose(input_dataset.read([1, 2, 3]), (1, 2, 0)) generic_path = run_dir / "generic_difference_mask.tif" generic_mask = None if generic_path.is_file(): with rasterio.open(generic_path) as generic_dataset: generic_mask = generic_dataset.read(1) > 0 output_dir.mkdir(parents=True, exist_ok=False) rows: list[dict] = [] contact_items: list[tuple[str, Image.Image]] = [] for threshold in thresholds: for area in areas: label = f"threshold-{threshold:.2f}_area-{area}" item_dir = output_dir / label item_dir.mkdir() raw = (probability >= threshold).astype(np.uint8) * 255 cleaned, raw_components = clean_mask(raw, area) cleaned_components = int(cv2.connectedComponents((cleaned > 0).astype(np.uint8), connectivity=8)[0] - 1) mask_path = item_dir / "change_mask.tif" with rasterio.open(mask_path, "w", driver="GTiff", height=cleaned.shape[0], width=cleaned.shape[1], count=1, dtype="uint8", transform=transform, crs=crs, compress="lzw") as destination: destination.write(cleaned, 1) # Vectorization is intentionally deferred: full-resolution GeoTIFF # polygons are expensive and are not needed to compare thresholds. features = demo._feature_summary(cleaned, probability, transform, crs is not None) (item_dir / "regions.json").write_text(json.dumps({"features": features}, ensure_ascii=False, indent=2), encoding="utf-8") overlay = preview(image, cleaned, generic_mask) overlay.save(item_dir / "overlay_preview.jpg", quality=90) row = { "threshold": threshold, "minimum_area_pixels": area, "raw_components": raw_components, "cleaned_components": cleaned_components, "raw_changed_pixels": int((raw > 0).sum()), "changed_pixels": int((cleaned > 0).sum()), "changed_pixel_ratio": round(float((cleaned > 0).mean()), 6), "vector_feature_count": len(features), "directory": label, } (item_dir / "summary.json").write_text(json.dumps(row, ensure_ascii=False, indent=2), encoding="utf-8") rows.append(row) contact_items.append((f"T={threshold:.2f} / A={area} | {len(features)} regions", overlay)) (output_dir / "scan_summary.json").write_text(json.dumps({"source_run": run_dir.name, "results": rows}, ensure_ascii=False, indent=2), encoding="utf-8") with (output_dir / "scan_summary.csv").open("w", newline="", encoding="utf-8-sig") as stream: writer = csv.DictWriter(stream, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) tile_width = max(image.width for _, image in contact_items) tile_height = max(image.height for _, image in contact_items) + 30 sheet = Image.new("RGB", (tile_width * 3, tile_height * 4), "#202b25") draw = ImageDraw.Draw(sheet) for index, (label, item) in enumerate(contact_items): x = (index % 3) * tile_width y = (index // 3) * tile_height sheet.paste(item, (x, y + 28)) draw.text((x + 6, y + 6), label, fill="white") sheet.save(output_dir / "parameter_scan_contact_sheet.jpg", quality=92) return {"source_run": run_dir.name, "output_dir": str(output_dir), "results": rows} def main() -> None: parser = argparse.ArgumentParser() parser.add_argument("--run-dir", type=Path, required=True) parser.add_argument("--output", type=Path, required=True) parser.add_argument("--threshold", type=float, action="append", dest="thresholds") parser.add_argument("--minimum-area", type=int, action="append", dest="areas") args = parser.parse_args() thresholds = args.thresholds or [0.3, 0.4, 0.5] areas = args.areas or [16, 64, 256, 686] result = run_scan(args.run_dir, args.output, thresholds, areas) print(json.dumps(result, ensure_ascii=False, indent=2)) if __name__ == "__main__": main()