src.acoustools.Mesh

  1from acoustools.Utilities import device, DTYPE, BOARD_POSITIONS
  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
 10from typing import Literal
 11
 12
 13def board_name(board:Tensor) -> str:
 14    '''
 15    Returns the name for a board, TOP and/or BOTTOM, used in cache system
 16    :param board: The board to use
 17    :return: name of board as `<'TOP'><'BOTTOM'><M>` for `M` transducers in the boards 
 18    '''
 19    M = board.shape[0]
 20
 21    top = "TOP" if 1 in torch.sign(board[:,2]) else ""
 22    bottom = "BOTTOM" if -1 in torch.sign(board[:,2]) else ""
 23    return top+bottom+str(M)
 24
 25def scatterer_file_name(scatterer:Mesh) ->str:
 26    '''
 27    Get a unique name to describe a scatterer position, calls `str(scatterer.coordinates)`
 28    ONLY USE TO SET FILENAME, USE `scatterer.filename` TO GET
 29    :param scatterer: The Mesh to use
 30    :return: Scatterer name
 31    
 32    '''
 33
 34    f_name = str(list(scatterer.coordinates)) + str(scatterer.cell_normals)
 35    return f_name
 36
 37def load_scatterer(path:str, compute_areas:bool = True, compute_normals:bool=True, dx:float=0,
 38                   dy:float=0,dz:float=0, rotx:float=0, roty:float=0, rotz:float=0, root_path:str="", force:bool=False, flip_normals=False) -> Mesh:
 39    '''
 40    Loads a scatterer as a `vedo` `Mesh` and applies translations as needed
 41    :param path: The name of the scatterer to load
 42    :param compute_areas: if `True` will call `scatterer.compute_cell_size()`. Default `True`
 43    :param compute_normals: if `True` will call `scatterer.compute_normals()`. Default `True`
 44    :param dx: Translation in the x direction to apply
 45    :param dy: Translation in the y direction to apply
 46    :param dz: Translation in the z direction to apply
 47    :param rotx: Rotation around the x axis to apply
 48    :param roty: Rotation around the y axis to apply
 49    :param rotz: Rotation around the z axis to apply
 50    :param root_path: The folder containing the file, the scatterer to be loaded will be loaded from `root_path+path`
 51    :return: The `vedo` `Mesh` of the scatterer
 52    '''
 53    scatterer = vedo.load(root_path+path, force=force)
 54    
 55    if scatterer is not None:
 56        if compute_areas: scatterer.compute_cell_size()
 57        if compute_normals: 
 58            scatterer.compute_normals()
 59            if flip_normals: scatterer.flip_normals()
 60
 61        scatterer.metadata["rotX"] = 0
 62        scatterer.metadata["rotY"] = 0
 63        scatterer.metadata["rotZ"] = 0
 64
 65        # scatterer.filename = scatterer.filename.split("/")[-1]
 66        scatterer.filename = scatterer_file_name(scatterer)
 67
 68        scatterer.metadata["FILE"] = scatterer.filename.split(".")[0]
 69
 70
 71        rotate(scatterer,(1,0,0),rotx)
 72        rotate(scatterer,(0,1,0),roty)
 73        rotate(scatterer,(0,0,1),rotz)
 74
 75        translate(scatterer,dx,dy,dz)
 76    else:
 77        raise ValueError(f"File not found at {path} - please check the path")
 78
 79    return scatterer
 80
 81def mesh_to_board(path:str, compute_areas:bool = True, compute_normals:bool=True, dx:float=0,
 82                   dy:float=0,dz:float=0, rotx:float=0, roty:float=0, rotz:float=0, root_path:str="", force:bool=False, flip_normals = True, diameter = 2*BOARD_POSITIONS, centre=True):
 83    '''
 84    Loads a scatterer as a `vedo` `Mesh` and interprets it as a transducer board with a transducer at each mesh centre
 85    :param path: The name of the scatterer to load
 86    :param compute_areas: if `True` will call `scatterer.compute_cell_size()`. Default `True`
 87    :param compute_normals: if `True` will call `scatterer.compute_normals()`. Default `True`
 88    :param dx: Translation in the x direction to apply
 89    :param dy: Translation in the y direction to apply
 90    :param dz: Translation in the z direction to apply
 91    :param rotx: Rotation around the x axis to apply
 92    :param roty: Rotation around the y axis to apply
 93    :param rotz: Rotation around the z axis to apply
 94    :param root_path: The folder containing the file, the scatterer to be loaded will be loaded from `root_path+path`
 95    :param flip_normals: If True will flip the normals 
 96    :param diameter: Size to scale the mesh to in thr x-axis
 97    :return: The `vedo` `Mesh` of the scatterer
 98    '''
 99    
