Spaces:
Running on Zero
Running on Zero
| """ | |
| 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) | |