shuishen
9 hours ago 2ae460fc4a4c2419cf44329783d49a739e2a04ea
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
"""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, 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
    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))
    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)
            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()