LandscapeUtility
1# coding: utf-8 2# Date: 2021/11/27 3# Author: Jianbo QI 4import os 5import sys 6import trimesh 7import numpy as np 8from LSBoundingBox import LSBoundingBox 9from PostProcessing import RasterHelper 10import alphashape 11 12 13class LandscapeUtility: 14 def __init__(self, land_scape): 15 self.land_scape = land_scape 16 self.width_extent = self.land_scape.get_terrain().extent_width 17 self.height_extent = self.land_scape.get_terrain().extent_height 18 self.tot_area = 0 19 20 # 0 1 2 3 for left lower right upper 21 self.__slice_planes = {0: [[1, 0, 0], [0, 0, 0]], 22 1: [[0, -1, 0], [0, self.height_extent, 0]], 23 2: [[-1, 0, 0], [self.width_extent, 0, 0]], 24 3: [[0, 1, 0], [0, 0, 0]]} 25 26 def __get_rotated_scaled_meshes(self, ignore_list=None): 27 if ignore_list is None: 28 ignore_list = [] 29 objects = self.land_scape.get_objects() 30 instances = self.land_scape.get_instances() 31 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 32 obj_areas = {} 33 obj_bounds = {} 34 obj_meshes = {} 35 for scene_obj in objects: 36 scene_obj_area = 0 37 obj_bounding_box = LSBoundingBox() 38 meshes = [] 39 for obj_component in objects[scene_obj]: 40 if obj_component not in ignore_list: 41 obj_component_path = os.path.join(parameter_dir, obj_component) 42 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 43 scene_obj_area += obj_comp_mesh.area 44 bbox = obj_comp_mesh.bounding_box 45 comp_bounding_box = LSBoundingBox() 46 comp_bounding_box.init_from_bounds(bbox.bounds) 47 obj_bounding_box.add_child(comp_bounding_box) 48 meshes.append(obj_comp_mesh) 49 obj_areas[scene_obj] = scene_obj_area 50 obj_bounds[scene_obj] = obj_bounding_box 51 obj_meshes[scene_obj] = meshes 52 53 instance_dict = {} 54 for instance in instances: 55 instance_dict[instance] = {} 56 pos_index = 0 57 for pos in instances[instance]: 58 instance_dict[instance][pos_index] = {} 59 obj_bound = obj_bounds[instance] 60 obj_mesh = obj_meshes[instance] 61 obj_area = obj_areas[instance] 62 if len(pos) == 10: # scale 63 x_scale_size = pos[7] 64 y_scale_size = pos[9] # vertical direction 65 z_scale_size = pos[8] 66 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 67 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 68 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 69 obj_mesh = [] 70 obj_bound = LSBoundingBox() 71 obj_area = 0 72 for mesh in obj_meshes[instance]: 73 new_mesh = mesh.copy() 74 new_mesh.vertices[:, 0] *= x_scale_factor 75 new_mesh.vertices[:, 1] *= y_scale_factor 76 new_mesh.vertices[:, 2] *= z_scale_factor 77 obj_mesh.append(new_mesh) 78 obj_area += new_mesh.area 79 comp_bounding_box = LSBoundingBox() 80 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 81 obj_bound.add_child(comp_bounding_box) 82 83 if pos[3] != 0: 84 obj_bound = LSBoundingBox() 85 obj_mesh1 = [] 86 obj_area = 0 87 for mesh in obj_mesh: 88 new_mesh = mesh.copy() 89 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 90 axis = [0, 1, 0] 91 if len(pos) >= 7: 92 axis = [-pos[4], pos[6], -pos[5]] 93 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 94 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 95 obj_mesh1.append(rotated_mesh) 96 obj_area += rotated_mesh.area 97 comp_bounding_box = LSBoundingBox() 98 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 99 obj_bound.add_child(comp_bounding_box) 100 obj_mesh = obj_mesh1 101 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 102 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 103 instance_dict[instance][pos_index]["obj_area"] = obj_area 104 pos_index += 1 105 return instance_dict 106 107 def compute_crown_radius(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5): 108 instance_dict = self.__get_rotated_scaled_meshes() 109 instances = self.land_scape.get_instances() 110 f_out = open(out_file, "w") 111 f_out.write("OBJ_Name X Y Z Radius_X Radius_Y\n") 112 for obj_name in instance_dict: 113 positions = instances[obj_name] 114 for pos_index in instance_dict[obj_name]: 115 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 116 pos = positions[pos_index] 117 radius_x = obj_bound.get_x_extents()*0.5 118 radius_z = obj_bound.get_z_extents()*0.5 119 f_out.write(f"{obj_name} %.4f %.4f %.4f %.4f %.4f\n" % (pos[0], pos[1], pos[2], radius_x, radius_z)) 120 f_out.close() 121 122 # write image 123 if out_mask_img_file is not None and out_mask_img_file != "": 124 from PIL import Image, ImageDraw 125 image_width = int(np.ceil(self.width_extent / mask_img_resolution)) 126 image_height = int(np.ceil(self.height_extent / mask_img_resolution)) 127 img = Image.new('L', (image_width, image_height), 0) 128 obj_index = 1 129 for obj_name in instance_dict: 130 positions = instances[obj_name] 131 for pos_index in instance_dict[obj_name]: 132 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 133 pos = positions[pos_index] 134 shape = [(int((pos[0]-obj_bound.maxX)/mask_img_resolution), int((pos[1]-obj_bound.maxZ)/mask_img_resolution)), 135 (int((pos[0]-obj_bound.minX)/mask_img_resolution), int((pos[1]-obj_bound.minZ)/mask_img_resolution))] 136 ImageDraw.Draw(img).ellipse(shape, outline=obj_index, fill=obj_index) 137 obj_index += 1 138 mask = np.array(img).astype(int) 139 out_path = out_mask_img_file 140 if out_path.endswith(".tif"): 141 out_path = out_path[0:-4] 142 RasterHelper.saveToHdr_no_transform(mask, out_path, [], "GTiff") 143 144 def __compute_alphashape_2d(self, point_2d): 145 if len(point_2d) > 5000: 146 indexes = np.random.randint(0, point_2d.shape[0], size=5000) 147 point_2d = point_2d[indexes] 148 radius = 5 149 alpha_shape = alphashape.alphashape(point_2d, radius) 150 while 1: 151 try: 152 coords = np.array(alpha_shape.boundary.coords) 153 except NotImplementedError: 154 radius += -1 155 alpha_shape = alphashape.alphashape(point_2d, radius) 156 else: 157 break 158 return coords 159 160 def compute_crown_boundary(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5): 161 objects = self.land_scape.get_objects() 162 instances = self.land_scape.get_instances() 163 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 164 obj_areas = {} 165 obj_bounds = {} 166 obj_meshes = {} 167 obj_boundary2ds = {} 168 for scene_obj in objects: 169 scene_obj_area = 0 170 obj_bounding_box = LSBoundingBox() 171 meshes = [] 172 obj_points_2d = None 173 for obj_component in objects[scene_obj]: 174 obj_component_path = os.path.join(parameter_dir, obj_component) 175 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 176 scene_obj_area += obj_comp_mesh.area 177 bbox = obj_comp_mesh.bounding_box 178 comp_bounding_box = LSBoundingBox() 179 comp_bounding_box.init_from_bounds(bbox.bounds) 180 obj_bounding_box.add_child(comp_bounding_box) 181 meshes.append(obj_comp_mesh) 182 if obj_points_2d is None: 183 obj_points_2d = obj_comp_mesh.triangles_center[:, [0,2]] 184 else: 185 obj_points_2d = trimesh.util.vstack_empty((obj_points_2d, obj_comp_mesh.triangles_center[:, [0,2]])) 186 # catch exception if alphashape has muilt-polygons 187 obj_boundary2ds[scene_obj] = self.__compute_alphashape_2d(obj_points_2d) 188 obj_areas[scene_obj] = scene_obj_area 189 obj_bounds[scene_obj] = obj_bounding_box 190 obj_meshes[scene_obj] = meshes 191 192 instance_dict = {} 193 for instance in instances: 194 instance_dict[instance] = {} 195 pos_index = 0 196 for pos in instances[instance]: 197 instance_dict[instance][pos_index] = {} 198 obj_bound = obj_bounds[instance] 199 obj_mesh = obj_meshes[instance] 200 obj_area = obj_areas[instance] 201 obj_boundary2d = obj_boundary2ds[instance] 202 if len(pos) == 10: # scale 203 x_scale_size = pos[7] 204 y_scale_size = pos[9] # vertical direction 205 z_scale_size = pos[8] 206 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 207 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 208 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 209 obj_mesh = [] 210 obj_bound = LSBoundingBox() 211 obj_area = 0 212 obj_points_2d = None 213 for mesh in obj_meshes[instance]: 214 new_mesh = mesh.copy() 215 new_mesh.vertices[:, 0] *= x_scale_factor 216 new_mesh.vertices[:, 1] *= y_scale_factor 217 new_mesh.vertices[:, 2] *= z_scale_factor 218 obj_mesh.append(new_mesh) 219 obj_area += new_mesh.area 220 comp_bounding_box = LSBoundingBox() 221 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 222 obj_bound.add_child(comp_bounding_box) 223 if obj_points_2d is None: 224 obj_points_2d = new_mesh.triangles_center[:, [0,2]] 225 else: 226 obj_points_2d = trimesh.util.vstack_empty( 227 (obj_points_2d, new_mesh.triangles_center[:, [0,2]])) 228 obj_boundary2d = self.__compute_alphashape_2d(obj_points_2d) 229 230 if pos[3] != 0: 231 obj_bound = LSBoundingBox() 232 obj_mesh1 = [] 233 obj_area = 0 234 obj_points_2d = None 235 for mesh in obj_mesh: 236 new_mesh = mesh.copy() 237 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 238 axis = [0, 1, 0] 239 if len(pos) >= 7: 240 axis = [-pos[4], pos[6], -pos[5]] 241 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 242 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 243 obj_mesh1.append(rotated_mesh) 244 obj_area += rotated_mesh.area 245 comp_bounding_box = LSBoundingBox() 246 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 247 obj_bound.add_child(comp_bounding_box) 248 if obj_points_2d is None: 249 obj_points_2d = new_mesh.triangles_center[:, [0,2]] 250 else: 251 obj_points_2d = trimesh.util.vstack_empty( 252 (obj_points_2d, new_mesh.triangles_center[:, [0,2]])) 253 obj_mesh = obj_mesh1 254 obj_boundary2d = self.__compute_alphashape_2d(obj_points_2d) 255 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 256 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 257 instance_dict[instance][pos_index]["obj_area"] = obj_area 258 instance_dict[instance][pos_index]["obj_boundary2d"] = obj_boundary2d 259 pos_index += 1 260 # output 261 f_out = open(out_file, "w") 262 f_out.write("OBJ_Name X Y Z Boundary_Point1_X Boundary_Point1_Y, Boundary_Point2_X, Boundary_Point2_Y,...\n") 263 polygons = dict() 264 for obj_name in instance_dict: 265 positions = instances[obj_name] 266 polygons[obj_name] = [] 267 for pos_index in instance_dict[obj_name]: 268 pos_polygon = [] 269 pos = positions[pos_index] 270 f_out.write(f"{obj_name} %.4f %.4f %.4f " % (pos[0], pos[1], pos[2])) 271 for bound_p in instance_dict[obj_name][pos_index]["obj_boundary2d"]: 272 coord = (pos[0]-bound_p[0], pos[1]-bound_p[1]) 273 f_out.write("%.4f %.4f " % coord) 274 pos_polygon.append((int(coord[0]/mask_img_resolution), int(coord[1]/mask_img_resolution))) 275 f_out.write("\n") 276 polygons[obj_name].append(pos_polygon) 277 f_out.close() 278 279 if out_mask_img_file is not None and out_mask_img_file != "": 280 from PIL import Image, ImageDraw 281 image_width = int(np.ceil(self.width_extent / mask_img_resolution)) 282 image_height = int(np.ceil(self.height_extent / mask_img_resolution)) 283 img = Image.new('L', (image_width, image_height), 0) 284 obj_index = 1 285 for obj_name in polygons: 286 for polygon in polygons[obj_name]: 287 ImageDraw.Draw(img).polygon(polygon, outline=obj_index, fill=obj_index) 288 obj_index += 1 289 mask = np.array(img).astype(int) 290 out_path = out_mask_img_file 291 if out_path.endswith(".tif"): 292 out_path = out_path[0:-4] 293 RasterHelper.saveToHdr_no_transform(mask, out_path, [], "GTiff") 294 295 def parse_meshes(self, pos, meshes, slice_planes): 296 for mesh in meshes: 297 new_mesh = mesh.copy() 298 new_mesh.vertices[:, 0] *= -1 299 new_mesh.vertices[:, [1, 2]] = new_mesh.vertices[:, [2, 1]] 300 new_mesh.vertices[:, 1] *= -1 301 new_mesh.apply_translation([pos[0], pos[1], pos[2]]) 302 slice_mesh = new_mesh 303 for slice_id in slice_planes: 304 slice_mesh = trimesh.intersections.slice_mesh_plane(slice_mesh, self.__slice_planes[slice_id][0] 305 , self.__slice_planes[slice_id][1]) 306 self.tot_area += slice_mesh.area 307 308 def compute_scene_lai(self, exclude_object_outside=False, ignore_list=None): 309 if ignore_list is None: 310 ignore_list = [] 311 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 312 objects = self.land_scape.get_objects() 313 instances = self.land_scape.get_instances() 314 obj_areas = {} 315 obj_bounds = {} 316 obj_meshes = {} 317 tot_area = 0 318 for scene_obj in objects: 319 scene_obj_area = 0 320 obj_bounding_box = LSBoundingBox() 321 meshes = [] 322 for obj_component in objects[scene_obj]: 323 if obj_component not in ignore_list: 324 obj_component_path = os.path.join(parameter_dir, obj_component) 325 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 326 scene_obj_area += obj_comp_mesh.area 327 bbox = obj_comp_mesh.bounding_box 328 comp_bounding_box = LSBoundingBox() 329 comp_bounding_box.init_from_bounds(bbox.bounds) 330 obj_bounding_box.add_child(comp_bounding_box) 331 meshes.append(obj_comp_mesh) 332 tot_area += len(instances[scene_obj]) * scene_obj_area 333 obj_areas[scene_obj] = scene_obj_area 334 obj_bounds[scene_obj] = obj_bounding_box 335 obj_meshes[scene_obj] = meshes 336 if not exclude_object_outside: 337 return tot_area / (self.width_extent * self.height_extent) 338 else: 339 for instance in instances: 340 for pos in instances[instance]: 341 obj_bound = obj_bounds[instance] 342 obj_mesh = obj_meshes[instance] 343 obj_area = obj_areas[instance] 344 if len(pos) == 10: # scale 345 x_scale_size = pos[7] 346 y_scale_size = pos[9] # vertical direction 347 z_scale_size = pos[8] 348 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 349 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 350 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 351 obj_mesh = [] 352 obj_bound = LSBoundingBox() 353 obj_area = 0 354 for mesh in obj_meshes[instance]: 355 new_mesh = mesh.copy() 356 new_mesh.vertices[:, 0] *= x_scale_factor 357 new_mesh.vertices[:, 1] *= y_scale_factor 358 new_mesh.vertices[:, 2] *= z_scale_factor 359 obj_mesh.append(new_mesh) 360 obj_area += new_mesh.area 361 comp_bounding_box = LSBoundingBox() 362 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 363 obj_bound.add_child(comp_bounding_box) 364 if pos[3] != 0: 365 # recalculate bounds 366 obj_bound = LSBoundingBox() 367 obj_mesh1 = [] 368 obj_area = 0 369 for mesh in obj_mesh: 370 new_mesh = mesh.copy() 371 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 372 axis = [0, 1, 0] 373 if len(pos) >= 7: 374 axis = [-pos[4], pos[6], -pos[5]] 375 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 376 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 377 obj_mesh1.append(rotated_mesh) 378 obj_area += rotated_mesh.area 379 comp_bounding_box = LSBoundingBox() 380 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 381 obj_bound.add_child(comp_bounding_box) 382 obj_mesh = obj_mesh1 383 left_x = pos[0] - obj_bound.maxX 384 right_x = pos[0] - obj_bound.minX 385 upper_y = pos[1] - obj_bound.maxZ 386 lower_y = pos[1] - obj_bound.minZ 387 if left_x >= 0 and right_x <= self.width_extent \ 388 and upper_y >= 0 and lower_y <= self.height_extent: 389 self.tot_area += obj_area 390 else: 391 slice_planes = [] # 0 1 2 3 for left lower right upper 392 if left_x < 0: 393 slice_planes.append(0) 394 if right_x > self.width_extent: 395 slice_planes.append(2) 396 if upper_y < 0: 397 slice_planes.append(3) 398 if lower_y > self.height_extent: 399 slice_planes.append(1) 400 self.parse_meshes(pos, obj_mesh, slice_planes) 401 return self.tot_area / (self.width_extent * self.height_extent) 402 403 def compute_lai3d(self, out_path, rows, cols, layers, compute_width=-1, compute_height=-1, ignore_list=None): 404 if ignore_list is None: 405 ignore_list = [] 406 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 407 if compute_width <= 0: 408 compute_width = self.width_extent 409 if compute_height <= 0: 410 compute_height = self.height_extent 411 c_left_x = 0.5 * (self.width_extent - compute_width) 412 c_left_y = 0.5 * (self.height_extent - compute_height) 413 lai3d = np.zeros((rows, cols, layers)) 414 step_x = compute_width / cols 415 step_y = compute_height / rows 416 objects = self.land_scape.get_objects() 417 instances = self.land_scape.get_instances() 418 obj_areas = {} 419 obj_bounds = {} 420 obj_meshes = {} 421 for scene_obj in objects: 422 scene_obj_area = 0 423 obj_bounding_box = LSBoundingBox() 424 meshes = [] 425 for obj_component in objects[scene_obj]: 426 if obj_component not in ignore_list: 427 obj_component_path = os.path.join(parameter_dir, obj_component) 428 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 429 scene_obj_area += obj_comp_mesh.area 430 bbox = obj_comp_mesh.bounding_box 431 comp_bounding_box = LSBoundingBox() 432 comp_bounding_box.init_from_bounds(bbox.bounds) 433 obj_bounding_box.add_child(comp_bounding_box) 434 meshes.append(obj_comp_mesh) 435 obj_areas[scene_obj] = scene_obj_area 436 obj_bounds[scene_obj] = obj_bounding_box 437 obj_meshes[scene_obj] = meshes 438 439 instance_dict = {} 440 depth = 0 441 for instance in instances: 442 instance_dict[instance] = {} 443 pos_index = 0 444 for pos in instances[instance]: 445 instance_dict[instance][pos_index] = {} 446 obj_bound = obj_bounds[instance] 447 obj_mesh = obj_meshes[instance] 448 obj_area = obj_areas[instance] 449 if len(pos) == 10: # scale 450 x_scale_size = pos[7] 451 y_scale_size = pos[9] # vertical direction 452 z_scale_size = pos[8] 453 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 454 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 455 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 456 obj_mesh = [] 457 obj_bound = LSBoundingBox() 458 obj_area = 0 459 for mesh in obj_meshes[instance]: 460 new_mesh = mesh.copy() 461 new_mesh.vertices[:, 0] *= x_scale_factor 462 new_mesh.vertices[:, 1] *= y_scale_factor 463 new_mesh.vertices[:, 2] *= z_scale_factor 464 obj_mesh.append(new_mesh) 465 obj_area += new_mesh.area 466 comp_bounding_box = LSBoundingBox() 467 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 468 obj_bound.add_child(comp_bounding_box) 469 470 if pos[3] != 0: 471 obj_bound = LSBoundingBox() 472 obj_mesh1 = [] 473 obj_area = 0 474 for mesh in obj_mesh: 475 new_mesh = mesh.copy() 476 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 477 axis = [0, 1, 0] 478 if len(pos) >= 7: 479 axis = [-pos[4], pos[6], -pos[5]] 480 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 481 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 482 obj_mesh1.append(rotated_mesh) 483 obj_area += rotated_mesh.area 484 comp_bounding_box = LSBoundingBox() 485 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 486 obj_bound.add_child(comp_bounding_box) 487 obj_mesh = obj_mesh1 488 tmp_height = obj_bound.maxY + pos[2] 489 if tmp_height > depth: 490 depth = tmp_height 491 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 492 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 493 instance_dict[instance][pos_index]["obj_area"] = obj_area 494 pos_index += 1 495 depth += 1e-5 # a small offset 496 step_z = depth / layers 497 for obj_name in instance_dict: 498 positions = instances[obj_name] 499 for pos_index in instance_dict[obj_name]: 500 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 501 obj_mesh = instance_dict[obj_name][pos_index]["obj_mesh"] 502 obj_area = instance_dict[obj_name][pos_index]["obj_area"] 503 pos = positions[pos_index] 504 left_x = pos[0] - obj_bound.maxX 505 right_x = pos[0] - obj_bound.minX 506 upper_y = pos[1] - obj_bound.maxZ 507 lower_y = pos[1] - obj_bound.minZ 508 lower_z = pos[2] + obj_bound.minY 509 upper_z = pos[2] + obj_bound.maxY 510 left_index = int(np.floor((left_x - c_left_x) / step_x)) 511 right_index = int(np.floor((right_x - c_left_x) / step_x)) 512 upper_index = int(np.floor((upper_y - c_left_y) / step_y)) 513 lower_index = int(np.floor((lower_y - c_left_y) / step_y)) 514 bottom_index = int(np.floor(lower_z / step_z)) 515 above_index = int(np.floor(upper_z / step_z)) 516 if left_index == right_index and upper_index == lower_index and bottom_index == above_index: 517 if 0 <= left_index < cols and 0 <= lower_index < rows and 0 <= bottom_index < layers: 518 lai3d[lower_index, left_index, bottom_index] += obj_area 519 else: 520 for mesh in obj_mesh: 521 tri_centers = mesh.triangles_center 522 tri_areas = mesh.area_faces 523 vert_x = -tri_centers[:, 0] 524 vert_y = -tri_centers[:, 2] 525 vert_z = tri_centers[:, 1] 526 vert_x_trans = vert_x + pos[0] 527 vert_y_trans = vert_y + pos[1] 528 vert_z_trans = vert_z + pos[2] 529 tri_rows = np.floor((vert_y_trans - c_left_y) / step_y).astype(int) 530 tri_cols = np.floor((vert_x_trans - c_left_x) / step_x).astype(int) 531 tri_layers = np.floor(vert_z_trans / step_z).astype(int) 532 unique_rows = set(tri_rows) 533 unique_cols = set(tri_cols) 534 unique_layers = set(tri_layers) 535 for row in unique_rows: 536 for col in unique_cols: 537 for layer in unique_layers: 538 query = (tri_rows == row) & (tri_cols == col) & (tri_layers == layer) 539 area_values = tri_areas[query] 540 if len(area_values) > 0: 541 if 0 <= row < rows and 0 <= col < cols and 0 <= layer < layers: 542 lai3d[row, col, layer] += sum(area_values) 543 if out_path == "": 544 out_path = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Results", "LAI.txt") 545 print("Save file to:", out_path) 546 f = open(out_path, "w") 547 f.write("Scene Size [X, Y, Z]: %.4f, %.4f, %.4f\n" % (compute_width, compute_height, depth)) 548 f.write("LAI Dimension [Rows, Cols, layers]: %d, %d, %d\n" % (rows, cols, layers)) 549 for d in range(0, layers): 550 layer_lai = lai3d[:, :, d] 551 for row in range(rows): 552 for col in range(cols): 553 f.write("%.6f " % (layer_lai[row, col] / (step_x * step_y))) 554 f.write("\n") 555 f.write("\n") 556 f.close() 557 558 def compute_chm(self, out_path, resolution, ignore_list=None): 559 if ignore_list is None: 560 ignore_list = [] 561 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 562 objects = self.land_scape.get_objects() 563 instances = self.land_scape.get_instances() 564 tot_cols = int(np.ceil(self.width_extent / resolution)) 565 tot_rows = int(np.ceil(self.height_extent / resolution)) 566 dsm = np.zeros((tot_rows, tot_cols)) 567 obj_bounds = {} 568 obj_meshes = {} 569 for scene_obj in objects: 570 obj_bounding_box = LSBoundingBox() 571 meshes = [] 572 for obj_component in objects[scene_obj]: 573 if obj_component not in ignore_list: 574 obj_component_path = os.path.join(parameter_dir, obj_component) 575 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 576 bbox = obj_comp_mesh.bounding_box 577 comp_bounding_box = LSBoundingBox() 578 comp_bounding_box.init_from_bounds(bbox.bounds) 579 obj_bounding_box.add_child(comp_bounding_box) 580 meshes.append(obj_comp_mesh) 581 obj_bounds[scene_obj] = obj_bounding_box 582 obj_meshes[scene_obj] = meshes 583 584 for instance in instances: 585 for pos in instances[instance]: 586 obj_bound = obj_bounds[instance] 587 obj_mesh = obj_meshes[instance] 588 if len(pos) == 10: # scale 589 x_scale_size = pos[7] 590 y_scale_size = pos[9] # vertical direction 591 z_scale_size = pos[8] 592 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 593 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 594 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 595 obj_mesh = [] 596 for mesh in obj_meshes[instance]: 597 new_mesh = mesh.copy() 598 new_mesh.vertices[:, 0] *= x_scale_factor 599 new_mesh.vertices[:, 1] *= y_scale_factor 600 new_mesh.vertices[:, 2] *= z_scale_factor 601 obj_mesh.append(new_mesh) 602 if pos[3] != 0: 603 obj_mesh1 = [] 604 for mesh in obj_mesh: 605 new_mesh = mesh.copy() 606 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 607 axis = [0, 1, 0] 608 if len(pos) >= 7: 609 axis = [-pos[4], pos[6], -pos[5]] 610 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 611 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 612 obj_mesh1.append(rotated_mesh) 613 obj_mesh = obj_mesh1 614 for mesh in obj_mesh: 615 vertices = mesh.vertices 616 vert_x = -vertices[:, 0] 617 vert_y = -vertices[:, 2] 618 vert_z = vertices[:, 1] 619 vert_x_trans = vert_x + pos[0] 620 vert_y_trans = vert_y + pos[1] 621 vert_z_trans = vert_z + pos[2] 622 rows = (vert_y_trans / resolution).astype(int) 623 cols = (vert_x_trans / resolution).astype(int) 624 unique_rows = set(rows) 625 unique_cols = set(cols) 626 for row in unique_rows: 627 for col in unique_cols: 628 query = (rows == row) & (cols == col) 629 z_values = vert_z_trans[query] 630 if len(z_values) > 0: 631 max_value = z_values.max() 632 if 0 <= row < tot_rows and 0 <= col < tot_cols and max_value > dsm[row, col]: 633 dsm[row, col] = max_value 634 if out_path.endswith(".tif"): 635 out_path = out_path[0:-4] 636 RasterHelper.saveToHdr_no_transform(dsm, out_path, [], "GTiff")
class
LandscapeUtility:
14class LandscapeUtility: 15 def __init__(self, land_scape): 16 self.land_scape = land_scape 17 self.width_extent = self.land_scape.get_terrain().extent_width 18 self.height_extent = self.land_scape.get_terrain().extent_height 19 self.tot_area = 0 20 21 # 0 1 2 3 for left lower right upper 22 self.__slice_planes = {0: [[1, 0, 0], [0, 0, 0]], 23 1: [[0, -1, 0], [0, self.height_extent, 0]], 24 2: [[-1, 0, 0], [self.width_extent, 0, 0]], 25 3: [[0, 1, 0], [0, 0, 0]]} 26 27 def __get_rotated_scaled_meshes(self, ignore_list=None): 28 if ignore_list is None: 29 ignore_list = [] 30 objects = self.land_scape.get_objects() 31 instances = self.land_scape.get_instances() 32 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 33 obj_areas = {} 34 obj_bounds = {} 35 obj_meshes = {} 36 for scene_obj in objects: 37 scene_obj_area = 0 38 obj_bounding_box = LSBoundingBox() 39 meshes = [] 40 for obj_component in objects[scene_obj]: 41 if obj_component not in ignore_list: 42 obj_component_path = os.path.join(parameter_dir, obj_component) 43 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 44 scene_obj_area += obj_comp_mesh.area 45 bbox = obj_comp_mesh.bounding_box 46 comp_bounding_box = LSBoundingBox() 47 comp_bounding_box.init_from_bounds(bbox.bounds) 48 obj_bounding_box.add_child(comp_bounding_box) 49 meshes.append(obj_comp_mesh) 50 obj_areas[scene_obj] = scene_obj_area 51 obj_bounds[scene_obj] = obj_bounding_box 52 obj_meshes[scene_obj] = meshes 53 54 instance_dict = {} 55 for instance in instances: 56 instance_dict[instance] = {} 57 pos_index = 0 58 for pos in instances[instance]: 59 instance_dict[instance][pos_index] = {} 60 obj_bound = obj_bounds[instance] 61 obj_mesh = obj_meshes[instance] 62 obj_area = obj_areas[instance] 63 if len(pos) == 10: # scale 64 x_scale_size = pos[7] 65 y_scale_size = pos[9] # vertical direction 66 z_scale_size = pos[8] 67 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 68 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 69 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 70 obj_mesh = [] 71 obj_bound = LSBoundingBox() 72 obj_area = 0 73 for mesh in obj_meshes[instance]: 74 new_mesh = mesh.copy() 75 new_mesh.vertices[:, 0] *= x_scale_factor 76 new_mesh.vertices[:, 1] *= y_scale_factor 77 new_mesh.vertices[:, 2] *= z_scale_factor 78 obj_mesh.append(new_mesh) 79 obj_area += new_mesh.area 80 comp_bounding_box = LSBoundingBox() 81 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 82 obj_bound.add_child(comp_bounding_box) 83 84 if pos[3] != 0: 85 obj_bound = LSBoundingBox() 86 obj_mesh1 = [] 87 obj_area = 0 88 for mesh in obj_mesh: 89 new_mesh = mesh.copy() 90 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 91 axis = [0, 1, 0] 92 if len(pos) >= 7: 93 axis = [-pos[4], pos[6], -pos[5]] 94 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 95 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 96 obj_mesh1.append(rotated_mesh) 97 obj_area += rotated_mesh.area 98 comp_bounding_box = LSBoundingBox() 99 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 100 obj_bound.add_child(comp_bounding_box) 101 obj_mesh = obj_mesh1 102 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 103 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 104 instance_dict[instance][pos_index]["obj_area"] = obj_area 105 pos_index += 1 106 return instance_dict 107 108 def compute_crown_radius(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5): 109 instance_dict = self.__get_rotated_scaled_meshes() 110 instances = self.land_scape.get_instances() 111 f_out = open(out_file, "w") 112 f_out.write("OBJ_Name X Y Z Radius_X Radius_Y\n") 113 for obj_name in instance_dict: 114 positions = instances[obj_name] 115 for pos_index in instance_dict[obj_name]: 116 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 117 pos = positions[pos_index] 118 radius_x = obj_bound.get_x_extents()*0.5 119 radius_z = obj_bound.get_z_extents()*0.5 120 f_out.write(f"{obj_name} %.4f %.4f %.4f %.4f %.4f\n" % (pos[0], pos[1], pos[2], radius_x, radius_z)) 121 f_out.close() 122 123 # write image 124 if out_mask_img_file is not None and out_mask_img_file != "": 125 from PIL import Image, ImageDraw 126 image_width = int(np.ceil(self.width_extent / mask_img_resolution)) 127 image_height = int(np.ceil(self.height_extent / mask_img_resolution)) 128 img = Image.new('L', (image_width, image_height), 0) 129 obj_index = 1 130 for obj_name in instance_dict: 131 positions = instances[obj_name] 132 for pos_index in instance_dict[obj_name]: 133 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 134 pos = positions[pos_index] 135 shape = [(int((pos[0]-obj_bound.maxX)/mask_img_resolution), int((pos[1]-obj_bound.maxZ)/mask_img_resolution)), 136 (int((pos[0]-obj_bound.minX)/mask_img_resolution), int((pos[1]-obj_bound.minZ)/mask_img_resolution))] 137 ImageDraw.Draw(img).ellipse(shape, outline=obj_index, fill=obj_index) 138 obj_index += 1 139 mask = np.array(img).astype(int) 140 out_path = out_mask_img_file 141 if out_path.endswith(".tif"): 142 out_path = out_path[0:-4] 143 RasterHelper.saveToHdr_no_transform(mask, out_path, [], "GTiff") 144 145 def __compute_alphashape_2d(self, point_2d): 146 if len(point_2d) > 5000: 147 indexes = np.random.randint(0, point_2d.shape[0], size=5000) 148 point_2d = point_2d[indexes] 149 radius = 5 150 alpha_shape = alphashape.alphashape(point_2d, radius) 151 while 1: 152 try: 153 coords = np.array(alpha_shape.boundary.coords) 154 except NotImplementedError: 155 radius += -1 156 alpha_shape = alphashape.alphashape(point_2d, radius) 157 else: 158 break 159 return coords 160 161 def compute_crown_boundary(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5): 162 objects = self.land_scape.get_objects() 163 instances = self.land_scape.get_instances() 164 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 165 obj_areas = {} 166 obj_bounds = {} 167 obj_meshes = {} 168 obj_boundary2ds = {} 169 for scene_obj in objects: 170 scene_obj_area = 0 171 obj_bounding_box = LSBoundingBox() 172 meshes = [] 173 obj_points_2d = None 174 for obj_component in objects[scene_obj]: 175 obj_component_path = os.path.join(parameter_dir, obj_component) 176 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 177 scene_obj_area += obj_comp_mesh.area 178 bbox = obj_comp_mesh.bounding_box 179 comp_bounding_box = LSBoundingBox() 180 comp_bounding_box.init_from_bounds(bbox.bounds) 181 obj_bounding_box.add_child(comp_bounding_box) 182 meshes.append(obj_comp_mesh) 183 if obj_points_2d is None: 184 obj_points_2d = obj_comp_mesh.triangles_center[:, [0,2]] 185 else: 186 obj_points_2d = trimesh.util.vstack_empty((obj_points_2d, obj_comp_mesh.triangles_center[:, [0,2]])) 187 # catch exception if alphashape has muilt-polygons 188 obj_boundary2ds[scene_obj] = self.__compute_alphashape_2d(obj_points_2d) 189 obj_areas[scene_obj] = scene_obj_area 190 obj_bounds[scene_obj] = obj_bounding_box 191 obj_meshes[scene_obj] = meshes 192 193 instance_dict = {} 194 for instance in instances: 195 instance_dict[instance] = {} 196 pos_index = 0 197 for pos in instances[instance]: 198 instance_dict[instance][pos_index] = {} 199 obj_bound = obj_bounds[instance] 200 obj_mesh = obj_meshes[instance] 201 obj_area = obj_areas[instance] 202 obj_boundary2d = obj_boundary2ds[instance] 203 if len(pos) == 10: # scale 204 x_scale_size = pos[7] 205 y_scale_size = pos[9] # vertical direction 206 z_scale_size = pos[8] 207 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 208 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 209 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 210 obj_mesh = [] 211 obj_bound = LSBoundingBox() 212 obj_area = 0 213 obj_points_2d = None 214 for mesh in obj_meshes[instance]: 215 new_mesh = mesh.copy() 216 new_mesh.vertices[:, 0] *= x_scale_factor 217 new_mesh.vertices[:, 1] *= y_scale_factor 218 new_mesh.vertices[:, 2] *= z_scale_factor 219 obj_mesh.append(new_mesh) 220 obj_area += new_mesh.area 221 comp_bounding_box = LSBoundingBox() 222 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 223 obj_bound.add_child(comp_bounding_box) 224 if obj_points_2d is None: 225 obj_points_2d = new_mesh.triangles_center[:, [0,2]] 226 else: 227 obj_points_2d = trimesh.util.vstack_empty( 228 (obj_points_2d, new_mesh.triangles_center[:, [0,2]])) 229 obj_boundary2d = self.__compute_alphashape_2d(obj_points_2d) 230 231 if pos[3] != 0: 232 obj_bound = LSBoundingBox() 233 obj_mesh1 = [] 234 obj_area = 0 235 obj_points_2d = None 236 for mesh in obj_mesh: 237 new_mesh = mesh.copy() 238 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 239 axis = [0, 1, 0] 240 if len(pos) >= 7: 241 axis = [-pos[4], pos[6], -pos[5]] 242 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 243 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 244 obj_mesh1.append(rotated_mesh) 245 obj_area += rotated_mesh.area 246 comp_bounding_box = LSBoundingBox() 247 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 248 obj_bound.add_child(comp_bounding_box) 249 if obj_points_2d is None: 250 obj_points_2d = new_mesh.triangles_center[:, [0,2]] 251 else: 252 obj_points_2d = trimesh.util.vstack_empty( 253 (obj_points_2d, new_mesh.triangles_center[:, [0,2]])) 254 obj_mesh = obj_mesh1 255 obj_boundary2d = self.__compute_alphashape_2d(obj_points_2d) 256 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 257 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 258 instance_dict[instance][pos_index]["obj_area"] = obj_area 259 instance_dict[instance][pos_index]["obj_boundary2d"] = obj_boundary2d 260 pos_index += 1 261 # output 262 f_out = open(out_file, "w") 263 f_out.write("OBJ_Name X Y Z Boundary_Point1_X Boundary_Point1_Y, Boundary_Point2_X, Boundary_Point2_Y,...\n") 264 polygons = dict() 265 for obj_name in instance_dict: 266 positions = instances[obj_name] 267 polygons[obj_name] = [] 268 for pos_index in instance_dict[obj_name]: 269 pos_polygon = [] 270 pos = positions[pos_index] 271 f_out.write(f"{obj_name} %.4f %.4f %.4f " % (pos[0], pos[1], pos[2])) 272 for bound_p in instance_dict[obj_name][pos_index]["obj_boundary2d"]: 273 coord = (pos[0]-bound_p[0], pos[1]-bound_p[1]) 274 f_out.write("%.4f %.4f " % coord) 275 pos_polygon.append((int(coord[0]/mask_img_resolution), int(coord[1]/mask_img_resolution))) 276 f_out.write("\n") 277 polygons[obj_name].append(pos_polygon) 278 f_out.close() 279 280 if out_mask_img_file is not None and out_mask_img_file != "": 281 from PIL import Image, ImageDraw 282 image_width = int(np.ceil(self.width_extent / mask_img_resolution)) 283 image_height = int(np.ceil(self.height_extent / mask_img_resolution)) 284 img = Image.new('L', (image_width, image_height), 0) 285 obj_index = 1 286 for obj_name in polygons: 287 for polygon in polygons[obj_name]: 288 ImageDraw.Draw(img).polygon(polygon, outline=obj_index, fill=obj_index) 289 obj_index += 1 290 mask = np.array(img).astype(int) 291 out_path = out_mask_img_file 292 if out_path.endswith(".tif"): 293 out_path = out_path[0:-4] 294 RasterHelper.saveToHdr_no_transform(mask, out_path, [], "GTiff") 295 296 def parse_meshes(self, pos, meshes, slice_planes): 297 for mesh in meshes: 298 new_mesh = mesh.copy() 299 new_mesh.vertices[:, 0] *= -1 300 new_mesh.vertices[:, [1, 2]] = new_mesh.vertices[:, [2, 1]] 301 new_mesh.vertices[:, 1] *= -1 302 new_mesh.apply_translation([pos[0], pos[1], pos[2]]) 303 slice_mesh = new_mesh 304 for slice_id in slice_planes: 305 slice_mesh = trimesh.intersections.slice_mesh_plane(slice_mesh, self.__slice_planes[slice_id][0] 306 , self.__slice_planes[slice_id][1]) 307 self.tot_area += slice_mesh.area 308 309 def compute_scene_lai(self, exclude_object_outside=False, ignore_list=None): 310 if ignore_list is None: 311 ignore_list = [] 312 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 313 objects = self.land_scape.get_objects() 314 instances = self.land_scape.get_instances() 315 obj_areas = {} 316 obj_bounds = {} 317 obj_meshes = {} 318 tot_area = 0 319 for scene_obj in objects: 320 scene_obj_area = 0 321 obj_bounding_box = LSBoundingBox() 322 meshes = [] 323 for obj_component in objects[scene_obj]: 324 if obj_component not in ignore_list: 325 obj_component_path = os.path.join(parameter_dir, obj_component) 326 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 327 scene_obj_area += obj_comp_mesh.area 328 bbox = obj_comp_mesh.bounding_box 329 comp_bounding_box = LSBoundingBox() 330 comp_bounding_box.init_from_bounds(bbox.bounds) 331 obj_bounding_box.add_child(comp_bounding_box) 332 meshes.append(obj_comp_mesh) 333 tot_area += len(instances[scene_obj]) * scene_obj_area 334 obj_areas[scene_obj] = scene_obj_area 335 obj_bounds[scene_obj] = obj_bounding_box 336 obj_meshes[scene_obj] = meshes 337 if not exclude_object_outside: 338 return tot_area / (self.width_extent * self.height_extent) 339 else: 340 for instance in instances: 341 for pos in instances[instance]: 342 obj_bound = obj_bounds[instance] 343 obj_mesh = obj_meshes[instance] 344 obj_area = obj_areas[instance] 345 if len(pos) == 10: # scale 346 x_scale_size = pos[7] 347 y_scale_size = pos[9] # vertical direction 348 z_scale_size = pos[8] 349 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 350 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 351 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 352 obj_mesh = [] 353 obj_bound = LSBoundingBox() 354 obj_area = 0 355 for mesh in obj_meshes[instance]: 356 new_mesh = mesh.copy() 357 new_mesh.vertices[:, 0] *= x_scale_factor 358 new_mesh.vertices[:, 1] *= y_scale_factor 359 new_mesh.vertices[:, 2] *= z_scale_factor 360 obj_mesh.append(new_mesh) 361 obj_area += new_mesh.area 362 comp_bounding_box = LSBoundingBox() 363 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 364 obj_bound.add_child(comp_bounding_box) 365 if pos[3] != 0: 366 # recalculate bounds 367 obj_bound = LSBoundingBox() 368 obj_mesh1 = [] 369 obj_area = 0 370 for mesh in obj_mesh: 371 new_mesh = mesh.copy() 372 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 373 axis = [0, 1, 0] 374 if len(pos) >= 7: 375 axis = [-pos[4], pos[6], -pos[5]] 376 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 377 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 378 obj_mesh1.append(rotated_mesh) 379 obj_area += rotated_mesh.area 380 comp_bounding_box = LSBoundingBox() 381 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 382 obj_bound.add_child(comp_bounding_box) 383 obj_mesh = obj_mesh1 384 left_x = pos[0] - obj_bound.maxX 385 right_x = pos[0] - obj_bound.minX 386 upper_y = pos[1] - obj_bound.maxZ 387 lower_y = pos[1] - obj_bound.minZ 388 if left_x >= 0 and right_x <= self.width_extent \ 389 and upper_y >= 0 and lower_y <= self.height_extent: 390 self.tot_area += obj_area 391 else: 392 slice_planes = [] # 0 1 2 3 for left lower right upper 393 if left_x < 0: 394 slice_planes.append(0) 395 if right_x > self.width_extent: 396 slice_planes.append(2) 397 if upper_y < 0: 398 slice_planes.append(3) 399 if lower_y > self.height_extent: 400 slice_planes.append(1) 401 self.parse_meshes(pos, obj_mesh, slice_planes) 402 return self.tot_area / (self.width_extent * self.height_extent) 403 404 def compute_lai3d(self, out_path, rows, cols, layers, compute_width=-1, compute_height=-1, ignore_list=None): 405 if ignore_list is None: 406 ignore_list = [] 407 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 408 if compute_width <= 0: 409 compute_width = self.width_extent 410 if compute_height <= 0: 411 compute_height = self.height_extent 412 c_left_x = 0.5 * (self.width_extent - compute_width) 413 c_left_y = 0.5 * (self.height_extent - compute_height) 414 lai3d = np.zeros((rows, cols, layers)) 415 step_x = compute_width / cols 416 step_y = compute_height / rows 417 objects = self.land_scape.get_objects() 418 instances = self.land_scape.get_instances() 419 obj_areas = {} 420 obj_bounds = {} 421 obj_meshes = {} 422 for scene_obj in objects: 423 scene_obj_area = 0 424 obj_bounding_box = LSBoundingBox() 425 meshes = [] 426 for obj_component in objects[scene_obj]: 427 if obj_component not in ignore_list: 428 obj_component_path = os.path.join(parameter_dir, obj_component) 429 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 430 scene_obj_area += obj_comp_mesh.area 431 bbox = obj_comp_mesh.bounding_box 432 comp_bounding_box = LSBoundingBox() 433 comp_bounding_box.init_from_bounds(bbox.bounds) 434 obj_bounding_box.add_child(comp_bounding_box) 435 meshes.append(obj_comp_mesh) 436 obj_areas[scene_obj] = scene_obj_area 437 obj_bounds[scene_obj] = obj_bounding_box 438 obj_meshes[scene_obj] = meshes 439 440 instance_dict = {} 441 depth = 0 442 for instance in instances: 443 instance_dict[instance] = {} 444 pos_index = 0 445 for pos in instances[instance]: 446 instance_dict[instance][pos_index] = {} 447 obj_bound = obj_bounds[instance] 448 obj_mesh = obj_meshes[instance] 449 obj_area = obj_areas[instance] 450 if len(pos) == 10: # scale 451 x_scale_size = pos[7] 452 y_scale_size = pos[9] # vertical direction 453 z_scale_size = pos[8] 454 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 455 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 456 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 457 obj_mesh = [] 458 obj_bound = LSBoundingBox() 459 obj_area = 0 460 for mesh in obj_meshes[instance]: 461 new_mesh = mesh.copy() 462 new_mesh.vertices[:, 0] *= x_scale_factor 463 new_mesh.vertices[:, 1] *= y_scale_factor 464 new_mesh.vertices[:, 2] *= z_scale_factor 465 obj_mesh.append(new_mesh) 466 obj_area += new_mesh.area 467 comp_bounding_box = LSBoundingBox() 468 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 469 obj_bound.add_child(comp_bounding_box) 470 471 if pos[3] != 0: 472 obj_bound = LSBoundingBox() 473 obj_mesh1 = [] 474 obj_area = 0 475 for mesh in obj_mesh: 476 new_mesh = mesh.copy() 477 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 478 axis = [0, 1, 0] 479 if len(pos) >= 7: 480 axis = [-pos[4], pos[6], -pos[5]] 481 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 482 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 483 obj_mesh1.append(rotated_mesh) 484 obj_area += rotated_mesh.area 485 comp_bounding_box = LSBoundingBox() 486 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 487 obj_bound.add_child(comp_bounding_box) 488 obj_mesh = obj_mesh1 489 tmp_height = obj_bound.maxY + pos[2] 490 if tmp_height > depth: 491 depth = tmp_height 492 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 493 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 494 instance_dict[instance][pos_index]["obj_area"] = obj_area 495 pos_index += 1 496 depth += 1e-5 # a small offset 497 step_z = depth / layers 498 for obj_name in instance_dict: 499 positions = instances[obj_name] 500 for pos_index in instance_dict[obj_name]: 501 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 502 obj_mesh = instance_dict[obj_name][pos_index]["obj_mesh"] 503 obj_area = instance_dict[obj_name][pos_index]["obj_area"] 504 pos = positions[pos_index] 505 left_x = pos[0] - obj_bound.maxX 506 right_x = pos[0] - obj_bound.minX 507 upper_y = pos[1] - obj_bound.maxZ 508 lower_y = pos[1] - obj_bound.minZ 509 lower_z = pos[2] + obj_bound.minY 510 upper_z = pos[2] + obj_bound.maxY 511 left_index = int(np.floor((left_x - c_left_x) / step_x)) 512 right_index = int(np.floor((right_x - c_left_x) / step_x)) 513 upper_index = int(np.floor((upper_y - c_left_y) / step_y)) 514 lower_index = int(np.floor((lower_y - c_left_y) / step_y)) 515 bottom_index = int(np.floor(lower_z / step_z)) 516 above_index = int(np.floor(upper_z / step_z)) 517 if left_index == right_index and upper_index == lower_index and bottom_index == above_index: 518 if 0 <= left_index < cols and 0 <= lower_index < rows and 0 <= bottom_index < layers: 519 lai3d[lower_index, left_index, bottom_index] += obj_area 520 else: 521 for mesh in obj_mesh: 522 tri_centers = mesh.triangles_center 523 tri_areas = mesh.area_faces 524 vert_x = -tri_centers[:, 0] 525 vert_y = -tri_centers[:, 2] 526 vert_z = tri_centers[:, 1] 527 vert_x_trans = vert_x + pos[0] 528 vert_y_trans = vert_y + pos[1] 529 vert_z_trans = vert_z + pos[2] 530 tri_rows = np.floor((vert_y_trans - c_left_y) / step_y).astype(int) 531 tri_cols = np.floor((vert_x_trans - c_left_x) / step_x).astype(int) 532 tri_layers = np.floor(vert_z_trans / step_z).astype(int) 533 unique_rows = set(tri_rows) 534 unique_cols = set(tri_cols) 535 unique_layers = set(tri_layers) 536 for row in unique_rows: 537 for col in unique_cols: 538 for layer in unique_layers: 539 query = (tri_rows == row) & (tri_cols == col) & (tri_layers == layer) 540 area_values = tri_areas[query] 541 if len(area_values) > 0: 542 if 0 <= row < rows and 0 <= col < cols and 0 <= layer < layers: 543 lai3d[row, col, layer] += sum(area_values) 544 if out_path == "": 545 out_path = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Results", "LAI.txt") 546 print("Save file to:", out_path) 547 f = open(out_path, "w") 548 f.write("Scene Size [X, Y, Z]: %.4f, %.4f, %.4f\n" % (compute_width, compute_height, depth)) 549 f.write("LAI Dimension [Rows, Cols, layers]: %d, %d, %d\n" % (rows, cols, layers)) 550 for d in range(0, layers): 551 layer_lai = lai3d[:, :, d] 552 for row in range(rows): 553 for col in range(cols): 554 f.write("%.6f " % (layer_lai[row, col] / (step_x * step_y))) 555 f.write("\n") 556 f.write("\n") 557 f.close() 558 559 def compute_chm(self, out_path, resolution, ignore_list=None): 560 if ignore_list is None: 561 ignore_list = [] 562 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 563 objects = self.land_scape.get_objects() 564 instances = self.land_scape.get_instances() 565 tot_cols = int(np.ceil(self.width_extent / resolution)) 566 tot_rows = int(np.ceil(self.height_extent / resolution)) 567 dsm = np.zeros((tot_rows, tot_cols)) 568 obj_bounds = {} 569 obj_meshes = {} 570 for scene_obj in objects: 571 obj_bounding_box = LSBoundingBox() 572 meshes = [] 573 for obj_component in objects[scene_obj]: 574 if obj_component not in ignore_list: 575 obj_component_path = os.path.join(parameter_dir, obj_component) 576 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 577 bbox = obj_comp_mesh.bounding_box 578 comp_bounding_box = LSBoundingBox() 579 comp_bounding_box.init_from_bounds(bbox.bounds) 580 obj_bounding_box.add_child(comp_bounding_box) 581 meshes.append(obj_comp_mesh) 582 obj_bounds[scene_obj] = obj_bounding_box 583 obj_meshes[scene_obj] = meshes 584 585 for instance in instances: 586 for pos in instances[instance]: 587 obj_bound = obj_bounds[instance] 588 obj_mesh = obj_meshes[instance] 589 if len(pos) == 10: # scale 590 x_scale_size = pos[7] 591 y_scale_size = pos[9] # vertical direction 592 z_scale_size = pos[8] 593 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 594 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 595 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 596 obj_mesh = [] 597 for mesh in obj_meshes[instance]: 598 new_mesh = mesh.copy() 599 new_mesh.vertices[:, 0] *= x_scale_factor 600 new_mesh.vertices[:, 1] *= y_scale_factor 601 new_mesh.vertices[:, 2] *= z_scale_factor 602 obj_mesh.append(new_mesh) 603 if pos[3] != 0: 604 obj_mesh1 = [] 605 for mesh in obj_mesh: 606 new_mesh = mesh.copy() 607 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 608 axis = [0, 1, 0] 609 if len(pos) >= 7: 610 axis = [-pos[4], pos[6], -pos[5]] 611 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 612 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 613 obj_mesh1.append(rotated_mesh) 614 obj_mesh = obj_mesh1 615 for mesh in obj_mesh: 616 vertices = mesh.vertices 617 vert_x = -vertices[:, 0] 618 vert_y = -vertices[:, 2] 619 vert_z = vertices[:, 1] 620 vert_x_trans = vert_x + pos[0] 621 vert_y_trans = vert_y + pos[1] 622 vert_z_trans = vert_z + pos[2] 623 rows = (vert_y_trans / resolution).astype(int) 624 cols = (vert_x_trans / resolution).astype(int) 625 unique_rows = set(rows) 626 unique_cols = set(cols) 627 for row in unique_rows: 628 for col in unique_cols: 629 query = (rows == row) & (cols == col) 630 z_values = vert_z_trans[query] 631 if len(z_values) > 0: 632 max_value = z_values.max() 633 if 0 <= row < tot_rows and 0 <= col < tot_cols and max_value > dsm[row, col]: 634 dsm[row, col] = max_value 635 if out_path.endswith(".tif"): 636 out_path = out_path[0:-4] 637 RasterHelper.saveToHdr_no_transform(dsm, out_path, [], "GTiff")
LandscapeUtility(land_scape)
15 def __init__(self, land_scape): 16 self.land_scape = land_scape 17 self.width_extent = self.land_scape.get_terrain().extent_width 18 self.height_extent = self.land_scape.get_terrain().extent_height 19 self.tot_area = 0 20 21 # 0 1 2 3 for left lower right upper 22 self.__slice_planes = {0: [[1, 0, 0], [0, 0, 0]], 23 1: [[0, -1, 0], [0, self.height_extent, 0]], 24 2: [[-1, 0, 0], [self.width_extent, 0, 0]], 25 3: [[0, 1, 0], [0, 0, 0]]}
def
compute_crown_radius(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5):
108 def compute_crown_radius(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5): 109 instance_dict = self.__get_rotated_scaled_meshes() 110 instances = self.land_scape.get_instances() 111 f_out = open(out_file, "w") 112 f_out.write("OBJ_Name X Y Z Radius_X Radius_Y\n") 113 for obj_name in instance_dict: 114 positions = instances[obj_name] 115 for pos_index in instance_dict[obj_name]: 116 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 117 pos = positions[pos_index] 118 radius_x = obj_bound.get_x_extents()*0.5 119 radius_z = obj_bound.get_z_extents()*0.5 120 f_out.write(f"{obj_name} %.4f %.4f %.4f %.4f %.4f\n" % (pos[0], pos[1], pos[2], radius_x, radius_z)) 121 f_out.close() 122 123 # write image 124 if out_mask_img_file is not None and out_mask_img_file != "": 125 from PIL import Image, ImageDraw 126 image_width = int(np.ceil(self.width_extent / mask_img_resolution)) 127 image_height = int(np.ceil(self.height_extent / mask_img_resolution)) 128 img = Image.new('L', (image_width, image_height), 0) 129 obj_index = 1 130 for obj_name in instance_dict: 131 positions = instances[obj_name] 132 for pos_index in instance_dict[obj_name]: 133 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 134 pos = positions[pos_index] 135 shape = [(int((pos[0]-obj_bound.maxX)/mask_img_resolution), int((pos[1]-obj_bound.maxZ)/mask_img_resolution)), 136 (int((pos[0]-obj_bound.minX)/mask_img_resolution), int((pos[1]-obj_bound.minZ)/mask_img_resolution))] 137 ImageDraw.Draw(img).ellipse(shape, outline=obj_index, fill=obj_index) 138 obj_index += 1 139 mask = np.array(img).astype(int) 140 out_path = out_mask_img_file 141 if out_path.endswith(".tif"): 142 out_path = out_path[0:-4] 143 RasterHelper.saveToHdr_no_transform(mask, out_path, [], "GTiff")
def
compute_crown_boundary(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5):
161 def compute_crown_boundary(self, out_file, out_mask_img_file=None, mask_img_resolution=0.5): 162 objects = self.land_scape.get_objects() 163 instances = self.land_scape.get_instances() 164 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 165 obj_areas = {} 166 obj_bounds = {} 167 obj_meshes = {} 168 obj_boundary2ds = {} 169 for scene_obj in objects: 170 scene_obj_area = 0 171 obj_bounding_box = LSBoundingBox() 172 meshes = [] 173 obj_points_2d = None 174 for obj_component in objects[scene_obj]: 175 obj_component_path = os.path.join(parameter_dir, obj_component) 176 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 177 scene_obj_area += obj_comp_mesh.area 178 bbox = obj_comp_mesh.bounding_box 179 comp_bounding_box = LSBoundingBox() 180 comp_bounding_box.init_from_bounds(bbox.bounds) 181 obj_bounding_box.add_child(comp_bounding_box) 182 meshes.append(obj_comp_mesh) 183 if obj_points_2d is None: 184 obj_points_2d = obj_comp_mesh.triangles_center[:, [0,2]] 185 else: 186 obj_points_2d = trimesh.util.vstack_empty((obj_points_2d, obj_comp_mesh.triangles_center[:, [0,2]])) 187 # catch exception if alphashape has muilt-polygons 188 obj_boundary2ds[scene_obj] = self.__compute_alphashape_2d(obj_points_2d) 189 obj_areas[scene_obj] = scene_obj_area 190 obj_bounds[scene_obj] = obj_bounding_box 191 obj_meshes[scene_obj] = meshes 192 193 instance_dict = {} 194 for instance in instances: 195 instance_dict[instance] = {} 196 pos_index = 0 197 for pos in instances[instance]: 198 instance_dict[instance][pos_index] = {} 199 obj_bound = obj_bounds[instance] 200 obj_mesh = obj_meshes[instance] 201 obj_area = obj_areas[instance] 202 obj_boundary2d = obj_boundary2ds[instance] 203 if len(pos) == 10: # scale 204 x_scale_size = pos[7] 205 y_scale_size = pos[9] # vertical direction 206 z_scale_size = pos[8] 207 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 208 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 209 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 210 obj_mesh = [] 211 obj_bound = LSBoundingBox() 212 obj_area = 0 213 obj_points_2d = None 214 for mesh in obj_meshes[instance]: 215 new_mesh = mesh.copy() 216 new_mesh.vertices[:, 0] *= x_scale_factor 217 new_mesh.vertices[:, 1] *= y_scale_factor 218 new_mesh.vertices[:, 2] *= z_scale_factor 219 obj_mesh.append(new_mesh) 220 obj_area += new_mesh.area 221 comp_bounding_box = LSBoundingBox() 222 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 223 obj_bound.add_child(comp_bounding_box) 224 if obj_points_2d is None: 225 obj_points_2d = new_mesh.triangles_center[:, [0,2]] 226 else: 227 obj_points_2d = trimesh.util.vstack_empty( 228 (obj_points_2d, new_mesh.triangles_center[:, [0,2]])) 229 obj_boundary2d = self.__compute_alphashape_2d(obj_points_2d) 230 231 if pos[3] != 0: 232 obj_bound = LSBoundingBox() 233 obj_mesh1 = [] 234 obj_area = 0 235 obj_points_2d = None 236 for mesh in obj_mesh: 237 new_mesh = mesh.copy() 238 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 239 axis = [0, 1, 0] 240 if len(pos) >= 7: 241 axis = [-pos[4], pos[6], -pos[5]] 242 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 243 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 244 obj_mesh1.append(rotated_mesh) 245 obj_area += rotated_mesh.area 246 comp_bounding_box = LSBoundingBox() 247 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 248 obj_bound.add_child(comp_bounding_box) 249 if obj_points_2d is None: 250 obj_points_2d = new_mesh.triangles_center[:, [0,2]] 251 else: 252 obj_points_2d = trimesh.util.vstack_empty( 253 (obj_points_2d, new_mesh.triangles_center[:, [0,2]])) 254 obj_mesh = obj_mesh1 255 obj_boundary2d = self.__compute_alphashape_2d(obj_points_2d) 256 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 257 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 258 instance_dict[instance][pos_index]["obj_area"] = obj_area 259 instance_dict[instance][pos_index]["obj_boundary2d"] = obj_boundary2d 260 pos_index += 1 261 # output 262 f_out = open(out_file, "w") 263 f_out.write("OBJ_Name X Y Z Boundary_Point1_X Boundary_Point1_Y, Boundary_Point2_X, Boundary_Point2_Y,...\n") 264 polygons = dict() 265 for obj_name in instance_dict: 266 positions = instances[obj_name] 267 polygons[obj_name] = [] 268 for pos_index in instance_dict[obj_name]: 269 pos_polygon = [] 270 pos = positions[pos_index] 271 f_out.write(f"{obj_name} %.4f %.4f %.4f " % (pos[0], pos[1], pos[2])) 272 for bound_p in instance_dict[obj_name][pos_index]["obj_boundary2d"]: 273 coord = (pos[0]-bound_p[0], pos[1]-bound_p[1]) 274 f_out.write("%.4f %.4f " % coord) 275 pos_polygon.append((int(coord[0]/mask_img_resolution), int(coord[1]/mask_img_resolution))) 276 f_out.write("\n") 277 polygons[obj_name].append(pos_polygon) 278 f_out.close() 279 280 if out_mask_img_file is not None and out_mask_img_file != "": 281 from PIL import Image, ImageDraw 282 image_width = int(np.ceil(self.width_extent / mask_img_resolution)) 283 image_height = int(np.ceil(self.height_extent / mask_img_resolution)) 284 img = Image.new('L', (image_width, image_height), 0) 285 obj_index = 1 286 for obj_name in polygons: 287 for polygon in polygons[obj_name]: 288 ImageDraw.Draw(img).polygon(polygon, outline=obj_index, fill=obj_index) 289 obj_index += 1 290 mask = np.array(img).astype(int) 291 out_path = out_mask_img_file 292 if out_path.endswith(".tif"): 293 out_path = out_path[0:-4] 294 RasterHelper.saveToHdr_no_transform(mask, out_path, [], "GTiff")
def
parse_meshes(self, pos, meshes, slice_planes):
296 def parse_meshes(self, pos, meshes, slice_planes): 297 for mesh in meshes: 298 new_mesh = mesh.copy() 299 new_mesh.vertices[:, 0] *= -1 300 new_mesh.vertices[:, [1, 2]] = new_mesh.vertices[:, [2, 1]] 301 new_mesh.vertices[:, 1] *= -1 302 new_mesh.apply_translation([pos[0], pos[1], pos[2]]) 303 slice_mesh = new_mesh 304 for slice_id in slice_planes: 305 slice_mesh = trimesh.intersections.slice_mesh_plane(slice_mesh, self.__slice_planes[slice_id][0] 306 , self.__slice_planes[slice_id][1]) 307 self.tot_area += slice_mesh.area
def
compute_scene_lai(self, exclude_object_outside=False, ignore_list=None):
309 def compute_scene_lai(self, exclude_object_outside=False, ignore_list=None): 310 if ignore_list is None: 311 ignore_list = [] 312 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 313 objects = self.land_scape.get_objects() 314 instances = self.land_scape.get_instances() 315 obj_areas = {} 316 obj_bounds = {} 317 obj_meshes = {} 318 tot_area = 0 319 for scene_obj in objects: 320 scene_obj_area = 0 321 obj_bounding_box = LSBoundingBox() 322 meshes = [] 323 for obj_component in objects[scene_obj]: 324 if obj_component not in ignore_list: 325 obj_component_path = os.path.join(parameter_dir, obj_component) 326 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 327 scene_obj_area += obj_comp_mesh.area 328 bbox = obj_comp_mesh.bounding_box 329 comp_bounding_box = LSBoundingBox() 330 comp_bounding_box.init_from_bounds(bbox.bounds) 331 obj_bounding_box.add_child(comp_bounding_box) 332 meshes.append(obj_comp_mesh) 333 tot_area += len(instances[scene_obj]) * scene_obj_area 334 obj_areas[scene_obj] = scene_obj_area 335 obj_bounds[scene_obj] = obj_bounding_box 336 obj_meshes[scene_obj] = meshes 337 if not exclude_object_outside: 338 return tot_area / (self.width_extent * self.height_extent) 339 else: 340 for instance in instances: 341 for pos in instances[instance]: 342 obj_bound = obj_bounds[instance] 343 obj_mesh = obj_meshes[instance] 344 obj_area = obj_areas[instance] 345 if len(pos) == 10: # scale 346 x_scale_size = pos[7] 347 y_scale_size = pos[9] # vertical direction 348 z_scale_size = pos[8] 349 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 350 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 351 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 352 obj_mesh = [] 353 obj_bound = LSBoundingBox() 354 obj_area = 0 355 for mesh in obj_meshes[instance]: 356 new_mesh = mesh.copy() 357 new_mesh.vertices[:, 0] *= x_scale_factor 358 new_mesh.vertices[:, 1] *= y_scale_factor 359 new_mesh.vertices[:, 2] *= z_scale_factor 360 obj_mesh.append(new_mesh) 361 obj_area += new_mesh.area 362 comp_bounding_box = LSBoundingBox() 363 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 364 obj_bound.add_child(comp_bounding_box) 365 if pos[3] != 0: 366 # recalculate bounds 367 obj_bound = LSBoundingBox() 368 obj_mesh1 = [] 369 obj_area = 0 370 for mesh in obj_mesh: 371 new_mesh = mesh.copy() 372 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 373 axis = [0, 1, 0] 374 if len(pos) >= 7: 375 axis = [-pos[4], pos[6], -pos[5]] 376 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 377 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 378 obj_mesh1.append(rotated_mesh) 379 obj_area += rotated_mesh.area 380 comp_bounding_box = LSBoundingBox() 381 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 382 obj_bound.add_child(comp_bounding_box) 383 obj_mesh = obj_mesh1 384 left_x = pos[0] - obj_bound.maxX 385 right_x = pos[0] - obj_bound.minX 386 upper_y = pos[1] - obj_bound.maxZ 387 lower_y = pos[1] - obj_bound.minZ 388 if left_x >= 0 and right_x <= self.width_extent \ 389 and upper_y >= 0 and lower_y <= self.height_extent: 390 self.tot_area += obj_area 391 else: 392 slice_planes = [] # 0 1 2 3 for left lower right upper 393 if left_x < 0: 394 slice_planes.append(0) 395 if right_x > self.width_extent: 396 slice_planes.append(2) 397 if upper_y < 0: 398 slice_planes.append(3) 399 if lower_y > self.height_extent: 400 slice_planes.append(1) 401 self.parse_meshes(pos, obj_mesh, slice_planes) 402 return self.tot_area / (self.width_extent * self.height_extent)
def
compute_lai3d( self, out_path, rows, cols, layers, compute_width=-1, compute_height=-1, ignore_list=None):
404 def compute_lai3d(self, out_path, rows, cols, layers, compute_width=-1, compute_height=-1, ignore_list=None): 405 if ignore_list is None: 406 ignore_list = [] 407 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 408 if compute_width <= 0: 409 compute_width = self.width_extent 410 if compute_height <= 0: 411 compute_height = self.height_extent 412 c_left_x = 0.5 * (self.width_extent - compute_width) 413 c_left_y = 0.5 * (self.height_extent - compute_height) 414 lai3d = np.zeros((rows, cols, layers)) 415 step_x = compute_width / cols 416 step_y = compute_height / rows 417 objects = self.land_scape.get_objects() 418 instances = self.land_scape.get_instances() 419 obj_areas = {} 420 obj_bounds = {} 421 obj_meshes = {} 422 for scene_obj in objects: 423 scene_obj_area = 0 424 obj_bounding_box = LSBoundingBox() 425 meshes = [] 426 for obj_component in objects[scene_obj]: 427 if obj_component not in ignore_list: 428 obj_component_path = os.path.join(parameter_dir, obj_component) 429 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 430 scene_obj_area += obj_comp_mesh.area 431 bbox = obj_comp_mesh.bounding_box 432 comp_bounding_box = LSBoundingBox() 433 comp_bounding_box.init_from_bounds(bbox.bounds) 434 obj_bounding_box.add_child(comp_bounding_box) 435 meshes.append(obj_comp_mesh) 436 obj_areas[scene_obj] = scene_obj_area 437 obj_bounds[scene_obj] = obj_bounding_box 438 obj_meshes[scene_obj] = meshes 439 440 instance_dict = {} 441 depth = 0 442 for instance in instances: 443 instance_dict[instance] = {} 444 pos_index = 0 445 for pos in instances[instance]: 446 instance_dict[instance][pos_index] = {} 447 obj_bound = obj_bounds[instance] 448 obj_mesh = obj_meshes[instance] 449 obj_area = obj_areas[instance] 450 if len(pos) == 10: # scale 451 x_scale_size = pos[7] 452 y_scale_size = pos[9] # vertical direction 453 z_scale_size = pos[8] 454 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 455 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 456 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 457 obj_mesh = [] 458 obj_bound = LSBoundingBox() 459 obj_area = 0 460 for mesh in obj_meshes[instance]: 461 new_mesh = mesh.copy() 462 new_mesh.vertices[:, 0] *= x_scale_factor 463 new_mesh.vertices[:, 1] *= y_scale_factor 464 new_mesh.vertices[:, 2] *= z_scale_factor 465 obj_mesh.append(new_mesh) 466 obj_area += new_mesh.area 467 comp_bounding_box = LSBoundingBox() 468 comp_bounding_box.init_from_bounds(new_mesh.bounding_box.bounds) 469 obj_bound.add_child(comp_bounding_box) 470 471 if pos[3] != 0: 472 obj_bound = LSBoundingBox() 473 obj_mesh1 = [] 474 obj_area = 0 475 for mesh in obj_mesh: 476 new_mesh = mesh.copy() 477 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 478 axis = [0, 1, 0] 479 if len(pos) >= 7: 480 axis = [-pos[4], pos[6], -pos[5]] 481 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 482 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 483 obj_mesh1.append(rotated_mesh) 484 obj_area += rotated_mesh.area 485 comp_bounding_box = LSBoundingBox() 486 comp_bounding_box.init_from_bounds(rotated_mesh.bounding_box.bounds) 487 obj_bound.add_child(comp_bounding_box) 488 obj_mesh = obj_mesh1 489 tmp_height = obj_bound.maxY + pos[2] 490 if tmp_height > depth: 491 depth = tmp_height 492 instance_dict[instance][pos_index]["obj_bounds"] = obj_bound 493 instance_dict[instance][pos_index]["obj_mesh"] = obj_mesh 494 instance_dict[instance][pos_index]["obj_area"] = obj_area 495 pos_index += 1 496 depth += 1e-5 # a small offset 497 step_z = depth / layers 498 for obj_name in instance_dict: 499 positions = instances[obj_name] 500 for pos_index in instance_dict[obj_name]: 501 obj_bound = instance_dict[obj_name][pos_index]["obj_bounds"] 502 obj_mesh = instance_dict[obj_name][pos_index]["obj_mesh"] 503 obj_area = instance_dict[obj_name][pos_index]["obj_area"] 504 pos = positions[pos_index] 505 left_x = pos[0] - obj_bound.maxX 506 right_x = pos[0] - obj_bound.minX 507 upper_y = pos[1] - obj_bound.maxZ 508 lower_y = pos[1] - obj_bound.minZ 509 lower_z = pos[2] + obj_bound.minY 510 upper_z = pos[2] + obj_bound.maxY 511 left_index = int(np.floor((left_x - c_left_x) / step_x)) 512 right_index = int(np.floor((right_x - c_left_x) / step_x)) 513 upper_index = int(np.floor((upper_y - c_left_y) / step_y)) 514 lower_index = int(np.floor((lower_y - c_left_y) / step_y)) 515 bottom_index = int(np.floor(lower_z / step_z)) 516 above_index = int(np.floor(upper_z / step_z)) 517 if left_index == right_index and upper_index == lower_index and bottom_index == above_index: 518 if 0 <= left_index < cols and 0 <= lower_index < rows and 0 <= bottom_index < layers: 519 lai3d[lower_index, left_index, bottom_index] += obj_area 520 else: 521 for mesh in obj_mesh: 522 tri_centers = mesh.triangles_center 523 tri_areas = mesh.area_faces 524 vert_x = -tri_centers[:, 0] 525 vert_y = -tri_centers[:, 2] 526 vert_z = tri_centers[:, 1] 527 vert_x_trans = vert_x + pos[0] 528 vert_y_trans = vert_y + pos[1] 529 vert_z_trans = vert_z + pos[2] 530 tri_rows = np.floor((vert_y_trans - c_left_y) / step_y).astype(int) 531 tri_cols = np.floor((vert_x_trans - c_left_x) / step_x).astype(int) 532 tri_layers = np.floor(vert_z_trans / step_z).astype(int) 533 unique_rows = set(tri_rows) 534 unique_cols = set(tri_cols) 535 unique_layers = set(tri_layers) 536 for row in unique_rows: 537 for col in unique_cols: 538 for layer in unique_layers: 539 query = (tri_rows == row) & (tri_cols == col) & (tri_layers == layer) 540 area_values = tri_areas[query] 541 if len(area_values) > 0: 542 if 0 <= row < rows and 0 <= col < cols and 0 <= layer < layers: 543 lai3d[row, col, layer] += sum(area_values) 544 if out_path == "": 545 out_path = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Results", "LAI.txt") 546 print("Save file to:", out_path) 547 f = open(out_path, "w") 548 f.write("Scene Size [X, Y, Z]: %.4f, %.4f, %.4f\n" % (compute_width, compute_height, depth)) 549 f.write("LAI Dimension [Rows, Cols, layers]: %d, %d, %d\n" % (rows, cols, layers)) 550 for d in range(0, layers): 551 layer_lai = lai3d[:, :, d] 552 for row in range(rows): 553 for col in range(cols): 554 f.write("%.6f " % (layer_lai[row, col] / (step_x * step_y))) 555 f.write("\n") 556 f.write("\n") 557 f.close()
def
compute_chm(self, out_path, resolution, ignore_list=None):
559 def compute_chm(self, out_path, resolution, ignore_list=None): 560 if ignore_list is None: 561 ignore_list = [] 562 parameter_dir = os.path.join(self.land_scape.get_sim().get_sim_dir(), "Parameters") 563 objects = self.land_scape.get_objects() 564 instances = self.land_scape.get_instances() 565 tot_cols = int(np.ceil(self.width_extent / resolution)) 566 tot_rows = int(np.ceil(self.height_extent / resolution)) 567 dsm = np.zeros((tot_rows, tot_cols)) 568 obj_bounds = {} 569 obj_meshes = {} 570 for scene_obj in objects: 571 obj_bounding_box = LSBoundingBox() 572 meshes = [] 573 for obj_component in objects[scene_obj]: 574 if obj_component not in ignore_list: 575 obj_component_path = os.path.join(parameter_dir, obj_component) 576 obj_comp_mesh = trimesh.load_mesh(obj_component_path) 577 bbox = obj_comp_mesh.bounding_box 578 comp_bounding_box = LSBoundingBox() 579 comp_bounding_box.init_from_bounds(bbox.bounds) 580 obj_bounding_box.add_child(comp_bounding_box) 581 meshes.append(obj_comp_mesh) 582 obj_bounds[scene_obj] = obj_bounding_box 583 obj_meshes[scene_obj] = meshes 584 585 for instance in instances: 586 for pos in instances[instance]: 587 obj_bound = obj_bounds[instance] 588 obj_mesh = obj_meshes[instance] 589 if len(pos) == 10: # scale 590 x_scale_size = pos[7] 591 y_scale_size = pos[9] # vertical direction 592 z_scale_size = pos[8] 593 x_scale_factor = x_scale_size / obj_bound.get_x_extents() 594 y_scale_factor = y_scale_size / obj_bound.get_y_extents() 595 z_scale_factor = z_scale_size / obj_bound.get_z_extents() 596 obj_mesh = [] 597 for mesh in obj_meshes[instance]: 598 new_mesh = mesh.copy() 599 new_mesh.vertices[:, 0] *= x_scale_factor 600 new_mesh.vertices[:, 1] *= y_scale_factor 601 new_mesh.vertices[:, 2] *= z_scale_factor 602 obj_mesh.append(new_mesh) 603 if pos[3] != 0: 604 obj_mesh1 = [] 605 for mesh in obj_mesh: 606 new_mesh = mesh.copy() 607 angle = np.radians(pos[3]) # 旋转角度（以弧度为单位） 608 axis = [0, 1, 0] 609 if len(pos) >= 7: 610 axis = [-pos[4], pos[6], -pos[5]] 611 rotation_matrix = trimesh.transformations.rotation_matrix(angle, axis) 612 rotated_mesh = new_mesh.apply_transform(rotation_matrix) 613 obj_mesh1.append(rotated_mesh) 614 obj_mesh = obj_mesh1 615 for mesh in obj_mesh: 616 vertices = mesh.vertices 617 vert_x = -vertices[:, 0] 618 vert_y = -vertices[:, 2] 619 vert_z = vertices[:, 1] 620 vert_x_trans = vert_x + pos[0] 621 vert_y_trans = vert_y + pos[1] 622 vert_z_trans = vert_z + pos[2] 623 rows = (vert_y_trans / resolution).astype(int) 624 cols = (vert_x_trans / resolution).astype(int) 625 unique_rows = set(rows) 626 unique_cols = set(cols) 627 for row in unique_rows: 628 for col in unique_cols: 629 query = (rows == row) & (cols == col) 630 z_values = vert_z_trans[query] 631 if len(z_values) > 0: 632 max_value = z_values.max() 633 if 0 <= row < tot_rows and 0 <= col < tot_cols and max_value > dsm[row, col]: 634 dsm[row, col] = max_value 635 if out_path.endswith(".tif"): 636 out_path = out_path[0:-4] 637 RasterHelper.saveToHdr_no_transform(dsm, out_path, [], "GTiff")