src.acoustools.Mesh

  1from acoustools.Utilities import device, DTYPE
  2import acoustools.Constants as Constants
  3
  4import vedo, torch
  5import matplotlib.pyplot as plt
  6import numpy as np
  7
  8from torch import Tensor
  9from vedo import Mesh
 10
 11
 12def board_name(board:Tensor) -> str:
 13    '''
 14    Returns the name for a board, TOP and/or BOTTOM, used in cache system
 15    :param board: The board to use
 16    :return: name of board as `<'TOP'><'BOTTOM'><M>` for `M` transducers in the boards 
 17    '''
 18    M = board.shape[0]
 19
 20    top = "TOP" if 1 in torch.sign(board[:,2]) else ""
 21    bottom = "BOTTOM" if -1 in torch.sign(board[:,2]) else ""
 22    return top+bottom+str(M)
 23
 24def scatterer_file_name(scatterer:Mesh) ->str:
 25    '''
 26    Get a unique name to describe a scatterer position, calls `str(scatterer.coordinates)`
 27    ONLY USE TO SET FILENAME, USE `scatterer.filename` TO GET
 28    :param scatterer: The Mesh to use
 29    :return: Scatterer name
 30    
 31    '''
 32
 33    f_name = str(list(scatterer.coordinates))
 34    return f_name
 35
 36def load_scatterer(path:str, compute_areas:bool = True, compute_normals:bool=True, dx:float=0,
 37                   dy:float=0,dz:float=0, rotx:float=0, roty:float=0, rotz:float=0, root_path:str="", force:bool=False) -> Mesh:
 38    '''
 39    Loads a scatterer as a `vedo` `Mesh` and applies translations as needed
 40    :param path: The name of the scatterer to load
 41    :param compute_areas: if `True` will call `scatterer.compute_cell_size()`. Default `True`
 42    :param compute_normals: if `True` will call `scatterer.compute_normals()`. Default `True`
 43    :param dx: Translation in the x direction to apply
 44    :param dy: Translation in the y direction to apply
 45    :param dz: Translation in the z direction to apply
 46    :param rotx: Rotation around the x axis to apply
 47    :param roty: Rotation around the y axis to apply
 48    :param rotz: Rotation around the z axis to apply
 49    :param root_path: The folder containing the file, the scatterer to be loaded will be loaded from `root_path+path`
 50    :return: The `vedo` `Mesh` of the scatterer
 51    '''
 52    scatterer = vedo.load(root_path+path, force=force)
 53    
 54    if scatterer is not None:
 55        if compute_areas: scatterer.compute_cell_size()
 56        if compute_normals: 
 57            scatterer.compute_normals()
 58
 59        scatterer.metadata["rotX"] = 0
 60        scatterer.metadata["rotY"] = 0
 61        scatterer.metadata["rotZ"] = 0
 62
 63        # scatterer.filename = scatterer.filename.split("/")[-1]
 64        scatterer.filename = scatterer_file_name(scatterer)
 65
 66        scatterer.metadata["FILE"] = scatterer.filename.split(".")[0]
 67
 68
 69        rotate(scatterer,(1,0,0),rotx)
 70        rotate(scatterer,(0,1,0),roty)
 71        rotate(scatterer,(0,0,1),rotz)
 72
 73        translate(scatterer,dx,dy,dz)
 74
 75
 76    return scatterer
 77
 78def calculate_features(scatterer:Mesh, compute_areas:bool = True, compute_normals:bool=True):
 79    '''
 80    @private
 81    '''
 82    if compute_areas: scatterer.compute_cell_size()
 83    if compute_normals: scatterer.compute_normals()
 84
 85    scatterer.filename = scatterer_file_name(scatterer)
 86    scatterer.metadata["FILE"] = scatterer.filename.split(".")[0]
 87
 88
 89def load_multiple_scatterers(paths:list[str],  compute_areas:bool = True, compute_normals:bool=True, 
 90                             dxs:list[int]=[],dys:list[int]=[],dzs:list[int]=[], rotxs:list[int]=[], rotys:list[int]=[], rotzs:list[int]=[], root_path:str="") -> Mesh:
 91    '''
 92    Loads multiple scatterers and combines them into a single scatterer object
 93    :param path: The name of the scatterers to load
 94    :param compute_areas: if true will call `scatterer.compute_cell_size()`. Default True
 95    :param compute_normals: if true will call `scatterer.compute_normals()`. Default True
 96    :param dxs: List of translations in the x direction to apply to each scatterer
 97    :param dys: List of translations in the y direction to apply to each scatterer
 98    :param dzs: List of translations in the z direction to apply to each scatterer
 99    :param rotxs: List pf rotations around the x axis to apply to each scatterer
100    :param rotys: List pf rotations around the y axis to apply to each scatterer
101    :param rotzs: List pf rotations around the z axis to apply to each scatterer
102    :param root_path: The folder containing the file, the scatterer to be loaded will be loaded from `root_path+path`
103    :return: A merged mesh from all of the paths provided
104    '''
105    dxs += [0] * (len(paths) - len(dxs))
106    dys += [0] * (len(paths) - len(dys))
107    dzs += [0] * (len(paths) - len(dzs))
108
109    rotxs += [0] * (len(paths) - len(rotxs))
110    rotys += [0] * (len(paths) - len(rotys))
111    rotzs += [0] * (len(paths) - len(rotzs))
112
113    scatterers = []
114    for i,path in enumerate(paths):
115        scatterer = load_scatterer(path, compute_areas, compute_normals, dxs[i],dys[i],dzs[i],rotxs[i],rotys[i],rotzs[i],root_path)
116        scatterers.append(scatterer)
117    combined = merge_scatterers(*scatterers)
118    return combined
119
120def merge_scatterers(*scatterers:Mesh, flag:bool=False) ->Mesh:
121    '''
122    Combines any number of scatterers into a single scatterer\n
123    :param scatterers: any number of scatterers to combine
124    :param flag: Value will be passed to `vedo.merge`
125    :return: the combined scatterer
126    '''
127    names = []
128    Fnames = []
129    for scatterer in scatterers:
130        names.append(scatterer_file_name(scatterer))
131        Fnames.append(scatterer.metadata["FILE"][0])
132    
133    if flag:
134        combined = vedo.merge(scatterers, flag=True)
135    else:
136        combined = vedo.merge(scatterers)
137    combined.filename = "".join(names)
138    combined.metadata["FILE"] = "".join(Fnames)
139    return combined
140
141
142def scale_to_diameter(scatterer:Mesh , diameter: float, reset:bool=True, origin:bool=True) -> None:
143    '''
144    Scale a mesh to a given diameter in the x-axis and recomputes normals and areas \n
145    Modifies scatterer in place so does not return anything.\n
146
147    :param scatterer: The scatterer to scale
148    :param diameter: The diameter target
149    '''
150    x1,x2,y1,y2,z1,z2 = scatterer.bounds()
151    diameter_sphere = x2 - x1
152    scatterer.scale(diameter/diameter_sphere,reset=reset, origin=origin)
153    scatterer.compute_cell_size()
154    scatterer.compute_normals()
155    scatterer.filename = scatterer_file_name(scatterer)
156    
157
158def get_plane(scatterer: Mesh, origin:tuple[int]=(0,0,0), normal:tuple[int]=(1,0,0)) -> Mesh:
159    '''
160    Get intersection of a scatterer and a plane\n
161    :param scatterer: The scatterer to intersect
162    :param origin: A point on the plane as a tuple `(x,y,z)`. Default `(0,0,0)`
163    :param normal: The normal to the plane at `point` as a tuple (x,y,z). Default `(1,0,0)`
164    :return: new `Mesh` Containing the intersection of the plane and the scatterer
165    '''
166    intersection = scatterer.clone().intersect_with_plane(origin,normal)
167    intersection.filename = scatterer.filename + "plane" + str(origin)+str(normal)
168    return intersection
169
170def get_lines_from_plane(scatterer:Mesh, origin:tuple[int]=(0,0,0), normal:tuple[int]=(1,0,0)) -> list[int]:
171    '''
172    Gets the edges on a plane from the intersection between a scatterer and the plane\n
173    :param scatterer: The scatterer to intersect
174    :param origin: A point on the plane as a tuple `(x,y,z)`. Default `(0,0,0)`
175    :param normal: The normal to the plane at `point` as a tuple (x,y,z). Default `(1,0,0)`
176    :return: a list of edges in the plane 
177    '''
178
179    mask = [0,0,0]
180    for i in range(3):
181        mask[i] =not normal[i]
182    mask = np.array(mask)
183
184    intersection = get_plane(scatterer, origin, normal)
185    verticies = intersection.vertices
186    lines = intersection.lines
187
188    connections = []
189
190    for i in range(len(lines)):
191        connections.append([verticies[lines[i][0]][mask],verticies[lines[i][1]][mask]])
192
193    return connections
194
195def plot_plane(connections:list[int]) -> None:
196    '''
197    Plot a set of edges assuming they are co-planar\n
198    :param connections: list of connections to plot
199    '''
200    
201    for con in connections:
202        xs = [con[0][0], con[1][0]]
203        ys = [con[0][1], con[1][1]]
204        plt.plot(xs,ys,color = "blue")
205
206    plt.xlim((-0.06,0.06))
207    plt.ylim((-0.06,0.06))
208    plt.show()
209
210def get_normals_as_points(*scatterers:Mesh, permute_to_points:bool=True) -> Tensor:
211    '''
212    Returns the normal vectors to the surface of a scatterer as a `torch` `Tensor` as acoustools points\n
213    :param scatterers: The scatterer to use
214    :param permute_to_points: If true will permute the order of coordinates to agree with what acoustools expects.
215    :return: normals
216    '''
217    norm_list = []
218    for scatterer in scatterers:
219        scatterer.compute_normals()
220        norm =  torch.tensor(scatterer.cell_normals).to(device)
221
222        if permute_to_points:
223            norm = torch.permute(norm,(1,0))
224        
225        norm_list.append(norm.to(DTYPE))
226    
227    return torch.stack(norm_list)
228
229def get_centre_of_mass_as_points(*scatterers:Mesh, permute_to_points:bool=True) ->Tensor:
230    '''
231    Returns the centre of mass(es) of a scatterer(s) as a `torch` `Tensor` as acoustools points\n
232    :param scatterers: The scatterer(s) to use
233    :param permute_to_points: If true will permute the order of coordinates to agree with what acoustools expects.
234    :return: centre of mass(es)
235    '''
236    centres_list = []
237    for scatterer in scatterers:
238        centre_of_mass =  torch.tensor(scatterer.center_of_mass()).to(device)
239
240        if permute_to_points:
241            centre_of_mass = torch.unsqueeze(centre_of_mass,1)
242        
243        centres_list.append(centre_of_mass.to(DTYPE))
244    
245    return torch.real(torch.stack(centres_list))
246
247
248def get_centres_as_points(*scatterers:Mesh, permute_to_points:bool=True, add_normals:bool=False, normal_scale:float=0.001) ->Tensor:
249    '''
250    Returns the centre of scatterer faces as a `torch` `Tensor` as acoustools points\n
251    :param scatterers: The scatterer to use
252    :param permute_to_points: If `True` will permute the order of coordinates to agree with what acoustools expects.
253    :return: centres
254    '''
255    centre_list = []
256    for scatterer in scatterers:
257        centres =  torch.tensor(scatterer.cell_centers().points).to(device)
258
259        if permute_to_points:
260            centres = torch.permute(centres,(1,0)).unsqueeze_(0)
261        
262        if add_normals:
263            norms= get_normals_as_points(scatterer)
264            centres += norms.real * normal_scale
265        
266        centre_list.append(centres.to(DTYPE))
267    centres = torch.cat(centre_list,dim=0)
268    return centres
269
270def get_areas(*scatterers: Mesh) -> Tensor:
271    '''
272    Returns the areas of faces of any number of scatterers\n
273    :param scatterers: The scatterers to use.
274    :return: areas
275    '''
276    area_list = []
277    for scatterer in scatterers:
278        scatterer.compute_cell_size()
279        area_list.append(torch.Tensor(scatterer.celldata["Area"]).to(device))
280    
281    return torch.stack(area_list)
282
283def get_weight(scatterer:Mesh, density:float=Constants.p_p, g:float=9.81) -> float:
284    '''
285    Get the weight of a scatterer\\
286    :param scatterer: The scatterer to use\\
287    :param density: The density to use. Default density for EPS\\
288    :param g: value for g to use. Default 9.81\\
289    :return: weight
290    '''
291    mass = scatterer.volume() * density
292    return g * mass
293
294def translate(scatterer:Mesh, dx:float=0,dy:float=0,dz:float=0) -> None:
295    '''
296    Translates a scatterer by (dx,dy,dz) \n
297    Modifies inplace so does not return a value \n
298    :param scatterer: The scatterer to use
299    :param dx: Translation in the x direction
300    :param dy: Translation in the y direction
301    :param dz: Translation in the z direction
302    '''
303    scatterer.shift(np.array([dx,dy,dz]))
304    scatterer.filename = scatterer_file_name(scatterer)
305
306def rotate(scatterer:Mesh, axis:tuple[int], rot:float, centre:tuple[int]=(0, 0, 0), rotate_around_COM:bool=False):
307    '''
308    Rotates a scatterer in axis by rot\n
309    Modifies inplace so does not return a value\n
310    :param scatterer: The scatterer to use
311    :param axis: The axis to rotate in
312    :param rot: Angle to rotate in degrees
313    :param centre: point to rotate around
314    :param rotate_around_COM: If True will set `centre` to `scatterer`s centre of mass
315    '''
316    if rotate_around_COM:
317        centre = vedo.vector(get_centre_of_mass_as_points(scatterer).cpu().detach().squeeze())
318
319    if axis[0]:
320        scatterer.metadata["rotX"] = scatterer.metadata["rotX"] + rot
321    if axis[1]:
322        scatterer.metadata["rotY"] = scatterer.metadata["rotY"] + rot
323    if axis[2]:
324        scatterer.metadata["rotZ"] = scatterer.metadata["rotZ"] + rot
325    scatterer.rotate(rot, axis,point=centre)
326    scatterer.filename = scatterer_file_name(scatterer)
327
328 
329def downsample(scatterer:Mesh, factor:int=2, n:int|None=None, method:str='quadric', boundaries:bool=False, compute_areas:bool=True, compute_normals:bool=True) -> Mesh:
330    '''
331    Downsamples a mesh to have `factor` less elements\n
332    :param scatterer: The scatterer to use
333    :param factor: The factor to downsample by
334    :param n: The desired number of final points, passed to `Vedo.Mesh.decimate`
335    :param method:, `boundaries` - passed to `vedo.decimate`
336    :param compute_areas: if true will call `scatterer.compute_cell_size()`. Default `True`
337    :param compute_normals: if true will call `scatterer.compute_normals()`. Default `True`
338    :return: downsampled mesh
339    '''
340    scatterer_small =  scatterer.decimate(1/factor, n, method, boundaries)
341    
342    scatterer_small.metadata["rotX"] = scatterer.metadata["rotX"]
343    scatterer_small.metadata["rotY"] = scatterer.metadata["rotY"]
344    scatterer_small.metadata["rotZ"] = scatterer.metadata["rotZ"]
345
346    if compute_areas: scatterer_small.compute_cell_size()
347    if compute_normals: 
348        scatterer_small.compute_normals()
349
350    scatterer_small.filename = scatterer_file_name(scatterer_small)  + "-scale-" + str(factor)
351
352
353    return scatterer_small
354
355
356def centre_scatterer(scatterer:Mesh) -> list[int]:
357    '''
358    Translate scatterer so the centre of mass is at (0,0,0)\n
359    Modifies Mesh in place \n
360    :param scatterer: Scatterer to centre
361    :return: Returns the amount needed to move in each direction
362    '''
363    com = get_centre_of_mass_as_points(scatterer).cpu()
364    correction = [-1*com[:,0].item(), -1*com[:,1].item(), -1*com[:,2].item()]
365    translate(scatterer, dx = correction[0], dy = correction[1], dz=  correction[2])
366
367    return correction
368
369
370def get_edge_data(scatterer:Mesh, wavelength:float=Constants.wavelength, print_output:bool=True, break_down_average:bool=False) -> None|tuple[float]:
371    '''
372    Get the maximum, minimum and average size of edges in a mesh. Optionally prints or returns the result.\n
373    :param scatterer: Mesh of interest
374    :param wavelength: Wavenelgth size for printing results as multiple of some wavelength
375    :param print_output: If True, prints results else returns values
376    :break_down_average: If True will also return (distance_sum, N)
377    :return: None if `print_outputs` is `True` else returns `(max_distance, min_distance, average_distance)` and optionally  (distance_sum, N)
378
379    '''
380    points = scatterer.vertices
381
382    distance_sum = 0
383    N = 0
384
385    max_distance = 0
386    min_distance = 100000000
387
388
389    for (start,end) in scatterer.edges:
390        start_point = points[start]
391        end_point = points[end]
392        sqvec = torch.Tensor((start_point-end_point)**2)
393        # print(sqvec, torch.sum(sqvec)**0.5)
394        distance = torch.sum(sqvec)**0.5
395        distance_sum += distance
396        N += 1
397        if distance < min_distance:
398            min_distance = distance
399        if distance > max_distance:
400            max_distance = distance
401
402    average_distance = distance_sum/N
403
404    if print_output:
405        print('Max Distance', max_distance.item(),'=' ,max_distance.item()/wavelength, 'lambda')
406        print('Min Distance', min_distance.item(),'=', min_distance.item()/wavelength, 'lambda')
407        print('Ave Distance', average_distance.item(),'=', average_distance.item()/wavelength, 'lambda')
408    else:
409        if break_down_average:
410            return (max_distance, min_distance, average_distance), (distance_sum, N)
411        else:
412            return (max_distance, min_distance, average_distance)
413
414
415def cut_mesh_to_walls(scatterer:Mesh, layer_z:float, layer_normal:tuple[float] = (0,0,-1.0), wall_thickness = 0.001) -> Mesh:
416    '''
417    Cuts a mesh with a given plane and then converts the result to have walls of a certain thickness \n
418    :param scatterer: Mesh to use
419    :param layer_z: coordinate of layer
420    :param layer_normal: Normal to layer (if not +- (0,0,1) then layer_z will not refer to a z coordinate)
421    :param wall_thickness: Thickness of the walls to returns
422    :return: Cut mesh with walls
423    '''
424
425    xmin,xmax, ymin,ymax, zmin,zmax = scatterer.bounds()
426    dx = xmax-xmin
427    dy = ymax-ymin
428
429    scale_x = (dx-2*wall_thickness) / dx
430    scale_y = (dy-2*wall_thickness) / dy
431
432    outler_layer = scatterer.cut_with_plane((0,0,layer_z),layer_normal)
433    inner_layer = outler_layer.clone()
434    inner_layer.scale((scale_x,scale_y,1), origin=False)
435
436    com_outer = get_centre_of_mass_as_points(outler_layer)
437    com_inner = get_centre_of_mass_as_points(inner_layer)
438
439    d_com = (com_outer - com_inner).squeeze()
440
441    translate(inner_layer, *d_com)
442
443    walls = vedo.merge(outler_layer, inner_layer)
444
445
446    boundaries_outer = outler_layer.boundaries()
447    boundaries_inner = inner_layer.boundaries()
448
449    strips = boundaries_outer.join_with_strips(boundaries_inner).triangulate()
450    
451
452    walls = vedo.merge(walls,strips)
453    
454    calculate_features(walls)
455    scatterer_file_name(walls)
456
457    return walls.clean()
458
459def cut_closed_scatterer(scatterer:Mesh,layer_z:float, normals=[(0,0,1)]):
460    origins=[(0,0,layer_z)]
461    closed_scatterer = scatterer.cut_closed_surface(origins=origins, normals=normals)
462    return closed_scatterer
463
464def get_volume(scatterer:Mesh):
465    '''
466    Returns the volume of a mesh
467    '''
468    return scatterer.volume()