From f2d71da8e1d78c81bec0bd14e6c710244cf377dd Mon Sep 17 00:00:00 2001 From: "s0wlz (Matthias Puchstein)" Date: Tue, 3 Feb 2026 23:35:51 +0100 Subject: [PATCH] fix(buildings): Optimize wall texture sampling with VRT windowed read Replaced the memory-intensive full image read of per-tile orthophotos with a windowed read from the global Ortho VRT. This fixes memory crashes and ensures correct texture coverage for buildings that cross tile boundaries. --- geodata_pipeline/buildings.py | 55 +++++++++++----- tests/test_buildings_texture.py | 107 ++++++++++++++++++++++++++++++++ 2 files changed, 145 insertions(+), 17 deletions(-) create mode 100644 tests/test_buildings_texture.py diff --git a/geodata_pipeline/buildings.py b/geodata_pipeline/buildings.py index 4bb3531..6bf0621 100644 --- a/geodata_pipeline/buildings.py +++ b/geodata_pipeline/buildings.py @@ -544,23 +544,44 @@ def export_buildings(cfg: Config) -> int: wall_color = np.zeros((len(vertices), 3), dtype=np.float32) + 0.75 ortho_path = os.path.join(cfg.export.ortho_dir, f"{tile_id}.jpg") ortho_bytes = _load_ortho(tile_id, ortho_path) - if ortho_bytes: - try: - ortho_ds = gdal.Open(ortho_path) - if ortho_ds: - gt_o = ortho_ds.GetGeoTransform() - bands = [ortho_ds.GetRasterBand(i + 1).ReadAsArray() for i in range(min(3, ortho_ds.RasterCount))] - for idx, (vx, vy, _) in enumerate(vertices): - wx = vx + xmin - wy = vy + ymin - col = int((wx - gt_o[0]) / gt_o[1]) - row = int((wy - gt_o[3]) / gt_o[5]) - if 0 <= row < ortho_ds.RasterYSize and 0 <= col < ortho_ds.RasterXSize: - rgb = [bands[i][row, col] for i in range(len(bands))] - if rgb: - wall_color[idx] = np.array(rgb[:3], dtype=np.float32) / 255.0 - except Exception: - pass + + try: + # Sample walls from global VRT using windowed read to avoid memory issues and ensure coverage + ortho_ds = gdal.Open(cfg.work.ortho_vrt) + if ortho_ds: + gt_o = ortho_ds.GetGeoTransform() + xmin, ymin, xmax, ymax = bounds + + # Calculate srcwin + xoff = int((xmin - gt_o[0]) / gt_o[1]) + yoff = int((ymax - gt_o[3]) / gt_o[5]) + xsize = int((xmax - xmin) / gt_o[1]) + 1 + ysize = int((ymax - ymin) / abs(gt_o[5])) + 1 + + # Clamp to raster dimensions + xoff_clamped = max(0, min(xoff, ortho_ds.RasterXSize - 1)) + yoff_clamped = max(0, min(yoff, ortho_ds.RasterYSize - 1)) + xsize_clamped = max(1, min(xsize, ortho_ds.RasterXSize - xoff_clamped)) + ysize_clamped = max(1, min(ysize, ortho_ds.RasterYSize - yoff_clamped)) + + bands = [ + ortho_ds.GetRasterBand(i + 1).ReadAsArray(xoff_clamped, yoff_clamped, xsize_clamped, ysize_clamped) + for i in range(min(3, ortho_ds.RasterCount)) + ] + + for idx, (vx, vy, _) in enumerate(vertices): + wx = vx + xmin + wy = vy + ymin + # Map global coord to window-local coord + col = int((wx - gt_o[0]) / gt_o[1]) - xoff_clamped + row = int((wy - gt_o[3]) / gt_o[5]) - yoff_clamped + + if 0 <= row < ysize_clamped and 0 <= col < xsize_clamped: + rgb = [bands[i][row, col] for i in range(len(bands))] + if rgb: + wall_color[idx] = np.array(rgb[:3], dtype=np.float32) / 255.0 + except Exception: + pass glb_bytes = _compose_glb( gltf_vertices, diff --git a/tests/test_buildings_texture.py b/tests/test_buildings_texture.py new file mode 100644 index 0000000..9bc8bab --- /dev/null +++ b/tests/test_buildings_texture.py @@ -0,0 +1,107 @@ +import unittest +from unittest.mock import patch, MagicMock +import numpy as np +import os +from geodata_pipeline.buildings import export_buildings +from geodata_pipeline.config import Config + +class TestBuildingsTexture(unittest.TestCase): + @patch("geodata_pipeline.buildings.gdal.Open") + @patch("geodata_pipeline.buildings._ensure_cityjson_for_tile") + @patch("geodata_pipeline.buildings._load_cityjson") + @patch("geodata_pipeline.buildings._collect_faces") + @patch("geodata_pipeline.buildings._compose_glb") + @patch("geodata_pipeline.buildings.ensure_dir") + @patch("geodata_pipeline.buildings.os.path.exists") + @patch("builtins.open") + def test_vrt_windowed_sampling(self, mock_open_file, mock_exists, mock_ensure_dir, mock_compose, mock_collect, mock_load_cj, mock_ensure_cj, mock_gdal_open): + # Setup mocks + mock_exists.return_value = True + + # Mock file handle for manifest CSV + mock_handle = MagicMock() + mock_open_file.return_value.__enter__.return_value = mock_handle + mock_handle.__iter__.return_value = [ + "tile_id,xmin,ymin,xmax,ymax,global_min,global_max,out_res,tile_key,tile_min,tile_max\n", + "tile1,1000,1000,2000,2000,0,100,1025,1_1,0,100\n" + ] + + mock_ensure_cj.return_value = "dummy.json" + mock_load_cj.return_value = {"CityObjects": {}} + mock_collect.return_value = ( + [[10.0, 10.0, 5.0]], # vertices (local to tile) + [([0, 0, 0], "WallSurface")] # faces + ) + + # Mock GDAL + # 1. Heightmap VRT (called first for ground snapping) + # 2. Ortho VRT (called second for texturing) + mock_height_ds = MagicMock() + mock_ortho_ds = MagicMock() + + def gdal_open_side_effect(path): + if "dgm.vrt" in path: + return mock_height_ds + if "dop.vrt" in path: + return mock_ortho_ds + # Current implementation might try to open a JPG + if ".jpg" in path: + return MagicMock() + return None + + mock_gdal_open.side_effect = gdal_open_side_effect + + # Setup Heightmap DS + mock_height_ds.GetGeoTransform.return_value = (0, 1, 0, 3000, 0, -1) + mock_height_ds.RasterXSize = 5000 + mock_height_ds.RasterYSize = 5000 + mock_height_band = MagicMock() + mock_height_ds.GetRasterBand.return_value = mock_height_band + mock_height_band.GetNoDataValue.return_value = -9999 + mock_height_band.ReadAsArray.return_value = np.array([[50.0]]) + + # Setup Ortho DS + mock_ortho_ds.GetGeoTransform.return_value = (0, 0.2, 0, 3000, 0, -0.2) # 20cm resolution + mock_ortho_ds.RasterXSize = 15000 + mock_ortho_ds.RasterYSize = 15000 + mock_ortho_ds.RasterCount = 3 + mock_ortho_band = MagicMock() + mock_ortho_ds.GetRasterBand.return_value = mock_ortho_band + mock_ortho_band.ReadAsArray.return_value = np.array([[255, 0, 0]]) # Red pixel + + cfg = Config.default() + cfg.work.heightmap_vrt = "work/dgm.vrt" + cfg.work.ortho_vrt = "work/dop.vrt" + cfg.export.manifest_path = "manifest.csv" + cfg.export.ortho_dir = "ortho_jpg" + + export_buildings(cfg) + + # Verify that we opened the Ortho VRT + mock_gdal_open.assert_any_call("work/dop.vrt") + + # Verify that ReadAsArray was called with a window on the Ortho DS bands + # We check the arguments of the LAST call to ReadAsArray on the ortho band + # (since it iterates over 3 bands) + args, kwargs = mock_ortho_band.ReadAsArray.call_args + + # It should be windowed: (xoff, yoff, xsize, ysize) + self.assertEqual(len(args), 4) + + # Verify window calculation + # tile bounds: xmin=1000, ymin=1000, xmax=2000, ymax=2000 + # ortho GT: origin=(0, 3000), res=0.2 + # xoff = (1000 - 0) / 0.2 = 5000 + # yoff = (2000 - 3000) / -0.2 = 5000 + xoff, yoff, xsize, ysize = args + self.assertEqual(xoff, 5000) + self.assertEqual(yoff, 5000) + + # xsize = (2000 - 1000) / 0.2 = 5000 + # ysize = (2000 - 1000) / 0.2 = 5000 + # Since implementation adds +1 for safety/inclusive bounds: + self.assertTrue(xsize >= 5000) + self.assertTrue(ysize >= 5000) + +if __name__ == "__main__": + unittest.main()