From 1be4d5acc89fe255660242a6ed05cc0490f8d5a4 Mon Sep 17 00:00:00 2001 From: pilar Date: Tue, 9 Jun 2020 11:34:12 -0400 Subject: [PATCH] Geometrical methods for dividing buildings by plane 2 --- city_model_structure/plying_with_geometry.py | 85 +++---------------- helpers/geometry.py | 84 ++++++++++++++++++ .../us_base_physics_parameters.py | 2 - 3 files changed, 97 insertions(+), 74 deletions(-) diff --git a/city_model_structure/plying_with_geometry.py b/city_model_structure/plying_with_geometry.py index bd8e93ca..b3d34d0e 100644 --- a/city_model_structure/plying_with_geometry.py +++ b/city_model_structure/plying_with_geometry.py @@ -1,81 +1,22 @@ import trimesh -from trimesh import intersections as inter -import open3d as o3d +from helpers.geometry import Geometry as Geo +import stl import numpy as np -def extract_points(segment_list): - point_list = np.asarray(segment_list[0]) - for segment in segment_list: - found = False - for new_point in segment: - for point in point_list: - comparison = new_point == point - if comparison.all(): - found = True - break - if not found: - point_list = np.concatenate((point_list, [new_point])) - return point_list - - -def point_cloud_to_mesh(point_list, normal_list): - pcd = o3d.geometry.PointCloud() - pcd.points = o3d.utility.Vector3dVector(point_list) - pcd.normals = o3d.utility.Vector3dVector(normal_list) - distances = pcd.compute_nearest_neighbor_distance() - avg_dist = np.mean(distances) - radius = 3 * avg_dist - bpa_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting( - pcd, - o3d.utility.DoubleVector([radius, radius * 2])) - mesh_result = trimesh.Trimesh(vertices=np.asarray(bpa_mesh.vertices), faces=np.asarray(bpa_mesh.triangles)) - return mesh_result - - -def merge_meshes(mesh1, mesh2): - v_1 = mesh1.vertices - f_1 = mesh1.faces - v_2 = mesh2.vertices - f_2 = mesh2.faces - length = len(v_1) - v_merge = np.concatenate((v_1, v_2)) - f_merge = np.asarray(f_1) - - for item in f_2: - point1 = item.item(0) + length - point2 = item.item(1) + length - point3 = item.item(2) + length - surface = np.asarray([point1, point2, point3]) - f_merge = np.concatenate((f_merge, [surface])) - - mesh_merge = trimesh.Trimesh(vertices=v_merge, faces=f_merge) - - return mesh_merge - - -# example +# example triangle vertices = [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]] -faces = [[0, 1, 2], [1, 3, 2], [1, 5, 3], [5, 3, 7], [0, 4, 1], [4, 5, 1], [4, 6, 5], [5, 6, 7], [0, 2, 4], [2, 6, 4], - [3, 7, 6], [2, 3, 6]] +#faces = [[0, 1, 2], [1, 3, 2], [1, 5, 3], [5, 7, 3], [0, 4, 1], [4, 5, 1], [4, 6, 5], [5, 6, 7], [0, 2, 4], [2, 6, 4], +# [3, 7, 6], [2, 3, 6]] + +# example rectangle +faces = [[0, 1, 3, 2], [1, 5, 7, 3], [0, 4, 5, 1], [4, 6, 7, 5], [0, 2, 6, 4], [3, 7, 6, 2]] mesh = trimesh.Trimesh(vertices=vertices, faces=faces) + normal_plane = [0, 0, 1] -point_plane = [0, 0, 0.5] -normal_plane_opp = [0, 0, -1] +point_plane = [0, 0, 2] -# Step 1: division and extraction of segments -mesh_1 = inter.slice_mesh_plane(mesh, normal_plane, point_plane) -mesh_1_segments = inter.mesh_plane(mesh, normal_plane, point_plane) +mesh_result = Geo.divide_mesh_by_plane(mesh, normal_plane, point_plane) -# Step 2: create mesh from point cloud -points = extract_points(mesh_1_segments) -normals = np.empty((len(points), 3)) -for i in range(0, len(normals)): - normals[i] = normal_plane_opp -mesh_2 = point_cloud_to_mesh(points, normals) - -# Step 3: merge both meshes -mesh_final = merge_meshes(mesh_1, mesh_2) - - -print(mesh_final.vertices, mesh_final.faces) +for mesh in mesh_result: + print(len(mesh.vertices), len(mesh.faces)) diff --git a/helpers/geometry.py b/helpers/geometry.py index c99dabff..f94c7508 100644 --- a/helpers/geometry.py +++ b/helpers/geometry.py @@ -1,5 +1,8 @@ import math import numpy as np +from trimesh import Trimesh +from trimesh import intersections +import open3d as o3d class Geometry: @@ -46,3 +49,84 @@ class Geometry: if remove_last: points = np.delete(points, rows - 1, 0) return points + + @staticmethod + def _segment_list_to_point_cloud(segment_list): + point_list = np.asarray(segment_list[0]) + for segment in segment_list: + for new_point in segment: + found = False + for point in point_list: + same_point = np.allclose(new_point, point) + if same_point: + found = True + break + if not found: + point_list = np.concatenate((point_list, [new_point])) + return point_list + + @staticmethod + def _point_cloud_to_mesh(point_list, normal_list): + # Return a mesh composed only by triangles + pcd = o3d.geometry.PointCloud() + pcd.points = o3d.utility.Vector3dVector(point_list) + pcd.normals = o3d.utility.Vector3dVector(normal_list) + distances = pcd.compute_nearest_neighbor_distance() + avg_dist = np.mean(distances) + radius = 3 * avg_dist + bpa_mesh = o3d.geometry.TriangleMesh().create_from_point_cloud_ball_pivoting( + pcd, o3d.utility.DoubleVector([radius, radius * 2])) + mesh_result = Trimesh(vertices=np.asarray(bpa_mesh.vertices), faces=np.asarray(bpa_mesh.triangles)) + return mesh_result + + @staticmethod + def _merge_meshes(mesh1, mesh2): + v_1 = mesh1.vertices + f_1 = mesh1.faces + v_2 = mesh2.vertices + f_2 = mesh2.faces + length = len(v_1) + v_merge = np.concatenate((v_1, v_2)) + f_merge = np.asarray(f_1) + + for item in f_2: + point1 = item.item(0) + length + point2 = item.item(1) + length + point3 = item.item(2) + length + surface = np.asarray([point1, point2, point3]) + f_merge = np.concatenate((f_merge, [surface])) + + mesh_merge = Trimesh(vertices=v_merge, faces=f_merge) + + return mesh_merge + + @staticmethod + def divide_mesh_by_plane(mesh, normal_plane, point_plane): + # The first mesh returns the positive side of the plane and the second the negative side. + # If the plane does not divide the mesh (i.e. it does not touch it or it is coplanar with one or more faces), + # then it returns only the original mesh. + normal_plane_opp = [None] * len(normal_plane) + for i in range(0, len(normal_plane)): + normal_plane_opp[i] = - normal_plane[i] + + normal = [normal_plane, normal_plane_opp] + normal_opp = [normal_plane_opp, normal_plane] + mesh_final = [] + for i in range(0, 2): + mesh_1 = intersections.slice_mesh_plane(mesh, normal[i], point_plane) + mesh_1_segments = intersections.mesh_plane(mesh, normal[i], point_plane) + boo = mesh.difference(mesh_1, engine='blender') + print(boo) + quit() + if len(mesh_1_segments) <= 0 or len(mesh_1.faces) == len(mesh.faces): + mesh_final.append(mesh) + break + else: + points = Geometry._segment_list_to_point_cloud(mesh_1_segments) + points_normals = [[None] * 3] * len(points) + for j in range(0, len(points_normals)): + points_normals[j] = normal_opp[i] + mesh_2 = Geometry._point_cloud_to_mesh(points, points_normals) + mesh_final.append(Geometry._merge_meshes(mesh_1, mesh_2)) + + return mesh_final diff --git a/physics/physics_feeders/us_base_physics_parameters.py b/physics/physics_feeders/us_base_physics_parameters.py index b261e4f4..8049a47f 100644 --- a/physics/physics_feeders/us_base_physics_parameters.py +++ b/physics/physics_feeders/us_base_physics_parameters.py @@ -69,8 +69,6 @@ class UsBasePhysicsParameters: material.no_mass = 'no_mass' in material_lib if 'thermal_resistance' in material_lib: material.thermal_resistance_m2k_w = material_lib['thermal_resistance']['#text'] - else: - material.thermal_resistance_m2k_w = None layer.material = material thermal_boundary.layers.append(layer) for opening in thermal_boundary.thermal_openings: