"""Run an explainable CPU image-region anomaly baseline. The anomaly scores come from lightweight image statistics and Isolation Forest. GeoAI is used for anomaly-mask vectorization, so this remains a B capability: GeoAI supplies the geospatial output workflow while ecosystem code supplies the anomaly model. A high score is only a review candidate, not a business event. """ from __future__ import annotations import argparse import json import math import time from dataclasses import dataclass from datetime import UTC, datetime from importlib import metadata as importlib_metadata from pathlib import Path from typing import Any, Iterable import cv2 import geopandas as gpd import numpy as np import pandas as pd import rasterio from PIL import Image from rasterio.features import geometry_mask, shapes from rasterio.transform import Affine from shapely.geometry import Polygon, mapping, shape from sklearn.ensemble import IsolationForest from sklearn.preprocessing import StandardScaler SUPPORTED_SUFFIXES = {".jpg", ".jpeg", ".png", ".tif", ".tiff"} MASK_LABELS = {1: "rule_only", 2: "isolation_only", 3: "agreement"} MASK_COLORS = { 0: np.array([0, 0, 0], dtype=np.uint8), 1: np.array([255, 158, 44], dtype=np.uint8), 2: np.array([47, 128, 237], dtype=np.uint8), 3: np.array([224, 49, 49], dtype=np.uint8), } FEATURE_NAMES = [ "rgb_mean_r", "rgb_mean_g", "rgb_mean_b", "rgb_std_r", "rgb_std_g", "rgb_std_b", "hsv_mean_h", "hsv_mean_s", "hsv_mean_v", "hsv_std_h", "hsv_std_s", "hsv_std_v", "gray_entropy", "edge_density", "laplacian_variance", "dark_ratio", "bright_ratio", ] @dataclass(frozen=True) class RasterInfo: transform: Affine crs: str | None georeferenced: bool source_bands: list[int] normalization: str def axis_positions(length: int, tile_size: int, stride: int) -> list[int]: """Return deterministic starts and always anchor the final tile to the edge.""" if length < tile_size: raise ValueError(f"Image dimension {length} is smaller than tile size {tile_size}.") positions = list(range(0, length - tile_size + 1, stride)) last = length - tile_size if positions[-1] != last: positions.append(last) return positions def tile_windows(width: int, height: int, tile_size: int, stride: int) -> list[tuple[int, int, int, int]]: if tile_size <= 0 or stride <= 0 or stride > tile_size: raise ValueError("tile_size must be positive and stride must be in 1..tile_size.") return [ (x, y, x + tile_size, y + tile_size) for y in axis_positions(height, tile_size, stride) for x in axis_positions(width, tile_size, stride) ] def _normalize_band(band: np.ndarray) -> np.ndarray: values = band.astype(np.float32) finite = values[np.isfinite(values)] if finite.size == 0: return np.zeros(band.shape, dtype=np.uint8) low, high = np.percentile(finite, [2.0, 98.0]) if high <= low: high = low + 1.0 return np.clip((values - low) * 255.0 / (high - low), 0, 255).astype(np.uint8) def read_rgb(path: Path) -> tuple[np.ndarray, RasterInfo]: """Read common imagery and retain spatial metadata only when it is real.""" if path.suffix.lower() in {".tif", ".tiff"}: with rasterio.open(path) as src: if src.count < 3: raise ValueError(f"GeoTIFF requires at least three bands: {path.name}") raw = src.read([1, 2, 3]) if raw.dtype == np.uint8: rgb = np.moveaxis(raw, 0, 2) normalization = "uint8_identity" else: rgb = np.stack([_normalize_band(raw[index]) for index in range(3)], axis=2) normalization = "per_band_percentile_2_98" transform = src.transform crs = src.crs.to_string() if src.crs else None georeferenced = bool(src.crs and src.transform != Affine.identity()) return rgb, RasterInfo(transform, crs, georeferenced, [1, 2, 3], normalization) with Image.open(path) as image: rgb = np.asarray(image.convert("RGB")) return rgb, RasterInfo(Affine.identity(), None, False, [1, 2, 3], "pillow_rgb") def extract_features(tile: np.ndarray) -> np.ndarray: """Extract explainable colour, texture, edge, focus, and exposure features.""" if tile.ndim != 3 or tile.shape[2] != 3: raise ValueError("Expected an RGB tile with shape (height, width, 3).") rgb = tile.astype(np.float32) hsv = cv2.cvtColor(tile, cv2.COLOR_RGB2HSV).astype(np.float32) gray = cv2.cvtColor(tile, cv2.COLOR_RGB2GRAY) histogram = cv2.calcHist([gray], [0], None, [32], [0, 256]).ravel().astype(np.float64) probabilities = histogram / max(float(histogram.sum()), 1.0) probabilities = probabilities[probabilities > 0] entropy = float(-(probabilities * np.log2(probabilities)).sum()) edges = cv2.Canny(gray, 80, 160) laplacian_variance = float(cv2.Laplacian(gray, cv2.CV_64F).var()) values = [ *rgb.mean(axis=(0, 1)).tolist(), *rgb.std(axis=(0, 1)).tolist(), *hsv.mean(axis=(0, 1)).tolist(), *hsv.std(axis=(0, 1)).tolist(), entropy, float(np.mean(edges > 0)), laplacian_variance, float(np.mean(gray < 32)), float(np.mean(gray > 223)), ] return np.asarray(values, dtype=np.float64) def image_features( rgb: np.ndarray, tile_size: int, stride: int ) -> tuple[np.ndarray, list[tuple[int, int, int, int]]]: height, width = rgb.shape[:2] windows = tile_windows(width, height, tile_size, stride) features = np.vstack([extract_features(rgb[y0:y1, x0:x1]) for x0, y0, x1, y1 in windows]) return features, windows def robust_rule_scores( features: np.ndarray, center: np.ndarray, scale: np.ndarray, feature_cutoffs: np.ndarray | None = None, ) -> tuple[np.ndarray, np.ndarray]: z_scores = np.abs((features - center) / scale) normalized = z_scores if feature_cutoffs is None else z_scores / feature_cutoffs return normalized.max(axis=1), normalized.argmax(axis=1) def _spatial_profile(features: np.ndarray) -> np.ndarray: """Remove image-wide colour/lighting shifts while retaining local spatial patterns.""" center = np.median(features, axis=0) mad = 1.4826 * np.median(np.abs(features - center), axis=0) fallback = features.std(axis=0) scale = np.where(mad > 1e-9, mad, np.where(fallback > 1e-9, fallback, 1.0)) return (features - center) / scale def aligned_rule_scores( features: np.ndarray, center: np.ndarray, scale: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: """Compare every target tile only with the same tile position in the references.""" profile = _spatial_profile(features) z_scores = np.abs((profile - center) / scale) top_count = min(3, z_scores.shape[1]) scores = np.sort(z_scores, axis=1)[:, -top_count:].mean(axis=1) return scores, z_scores.argmax(axis=1) def fit_aligned_rule_model( reference_feature_sets: list[np.ndarray], threshold_quantile: float ) -> dict[str, Any]: """Fit a fixed-camera rule model and calibrate it by leave-one-reference-out scores.""" if len(reference_feature_sets) < 3: raise ValueError("Aligned spatial mode requires at least three reference images.") if len({values.shape for values in reference_feature_sets}) != 1: raise ValueError("Aligned spatial mode requires identical reference image dimensions.") profiles = np.stack([_spatial_profile(values) for values in reference_feature_sets]) feature_floor = np.maximum( profiles.reshape(-1, profiles.shape[-1]).std(axis=0) * 0.15, 0.1 ) def spatial_stats(values: np.ndarray) -> tuple[np.ndarray, np.ndarray]: center = np.median(values, axis=0) mad = 1.4826 * np.median(np.abs(values - center), axis=0) return center, np.maximum(mad, feature_floor) calibration_scores: list[np.ndarray] = [] for index in range(profiles.shape[0]): center, scale = spatial_stats(np.delete(profiles, index, axis=0)) z_scores = np.abs((profiles[index] - center) / scale) top_count = min(3, z_scores.shape[1]) calibration_scores.append(np.sort(z_scores, axis=1)[:, -top_count:].mean(axis=1)) center, scale = spatial_stats(profiles) reference_scores = np.concatenate(calibration_scores) return { "rule_mode": "aligned", "aligned_center": center, "aligned_scale": scale, "rule_threshold": float(np.quantile(reference_scores, threshold_quantile)), "reference_rule_scores": reference_scores, } def _robust_image_channel(values: np.ndarray) -> np.ndarray: low, high = np.percentile(values, [10.0, 90.0]) scale = max(float(high - low), 1e-3) return (values - np.median(values)) / scale def _local_appearance_profile(rgb: np.ndarray) -> np.ndarray: """Represent generic local colour and luminance without any class/colour rules.""" lab = cv2.cvtColor(rgb, cv2.COLOR_RGB2LAB).astype(np.float32) normalized = np.stack( [_robust_image_channel(lab[:, :, index]) for index in range(3)], axis=2 ).astype(np.float32) return cv2.GaussianBlur(normalized, (0, 0), 2.0) def _local_structure_profile(rgb: np.ndarray) -> np.ndarray: """Represent local intensity, edges, and contrast for class-agnostic changes.""" gray = cv2.cvtColor(rgb, cv2.COLOR_RGB2GRAY).astype(np.float32) / 255.0 normalized = _robust_image_channel(gray).astype(np.float32) smooth = cv2.GaussianBlur(normalized, (0, 0), 2.0) gradient_x = cv2.Sobel(smooth, cv2.CV_32F, 1, 0, ksize=3) gradient_y = cv2.Sobel(smooth, cv2.CV_32F, 0, 1, ksize=3) gradient = np.sqrt(gradient_x * gradient_x + gradient_y * gradient_y) local_mean = cv2.GaussianBlur(normalized, (0, 0), 4.0) local_square_mean = cv2.GaussianBlur(normalized * normalized, (0, 0), 4.0) local_contrast = np.sqrt(np.maximum(local_square_mean - local_mean * local_mean, 0.0)) return np.stack([smooth, gradient, local_contrast], axis=2).astype(np.float32) def _nearest_profile_distance(profile: np.ndarray, references: list[np.ndarray]) -> np.ndarray: distances = [np.sqrt(np.sum((profile - reference) ** 2, axis=2)) for reference in references] return np.min(np.stack(distances), axis=0).astype(np.float32) def _local_processing_shape(shape_: tuple[int, int], max_dimension: int = 2048) -> tuple[int, int]: height, width = shape_ scale = min(1.0, max_dimension / max(height, width)) return max(1, round(height * scale)), max(1, round(width * scale)) def _resize_rgb(rgb: np.ndarray, shape_: tuple[int, int]) -> np.ndarray: if rgb.shape[:2] == shape_: return rgb return cv2.resize(rgb, (shape_[1], shape_[0]), interpolation=cv2.INTER_AREA) def fit_local_change_model( reference_images: list[np.ndarray], threshold_quantile: float ) -> dict[str, Any]: """Calibrate generic pixel-local appearance/structure change from normal references.""" if len(reference_images) < 3 or len({image.shape for image in reference_images}) != 1: raise ValueError("Local change scoring requires at least three equal-sized RGB references.") original_shape = reference_images[0].shape[:2] processing_shape = _local_processing_shape(original_shape) resized = [_resize_rgb(image, processing_shape) for image in reference_images] appearance_profiles = [_local_appearance_profile(image) for image in resized] structure_profiles = [_local_structure_profile(image) for image in resized] def calibrate(profiles: list[np.ndarray]) -> tuple[float, list[float]]: scores: list[np.ndarray] = [] for index, profile in enumerate(profiles): others = [value for other_index, value in enumerate(profiles) if other_index != index] scores.append(_nearest_profile_distance(profile, others)) merged = np.concatenate([values.ravel() for values in scores]) threshold = float(np.quantile(merged, threshold_quantile)) coverages = [float(np.mean(values > threshold)) for values in scores] return max(threshold, 1e-6), coverages appearance_threshold, appearance_coverages = calibrate(appearance_profiles) structure_threshold, structure_coverages = calibrate(structure_profiles) return { "local_change_enabled": True, "local_original_shape": original_shape, "local_processing_shape": processing_shape, "local_appearance_profiles": appearance_profiles, "local_structure_profiles": structure_profiles, "local_appearance_threshold": appearance_threshold, "local_structure_threshold": structure_threshold, "local_reference_appearance_coverages": appearance_coverages, "local_reference_structure_coverages": structure_coverages, } def local_change_scores(rgb: np.ndarray, models: dict[str, Any]) -> tuple[np.ndarray, np.ndarray]: """Return threshold-normalized generic local scores and dominant-channel codes.""" if rgb.shape[:2] != tuple(models["local_original_shape"]): raise ValueError("Local change scoring requires the target to match reference dimensions.") resized = _resize_rgb(rgb, tuple(models["local_processing_shape"])) appearance = _nearest_profile_distance( _local_appearance_profile(resized), models["local_appearance_profiles"] ) / models["local_appearance_threshold"] structure = _nearest_profile_distance( _local_structure_profile(resized), models["local_structure_profiles"] ) / models["local_structure_threshold"] if resized.shape[:2] != rgb.shape[:2]: size = (rgb.shape[1], rgb.shape[0]) appearance = cv2.resize(appearance, size, interpolation=cv2.INTER_LINEAR) structure = cv2.resize(structure, size, interpolation=cv2.INTER_LINEAR) return np.maximum(appearance, structure), (structure > appearance).astype(np.uint8) def clean_local_change_flag(score_map: np.ndarray, minimum_area: int = 512) -> np.ndarray: """Discard isolated/thin registration noise while retaining coherent anomalies.""" count, labels, stats, _ = cv2.connectedComponentsWithStats( (score_map > 1.0).astype(np.uint8), connectivity=8 ) cleaned = np.zeros(score_map.shape, dtype=bool) for label in range(1, count): area = int(stats[label, cv2.CC_STAT_AREA]) width = int(stats[label, cv2.CC_STAT_WIDTH]) height = int(stats[label, cv2.CC_STAT_HEIGHT]) fill_ratio = area / max(width * height, 1) coherent_small_region = fill_ratio >= 0.1 and min(width, height) >= 16 if area >= minimum_area and (area >= 4096 or coherent_small_region): cleaned[labels == label] = True return cleaned def fit_models( reference_features: np.ndarray, threshold_quantile: float, random_state: int, aligned_reference_features: list[np.ndarray] | None = None, aligned_reference_images: list[np.ndarray] | None = None, ) -> dict[str, Any]: if reference_features.shape[0] < 50: raise ValueError( f"At least 50 reference tiles are required; found {reference_features.shape[0]}." ) if not 0.9 <= threshold_quantile < 1.0: raise ValueError("threshold_quantile must be in [0.9, 1.0).") center = np.median(reference_features, axis=0) mad = np.median(np.abs(reference_features - center), axis=0) fallback = reference_features.std(axis=0) scale = np.where(mad > 1e-9, 1.4826 * mad, np.where(fallback > 1e-9, fallback, 1.0)) reference_z = np.abs((reference_features - center) / scale) feature_cutoffs = np.quantile(reference_z, threshold_quantile, axis=0) feature_cutoffs = np.maximum(feature_cutoffs, 1.0) rule_scores, _ = robust_rule_scores(reference_features, center, scale, feature_cutoffs) scaler = StandardScaler().fit(reference_features) standardized = scaler.transform(reference_features) isolation = IsolationForest( n_estimators=200, contamination="auto", random_state=random_state, n_jobs=-1, ).fit(standardized) isolation_scores = -isolation.score_samples(standardized) result = { "rule_mode": "global", "center": center, "scale": scale, "feature_cutoffs": feature_cutoffs, "scaler": scaler, "isolation": isolation, "rule_threshold": float(np.quantile(rule_scores, threshold_quantile)), "isolation_threshold": float(np.quantile(isolation_scores, threshold_quantile)), "reference_rule_scores": rule_scores, "reference_isolation_scores": isolation_scores, } if aligned_reference_features is not None: result.update(fit_aligned_rule_model(aligned_reference_features, threshold_quantile)) if aligned_reference_images is not None: result.update(fit_local_change_model(aligned_reference_images, threshold_quantile)) return result def score_features(features: np.ndarray, models: dict[str, Any]) -> tuple[np.ndarray, np.ndarray, np.ndarray]: if models.get("rule_mode") == "aligned": if features.shape[0] != models["aligned_center"].shape[0]: raise ValueError( "Aligned spatial mode requires target images with the same dimensions as references." ) rule_scores, reason_indices = aligned_rule_scores( features, models["aligned_center"], models["aligned_scale"] ) else: rule_scores, reason_indices = robust_rule_scores( features, models["center"], models["scale"], models["feature_cutoffs"] ) isolation_scores = -models["isolation"].score_samples(models["scaler"].transform(features)) return rule_scores, isolation_scores, reason_indices def scores_to_raster( scores: np.ndarray, windows: list[tuple[int, int, int, int]], shape_: tuple[int, int], ) -> np.ndarray: sums = np.zeros(shape_, dtype=np.float32) counts = np.zeros(shape_, dtype=np.uint16) for score, (x0, y0, x1, y1) in zip(scores, windows, strict=True): sums[y0:y1, x0:x1] += float(score) counts[y0:y1, x0:x1] += 1 return sums / np.maximum(counts, 1) def _heatmap(score_map: np.ndarray, threshold: float) -> np.ndarray: scaled = np.clip(score_map / max(threshold, 1e-9), 0.0, 2.0) * 127.5 colored_bgr = cv2.applyColorMap(scaled.astype(np.uint8), cv2.COLORMAP_TURBO) return cv2.cvtColor(colored_bgr, cv2.COLOR_BGR2RGB) def _write_mask_raster(path: Path, mask: np.ndarray, info: RasterInfo) -> None: with rasterio.open( path, "w", driver="GTiff", height=mask.shape[0], width=mask.shape[1], count=1, dtype="uint8", transform=info.transform, crs=info.crs, nodata=0, compress="deflate", ) as dst: dst.write(mask, 1) def _fallback_geometries(binary: np.ndarray, transform: Affine, crs: str | None) -> gpd.GeoDataFrame: records: list[dict[str, Any]] = [] for geometry, value in shapes(binary.astype(np.uint8), mask=binary, transform=transform): if value != 1: continue polygon = shape(geometry) if polygon.is_empty: continue records.append({"geometry": polygon}) return gpd.GeoDataFrame(records, geometry="geometry", crs=crs) def _polygon_score_stats( geometry: Any, transform: Affine, shape_: tuple[int, int], rule_map: np.ndarray, isolation_map: np.ndarray, ) -> tuple[float, float, int]: inverse = ~transform corners = [ inverse * (geometry.bounds[0], geometry.bounds[1]), inverse * (geometry.bounds[0], geometry.bounds[3]), inverse * (geometry.bounds[2], geometry.bounds[1]), inverse * (geometry.bounds[2], geometry.bounds[3]), ] xs = [point[0] for point in corners] ys = [point[1] for point in corners] x0 = max(0, int(math.floor(min(xs)))) x1 = min(shape_[1], int(math.ceil(max(xs))) + 1) y0 = max(0, int(math.floor(min(ys)))) y1 = min(shape_[0], int(math.ceil(max(ys))) + 1) if x1 <= x0 or y1 <= y0: return 0.0, 0.0, 0 local_transform = transform * Affine.translation(x0, y0) inside = geometry_mask( [mapping(geometry)], out_shape=(y1 - y0, x1 - x0), transform=local_transform, invert=True, ) if not inside.any(): return 0.0, 0.0, 0 return ( float(rule_map[y0:y1, x0:x1][inside].max()), float(isolation_map[y0:y1, x0:x1][inside].max()), int(inside.sum()), ) def vectorize_candidates( mask: np.ndarray, rule_map: np.ndarray, isolation_map: np.ndarray, transform: Affine, crs: str | None, tile_rows: list[dict[str, Any]], output_path: Path, minimum_area: int, source_file: str, ) -> tuple[int, str]: frames: list[gpd.GeoDataFrame] = [] vectorizers: list[str] = [] for code in (1, 2, 3): binary = mask == code if not binary.any(): continue reference = _fallback_geometries(binary, transform, crs) temporary = output_path.with_name(f".{output_path.stem}-{code}.tif") _write_mask_raster( temporary, np.where(binary, 255, 0).astype(np.uint8), RasterInfo(transform, crs, bool(crs and transform != Affine.identity()), [1], "binary"), ) selected = reference source = "rasterio.features.shapes fallback" try: from geoai import masks_to_vector candidate = masks_to_vector( str(temporary), min_object_area=minimum_area, simplify_tolerance=2.0 ) pixel_area = abs(transform.a * transform.e - transform.b * transform.d) or 1.0 reference = reference[ reference.geometry.area / pixel_area >= minimum_area ].copy() reference_area = float(reference.geometry.area.sum()) candidate_area = float(candidate.geometry.area.sum()) if not candidate.empty else 0.0 ratio = candidate_area / reference_area if reference_area else 0.0 if len(candidate) == len(reference) and (not reference_area or 0.9 <= ratio <= 1.1): selected = candidate[["geometry"]].copy() source = "geoai.masks_to_vector" else: selected = reference source = "geoai.masks_to_vector + rasterio completeness repair" except Exception: pixel_area = abs(transform.a * transform.e - transform.b * transform.d) or 1.0 selected = reference[reference.geometry.area / pixel_area >= minimum_area].copy() finally: temporary.unlink(missing_ok=True) vectorizers.append(source) rows: list[dict[str, Any]] = [] for _, record in selected.iterrows(): geometry = record.geometry max_rule, max_isolation, area_pixels = _polygon_score_stats( geometry, transform, mask.shape, rule_map, isolation_map ) reason = "" best_rule = -math.inf for tile in tile_rows: corners = [ transform * (tile["x0"], tile["y0"]), transform * (tile["x1"], tile["y0"]), transform * (tile["x1"], tile["y1"]), transform * (tile["x0"], tile["y1"]), ] tile_geometry = Polygon(corners) if geometry.intersects(tile_geometry) and tile["rule_score"] > best_rule: best_rule = float(tile["rule_score"]) reason = str(tile["reason_feature"]) rows.append( { "geometry": geometry, "source_file": source_file, "mask_code": code, "method": MASK_LABELS[code], "rule_flag": code in (1, 3), "isolation_flag": code in (2, 3), "agreement": code == 3, "max_rule_score": round(max_rule, 6), "max_isolation_score": round(max_isolation, 6), "reason_feature": reason, "area_pixels": area_pixels, } ) if rows: frames.append(gpd.GeoDataFrame(rows, geometry="geometry", crs=crs)) if frames: merged = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), geometry="geometry", crs=crs) merged.insert(0, "feature_id", np.arange(1, len(merged) + 1)) else: merged = gpd.GeoDataFrame( { "feature_id": pd.Series(dtype="int64"), "mask_code": pd.Series(dtype="int64"), "method": pd.Series(dtype="str"), "geometry": gpd.GeoSeries([], crs=crs), }, geometry="geometry", crs=crs, ) merged.to_file(output_path, driver="GeoJSON") unique_vectorizers = set(vectorizers) if "geoai.masks_to_vector + rasterio completeness repair" in unique_vectorizers: vectorizer = "geoai.masks_to_vector + rasterio completeness repair" elif "geoai.masks_to_vector" in unique_vectorizers: vectorizer = "geoai.masks_to_vector" elif unique_vectorizers: vectorizer = "rasterio.features.shapes fallback" else: vectorizer = "no candidates" return len(merged), vectorizer def _truth_entry(manifest: dict[str, Any] | None, filename: str) -> dict[str, Any] | None: if not manifest: return None images = manifest.get("images", {}) entry = images.get(filename) return entry if isinstance(entry, dict) else None def _binary_metrics(predicted: np.ndarray, truth: np.ndarray) -> dict[str, float | int]: predicted = predicted.astype(bool) truth = truth.astype(bool) true_positive = int(np.count_nonzero(predicted & truth)) false_positive = int(np.count_nonzero(predicted & ~truth)) false_negative = int(np.count_nonzero(~predicted & truth)) precision = true_positive / (true_positive + false_positive) if true_positive + false_positive else 0.0 recall = true_positive / (true_positive + false_negative) if true_positive + false_negative else 0.0 union = true_positive + false_positive + false_negative return { "true_positive_pixels": true_positive, "false_positive_pixels": false_positive, "false_negative_pixels": false_negative, "precision": round(precision, 6), "recall": round(recall, 6), "iou": round(true_positive / union if union else 0.0, 6), } def evaluate_truth( mask: np.ndarray, entry: dict[str, Any], manifest_dir: Path, expected_shape: tuple[int, int], ) -> dict[str, Any]: truth_path = (manifest_dir / str(entry["mask"])).resolve() with Image.open(truth_path) as image: truth = np.asarray(image.convert("L")) > 0 if truth.shape != expected_shape: raise ValueError(f"Truth mask shape {truth.shape} does not match image {expected_shape}.") methods = { "rule": np.isin(mask, [1, 3]), "isolation": np.isin(mask, [2, 3]), "agreement": mask == 3, } result: dict[str, Any] = {name: _binary_metrics(values, truth) for name, values in methods.items()} region_results: dict[str, list[dict[str, Any]]] = {} for name, values in methods.items(): regions: list[dict[str, Any]] = [] for region in entry.get("regions", []): x0, y0, x1, y1 = map(int, region["bbox"]) region_truth = truth[y0:y1, x0:x1] detected = values[y0:y1, x0:x1] overlap = int(np.count_nonzero(region_truth & detected)) truth_pixels = int(np.count_nonzero(region_truth)) regions.append( { "name": region["name"], "overlap_pixels": overlap, "truth_pixels": truth_pixels, "overlap_ratio": round(overlap / truth_pixels if truth_pixels else 0.0, 6), "hit": overlap > 0, } ) region_results[name] = regions result["regions"] = region_results return result def process_image( input_path: Path, output_dir: Path, models: dict[str, Any], tile_size: int, stride: int, truth_manifest: dict[str, Any] | None, truth_manifest_dir: Path | None, ) -> dict[str, Any]: started = time.perf_counter() rgb, raster_info = read_rgb(input_path) features, windows = image_features(rgb, tile_size, stride) rule_scores, isolation_scores, reason_indices = score_features(features, models) coarse_rule_map = scores_to_raster(rule_scores, windows, rgb.shape[:2]) isolation_map = scores_to_raster(isolation_scores, windows, rgb.shape[:2]) rule_map = coarse_rule_map local_score_map: np.ndarray | None = None local_reason_map: np.ndarray | None = None local_flag = np.zeros(rgb.shape[:2], dtype=bool) if models.get("local_change_enabled"): local_score_map, local_reason_map = local_change_scores(rgb, models) local_flag = clean_local_change_flag(local_score_map) local_equivalent = local_score_map * models["rule_threshold"] rule_map = np.maximum(coarse_rule_map, local_equivalent) rule_flag = (coarse_rule_map > models["rule_threshold"]) | local_flag isolation_flag = isolation_map > models["isolation_threshold"] encoded = rule_flag.astype(np.uint8) + isolation_flag.astype(np.uint8) * 2 stem = input_path.stem rule_heatmap_path = output_dir / f"{stem}.rule.heatmap.png" isolation_heatmap_path = output_dir / f"{stem}.isolation.heatmap.png" overlay_path = output_dir / f"{stem}.comparison.overlay.png" mask_path = output_dir / f"{stem}.anomaly-mask.tif" vector_path = output_dir / f"{stem}.anomalies.geojson" csv_path = output_dir / f"{stem}.tiles.csv" Image.fromarray(_heatmap(rule_map, models["rule_threshold"])).save(rule_heatmap_path) Image.fromarray(_heatmap(isolation_map, models["isolation_threshold"])).save( isolation_heatmap_path ) color_mask = np.zeros_like(rgb) for code, color in MASK_COLORS.items(): color_mask[encoded == code] = color overlay = rgb.copy() candidates = encoded > 0 overlay[candidates] = ( rgb[candidates].astype(np.float32) * 0.42 + color_mask[candidates].astype(np.float32) * 0.58 ).clip(0, 255).astype(np.uint8) Image.fromarray(overlay).save(overlay_path, quality=92) _write_mask_raster(mask_path, encoded, raster_info) tile_rows: list[dict[str, Any]] = [] for index, ((x0, y0, x1, y1), feature_values) in enumerate( zip(windows, features, strict=True) ): output_rule_score = float(rule_scores[index]) output_reason = FEATURE_NAMES[int(reason_indices[index])] local_tile_score = 0.0 if local_score_map is not None and local_reason_map is not None: local_values = local_score_map[y0:y1, x0:x1] local_tile_score = float(local_values.max()) local_equivalent = local_tile_score * models["rule_threshold"] if local_equivalent > output_rule_score: output_rule_score = local_equivalent maximum_position = np.unravel_index(int(local_values.argmax()), local_values.shape) output_reason = ( "local_structure_change" if local_reason_map[y0:y1, x0:x1][maximum_position] == 1 else "local_appearance_change" ) row: dict[str, Any] = { "tile_id": index + 1, "x0": x0, "y0": y0, "x1": x1, "y1": y1, "rule_score": output_rule_score, "local_change_score": local_tile_score, "isolation_score": float(isolation_scores[index]), "rule_anomaly": bool(rule_flag[y0:y1, x0:x1].any()), "isolation_anomaly": bool(isolation_scores[index] > models["isolation_threshold"]), "reason_feature": output_reason, } row.update({name: float(value) for name, value in zip(FEATURE_NAMES, feature_values, strict=True)}) tile_rows.append(row) pd.DataFrame(tile_rows).to_csv(csv_path, index=False, encoding="utf-8-sig") candidate_count, vectorizer = vectorize_candidates( encoded, rule_map, isolation_map, raster_info.transform, raster_info.crs, tile_rows, vector_path, minimum_area=512 if models.get("local_change_enabled") else max(64, tile_size * tile_size // 16), source_file=input_path.name, ) truth_metrics = None entry = _truth_entry(truth_manifest, input_path.name) if entry and truth_manifest_dir: truth_metrics = evaluate_truth(encoded, entry, truth_manifest_dir, rgb.shape[:2]) pixel_count = int(encoded.size) return { "file": input_path.name, "width": int(rgb.shape[1]), "height": int(rgb.shape[0]), "tile_count": len(windows), "rule_heatmap_file": rule_heatmap_path.name, "isolation_heatmap_file": isolation_heatmap_path.name, "overlay_file": overlay_path.name, "mask_file": mask_path.name, "vector_file": vector_path.name, "tiles_file": csv_path.name, "candidate_count": candidate_count, "rule_anomaly_pixels": int(np.count_nonzero(rule_flag)), "local_change_pixels": int(np.count_nonzero(local_flag)), "isolation_anomaly_pixels": int(np.count_nonzero(isolation_flag)), "agreement_pixels": int(np.count_nonzero(encoded == 3)), "rule_anomaly_coverage": round(float(np.count_nonzero(rule_flag) / pixel_count), 6), "local_change_coverage": round(float(np.count_nonzero(local_flag) / pixel_count), 6), "isolation_anomaly_coverage": round( float(np.count_nonzero(isolation_flag) / pixel_count), 6 ), "agreement_coverage": round(float(np.count_nonzero(encoded == 3) / pixel_count), 6), "georeferenced": raster_info.georeferenced, "crs": raster_info.crs, "source_bands": raster_info.source_bands, "normalization": raster_info.normalization, "coordinate_basis": "map_coordinates" if raster_info.georeferenced else "pixel_coordinates", "vectorizer": vectorizer, "truth_metrics": truth_metrics, "elapsed_seconds": round(time.perf_counter() - started, 3), } def collect_inputs(path: Path) -> list[Path]: if path.is_file(): candidates: Iterable[Path] = [path] elif path.is_dir(): candidates = sorted(item for item in path.iterdir() if item.is_file()) else: raise ValueError(f"Input path does not exist: {path}") inputs = [item for item in candidates if item.suffix.lower() in SUPPORTED_SUFFIXES] if not inputs: raise ValueError(f"No JPG, PNG, or GeoTIFF inputs were found in: {path}") return inputs def load_truth_manifest(path: Path | None) -> tuple[dict[str, Any] | None, Path | None]: if path is None: return None, None payload = json.loads(path.read_text(encoding="utf-8")) if not isinstance(payload.get("images"), dict): raise ValueError("Truth manifest must contain an 'images' object.") return payload, path.parent def run(args: argparse.Namespace) -> dict[str, Any]: if args.output.exists() and any(args.output.iterdir()): raise ValueError(f"Output directory is not empty: {args.output}. Use a new run directory.") args.output.mkdir(parents=True, exist_ok=True) reference_paths = collect_inputs(args.reference) input_paths = collect_inputs(args.input) started = time.perf_counter() reference_feature_sets: list[np.ndarray] = [] reference_rgbs: list[np.ndarray] = [] reference_images: list[dict[str, Any]] = [] reference_shapes: list[tuple[int, int]] = [] for path in reference_paths: rgb, _ = read_rgb(path) features, windows = image_features(rgb, args.tile_size, args.stride) reference_feature_sets.append(features) reference_rgbs.append(rgb) reference_shapes.append(rgb.shape[:2]) reference_images.append( { "file": path.name, "tile_count": len(windows), "width": rgb.shape[1], "height": rgb.shape[0], } ) reference_features = np.vstack(reference_feature_sets) requested_spatial_mode = getattr(args, "spatial_mode", "auto") if requested_spatial_mode not in {"auto", "global", "aligned"}: raise ValueError("spatial_mode must be auto, global, or aligned.") input_shapes = [read_rgb(path)[0].shape[:2] for path in input_paths] aligned_compatible = ( len(reference_paths) >= 3 and len(set(reference_shapes)) == 1 and all(shape_ == reference_shapes[0] for shape_ in input_shapes) ) if requested_spatial_mode == "aligned" and not aligned_compatible: raise ValueError( "Aligned spatial mode requires at least three references and all reference/target images " "to have identical dimensions." ) selected_rule_mode = ( "aligned" if requested_spatial_mode == "aligned" or (requested_spatial_mode == "auto" and aligned_compatible) else "global" ) models = fit_models( reference_features, args.threshold_quantile, args.random_state, reference_feature_sets if selected_rule_mode == "aligned" else None, reference_rgbs if selected_rule_mode == "aligned" else None, ) truth_manifest, truth_manifest_dir = load_truth_manifest(args.truth_manifest) images = [ process_image( path, args.output, models, args.tile_size, args.stride, truth_manifest, truth_manifest_dir, ) for path in input_paths ] metadata = { "capability": "09-anomaly-detection", "classification": "B", "created_at": datetime.now(UTC).isoformat(), "versions": { "python": __import__("platform").python_version(), "geoai-py": importlib_metadata.version("geoai-py"), "scikit-learn": importlib_metadata.version("scikit-learn"), "rasterio": importlib_metadata.version("rasterio"), "opencv-python-headless": importlib_metadata.version("opencv-python-headless"), }, "method": ( "position-aligned robust feature and generic local-change rules compared with global Isolation Forest" if models["rule_mode"] == "aligned" else "global robust feature rules compared with Isolation Forest" ), "model": "IsolationForest(n_estimators=200, contamination='auto')", "device": "CPU", "reference_path": args.reference.as_posix(), "input_path": args.input.as_posix(), "reference_count": len(reference_paths), "reference_tile_count": int(reference_features.shape[0]), "input_count": len(input_paths), "feature_names": FEATURE_NAMES, "parameters": { "tile_size": args.tile_size, "stride": args.stride, "threshold_quantile": args.threshold_quantile, "random_state": args.random_state, "minimum_reference_tiles": 50, "spatial_mode_requested": requested_spatial_mode, "rule_mode_selected": models["rule_mode"], }, "thresholds": { "rule_score": round(float(models["rule_threshold"]), 8), "isolation_score": round(float(models["isolation_threshold"]), 8), "local_appearance_normalized": 1.0 if models.get("local_change_enabled") else None, "local_structure_normalized": 1.0 if models.get("local_change_enabled") else None, }, "local_change": { "enabled": bool(models.get("local_change_enabled")), "method": "nearest normal-reference distance in robust Lab appearance and grayscale structure profiles", "processing_shape": list(models["local_processing_shape"]) if models.get("local_change_enabled") else None, "appearance_threshold": round(float(models["local_appearance_threshold"]), 8) if models.get("local_change_enabled") else None, "structure_threshold": round(float(models["local_structure_threshold"]), 8) if models.get("local_change_enabled") else None, "minimum_component_pixels": 512, "small_component_minimum_fill_ratio": 0.1, "class_or_colour_rules": False, }, "reference_images": reference_images, "images": images, "truth_manifest": args.truth_manifest.as_posix() if args.truth_manifest else None, "elapsed_seconds": round(time.perf_counter() - started, 3), "limitations": [ "异常表示相对正常参考影像的统计离群，只是人工复核候选，不是堵塞、损坏、渗漏或告警结论。", "受控注入异常只能验证检测链路，不能代表真实矿山场景精度。", "自然材质稀有区域、视角、光照、季节和清晰度变化都可能产生误报。", "同位置规则模式要求固定机位和相同像素尺寸；auto 条件不满足时自动回退到全局规则。", "局部变化通道不使用颜色或物品类别规则；显著机位偏移、移动阴影和生成式影像细节漂移仍可能造成误报。", "普通 JPG/PNG 使用像素坐标；只有含有效 CRS 和仿射变换的 GeoTIFF 才保留地图坐标。", ], } (args.output / "run_metadata.json").write_text( json.dumps(metadata, ensure_ascii=False, indent=2), encoding="utf-8" ) return metadata def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser(description="Run the CPU image-region anomaly baseline.") parser.add_argument("--reference", type=Path, required=True, help="Normal reference image or directory.") parser.add_argument("--input", type=Path, required=True, help="Target image or directory.") parser.add_argument("--output", type=Path, required=True, help="A new, empty output directory.") parser.add_argument("--tile-size", type=int, default=256) parser.add_argument("--stride", type=int, default=128) parser.add_argument("--threshold-quantile", type=float, default=0.995) parser.add_argument("--random-state", type=int, default=42) parser.add_argument( "--spatial-mode", choices=("auto", "global", "aligned"), default="auto", help="auto uses same-position rules for 3+ aligned references, otherwise global rules.", ) parser.add_argument("--truth-manifest", type=Path) return parser def main() -> int: args = build_parser().parse_args() try: metadata = run(args) except (OSError, ValueError, json.JSONDecodeError) as error: print(f"ERROR: {error}", file=__import__("sys").stderr) return 2 print(json.dumps(metadata, ensure_ascii=False, indent=2)) return 0 if __name__ == "__main__": raise SystemExit(main())