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]]}
land_scape
width_extent
height_extent
tot_area
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")