"""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())
|