100    scatterer = load_scatterer(path=path,compute_areas=compute_areas,  compute_normals=compute_normals, 
101                               dx=dx, dy=dy, dz=dz, rotx=rotx, roty=roty, rotz=rotz, root_path=root_path,force=force)
102    
103    if centre: centre_scatterer(scatterer)
104    
105    if diameter is not None: scale_to_diameter(scatterer, diameter)
106    
107    centres = get_centres_as_points(scatterer).squeeze(0).permute(1,0)
108    norms = get_normals_as_points(scatterer).squeeze(0).permute(1,0)  
109    if flip_normals: norms = norms * -1
110
111    return centres, norms
112    
113
114def calculate_features(scatterer:Mesh, compute_areas:bool = True, compute_normals:bool=True):
115    '''
116    @private
117    '''
118    if compute_areas: scatterer.compute_cell_size()
119    if compute_normals: scatterer.compute_normals()
120
121    scatterer.filename = scatterer_file_name(scatterer)
122    scatterer.metadata["FILE"] = scatterer.filename.split(".")[0]
123
124
125def load_multiple_scatterers(paths:list[str],  compute_areas:bool = True, compute_normals:bool=True, 
126                             dxs:list[int]=[],dys:list[int]=[],dzs:list[int]=[], rotxs:list[int]=[], rotys:list[int]=[], rotzs:list[int]=[], root_path:str="") -> Mesh:
127    '''
128    Loads multiple scatterers and combines them into a single scatterer object
129    :param path: The name of the scatterers to load
130    :param compute_areas: if true will call `scatterer.compute_cell_size()`. Default True
131    :param compute_normals: if true will call `scatterer.compute_normals()`. Default True
132    :param dxs: List of translations in the x direction to apply to each scatterer
133    :param dys: List of translations in the y direction to apply to each scatterer
134    :param dzs: List of translations in the z direction to apply to each scatterer
135    :param rotxs: List pf rotations around the x axis to apply to each scatterer
136    :param rotys: List pf rotations around the y axis to apply to each scatterer
137    :param rotzs: List pf rotations around the z axis to apply to each scatterer
138    :param root_path: The folder containing the file, the scatterer to be loaded will be loaded from `root_path+path`
139    :return: A merged mesh from all of the paths provided
140    '''
141    dxs += [0] * (len(paths) - len(dxs))
142    dys += [0] * (len(paths) - len(dys))
143    dzs += [0] * (len(paths) - len(dzs))
144
145    rotxs += [0] * (len(paths) - len(rotxs))
146    rotys += [0] * (len(paths) - len(rotys))
147    rotzs += [0] * (len(paths) - len(rotzs))
148
149    scatterers = []
150    for i,path in enumerate(paths):
151        scatterer = load_scatterer(path, compute_areas, compute_normals, dxs[i],dys[i],dzs[i],rotxs[i],rotys[i],rotzs[i],root_path)
152        scatterers.append(scatterer)
153    combined = merge_scatterers(*scatterers)
154    return combined
155
156def merge_scatterers(*scatterers:Mesh, flag:bool=False) ->Mesh:
157    '''
158    Combines any number of scatterers into a single scatterer\n
159    :param scatterers: any number of scatterers to combine
160    :param flag: Value will be passed to `vedo.merge`
161    :return: the combined scatterer
162    '''
163    names = []
164    Fnames = []
165    for scatterer in scatterers:
166        names.append(scatterer_file_name(scatterer))
167        Fnames.append(scatterer.metadata["FILE"][0])
168    
169    if flag:
170        combined = vedo.merge(scatterers, flag=True)
171    else:
172        combined = vedo.merge(scatterers)
173    combined.filename = "".join(names)
174    combined.metadata["FILE"] = "".join(Fnames)
175    return combined
176
177
178def scale_to_diameter(scatterer:Mesh , diameter: float, reset:bool=True, origin:bool=True) -> None:
179    '''
180    Scale a mesh to a given diameter in the x-axis and recomputes normals and areas \n
181    Modifies scatterer in place so does not return anything.\n
182
183    :param scatterer: The scatterer to scale
184    :param diameter: The diameter target
185    '''
186    x1,x2,y1,y2,z1,z2 = scatterer.bounds()
187    diameter_sphere = x2 - x1
188    scatterer.scale(diameter/diameter_sphere,reset=reset, origin=origin)
189    scatterer.compute_cell_size()
190    scatterer.compute_normals()
191    scatterer.filename = scatterer_file_name(scatterer)
192    
193def get_diameter(scatterer:Mesh):
194    x1,x2,y1,y2,z1,z2 = scatterer.bounds()
195    diameter_sphere = torch.norm(torch.Tensor([x2,]) - torch.Tensor([x1,]), p=2)
196    return diameter_sphere
197
198
199def get_plane(scatterer: Mesh, origin:tuple[int]=(0,0,0), normal:tuple[int]=(1,0,0)) -> Mesh:
200    '''
201    Get intersection of a scatterer and a plane\n
202    :param scatterer: The scatterer to intersect
203    :param origin: A point on the plane as a tuple `(x,y,z)`. Default `(0,0,0)`
204    :param normal: The normal to the plane at `point` as a tuple (x,y,z). Default `(1,0,0)`
205    :return: new `Mesh` Containing the intersection of the plane and the scatterer
206    '''
207    intersection = scatterer.clone().intersect_with_plane(origin,normal)
208    intersection.filename = scatterer.filename + "plane" + str(origin)+str(normal)
209    return intersection
210
211def get_lines_from_plane(scatterer:Mesh, origin:tuple[int]=(0,0,0), normal:tuple[int]=(1,0,0)) -> list[int]:
212    '''
213    Gets the edges on a plane from the intersection between a scatterer and the plane\n
214    :param scatterer: The scatterer to intersect
215    :param origin: A point on the plane as a tuple `(x,y,z)`. Default `(0,0,0)`
216    :param normal: The normal to the plane at `point` as a tuple (x,y,z). Default `(1,0,0)`
217    :return: a list of edges in the plane 
218    '''
219
220    mask = [0,0,0]
221    for i in range(3):
222        mask[i] =not normal[i]
223    mask = np.array(mask)
224
225    intersection = get_plane(scatterer, origin, normal)
226    verticies = intersection.vertices
227    lines = intersection.lines
228
229    connections = []
230
231    for i in range(len(lines)):
232        connections.append([verticies[lines[i][0]][mask],verticies[lines[i][1]][mask]])
233
234    return connections
235
236def plot_plane(connections:list[int]) -> None:
237    '''
238    Plot a set of edges assuming they are co-planar\n
239    :param connections: list of connections to plot
240    '''
241    
242    for con in connections:
243        xs = [con[0][0], con[1][0]]
244        ys = [con[0][1], con[1][1]]
245        plt.plot(xs,ys,color = "blue")
246
247    plt.xlim((-0.06,0.06))
248    plt.ylim((-0.06,0.06))
249    plt.show()
250
251def get_normals_as_points(*scatterers:Mesh, permute_to_points:bool=True) -> Tensor:
252    '''
253    Returns the normal vectors to the surface of a scatterer as a `torch` `Tensor` as acoustools points\n
254    :param scatterers: The scatterer to use
255    :param permute_to_points: If true will permute the order of coordinates to agree with what acoustools expects.
256    :return: normals
257    '''
258    norm_list = []
259    for scatterer in scatterers:
260        scatterer.compute_normals()
261        norm =  torch.tensor(scatterer.cell_normals).to(device)
262
263        if permute_to_points:
264            norm = torch.permute(norm,(1,0))
265        
266        norm_list.append(norm.to(DTYPE))
267    
268    return torch.stack(norm_list)
269
270def get_centre_of_mass_as_points(*scatterers:Mesh, permute_to_points:bool=True) ->Tensor:
271    '''
272    Returns the centre of mass(es) of a scatterer(s) as a `torch` `Tensor` as acoustools points\n
273    :param scatterers: The scatterer(s) to use
274    :param permute_to_points: If true will permute the order of coordinates to agree with what acoustools expects.
275    :return: centre of mass(es)
276    '''
277    centres_list = []
278    for scatterer in scatterers:
279        centre_of_mass =  torch.tensor(scatterer.center_of_mass()).to(DTYPE).to(device)
280
281        if permute_to_points:
282            centre_of_mass = torch.unsqueeze(centre_of_mass,1)
283        
284        centres_list.append(centre_of_mass.to(DTYPE))
285    
286    return torch.real(torch.stack(centres_list))
287
288
289def get_centres_as_points(*scatterers:Mesh, permute_to_points:bool=True, add_normals:bool=False, normal_scale:float=0.001) ->Tensor:
290    '''
291    Returns the centre of scatterer faces as a `torch` `Tensor` as acoustools points\n
292    :param scatterers: The scatterer to use
293    :param permute_to_points: If `True` will permute the order of coordinates to agree with what acoustools expects.
294    :return: centres
295    '''
296    centre_list = []
297    for scatterer in scatterers:
298        centres =  torch.tensor(scatterer.cell_centers().points).to(DTYPE).to(device)
299
300        if permute_to_points:
301            centres = torch.permute(centres,(1,0)).unsqueeze_(0)
302        
303        if add_normals:
304            norms= get_normals_as_points(scatterer)
305            centres += norms.real * normal_scale
306        
307        centre_list.append(centres)
308    centres = torch.cat(centre_list,dim=0)
309    return centres
310
311def get_verticies_as_points(*scatterers:Mesh):
312    '''
313    Gets the verticies of a mesh as a Tensor of AcousTools (B,3,N) points \n
314    :param Mesh: Mesh to use
315    :returns verticies: verticies as points
316    '''
317   
318    vert_list = []
319    for scatterer in scatterers:
320        vert =  torch.tensor(scatterer.vertices).to(DTYPE).to(device)
321        vert_list.append(vert)
322
323    verts = torch.cat(vert_list,dim=0).unsqueeze(0).permute(0,2,1)
324    return verts
325
326def get_cell_verticies(*scatterers:Mesh):
327    '''
328    Gets a tensor of (B,3,M,3) - batch x (xyz) x Faces x (vertex) \n
329    :param Mesh: Mesh to use
330    :returns verticies: verticies
331    '''
332    verts = get_verticies_as_points(*scatterers)
333    vert_list = []
334    for scatterer in scatterers:
335        cells = torch.tensor(scatterer.cells)
336        N = cells.shape[0]
337        cell_indexes = cells.flatten()
338        cell_verts = torch.index_select(verts, 2, cell_indexes)
339        cell_verts=cell_verts.reshape(1,3,N,3)
340
341
342        vert_list.append(cell_verts)
343    verts = torch.cat(vert_list,dim=0)
344    return verts
345
346
347def get_barycentric_points(*scatterers:Mesh, N=7, sum=True):
348    '''
349    @private
350    '''
351    
352
353    if N != 7: raise ValueError("Only N=7 is supported") #Allow for N as a parameter incase it it implemented in future
354
355    cell_verts = get_cell_verticies(*scatterers)
356
357    DUNAVANT_7 = torch.tensor([
358    [1/3, 1/3, 1/3, 0.225],
359    [0.0597158717, 0.4701420641, 0.4701420641, 0.1323941527],
360    [0.4701420641, 0.0597158717, 0.4701420641, 0.1323941527],
361    [0.4701420641, 0.4701420641, 0.0597158717, 0.1323941527],
362    [0.7974269853, 0.1012865073, 0.1012865073, 0.1259391805],
363    [0.1012865073, 0.7974269853, 0.1012865073, 0.1259391805],
364    [0.1012865073, 0.1012865073, 0.7974269853, 0.1259391805],
365    ])
366    DUNAVANT_7_abg = DUNAVANT_7[:,:3].permute(1,0).unsqueeze(0).unsqueeze(0).unsqueeze(0)
367
368
369    DUNAVANT_7_W = DUNAVANT_7[:,3]
370
371    cell_verts = cell_verts.unsqueeze(-1)
372    barycentric_verts = cell_verts * DUNAVANT_7_abg
373    if sum: barycentric_verts = barycentric_verts.sum(dim=3)
374
375    return barycentric_verts, DUNAVANT_7_W
376
377    
378
379def get_areas(*scatterers: Mesh) -> Tensor:
380    '''
381    Returns the areas of faces of any number of scatterers\n
382    :param scatterers: The scatterers to use.
383    :return: areas
384    '''
385    area_list = []
386    for scatterer in scatterers:
387        scatterer.compute_cell_size()
388        area_list.append(torch.Tensor(scatterer.celldata["Area"]).to(device))
389    
390    return torch.stack(area_list)
391
392def get_weight(scatterer:Mesh, density:float=Constants.p_p, g:float=9.81) -> float:
393    '''
394    Get the weight of a scatterer\\
395    :param scatterer: The scatterer to use\\
396    :param density: The density to use. Default density for EPS\\
397    :param g: value for g to use. Default 9.81\\
398    :return: weight
399    '''
400    mass = scatterer.volume() * density
401    return g * mass
402
403def translate(scatterer:Mesh, dx:float=0,dy:float=0,dz:float=0) -> None:
404    '''
405    Translates a scatterer by (dx,dy,dz) \n
406    Modifies inplace so does not return a value \n
407    :param scatterer: The scatterer to use
408    :param dx: Translation in the x direction
409    :param dy: Translation in the y direction
410    :param dz: Translation in the z direction
411    '''
412    scatterer.shift(np.array([dx,dy,dz]))
413    scatterer.filename = scatterer_file_name(scatterer)
414
415def rotate(scatterer:Mesh, axis:tuple[int], rot:float, centre:tuple[int]=(0, 0, 0), rotate_around_COM:bool=False):
416    '''
417    Rotates a scatterer in axis by rot\n
418    Modifies inplace so does not return a value\n
419    :param scatterer: The scatterer to use
420    :param axis: The axis to rotate in
421    :param rot: Angle to rotate in degrees
422    :param centre: point to rotate around
423    :param rotate_around_COM: If True will set `centre` to `scatterer`s centre of mass
424    '''
425    if rotate_around_COM:
426        centre = vedo.vector(get_centre_of_mass_as_points(scatterer).cpu().detach().squeeze())
427
428    if axis[0]:
429        scatterer.metadata["rotX"] = scatterer.metadata["rotX"] + rot
430    if axis[1]:
431        scatterer.metadata["rotY"] = scatterer.metadata["rotY"] + rot
432    if axis[2]:
433        scatterer.metadata["rotZ"] = scatterer.metadata["rotZ"] + rot
434    scatterer.rotate(rot, axis,point=centre)
435    scatterer.filename = scatterer_file_name(scatterer)
436
437 
438def 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:
439    '''
440    Downsamples a mesh to have `factor` less elements\n
441    :param scatterer: The scatterer to use
442    :param factor: The factor to downsample by
443    :param n: The desired number of final points, passed to `Vedo.Mesh.decimate`
444    :param method:, `boundaries` - passed to `vedo.decimate`
445    :param compute_areas: if true will call `scatterer.compute_cell_size()`. Default `True`
446    :param compute_normals: if true will call `scatterer.compute_normals()`. Default `True`
447    :return: downsampled mesh
448    '''
449    scatterer_small =  scatterer.decimate(1/factor, n, method, boundaries)
450    
451    scatterer_small.metadata["rotX"] = scatterer.metadata["rotX"]
452    scatterer_small.metadata["rotY"] = scatterer.metadata["rotY"]
453    scatterer_small.metadata["rotZ"] = scatterer.metadata["rotZ"]
454
455    if compute_areas: scatterer_small.compute_cell_size()
456    if compute_normals: 
457        scatterer_small.compute_normals()
458
459    scatterer_small.filename = scatterer_file_name(scatterer_small)  + "-scale-" + str(factor)
460
461
462    return scatterer_small
463
464
465def centre_scatterer(scatterer:Mesh) -> list[int]:
466    '''
467    Translate scatterer so the centre of mass is at (0,0,0)\n
468    Modifies Mesh in place \n
469    :param scatterer: Scatterer to centre
470    :return: Returns the amount needed to move in each direction
471    '''
472    com = get_centre_of_mass_as_points(scatterer).cpu()
473    correction = [-1*com[:,0].item(), -1*com[:,1].item(), -1*com[:,2].item()]
474    translate(scatterer, dx = correction[0], dy = correction[1], dz=  correction[2])
475
476    return correction
477
478
479def get_edge_data(scatterer:Mesh, wavelength:float=Constants.wavelength, print_output:bool=True, break_down_average:bool=False) -> None|tuple[float]:
480    '''
481    Get the maximum, minimum and average size of edges in a mesh. Optionally prints or returns the result.\n
482    :param scatterer: Mesh of interest
483    :param wavelength: Wavenelgth size for printing results as multiple of some wavelength
484    :param print_output: If True, prints results else returns values
485    :break_down_average: If True will also return (distance_sum, N)
486    :return: None if `print_outputs` is `True` else returns `(max_distance, min_distance, average_distance)` and optionally  (distance_sum, N)
487
488    '''
489    points = scatterer.vertices
490
491    distance_sum = 0
492    N = 0
493
494    max_distance = 0
495    min_distance = 100000000
496
497
498    for (start,end) in scatterer.edges:
499        start_point = points[start]
500        end_point = points[end]
501        sqvec = torch.Tensor((start_point-end_point)**2)
502        # print(sqvec, torch.sum(sqvec)**0.5)
503        distance = torch.sum(sqvec)**0.5
504        distance_sum += distance
505        N += 1
506        if distance < min_distance:
507            min_distance = distance
508        if distance > max_distance:
509            max_distance = distance
510
511    average_distance = distance_sum/N
512
513    if print_output:
514        print('Max Distance', max_distance.item(),'=' ,max_distance.item()/wavelength, 'lambda')
515        print('Min Distance', min_distance.item(),'=', min_distance.item()/wavelength, 'lambda')
516        print('Ave Distance', average_distance.item(),'=', average_distance.item()/wavelength, 'lambda')
517    else:
518        if break_down_average:
519            return (max_distance, min_distance, average_distance), (distance_sum, N)
520        else:
521            return (max_distance, min_distance, average_distance)
522
523
524def cut_mesh_to_walls(scatterer:Mesh, layer_z:float, layer_normal:tuple[float] = (0,0,-1.0), wall_thickness = 0.001) -> Mesh:
525    '''
526    Cuts a mesh with a given plane and then converts the result to have walls of a certain thickness \n
527    :param scatterer: Mesh to use
528    :param layer_z: coordinate of layer
529    :param layer_normal: Normal to layer (if not +- (0,0,1) then layer_z will not refer to a z coordinate)
530    :param wall_thickness: Thickness of the walls to returns
531    :return: Cut mesh with walls
532    '''
533
534    xmin,xmax, ymin,ymax, zmin,zmax = scatterer.bounds()
535    dx = xmax-xmin
536    dy = ymax-ymin
537
538    scale_x = (dx-2*wall_thickness) / dx
539    scale_y = (dy-2*wall_thickness) / dy
540
541    outler_layer = scatterer.cut_with_plane((0,0,layer_z),layer_normal)
542    inner_layer = outler_layer.clone()
543    inner_layer.scale((scale_x,scale_y,1), origin=False)
544
545    com_outer = get_centre_of_mass_as_points(outler_layer)
546    com_inner = get_centre_of_mass_as_points(inner_layer)
547
548    d_com = (com_outer - com_inner).squeeze()
549
550    translate(inner_layer, *d_com)
551
552    walls = vedo.merge(outler_layer, inner_layer)
553
554
555    boundaries_outer = outler_layer.boundaries()
556    boundaries_inner = inner_layer.boundaries()
557
558    strips = boundaries_outer.join_with_strips(boundaries_inner).triangulate()
559    
560
561    walls = vedo.merge(walls,strips)
562    
563    calculate_features(walls)
564    scatterer_file_name(walls)
565
566    return walls.clean()
567
568def cut_closed_scatterer(scatterer:Mesh,layer_z:float, normals=[(0,0,1)]):
569    '''
570    Cuts a scatterer across a z-plane\\
571    :param scatterer: Mesh
572    :param layer_z: height to cute
573    :param normals: Which way is up 
574    '''
575    origins=[(0,0,layer_z)]
576    closed_scatterer = scatterer.cut_closed_surface(origins=origins, normals=normals)
577    return closed_scatterer
578
579def get_volume(scatterer:Mesh):
580    '''
581    Returns the volume of a mesh
582    '''
583    return scatterer.volume()
584
585def insert_parasite(scatterer:Mesh, parasite_path:str = '/Sphere-lam1.stl', root_path:str="../BEMMedia", parasite_size:float=Constants.wavelength/4, parasite_offset:Tensor=None) -> Mesh:
586    '''
587    Inserts a parasitic body into an existing scatterer. Used to supress the resonance from BEM \n
588    See https://doi.org/10.1109/8.310000 \n
589    :param scatterer: The scatterer to insert parasite into
590    :param parasite_path: The path to the mesh to load and use as parasite
591    :param root_path: The folder to load the file from
592    :param parasite_size: The diameter to scale the parasite to
593    :param parasite_offset: Tensor of offsets for the parasite from the (0,0,0) point
594    :returns: Scatterer with parasite inserted
595    '''
596    parasite = load_scatterer(parasite_path, root_path=root_path)
597    centre_scatterer(parasite)
598    if parasite_offset is None:
599        parasite_offset = get_centre_of_mass_as_points(scatterer)
600
601    dx = parasite_offset[:,0].item()
602    dy = parasite_offset[:,1].item()
603    dz = parasite_offset[:,2].item()
604
605    translate(parasite, dx=dx, dy=dy, dz=dz)
606
607    scale_to_diameter(parasite, parasite_size)
608
609    infected_scatterer = merge_scatterers(scatterer, parasite)
610
611    return infected_scatterer
612
613def get_CHIEF_points(scatterer:Mesh, P=30, method:Literal['random', 'uniform', 'volume-random']='random', start:Literal['surface', 'centre']='surface', scale=0.001, scale_mode:Literal['abs','diameter-scale']='abs') -> Mesh:
614    '''
615    Generates internal points that can be used for the CHIEF BEM formulation (or any other reason)\n
616    :param scatterer: The scatterer to insert points into
617    :param P: Number of points. if P=-1 then P= number of mesh elements
618    :param method: The method used to generate points \n
619        - random: will move scale metres along each of P randomly selected normals \n
620        - uniform:  will move scale metres along each of P uniformly spaced normals (based on order coming from `Mesh.get_normals_as_points`) \n
621        - volume-random: will use `vedo.Mesh..generate_random_points` to generate P internal points
622    :param start: The point to use as the basis for generating points \n
623         - surface: Will step along normals from surface (will step in the -ve normal direction)
624         - centre: Will step along normal from centre of mass (will step in +ve normal direction)
625    :param scale: The distance in m to step 
626    :returns internal points:
627    '''
628
629    centre_norms = get_normals_as_points(scatterer, permute_to_points=False)
630
631    if scale_mode.lower() == 'diameter-scale':
632        d = get_diameter(scatterer)
633        scale = scale * d
634    
635
636    if start.lower() == 'centre':
637        centres = get_centre_of_mass_as_points(scatterer, permute_to_points=False).unsqueeze(1)
638        internal_points = centres + centre_norms * scale      
639
640    else:
641        centres = torch.tensor(scatterer.cell_centers().points, dtype=DTYPE, device=device)
642        internal_points = centres - centre_norms * scale
643
644    M = centre_norms.shape[1]
645    
646    if P == -1: P = M
647
648   
649
650
651    
652    if method.lower() == 'random':
653        indices = torch.randperm(M)[:P]
654        internal_points = internal_points[:, indices,:]
655
656    elif method.lower()== 'uniform':
657        idx = [i for i in range(M) if i%(int(M/P)) == 0]
658        internal_points = internal_points[:, idx,:]
659    elif method.lower() == 'volume-random':
660        internal_points = torch.Tensor(scatterer.generate_random_points(P).points).unsqueeze(0)
661
662    internal_points = internal_points.permute(0,2,1)
663
664
665    return internal_points