""" Module: process_image Description: Official Problem-Statement-Compliant Elevation Processing Engine for DepthWizard M2. Core Semantics: 1. MODE 1 — NON-GEOREFERENCED (PNG/JPG): - Relative Digital Surface Model (rDSM in [0, 1]) - mode = "relative", calibrated = False - crs = None, transform = None, gsd_x = None, gsd_y = None - Values are explicitly relative, not metric elevations. 2. MODE 2 — GEOREFERENCED (GeoTIFF): - Preserves valid CRS, affine transform, and GSD. - Spatial reference does NOT provide elevation. - If GCPs provided: calibrated = True, method = "gcp", metric DSM in metres. - If DEM provided: calibrated = True, method = "dem", metric DSM in metres. - If NO calibration: calibrated = False, method = "none", rDSM in [0, 1], NO hardcoded priors, NO invented elevation. """ import os import tempfile import rasterio import numpy as np from PIL import Image from typing import Optional, List, Dict, Any, Union, Tuple from rasterio.transform import Affine from rasterio.warp import reproject, Resampling from depth.depth_model import estimate_depth, load_model from depth.tiled_inference import estimate_depth_tiled from calibration.rdsm import make_rdsm from shadow_detection import compute_shadow_confidence from calibration.metric import ( fit_scale_offset, apply_scale_offset, fit_poly_calibration, apply_poly_calibration, fit_isotonic_calibration, apply_isotonic_calibration, calibrate_depth_to_dsm_nonlinear ) from calibration.gcp_calibration import ( fit_gcp_calibration_nonlinear, apply_gcp_calibration_nonlinear ) from depth.sat2map_refiner import ( regularize_roof_planes, extract_planar_segments, export_point_cloud_ply, export_mesh_obj ) # Integration hook for M7 Guided Filter Module try: from shadow.guided_filter import refine_depth_anything_map HAS_GUIDED_FILTER = True except Exception: HAS_GUIDED_FILTER = False # Integration hook for M4 Shadow-Cue Module import sys shadow_cue_path = os.path.join(os.path.dirname(__file__), "shadow-cue") if os.path.isdir(shadow_cue_path) and shadow_cue_path not in sys.path: sys.path.insert(0, shadow_cue_path) try: from shadow.run_full_pipeline import run_full_pipeline as run_m4_shadow_pipeline HAS_M4_SHADOW = True except Exception: HAS_M4_SHADOW = False def calculate_slope( height_map: np.ndarray, gsd_x: float = 1.0, gsd_y: float = 1.0 ) -> np.ndarray: """ Calculate surface slope angle in degrees [0, 90]. Parameters: - height_map: 2D float32 numpy array (metric elevation in m or relative rDSM) - gsd_x: Ground sampling distance along x-axis (meters/pixel) - gsd_y: Ground sampling distance along y-axis (meters/pixel) """ if height_map.ndim != 2: raise ValueError(f"height_map must be 2D, got shape {height_map.shape}") gsd_x = max(1e-4, float(gsd_x)) gsd_y = max(1e-4, float(gsd_y)) grad_y, grad_x = np.gradient(height_map.astype(np.float64), gsd_y, gsd_x) grad_mag = np.sqrt(grad_x ** 2 + grad_y ** 2) slope_rad = np.arctan(grad_mag) slope_deg = np.degrees(slope_rad).astype(np.float32) return np.clip(slope_deg, 0.0, 90.0) def calculate_hillshade( height_map: np.ndarray, azimuth_deg: float = 315.0, altitude_deg: float = 45.0, gsd_x: float = 1.0, gsd_y: float = 1.0, z_factor: float = 1.0 ) -> np.ndarray: """ Calculate 2D shaded relief (hillshade) array in [0, 255] uint8 using standard Horn's algorithm. """ if height_map.ndim != 2: raise ValueError(f"height_map must be 2D, got shape {height_map.shape}") gsd_x = max(1e-4, float(gsd_x)) gsd_y = max(1e-4, float(gsd_y)) # Convert geographic azimuth to mathematical radian angle azimuth_rad = np.radians(360.0 - azimuth_deg + 90.0) altitude_rad = np.radians(altitude_deg) dy, dx = np.gradient(height_map.astype(np.float64) * float(z_factor), gsd_y, gsd_x) slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) aspect_rad = np.arctan2(-dy, -dx) shaded = np.sin(altitude_rad) * np.cos(slope_rad) + \ np.cos(altitude_rad) * np.sin(slope_rad) * np.cos(azimuth_rad - aspect_rad) shaded = np.clip(shaded, 0.0, 1.0) return (shaded * 255.0).astype(np.uint8) def calculate_confidence_map( depth: np.ndarray, rgb: Optional[np.ndarray] = None, shadow_conf: Optional[np.ndarray] = None, border_margin: int = 15 ) -> np.ndarray: """ Compute heuristic confidence map C(x,y) in [0, 1]. Penalizes image boundaries, dark shadow zones, and high-frequency disparity noise. """ h, w = depth.shape conf = np.full((h, w), 0.90, dtype=np.float32) # 1. Boundary attenuation if border_margin > 0: y_idx, x_idx = np.ogrid[:h, :w] dist_left = x_idx dist_right = w - 1 - x_idx dist_top = y_idx dist_bottom = h - 1 - y_idx min_border_dist = np.minimum(np.minimum(dist_left, dist_right), np.minimum(dist_top, dist_bottom)) border_factor = np.clip(min_border_dist / float(border_margin), 0.35, 1.0).astype(np.float32) conf *= border_factor # 2. Shadow penalty if shadow_conf is not None: shadow_penalty = 1.0 - 0.35 * np.clip(shadow_conf, 0.0, 1.0) conf *= shadow_penalty # 3. High-gradient noise attenuation dy, dx = np.gradient(depth) grad_d = np.sqrt(dx ** 2 + dy ** 2) p95_grad = float(np.percentile(grad_d, 95.0)) if p95_grad > 0: noise_factor = np.clip(1.0 - 0.20 * (grad_d / (p95_grad * 2.0)), 0.5, 1.0).astype(np.float32) conf *= noise_factor return np.clip(conf, 0.0, 1.0).astype(np.float32) def fit_supplied_gcps( depth: np.ndarray, gcps: List[Dict[str, float]], transform: Optional[Affine] = None, a_prior: Optional[float] = None, lambda_prior: float = 0.0 ) -> Tuple[float, float, Dict[str, Any]]: r""" Fit metric linear calibration H = a*D + b strictly using user-supplied GCPs. Validation Rules: - Coordinates must lie within [0, W) and [0, H). - Duplicate pixel locations are detected and rejected. - Elevations must be finite. - At least 2 valid, non-collocated GCPs are required (otherwise ValueError). - Depth variation across GCPs must be non-zero (otherwise degenerate ValueError). Mathematical Formulation: If a_prior is not None and lambda_prior > 0: min_{a, b} \sum (a D_i + b - H_i)^2 + \lambda (a - a_prior)^2 Normal system: [[ \sum D_i^2 + \lambda, \sum D_i ], [ \sum D_i, K ]] [ a, b ]^T = [ \sum D_i H_i + \lambda a_prior, \sum H_i ]^T Else: Standard unconstrained OLS: A = [D, 1], params = (A^T A)^{-1} A^T H Residuals: residual_i = (a D_i + b) - H_i gcp_mae = mean(|residual_i|) gcp_rmse = sqrt(mean(residual_i^2)) """ if depth.ndim != 2: raise ValueError(f"Depth array must be 2D, got shape {depth.shape}") h, w = depth.shape if not gcps or len(gcps) == 0: raise ValueError("Cannot perform GCP calibration with empty GCP list.") seen_pixels = set() gcp_d_list = [] gcp_h_list = [] for idx, pt in enumerate(gcps): if not isinstance(pt, dict): raise ValueError(f"GCP at index {idx} must be a dictionary with 'x', 'y', 'elevation'.") if "x" not in pt or "y" not in pt or "elevation" not in pt: raise ValueError(f"GCP at index {idx} missing required key ('x', 'y', 'elevation').") px_x = float(pt["x"]) px_y = float(pt["y"]) elev = float(pt["elevation"]) if not np.isfinite(elev): raise ValueError(f"GCP at index {idx} has non-finite elevation: {elev}") # Check if coordinates are in geographic/projected units (convert if outside pixel bounds and transform exists) if transform is not None and (px_x >= w or px_y >= h or px_x < 0 or px_y < 0): col_f, row_f = ~transform * (px_x, px_y) r = int(round(row_f)) c_pt = int(round(col_f)) else: r = int(round(px_y)) c_pt = int(round(px_x)) if r < 0 or r >= h or c_pt < 0 or c_pt >= w: raise ValueError(f"GCP coordinate ({px_x}, {px_y}) maps to pixel ({c_pt}, {r}) which is out of raster bounds ({w}x{h}).") pixel_key = (r, c_pt) if pixel_key in seen_pixels: raise ValueError(f"Duplicate GCP pixel location detected at ({c_pt}, {r}). Each GCP must occupy a unique pixel.") seen_pixels.add(pixel_key) d_val = float(depth[r, c_pt]) if not np.isfinite(d_val): raise ValueError(f"Depth value at GCP pixel ({c_pt}, {r}) is non-finite.") gcp_d_list.append(d_val) gcp_h_list.append(elev) k = len(gcp_d_list) if k < 2: raise ValueError(f"GCP calibration requires at least 2 usable GCPs, got {k}.") d_arr = np.array(gcp_d_list, dtype=np.float64) h_arr = np.array(gcp_h_list, dtype=np.float64) # Check for depth variation (avoid degenerate fit) d_range = float(np.max(d_arr) - np.min(d_arr)) if d_range < 1e-7: raise ValueError(f"Degenerate GCP configuration: insufficient depth variation across GCP locations (range={d_range:.2e}).") if a_prior is not None and lambda_prior > 0: # Regularized / Prior OLS sum_d2 = float(np.sum(d_arr ** 2)) sum_d = float(np.sum(d_arr)) sum_dh = float(np.sum(d_arr * h_arr)) sum_h = float(np.sum(h_arr)) M = np.array([ [sum_d2 + float(lambda_prior), sum_d], [sum_d, float(k)] ], dtype=np.float64) v = np.array([ sum_dh + float(lambda_prior) * float(a_prior), sum_h ], dtype=np.float64) params = np.linalg.solve(M, v) scale_a = float(params[0]) offset_b = float(params[1]) else: # Standard unconstrained OLS A = np.column_stack([d_arr, np.ones_like(d_arr)]) params, _, _, _ = np.linalg.lstsq(A, h_arr, rcond=None) scale_a = float(params[0]) offset_b = float(params[1]) # Compute residuals and metrics h_pred = scale_a * d_arr + offset_b residuals = (h_pred - h_arr).tolist() gcp_mae = float(np.mean(np.abs(h_pred - h_arr))) gcp_rmse = float(np.sqrt(np.mean((h_pred - h_arr) ** 2))) gcp_info = { "method": "gcp", "gcp_count": int(k), "dem_source": None, "scale_a": float(scale_a), "offset_b": float(offset_b), "gcp_mae": float(gcp_mae), "gcp_rmse": float(gcp_rmse), "gcp_residuals": [float(r) for r in residuals], "a_prior": float(a_prior) if a_prior is not None else None, "lambda_prior": float(lambda_prior) } return scale_a, offset_b, gcp_info def fit_supplied_dem( depth: np.ndarray, dem_path: str, target_crs: Any, target_transform: Affine, target_shape: Tuple[int, int], terrain_percentile: float = 25.0 ) -> Tuple[float, float, Dict[str, Any]]: """ Align user-supplied lower-resolution DEM (e.g. SRTM) and fit ground elevation anchor. Principles: 1. DEM is reprojected and resampled to the RGB raster grid using bilinear interpolation. 2. Coarse terrain datum is distinguished from high-frequency monocular surface structure. 3. Low-relative-elevation pixels (bottom terrain_percentile %) are identified as terrain candidates. 4. Robust statistics on terrain anchors provide the vertical elevation datum anchor. """ if target_crs is None or target_transform is None: raise ValueError("DEM calibration requires a georeferenced target image with valid CRS and transform.") if not os.path.isfile(dem_path): raise FileNotFoundError(f"Supplied DEM file not found at: {dem_path}") dst_h, dst_w = target_shape reprojected_dem = np.zeros((dst_h, dst_w), dtype=np.float32) with rasterio.open(dem_path) as src_dem: if src_dem.crs is None: raise ValueError(f"Supplied DEM at {dem_path} is missing CRS metadata.") reproject( source=rasterio.band(src_dem, 1), destination=reprojected_dem, src_transform=src_dem.transform, src_crs=src_dem.crs, dst_transform=target_transform, dst_crs=target_crs, resampling=Resampling.bilinear ) # Identify valid spatial overlap valid_mask = np.isfinite(reprojected_dem) & (reprojected_dem > -1000.0) & np.isfinite(depth) valid_count = int(np.sum(valid_mask)) if valid_count < 20: raise ValueError(f"Insufficient spatial overlap between supplied DEM and target raster ({valid_count} pixels).") d_valid = depth[valid_mask].astype(np.float64) dem_valid = reprojected_dem[valid_mask].astype(np.float64) # Identify terrain candidates (lowest relative depth values) p_clamped = float(np.clip(float(terrain_percentile), 5.0, 95.0)) d_thresh = float(np.percentile(d_valid, p_clamped)) terrain_mask = valid_mask & (depth <= d_thresh) terrain_count = int(np.sum(terrain_mask)) if terrain_count < 10: terrain_mask = valid_mask terrain_count = valid_count d_terrain = depth[terrain_mask].astype(np.float64) dem_terrain = reprojected_dem[terrain_mask].astype(np.float64) terrain_anchor_elevation = float(np.median(dem_terrain)) # Fit robust linear relation on terrain anchors d_t_range = float(np.max(d_terrain) - np.min(d_terrain)) if d_t_range > 1e-4: A = np.column_stack([d_terrain, np.ones_like(d_terrain)]) params, _, _, _ = np.linalg.lstsq(A, dem_terrain, rcond=None) scale_a = float(params[0]) offset_b = float(params[1]) else: # Flat terrain fallback: anchor datum offset scale_a = 1.0 offset_b = float(terrain_anchor_elevation - np.median(d_terrain)) info = { "method": "dem", "gcp_count": 0, "dem_source": str(dem_path), "scale_a": float(scale_a), "offset_b": float(offset_b), "valid_overlap_pixels": int(valid_count), "terrain_anchor_count": int(terrain_count), "terrain_anchor_elevation": float(terrain_anchor_elevation), "terrain_percentile": float(p_clamped), "dem_min": float(np.nanmin(dem_valid)), "dem_max": float(np.nanmax(dem_valid)) } return scale_a, offset_b, info def process_image( path: str, gcps: Optional[List[Dict[str, float]]] = None, dem_path: Optional[str] = None, dem_file: Optional[str] = None, use_shadows: bool = True, use_guided_filter: bool = True, model=None, processor=None, device=None, refinement: Optional[Any] = None, external_depth: Optional[np.ndarray] = None, shadow_constraints: Optional[Dict[str, Any]] = None, a_prior: Optional[float] = None, lambda_prior: float = 0.0, terrain_percentile: float = 25.0, calibration_method: str = "linear", refine_roofs: bool = True ) -> Dict[str, Any]: """ Core M2 processing entry point. Parameters: - path: Filepath to RGB image (.tif, .png, .jpg, etc.) - gcps: Optional user-supplied Ground Control Points [{"x": px_x, "y": px_y, "elevation": z_m}, ...] - dem_path: Optional user-supplied lower-resolution DEM filepath (e.g. SRTM GeoTIFF) - dem_file: Optional alias for dem_path - use_shadows: Whether to compute shadow confidence - model, processor, device: Optional pre-loaded model instances - refinement: Optional M1/M3 spatial refinement configuration - external_depth: Optional externally supplied depth map from M1 (bypasses raw inference if present) - shadow_constraints: Optional M4 geometric shadow constraints dictionary - calibration_method: Calibration approach to use. Options: - "linear": Standard OLS linear calibration H = a*D + b (default) - "polynomial_deg2": 2nd-degree polynomial (requires >= 3 GCPs for GCP mode) - "polynomial_deg3": 3rd-degree polynomial (requires >= 4 GCPs for GCP mode) - "polynomial_deg4": 4th-degree polynomial (requires >= 5 GCPs for GCP mode) - "isotonic": Non-parametric monotonic regression (requires >= 5 GCPs for GCP mode) Returns: Problem-statement-compliant result dictionary: { "height_map": np.ndarray, # float32 (H, W) -> metric DSM (m) if calibrated, else rDSM in [0, 1] "width": int, "height": int, "mode": "relative" | "absolute", "calibrated": bool, # True ONLY if genuine GCP/DEM calibration source was supplied "georeferenced": bool, # True if image has valid CRS and affine transform "height_unit": "rel" | "m", # "m" ONLY if calibrated is True, otherwise "rel" "rgb": np.ndarray, # uint8 (H, W, 3) "slope_map": np.ndarray, # float32 (H, W) in degrees [0, 90] "confidence_map": np.ndarray, # float32 (H, W) in [0, 1] "crs": CRS or None, "transform": Affine or None, "metadata": dict } """ if not os.path.isfile(path): raise FileNotFoundError(f"Input image not found: {path}") ext = os.path.splitext(path)[1].lower() is_georeferenced = False crs = None transform = None gsd_x, gsd_y = None, None solar_elevation = None solar_azimuth = None rgb_tags = {} # Support either dem_path or dem_file effective_dem = dem_path if dem_path else dem_file # [1/5] Load Image & Inspect Georeferencing if ext in [".tif", ".tiff"]: try: with rasterio.open(path) as src: raw_crs = src.crs raw_transform = src.transform rgb_tags = src.tags() if raw_crs is not None and raw_transform is not None and not raw_transform.is_identity: is_georeferenced = True crs = raw_crs transform = raw_transform gsd_x = abs(float(transform.a)) gsd_y = abs(float(transform.e)) if "SOLAR_ELEVATION" in rgb_tags: solar_elevation = float(rgb_tags["SOLAR_ELEVATION"]) if "SOLAR_AZIMUTH" in rgb_tags: solar_azimuth = float(rgb_tags["SOLAR_AZIMUTH"]) data = src.read() if data.ndim == 3: rgb = data.transpose(1, 2, 0) if rgb.shape[2] > 3: rgb = rgb[:, :, :3] elif rgb.shape[2] == 1: rgb = np.repeat(rgb, 3, axis=2) elif data.ndim == 2: rgb = np.repeat(data[:, :, np.newaxis], 3, axis=2) except Exception: pil_img = Image.open(path).convert("RGB") rgb = np.array(pil_img, dtype=np.uint8) else: pil_img = Image.open(path).convert("RGB") rgb = np.array(pil_img, dtype=np.uint8) if rgb.dtype != np.uint8: if np.max(rgb) <= 1.0: rgb = (rgb * 255.0).astype(np.uint8) elif np.max(rgb) > 255.0: # Satellite 10-bit / 11-bit / 12-bit / 16-bit dynamic range (e.g. Cartosat, WorldView, Landsat) from data_io.cartosat_adapter import stretch_satellite_dynamic_range rgb = stretch_satellite_dynamic_range(rgb) else: rgb = np.clip(rgb, 0, 255).astype(np.uint8) h, w, c = rgb.shape # [2/5] Monocular Depth Inference / External Depth Hook tiled_inference_used = False if external_depth is not None and isinstance(external_depth, np.ndarray): depth = external_depth.astype(np.float32) if depth.shape != (h, w): depth_img = Image.fromarray(depth).resize((w, h), Image.Resampling.BILINEAR) depth = np.array(depth_img, dtype=np.float32) else: # Use tiled inference for high-res inputs (>768px in either dimension) if max(h, w) > 768: try: depth, _ = estimate_depth_tiled(rgb, tile_size=512, overlap=0.25, device=device) tiled_inference_used = True except Exception: if model is None: model, processor, device = load_model() depth = estimate_depth(rgb, model=model, processor=processor, device=device) else: if model is None: model, processor, device = load_model() depth = estimate_depth(rgb, model=model, processor=processor, device=device) # Robust percentile-based dynamic range normalization via calibration.rdsm try: rdsm = make_rdsm(depth, p_low=0.5, p_high=99.5) d_valid = depth[np.isfinite(depth)] d_min = float(np.percentile(d_valid, 0.5)) if len(d_valid) > 0 else 0.0 d_max = float(np.percentile(d_valid, 99.5)) if len(d_valid) > 0 else 1.0 d_span = (d_max - d_min) if (d_max - d_min) > 1e-6 else 1.0 except Exception: d_valid = depth[np.isfinite(depth)] if len(d_valid) > 0: d_p01 = float(np.percentile(d_valid, 0.5)) d_p99 = float(np.percentile(d_valid, 99.5)) d_span = (d_p99 - d_p01) if (d_p99 - d_p01) > 1e-6 else 1.0 rdsm = np.clip((depth - d_p01) / d_span, 0.0, 1.0).astype(np.float32) d_min = d_p01 d_max = d_p99 else: d_min, d_max, d_span = 0.0, 1.0, 1.0 rdsm = np.zeros_like(depth, dtype=np.float32) # M7 Guided Filter edge-preserving refinement guided_filter_applied = False if use_guided_filter and HAS_GUIDED_FILTER: try: rdsm_refined = refine_depth_anything_map(guide_image=rgb, raw_depth=rdsm, radius=16, eps=0.01) rdsm = np.clip(rdsm_refined, 0.0, 1.0).astype(np.float32) guided_filter_applied = True except Exception as gf_err: print(f"[M7 Guided Filter Warning] {gf_err}") # [3/5] Calibration: Relative vs Absolute Mode scale_a = None offset_b = None calibrated = False calibration_method = "none" calibration_info = { "method": "none", "gcp_count": 0, "dem_source": None, "scale_a": None, "offset_b": None } # Normalize calibration method if not calibration_method or calibration_method.lower() in ["none", "linear", "ols", "default"]: calibration_method = "linear" if is_georeferenced: # Path A: GCP Calibration if gcps is not None and len(gcps) > 0: # Determine minimum GCPs needed for non-linear method min_gcps_for_method = {"linear": 2, "polynomial_deg2": 3, "polynomial_deg3": 4, "polynomial_deg4": 5, "isotonic": 5} required_gcps = min_gcps_for_method.get(calibration_method, 2) use_nonlinear = calibration_method != "linear" and len(gcps) >= required_gcps if use_nonlinear: try: gcp_d_list = [] gcp_h_list = [] seen_pixels = set() for idx, pt in enumerate(gcps): px_x = float(pt["x"]) px_y = float(pt["y"]) elev = float(pt["elevation"]) if transform is not None and (px_x >= w or px_y >= h or px_x < 0 or px_y < 0): col_f, row_f = ~transform * (px_x, px_y) r = int(round(row_f)) c_pt = int(round(col_f)) else: r = int(round(px_y)) c_pt = int(round(px_x)) if r < 0 or r >= h or c_pt < 0 or c_pt >= w: continue pixel_key = (r, c_pt) if pixel_key in seen_pixels: continue seen_pixels.add(pixel_key) d_val = float(depth[r, c_pt]) if not np.isfinite(d_val): continue gcp_d_list.append(d_val) gcp_h_list.append(elev) gcp_d_arr = np.array(gcp_d_list, dtype=np.float64) gcp_h_arr = np.array(gcp_h_list, dtype=np.float64) if len(gcp_d_arr) >= required_gcps: model, nl_info = fit_gcp_calibration_nonlinear(gcp_d_arr, gcp_h_arr, method=calibration_method) height_map = apply_gcp_calibration_nonlinear(depth, model, method=calibration_method) calibrated = True mode = "absolute" height_unit = "m" calibration_method_used = nl_info["method"] calibration_info = nl_info calibration_info["a_prior"] = float(a_prior) if a_prior is not None else None calibration_info["lambda_prior"] = float(lambda_prior) else: raise ValueError(f"Insufficient valid GCPs ({len(gcp_d_arr)}) for {calibration_method}") except Exception as nl_err: print(f"[Non-linear GCP Warning] {nl_err}. Falling back to linear...") use_nonlinear = False if not use_nonlinear: try: scale_a, offset_b, gcp_info = fit_supplied_gcps( depth, gcps, transform=transform, a_prior=a_prior, lambda_prior=lambda_prior ) height_map = (scale_a * depth + offset_b).astype(np.float32) calibrated = True mode = "absolute" height_unit = "m" calibration_method_used = "gcp_linear" calibration_info = gcp_info except Exception as gcp_err: print(f"[GCP Calibration Warning] {gcp_err}. Falling back to DEM or relative mode...") if effective_dem is not None and os.path.isfile(effective_dem): scale_a, offset_b, dem_info = fit_supplied_dem( depth, effective_dem, target_crs=crs, target_transform=transform, target_shape=(h, w), terrain_percentile=terrain_percentile ) height_map = (scale_a * depth + offset_b).astype(np.float32) calibrated = True mode = "absolute" height_unit = "m" calibration_method_used = "dem" calibration_info = dem_info else: mode = "relative" calibrated = False height_unit = "rel" calibration_method_used = "none" calibration_method = "none" height_map = rdsm calibration_info = {"method": "none", "error": str(gcp_err), "gcp_count": len(gcps)} # Path B: User-supplied DEM Calibration elif effective_dem is not None and os.path.isfile(effective_dem): scale_a, offset_b, dem_info = fit_supplied_dem( depth, effective_dem, target_crs=crs, target_transform=transform, target_shape=(h, w), terrain_percentile=terrain_percentile ) # If non-linear method requested, refine DEM calibration with non-linear mapping if calibration_method != "linear": try: reprojected_dem = np.zeros((h, w), dtype=np.float32) with rasterio.open(effective_dem) as src_dem: if src_dem.crs is not None: from rasterio.warp import reproject as _reproject, Resampling as _Resampling _reproject( source=rasterio.band(src_dem, 1), destination=reprojected_dem, src_transform=src_dem.transform, src_crs=src_dem.crs, dst_transform=transform, dst_crs=crs, resampling=_Resampling.bilinear ) valid_mask = np.isfinite(reprojected_dem) & (reprojected_dem > -1000.0) & np.isfinite(depth) if np.sum(valid_mask) >= 50: nl_height, nl_method, nl_info = calibrate_depth_to_dsm_nonlinear( depth[valid_mask], reprojected_dem[valid_mask], method=calibration_method ) height_map = np.full((h, w), np.nan, dtype=np.float32) height_map[valid_mask] = nl_height # Fill NaN with linear fallback nan_mask = ~np.isfinite(height_map) if np.any(nan_mask): height_map[nan_mask] = (scale_a * depth[nan_mask] + offset_b).astype(np.float32) calibrated = True mode = "absolute" height_unit = "m" calibration_method_used = nl_info["method"] calibration_info = nl_info calibration_info["linear_fallback_applied"] = bool(np.any(nan_mask)) else: raise ValueError("Insufficient valid overlap for non-linear DEM calibration") except Exception as nl_err: print(f"[Non-linear DEM Warning] {nl_err}. Using linear DEM calibration.") height_map = (scale_a * depth + offset_b).astype(np.float32) calibrated = True mode = "absolute" height_unit = "m" calibration_method_used = "dem_linear" calibration_info = dem_info else: height_map = (scale_a * depth + offset_b).astype(np.float32) calibrated = True mode = "absolute" height_unit = "m" calibration_method_used = "dem_linear" calibration_info = dem_info # Path C: Single-View GeoTIFF without GCPs or DEM -> Pure relative rDSM [0, 1] else: mode = "relative" calibrated = False height_unit = "rel" calibration_method_used = "none" calibration_method = "none" height_map = rdsm calibration_info = { "method": "none", "gcp_count": 0, "dem_source": None, "scale_a": None, "offset_b": None } else: # Non-georeferenced JPG/PNG -> Pure relative rDSM [0, 1] mode = "relative" calibrated = False height_unit = "rel" calibration_method_used = "none" calibration_method = "none" height_map = rdsm calibration_info = { "method": "none", "gcp_count": 0, "dem_source": None, "scale_a": None, "offset_b": None } crs = None transform = None gsd_x, gsd_y = None, None # [4/5] Shadow Analysis & Shadow Constraints Hook (M4 Integration) shadow_conf = None m4_shadow_summary = None if use_shadows: if HAS_M4_SHADOW: temp_shadow_png = None try: # Ensure OpenCV can always read the image by saving a temporary standard 8-bit PNG with tempfile.NamedTemporaryFile(delete=False, suffix=".png") as tmp_png: Image.fromarray(rgb).save(tmp_png.name) temp_shadow_png = tmp_png.name m4_shadow_summary = run_m4_shadow_pipeline( image_path=str(temp_shadow_png), meters_per_pixel=gsd_x if is_georeferenced else None, sun_elevation_deg=solar_elevation, is_test_mode=False, generate_diagnostics=False ) except Exception as e: print(f"[M4 Shadow Warning] {e}") finally: if temp_shadow_png and os.path.isfile(temp_shadow_png): try: os.remove(temp_shadow_png) except Exception: pass shadow_conf = compute_shadow_confidence(rgb) shadow_mode = "metric" if solar_elevation is not None else "structural" else: shadow_mode = "disabled" # [4.5/5] sat2Map Roof Planar Regularization & Facet Extraction slope_gsd_x = gsd_x if (gsd_x is not None and calibrated) else 1.0 slope_gsd_y = gsd_y if (gsd_y is not None and calibrated) else 1.0 if refine_roofs: try: height_map = regularize_roof_planes(height_map, rgb=rgb) if not calibrated: h_m = float(np.nanmin(height_map)) h_M = float(np.nanmax(height_map)) if (h_M - h_m) > 1e-6: height_map = np.clip((height_map - h_m) / (h_M - h_m), 0.0, 1.0).astype(np.float32) except Exception as ref_err: print(f"[sat2Map Refiner Warning] {ref_err}") planar_info = extract_planar_segments(height_map, gsd=slope_gsd_x) # [5/5] Slope, Hillshade & Confidence Calculations slope_map = calculate_slope(height_map, gsd_x=slope_gsd_x, gsd_y=slope_gsd_y) confidence_map = calculate_confidence_map(depth, rgb=rgb, shadow_conf=shadow_conf, border_margin=15) # Calculate shaded relief (hillshade) for photorealistic 3D visualization hillshade_map = calculate_hillshade( height_map, azimuth_deg=315.0 if solar_azimuth is None else solar_azimuth, altitude_deg=45.0 if solar_elevation is None else max(15.0, min(75.0, solar_elevation)), gsd_x=slope_gsd_x, gsd_y=slope_gsd_y, z_factor=1.0 if calibrated else 2.5 ) # Error Map: If reference DEM supplied, compute residual error; else compute empirical uncertainty map if 'reprojected_dem' in locals() and reprojected_dem is not None: valid_dem = np.isfinite(reprojected_dem) & (reprojected_dem > -1000.0) error_map = np.where(valid_dem, height_map - reprojected_dem, np.nan).astype(np.float32) error_type = "ground_truth_residual" else: h_span = float(np.nanmax(height_map) - np.nanmin(height_map)) if np.any(np.isfinite(height_map)) else 1.0 scale_factor = (h_span * 0.12) if calibrated else 1.0 error_map = ((1.0 - confidence_map) * scale_factor).astype(np.float32) error_type = "estimated_uncertainty" metadata = { "input_path": str(path), "input_format": ext.lstrip("."), "mode": mode, "calibrated": calibrated, "georeferenced": bool(is_georeferenced), "height_unit": height_unit, "model": "Depth Anything V2 + sat2Map 3D Refiner", "calibration": calibration_info, "calibration_method_requested": str(calibration_method), "calibration_method_used": str(calibration_method_used), "shadow_mode": shadow_mode, "shadow_constraints": shadow_constraints if shadow_constraints is not None else m4_shadow_summary, "m4_active": HAS_M4_SHADOW, "width": int(w), "height": int(h), "crs": str(crs) if crs is not None else None, "gsd_x": float(gsd_x) if gsd_x is not None else None, "gsd_y": float(gsd_y) if gsd_y is not None else None, "guided_filter_applied": bool(guided_filter_applied), "tiled_inference_used": bool(tiled_inference_used), "refine_roofs_applied": bool(refine_roofs), "planar_facets": { "num_planes": int(planar_info["num_planes"]), "plane_tilts": planar_info["plane_tilts"], "mean_tilt_deg": float(planar_info["mean_tilt_deg"]) }, "error_type": error_type, "raw_depth_stats": { "min": float(d_min), "max": float(d_max), "span": float(d_span) } } return { "height_map": height_map.astype(np.float32), "depth_map": rdsm.astype(np.float32), # normalized [0,1] monocular depth "error_map": error_map.astype(np.float32), "error_type": error_type, "hillshade": hillshade_map.astype(np.uint8), "width": int(w), "height": int(h), "mode": mode, "calibrated": calibrated, "georeferenced": bool(is_georeferenced), "height_unit": height_unit, "rgb": rgb.astype(np.uint8), "slope_map": slope_map.astype(np.float32), "confidence_map": confidence_map.astype(np.float32), "planar_facets": planar_info, "crs": crs, "transform": transform, "metadata": metadata } def export_dsm(result: Dict[str, Any], output_path: str) -> str: """ Export DSM raster to disk. - If georeferenced GeoTIFF: writes float32 GeoTIFF with CRS and Affine transform. - If relative image: writes GeoTIFF or standard PNG. """ os.makedirs(os.path.dirname(os.path.abspath(output_path)), exist_ok=True) height_map = result["height_map"] h, w = height_map.shape crs = result.get("crs") transform = result.get("transform") ext = os.path.splitext(output_path)[1].lower() if crs is not None and transform is not None: with rasterio.open( output_path, "w", driver="GTiff", height=h, width=w, count=1, dtype=np.float32, crs=crs, transform=transform, nodata=-9999.0 ) as dst: dst.write(height_map.astype(np.float32), 1) else: if ext in [".tif", ".tiff"]: with rasterio.open( output_path, "w", driver="GTiff", height=h, width=w, count=1, dtype=np.float32, nodata=-9999.0 ) as dst: dst.write(height_map.astype(np.float32), 1) elif ext == ".npy": np.save(output_path, height_map) else: norm = (np.clip(height_map, 0.0, 1.0) * 255.0).astype(np.uint8) Image.fromarray(norm).save(output_path) return output_path def export_slope(result: Dict[str, Any], output_path: str) -> str: """ Export slope raster (in degrees) to disk. """ os.makedirs(os.path.dirname(os.path.abspath(output_path)), exist_ok=True) slope_map = result["slope_map"] h, w = slope_map.shape crs = result.get("crs") transform = result.get("transform") ext = os.path.splitext(output_path)[1].lower() if crs is not None and transform is not None: with rasterio.open( output_path, "w", driver="GTiff", height=h, width=w, count=1, dtype=np.float32, crs=crs, transform=transform, nodata=-9999.0 ) as dst: dst.write(slope_map.astype(np.float32), 1) else: if ext in [".tif", ".tiff"]: with rasterio.open( output_path, "w", driver="GTiff", height=h, width=w, count=1, dtype=np.float32, nodata=-9999.0 ) as dst: dst.write(slope_map.astype(np.float32), 1) elif ext == ".npy": np.save(output_path, slope_map) else: norm = (np.clip(slope_map / 90.0, 0.0, 1.0) * 255.0).astype(np.uint8) Image.fromarray(norm).save(output_path) return output_path def export_point_cloud(result: Dict[str, Any], output_path: str, step: int = 1) -> str: """ Export 3D calibrated point cloud in Stanford Binary PLY format. """ height_map = result["height_map"] rgb = result.get("rgb") gsd_x = result.get("metadata", {}).get("gsd_x") or 1.0 gsd_y = result.get("metadata", {}).get("gsd_y") or 1.0 return export_point_cloud_ply(height_map, rgb, output_path, gsd_x=gsd_x, gsd_y=gsd_y, step=step) def export_mesh(result: Dict[str, Any], output_path: str, step: int = 2) -> Tuple[str, str, str]: """ Export 3D textured mesh in Wavefront OBJ format. """ height_map = result["height_map"] rgb = result.get("rgb") gsd_x = result.get("metadata", {}).get("gsd_x") or 1.0 gsd_y = result.get("metadata", {}).get("gsd_y") or 1.0 return export_mesh_obj(height_map, rgb, output_path, gsd_x=gsd_x, gsd_y=gsd_y, step=step)