diff --git a/docs/examples/01_aperture_projection.py b/docs/examples/01_aperture_projection.py index ae2c812..9799548 100644 --- a/docs/examples/01_aperture_projection.py +++ b/docs/examples/01_aperture_projection.py @@ -28,7 +28,6 @@ import copy import numpy as np -import open3d as o3d # %% # Setting Farfield Resolution and Wavelength @@ -54,8 +53,7 @@ body, array, _ = data.exampleUAV(10e9) -# visualise UAV and Array -o3d.visualization.draw_geometries([body, array]) + # %% ## .. image:: ../_static/open3d_structure.png @@ -63,12 +61,9 @@ # crop the inner surface of the array trianglemesh (not strictly required, as the UAV main body provides blocking to # the hidden surfaces, but correctly an aperture will only have an outer face. surface_array = copy.deepcopy(array) -surface_array.triangles = o3d.utility.Vector3iVector( - np.asarray(array.triangles)[: len(array.triangles) // 2, :] -) -surface_array.triangle_normals = o3d.utility.Vector3dVector( - np.asarray(array.triangle_normals)[: len(array.triangle_normals) // 2, :] -) +surface_array.cells[0].data = np.asarray(array.cells[0].data)[: (array.cells[0].data).shape[0] // 2, :] + +surface_array.cell_data["normals"] = np.array(array.cell_data["normals"])[: (array.cells[0].data).shape[0] // 2] # %% # Structures @@ -104,7 +99,6 @@ # also as an open3d point cloud. This allows easy visualisation using :func:`open3d.visualization.draw_geometries`. # %% -o3d.visualization.draw_geometries([body, surface_array, pcd]) # %% # .. image:: ../_static/open3d_results_rendering.png diff --git a/docs/examples/02_coherently_polarised_array.py b/docs/examples/02_coherently_polarised_array.py index 3f9d45f..4a85de2 100644 --- a/docs/examples/02_coherently_polarised_array.py +++ b/docs/examples/02_coherently_polarised_array.py @@ -13,8 +13,7 @@ import copy import numpy as np -import open3d as o3d - +import meshio # %% # Setting Farfield Resolution and Wavelength # ------------------------------------------- @@ -39,29 +38,17 @@ body, array, source_coords = data.exampleUAV(10e9) -# %% -# Visualise the Resultant UAV and Array -# --------------------------------------- -# :func:`open3d.visualization.draw_geometries` can be used to visualise the open3d data -# structures :class:`open3d.geometry.PointCloud` and :class:`open3d.geometry.PointCloud` -mesh_frame = o3d.geometry.TriangleMesh.create_coordinate_frame( - size=0.5, origin=[0, 0, 0] -) -o3d.visualization.draw_geometries([body, array, source_coords, mesh_frame]) # %% -# .. image:: ../_static/UAVArraywithPoints.png +## .. image:: ../_static/open3d_structure.png # crop the inner surface of the array trianglemesh (not strictly required, as the UAV main body provides blocking to # the hidden surfaces, but correctly an aperture will only have an outer face. surface_array = copy.deepcopy(array) -surface_array.triangles = o3d.utility.Vector3iVector( - np.asarray(array.triangles)[: len(array.triangles) // 2, :] -) -surface_array.triangle_normals = o3d.utility.Vector3dVector( - np.asarray(array.triangle_normals)[: len(array.triangle_normals) // 2, :] -) +surface_array.cells[0].data = np.asarray(array.cells[0].data)[: (array.cells[0].data).shape[0] // 2, :] + +surface_array.cell_data["normals"] = np.array(array.cell_data["normals"])[: (array.cells[0].data).shape[0] // 2] from lyceanem.base_classes import structures @@ -69,11 +56,14 @@ from lyceanem.models.frequency_domain import calculate_farfield -from lyceanem.geometry.targets import source_cloud_from_shape -source_points, _ = source_cloud_from_shape(surface_array, wavelength * 0.5) -o3d.visualization.draw_geometries([body, array, source_points]) + +import pyvista as pv + + +source_points = surface_array.points + # %% # .. image:: ../_static/sourcecloudfromshapeuav.png diff --git a/docs/examples/03_frequency_domain_channel_modelling.py b/docs/examples/03_frequency_domain_channel_modelling.py index 771c7b4..47ab392 100644 --- a/docs/examples/03_frequency_domain_channel_modelling.py +++ b/docs/examples/03_frequency_domain_channel_modelling.py @@ -36,12 +36,7 @@ transmit_horn_structure, transmitting_antenna_surface_coords = TL.meshedHorn( 58e-3, 58e-3, 128e-3, 2e-3, 0.21, mesh_resolution ) -print("hi",transmit_horn_structure) -print("hii",transmit_horn_structure.points.shape) -print("hiiii",transmit_horn_structure.point_data) -print("len", len(transmit_horn_structure.cells)) -print("cells", transmit_horn_structure.cells[0]) -print("data", transmitting_antenna_surface_coords.cell_data) + receive_horn_structure, receiving_antenna_surface_coords = TL.meshedHorn( 58e-3, 58e-3, 128e-3, 2e-3, 0.21, mesh_resolution ) diff --git a/lyceanem/models/frequency_domain.py b/lyceanem/models/frequency_domain.py index c3b0bd6..920b218 100644 --- a/lyceanem/models/frequency_domain.py +++ b/lyceanem/models/frequency_domain.py @@ -8,6 +8,7 @@ from ..geometry import targets as TL from ..raycasting import rayfunctions as RF from ..utility.math_functions import calc_dv_norm +import meshio def aperture_projection( @@ -53,7 +54,7 @@ def aperture_projection( ) triangle_centroids, _ = GF.tri_centroids(aperture) triangle_areas = GF.tri_areas(aperture) - triangle_normals = np.asarray(aperture.triangle_normals) + triangle_normals = np.asarray(aperture.cell_data["normals"]) visible_patterns, pcd = RF.visiblespace( triangle_centroids, triangle_normals, @@ -142,17 +143,17 @@ def calculate_farfield( sinks, np.zeros((len(sinks), 3), dtype=np.float32), sink_normals, lengths ) sink_cloud = RF.points2pointcloud(sinks) - sink_cloud.normals = o3d.utility.Vector3dVector(sink_normals) + sink_cloud.point_data["normals"] = sink_normals num_sources = len(np.asarray(aperture_coords.points)) num_sinks = len(np.asarray(sink_cloud.points)) environment_triangles=GF.mesh_conversion(antenna_solid) if project_vectors: conformal_E_vectors = EM.calculate_conformalVectors( - desired_E_axis, np.asarray(aperture_coords.normals), antenna_axes + desired_E_axis, np.asarray(aperture_coords.point_data["normals"]), antenna_axes ) else: - if desired_E_axis.shape[0]==np.asarray(aperture_coords.normals).shape[0]: + if desired_E_axis.shape[0]==np.asarray(aperture_coords.point_data["normals"]).shape[0]: conformal_E_vectors=copy.deepcopy(desired_E_axis) else: conformal_E_vectors = np.repeat( @@ -161,15 +162,15 @@ def calculate_farfield( if scattering == 0: # only use the aperture point cloud, no scattering required. - scatter_points = o3d.geometry.PointCloud() + scatter_points = meshio.Mesh(points= np.empty((0, 3)), cells=[]) unified_model = np.append( np.asarray(aperture_coords.points).astype(np.float32), np.asarray(sink_cloud.points).astype(np.float32), axis=0, ) unified_normals = np.append( - np.asarray(aperture_coords.normals).astype(np.float32), - np.asarray(sink_cloud.normals).astype(np.float32), + np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), + np.asarray(sink_cloud.point_data["normals"]).astype(np.float32), axis=0, ) unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) @@ -198,15 +199,10 @@ def calculate_farfield( point_informationv2[0:num_sources]["pz"] = np.asarray( aperture_coords.points ).astype(np.float32)[:, 2] - point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 2] + normals = np.asarray(aperture_coords.point_data["normals"]) + point_informationv2[0:num_sources]["nx"] = normals[:, 0] + point_informationv2[0:num_sources]["ny"] = normals[:, 1] + point_informationv2[0:num_sources]["nz"] = normals[:, 2] # set position and velocity of sinks point_informationv2[num_sources : (num_sources + num_sinks)]["px"] = sinks[:, 0] point_informationv2[num_sources : (num_sources + num_sinks)]["py"] = sinks[:, 1] @@ -246,11 +242,11 @@ def calculate_farfield( ) unified_normals = np.append( np.append( - np.asarray(aperture_coords.normals).astype(np.float32), - np.asarray(sink_cloud.normals).astype(np.float32), + np.asarray(aperture_coords.point_data["normals"]).astype(np.float32), + np.asarray(sink_cloud.point_data["normals"]).astype(np.float32), axis=0, ), - np.asarray(scatter_points.normals).astype(np.float32), + np.asarray(scatter_points.point_data["normals"]).astype(np.float32), axis=0, ) unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) @@ -282,13 +278,13 @@ def calculate_farfield( aperture_coords.points ).astype(np.float32)[:, 2] point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.normals + aperture_coords.point_data["normals"] ).astype(np.float32)[:, 0] point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.normals + aperture_coords.point_data["normals"] ).astype(np.float32)[:, 1] point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.normals + aperture_coords.point_data["normals"] ).astype(np.float32)[:, 2] # point_informationv2[0:num_sources]['ex']=unified_weights[0:num_sources,0] # point_informationv2[0:num_sources]['ey']=unified_weights[0:num_sources,1] @@ -318,15 +314,10 @@ def calculate_farfield( point_informationv2[(num_sources + num_sinks) :]["pz"] = np.asarray( scatter_points.points ).astype(np.float32)[:, 2] - point_informationv2[(num_sources + num_sinks) :]["nx"] = np.asarray( - scatter_points.normals - ).astype(np.float32)[:, 0] - point_informationv2[(num_sources + num_sinks) :]["ny"] = np.asarray( - scatter_points.normals - ).astype(np.float32)[:, 1] - point_informationv2[(num_sources + num_sinks) :]["nz"] = np.asarray( - scatter_points.normals - ).astype(np.float32)[:, 2] + scatter_normals = np.asarray(scatter_points.point_data["normals"]) + point_informationv2[(num_sources + num_sinks) :]["nx"] = scatter_normals[:, 0] + point_informationv2[(num_sources + num_sinks) :]["ny"] = scatter_normals[:, 1] + point_informationv2[(num_sources + num_sinks) :]["nz"] = scatter_normals[:, 2] point_informationv2[:]["ex"] = unified_weights[:, 0] point_informationv2[:]["ey"] = unified_weights[:, 1] point_informationv2[:]["ez"] = unified_weights[:, 2] @@ -832,246 +823,3 @@ def calculate_scattering( return Ex, Ey, Ez -def calculate_scattering_isotropic( - aperture_coords, - sink_coords, - antenna_solid, - scatter_points=None, - wavelength=1.0, - scattering=0, - elements=False, - mesh_resolution=0.5, -): - """ - Based upon the aperture coordinates and solids, predict the farfield for the antenna. - """ - num_sources = len(np.asarray(aperture_coords.points)) - num_sinks = len(np.asarray(sink_coords.points)) - - if scattering == 0: - # only use the aperture point cloud, no scattering required. - scatter_points = o3d.geometry.PointCloud() - - unified_model = np.append( - np.asarray(aperture_coords.points).astype(np.float32), - np.asarray(sink_coords.points).astype(np.float32), - axis=0, - ) - unified_normals = np.append( - np.asarray(aperture_coords.normals).astype(np.float32), - np.asarray(sink_coords.normals).astype(np.float32), - axis=0, - ) - unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) - unified_weights[0:num_sources, :] = ( - 1.0 / num_sources - ) # set total amplitude to 1 for the aperture - unified_weights[num_sources : num_sources + num_sinks, :] = ( - 1 / num_sinks - ) # set total amplitude to 1 for the aperture - point_informationv2 = np.empty((len(unified_model)), dtype=scattering_t) - # set all sources as electric current sources, and permittivity and permeability as free space - point_informationv2[:]["Electric"] = True - point_informationv2[:]["permittivity"] = 8.8541878176e-12 - point_informationv2[:]["permeability"] = 1.25663706212e-6 - # set position, velocity, normal, and weight of sources - point_informationv2[0:num_sources]["px"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["py"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["pz"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 2] - point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 2] - # set position and velocity of sinks - point_informationv2[num_sources : (num_sources + num_sinks)]["px"] = np.asarray( - sink_coords.points - ).astype(np.float32)[:, 0] - point_informationv2[num_sources : (num_sources + num_sinks)]["py"] = np.asarray( - sink_coords.points - ).astype(np.float32)[:, 1] - point_informationv2[num_sources : (num_sources + num_sinks)]["pz"] = np.asarray( - sink_coords.points - ).astype(np.float32)[:, 2] - point_informationv2[num_sources : (num_sources + num_sinks)]["nx"] = np.asarray( - sink_coords.normals - ).astype(np.float32)[:, 0] - point_informationv2[num_sources : (num_sources + num_sinks)]["ny"] = np.asarray( - sink_coords.normals - ).astype(np.float32)[:, 1] - point_informationv2[num_sources : (num_sources + num_sinks)]["nz"] = np.asarray( - sink_coords.normals - ).astype(np.float32)[:, 2] - - point_informationv2[:]["ex"] = unified_weights[:, 0] - point_informationv2[:]["ey"] = unified_weights[:, 1] - point_informationv2[:]["ez"] = unified_weights[:, 2] - scatter_mask = np.zeros((point_informationv2.shape[0]), dtype=np.int32) - scatter_mask[0:num_sources] = 0 - scatter_mask[(num_sources + num_sinks) :] = 0 - - else: - # create scatter points on antenna solids based upon a half wavelength square - if scatter_points is None: - scatter_points, areas = TL.source_cloud_from_shape( - antenna_solid, 1e-6, (wavelength * mesh_resolution) ** 2 - ) - - unified_model = np.append( - np.append( - np.asarray(aperture_coords.points).astype(np.float32), - np.asarray(sink_coords.points).astype(np.float32), - axis=0, - ), - np.asarray(scatter_points.points).astype(np.float32), - axis=0, - ) - unified_normals = np.append( - np.append( - np.asarray(aperture_coords.normals).astype(np.float32), - np.asarray(sink_coords.normals).astype(np.float32), - axis=0, - ), - np.asarray(scatter_points.normals).astype(np.float32), - axis=0, - ) - unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) - unified_weights[0:num_sources, :] = ( - 1.0 / num_sources - ) # set total amplitude to 1 for the aperture - unified_weights[num_sources : num_sources + num_sinks, :] = ( - 1 / num_sinks - ) # set total amplitude to 1 for the aperture - unified_weights[ - num_sources + num_sinks :, : - ] = 1.0 # / len(np.asarray(scatter_points.points)) # set total amplitude to 1 for the aperture - point_informationv2 = np.empty((len(unified_model)), dtype=scattering_t) - # set all sources as magnetic current sources, and permittivity and permeability as free space - point_informationv2[:]["Electric"] = True - point_informationv2[:]["permittivity"] = 8.8541878176e-12 - point_informationv2[:]["permeability"] = 1.25663706212e-6 - # set position, velocity, normal, and weight of sources - point_informationv2[0:num_sources]["px"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["py"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["pz"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 2] - point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.normals - ).astype(np.float32)[:, 2] - # point_informationv2[0:num_sources]['ex']=unified_weights[0:num_sources,0] - # point_informationv2[0:num_sources]['ey']=unified_weights[0:num_sources,1] - # point_informationv2[0:num_sources]['ez']=unified_weights[0:num_sources,2] - # set position and velocity of sinks - point_informationv2[num_sources : (num_sources + num_sinks)]["px"] = np.asarray( - sink_coords.points - ).astype(np.float32)[:, 0] - point_informationv2[num_sources : (num_sources + num_sinks)]["py"] = np.asarray( - sink_coords.points - ).astype(np.float32)[:, 1] - point_informationv2[num_sources : (num_sources + num_sinks)]["pz"] = np.asarray( - sink_coords.points - ).astype(np.float32)[:, 2] - # point_informationv2[num_sources:(num_sources+num_sinks)]['vx']=0.0 - # point_informationv2[num_sources:(num_sources+num_sinks)]['vy']=0.0 - # point_informationv2[num_sources:(num_sources+num_sinks)]['vz']=0.0 - point_informationv2[num_sources : (num_sources + num_sinks)]["nx"] = np.asarray( - sink_coords.normals - ).astype(np.float32)[:, 0] - point_informationv2[num_sources : (num_sources + num_sinks)]["ny"] = np.asarray( - sink_coords.normals - ).astype(np.float32)[:, 1] - point_informationv2[num_sources : (num_sources + num_sinks)]["nz"] = np.asarray( - sink_coords.normals - ).astype(np.float32)[:, 2] - point_informationv2[(num_sources + num_sinks) :]["px"] = np.asarray( - scatter_points.points - ).astype(np.float32)[:, 0] - point_informationv2[(num_sources + num_sinks) :]["py"] = np.asarray( - scatter_points.points - ).astype(np.float32)[:, 1] - point_informationv2[(num_sources + num_sinks) :]["pz"] = np.asarray( - scatter_points.points - ).astype(np.float32)[:, 2] - point_informationv2[(num_sources + num_sinks) :]["nx"] = np.asarray( - scatter_points.normals - ).astype(np.float32)[:, 0] - point_informationv2[(num_sources + num_sinks) :]["ny"] = np.asarray( - scatter_points.normals - ).astype(np.float32)[:, 1] - point_informationv2[(num_sources + num_sinks) :]["nz"] = np.asarray( - scatter_points.normals - ).astype(np.float32)[:, 2] - point_informationv2[:]["ex"] = unified_weights[:, 0] - point_informationv2[:]["ey"] = unified_weights[:, 1] - point_informationv2[:]["ez"] = unified_weights[:, 2] - scatter_mask = np.zeros((point_informationv2.shape[0]), dtype=np.int32) - scatter_mask[0:num_sources] = 0 - scatter_mask[(num_sources + num_sinks) :] = scattering - - # full_index, initial_index = RF.integratedraycastersetup(num_sources, - # num_sinks, - # point_informationv2, - # RF.convertTriangles(antenna_solid), - # scatter_mask) - full_index, rays = RF.workchunkingv2( - np.asarray(aperture_coords.points).astype(np.float32), - np.asarray(sink_coords.points).astype(np.float32), - np.asarray(scatter_points.points).astype(np.float32), - antenna_solid.triangles_base_raycaster(), - scattering + 1, - ) - - if not elements: - # create efiles for model - scatter_map = EM.IsoGPUFreqDomain( - num_sources, num_sinks, full_index, point_informationv2, wavelength - ) - final_scattering_map = np.sum( - np.dot(np.ones((num_sources)), scatter_map[:, :, 0]) - ) - - # convert to etheta,ephi - - else: - # create efiles for model - - final_scattering_map = np.zeros((num_sources, num_sinks), dtype=np.complex64) - - for element in range(num_sources): - point_informationv2[0:num_sources]["ex"] = 0.0 - point_informationv2[0:num_sources]["ey"] = 0.0 - point_informationv2[0:num_sources]["ez"] = 0.0 - point_informationv2[element]["ex"] = 1.0 / num_sources - point_informationv2[element]["ey"] = 1.0 / num_sources - point_informationv2[element]["ez"] = 1.0 / num_sources - unified_weights[0:num_sources, :] = 0.0 - unified_weights[element, :] = 1.0 / num_sources - scatter_map = EM.EMGPUFreqDomain( - num_sources, num_sinks, full_index, point_informationv2, wavelength - ) - final_scattering_map[element, :] = np.dot( - np.ones((num_sources)), scatter_map[:, :, 0] - ) - - return final_scattering_map diff --git a/lyceanem/raycasting/rayfunctions.py b/lyceanem/raycasting/rayfunctions.py index 7262b38..3a9596a 100644 --- a/lyceanem/raycasting/rayfunctions.py +++ b/lyceanem/raycasting/rayfunctions.py @@ -1841,6 +1841,10 @@ def charge_rays_environment1Dv2(sources, sinks, environment_points, point_indexi ray payload to be sent to GPU """ + print("sources shape", sources.shape) + print("sinks shape", sinks.shape) + print("environment_points shape", environment_points.shape) + unified_model = np.append( np.append(sources, sinks, axis=0), environment_points, axis=0 ) diff --git a/lyceanem/tests/reflectordata.py b/lyceanem/tests/reflectordata.py index fc05a7e..628c377 100644 --- a/lyceanem/tests/reflectordata.py +++ b/lyceanem/tests/reflectordata.py @@ -5,9 +5,10 @@ import math import numpy as np -import open3d as o3d +import meshio import scipy.io as io from importlib_resources import files +import pyvista as pv import lyceanem.geometry.geometryfunctions as GF import lyceanem.geometry.targets as tl @@ -37,36 +38,40 @@ def exampleUAV(frequency): bodystream = files(lyceanem.tests.data).joinpath("UAV.stl") arraystream = files(lyceanem.tests.data).joinpath("UAVarray.stl") - body = o3d.io.read_triangle_mesh(str(bodystream)) - array = o3d.io.read_triangle_mesh(str(arraystream)) + body = meshio.read(str(bodystream)) + array = meshio.read(str(arraystream)) rotation_vector1 = np.asarray([0.0, np.deg2rad(90), 0.0]) rotation_vector2 = np.asarray([np.deg2rad(90), 0.0, 0.0]) - body = GF.open3drotate( - body, o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1) - ) - body = GF.open3drotate( - body, o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector2) - ) - array = GF.open3drotate( - array, o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1) - ) - array = GF.open3drotate( - array, o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector2) + body = GF.mesh_rotate(body, rotation_vector1) + body = GF.mesh_rotate(body, rotation_vector2) + array = GF.mesh_rotate(array, rotation_vector1) + array = GF.mesh_rotate(array, rotation_vector2) + + + body = GF.translate_mesh(body, np.array([0.25, 0, 0])+np.array([-0.18, 0, 0.0125])) + array = GF.translate_mesh(array, np.array([0.25, 0, 0])+np.array([-0.18, 0, 0.0125]) ) - body.translate(np.asarray([0.25, 0, 0]), relative=True) - array.translate(np.asarray([0.25, 0, 0]), relative=True) - array.translate(np.array([-0.18, 0, 0.0125]), relative=True) - body.translate(np.array([-0.18, 0, 0.0125]), relative=True) - array.compute_vertex_normals() - array.paint_uniform_color(np.array([0, 1.0, 1.0])) - body.compute_vertex_normals() - body.paint_uniform_color(np.array([0, 0.259, 0.145])) + + def structure_cells(array): + ## add collumn of 3s to beggining of each row + array = np.append(np.ones((array.shape[0], 1), dtype=np.int32) * 3, array, axis=1) + return array + pyvista_array = pv.PolyData(array.points, structure_cells(array.cells[0].data)) + pyvista_body = pv.PolyData(body.points, structure_cells(body.cells[0].data)) + pyvista_array.compute_normals(inplace=True) + pyvista_body.compute_normals(inplace=True) + + array.point_data['normals'] = pyvista_array.point_normals + body.point_data['normals'] = pyvista_body.point_normals + array.cell_data['normals'] = pyvista_array.cell_normals + body.cell_data['normals'] = pyvista_body.cell_normals + wavelength = 3e8 / frequency mesh_sep = wavelength * 0.5 # define UAV array, a line, a circle, then a line with reference to the array solid - array_vertices = np.asarray(array.vertices) + array_vertices = np.asarray(array.points) array_height = np.max(array_vertices[:, 2]) - np.min(array_vertices[:, 2]) top_index = ( array_vertices[:, 2] >= np.max(array_vertices[:, 2]) - 1e-4 @@ -190,215 +195,8 @@ def exampleUAV(frequency): axis=0, ) + np.array([0.1025, 0, -0.025]) source_pcd = RF.points2pointcloud(total_array) - source_pcd.normals = o3d.utility.Vector3dVector(total_array_normals) - source_pcd.translate(np.array([-0.18, 0, 0.0125]), relative=True) - - return body, array, source_pcd - - -def prototypescan(mesh_size, wavelength=3e8 / 24e9, importpcd=True): - if importpcd: - if mesh_size < (wavelength * 0.5): - # stream = pkg_resources.resource_stream(__name__, - # 'data/prototypescanpoints24GHztenthwavelength.ply') - stream = files("lyceanem.tests.data").joinpath( - "prototypescanpoints24GHztenthwavelength.ply" - ) - prototype_points = o3d.io.read_point_cloud(str(stream)) - else: - # stream = pkg_resources.resource_stream(__name__, - # 'data/prototypescanpoints24GHztenthwavelength.ply') - stream = files(lyceanem.tests.data).joinpath( - "prototypescanpoints24GHztenthwavelength.ply" - ) - prototype_points = o3d.io.read_point_cloud(str(stream)) - prototype_mesh, _ = tl.meshedReflector(0.3, 0.3, 6e-3, mesh_size) - else: - # stream = pkg_resources.resource_stream(__name__, - # 'data/plate_copper.mat') - stream = files("lyceanem.tests.data").joinpath("plate_copper.mat") - matdata = io.loadmat(str(stream))["Flat_prototype"] - copper_reference_point = ( - np.asarray( - [[np.min(matdata[:, 0]), np.min(matdata[:, 1]), np.min(matdata[:, 2])]] - ) - / 1000 - ) - prototype = o3d.geometry.PointCloud() - # convert from mm to si units - prototype.points = o3d.utility.Vector3dVector(matdata / 1000) - prototype.estimate_normals() - prototype.orient_normals_to_align_with_direction( - orientation_reference=[0.0, 0.0, 1.0] - ) - downsampled_prototype = prototype.voxel_down_sample(voxel_size=mesh_size * 1.1) - downsampled_prototype.translate(-copper_reference_point.ravel(), relative=True) - prototype.translate(-copper_reference_point.ravel(), relative=True) - downsampled_prototype.translate([-0.15, -0.15, 0], relative=True) - prototype.translate([-0.15, -0.15, 0], relative=True) - reflectorplate, scatter_points = tl.meshedReflector( - 0.3, 0.3, 6e-3, mesh_size, sides="all" - ) - temp_points = np.asarray(scatter_points.points) - temp_points_back = temp_points[temp_points[:, 2] < -1e-3, :] - temp_points_top = temp_points[temp_points[:, 1] > 1.49e-1, :] - top_normals = np.zeros((temp_points_top.shape[0], 3)) - top_normals[:, 1] = 1 - toppoints = o3d.geometry.PointCloud() - toppoints.points = o3d.utility.Vector3dVector(temp_points_top) - toppoints.normals = o3d.utility.Vector3dVector(top_normals) - sideapoints = copy.deepcopy(toppoints) - rotation_vector1 = np.radians(np.array([0, 0, 90])) - sideapoints = GF.open3drotate( - sideapoints, - o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1), - ) - bottompoints = copy.deepcopy(sideapoints) - bottompoints = GF.open3drotate( - bottompoints, - o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1), - ) - sidebpoints = copy.deepcopy(bottompoints) - sidebpoints = GF.open3drotate( - sidebpoints, - o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1), - ) - - back_normals = np.zeros((temp_points_back.shape[0], 3)) - back_normals[:, 2] = -1 - backpoints = o3d.geometry.PointCloud() - backpoints.points = o3d.utility.Vector3dVector(temp_points_back) - backpoints.normals = o3d.utility.Vector3dVector(back_normals) - # for now only interested in points on the front surface - temp = np.asarray(downsampled_prototype.points) - interest_reference = np.min(temp[:, 2]) - test_prototype = ( - downsampled_prototype - + backpoints - + toppoints - + bottompoints - + sideapoints - + sidebpoints - ) - ( - prototype_mesh, - output, - ) = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( - test_prototype, linear_fit=True - ) - prototype_mesh.compute_triangle_normals() - prototype_points = o3d.geometry.PointCloud() - vertices = np.asarray(prototype_mesh.vertices) - normals = np.asarray(prototype_mesh.vertex_normals) - prototype_points.points = o3d.utility.Vector3dVector( - vertices[vertices[:, 2] >= interest_reference, :] - ) - prototype_points.normals = o3d.utility.Vector3dVector( - normals[vertices[:, 2] >= interest_reference, :] - ) - return prototype_mesh, prototype_points - - -def mediumreferencescan(mesh_size, wavelength=3e8 / 24e9, importpcd=True): - if importpcd: - if mesh_size < (wavelength * 0.5): - # stream = pkg_resources.resource_stream(__name__, - # 'data/referencescanpoints24GHztenthwavelength.ply') - stream = files("lyceanem.tests.data").joinpath( - "referencescanpoints24GHztenthwavelength.ply" - ) - reflector_points = o3d.io.read_point_cloud(str(stream)) - else: - # stream = pkg_resources.resource_stream(__name__, - # 'data/referencescanpoints24GHztenthwavelength.ply') - stream = files("lyceanem.tests.data").joinpath( - "referencescanpoints24GHztenthwavelength.ply" - ) - reflector_points = o3d.io.read_point_cloud(str(stream)) - reference_mesh, _ = tl.meshedReflector(0.3, 0.3, 6e-3, mesh_size) - else: - # stream=pkg_resources.resource_stream(__name__,'data/Medium_Reference_Plate_Covered_Normal_1_plan_grid_0p5mm.txt') - stream = files("lyceanem.tests.data").joinpath( - "Medium_Reference_Plate_Covered_Normal_1_plan_grid_0p5mm.txt" - ) - temp = np.loadtxt(str(stream), delimiter=";", skiprows=2) - medium_reference_reflector = o3d.geometry.PointCloud() - medium_reference_reflector.points = o3d.utility.Vector3dVector(temp / 1000) - medium_reference_reflector.estimate_normals() - medium_reference_reflector.orient_normals_to_align_with_direction( - orientation_reference=[0.0, 0.0, 1.0] - ) - reference_point = ( - np.asarray([[np.min(temp[:, 0]), np.min(temp[:, 1]), np.min(temp[:, 2])]]) - / 1000 - ) - - downsampled_reflector = medium_reference_reflector.voxel_down_sample( - voxel_size=mesh_size * 1.1 - ) - downsampled_reflector.translate(-reference_point.ravel(), relative=True) - medium_reference_reflector.translate(-reference_point.ravel(), relative=True) - downsampled_reflector.translate([-0.15, -0.15, 0], relative=True) - medium_reference_reflector.translate([-0.15, -0.15, 0], relative=True) - reflectorplate, scatter_points = tl.meshedReflector( - 0.3, 0.3, 6e-3, mesh_size, sides="all" - ) - temp_points = np.asarray(scatter_points.points) - temp_points_back = temp_points[temp_points[:, 2] < -1e-3, :] - temp_points_top = temp_points[temp_points[:, 1] > 1.49e-1, :] - top_normals = np.zeros((temp_points_top.shape[0], 3)) - top_normals[:, 1] = 1 - toppoints = o3d.geometry.PointCloud() - toppoints.points = o3d.utility.Vector3dVector(temp_points_top) - toppoints.normals = o3d.utility.Vector3dVector(top_normals) - sideapoints = copy.deepcopy(toppoints) - rotation_vector1 = np.radians(np.array([0, 0, 90])) - sideapoints = GF.open3drotate( - sideapoints, - o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1), - ) + source_pcd.point_data['normals'] = total_array_normals + source_pcd = GF.translate_mesh(source_pcd, np.array([-0.18, 0, 0.0125])) - bottompoints = copy.deepcopy(sideapoints) - bottompoints = GF.open3drotate( - bottompoints, - o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1), - ) - - sidebpoints = copy.deepcopy(bottompoints) - sidebpoints = GF.open3drotate( - sidebpoints, - o3d.geometry.TriangleMesh.get_rotation_matrix_from_xyz(rotation_vector1), - ) - back_normals = np.zeros((temp_points_back.shape[0], 3)) - back_normals[:, 2] = -1 - backpoints = o3d.geometry.PointCloud() - backpoints.points = o3d.utility.Vector3dVector(temp_points_back) - backpoints.normals = o3d.utility.Vector3dVector(back_normals) - temp = np.asarray(downsampled_reflector.points) - interest_reference = np.min(temp[:, 2]) - test_prototype = ( - downsampled_reflector - + backpoints - + toppoints - + bottompoints - + sideapoints - + sidebpoints - ) - ( - reference_mesh, - output, - ) = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( - test_prototype, linear_fit=True - ) - reference_mesh.compute_triangle_normals() - reflector_points = o3d.geometry.PointCloud() - vertices = np.asarray(reference_mesh.vertices) - normals = np.asarray(reference_mesh.vertex_normals) - reflector_points.points = o3d.utility.Vector3dVector( - vertices[vertices[:, 2] >= interest_reference, :] - ) - reflector_points.normals = o3d.utility.Vector3dVector( - normals[vertices[:, 2] >= interest_reference, :] - ) - return reference_mesh, reflector_points + return body, array, source_pcd