src.acoustools.BEM.Propagator

  1import torch
  2from torch import Tensor
  3
  4from vedo import Mesh
  5
  6from acoustools.Utilities import TOP_BOARD
  7from acoustools.BEM.Forward_models import compute_E
  8from acoustools.BEM.Gradients import BEM_forward_model_grad, BEM_laplacian
  9import acoustools.Constants as Constants
 10
 11def propagate_BEM(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:Tensor|None=None,
 12                  E:Tensor|None=None,path:str="Media", use_cache_H: bool=True,print_lines:bool=False, 
 13                  p_ref=Constants.P_ref,k:float=Constants.k, betas:float|Tensor = 0, alphas:float|Tensor=1, a=None,c=None, 
 14                  internal_points=None,smooth_distance=0, h=None, BM_alpha=None,
 15                  transducer_radius = Constants.radius) ->Tensor:
 16    '''
 17    Propagates transducer phases to points using BEM\n
 18    :param activations: Transducer hologram
 19    :param points: Points to propagate to
 20    :param scatterer: The mesh used (as a `vedo` `mesh` object)
 21    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
 22    :param H: Precomputed H - if None H will be computed
 23    :param E: Precomputed E - if None E will be computed
 24    :param path: path to folder containing `BEMCache/ `
 25    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
 26    :param print_lines: if true prints messages detaling progress
 27    :param k: wavenumber
 28    :param alphas: Absorbance of each element, can be Tensor for element-wise attribution or a number for all elements
 29    :param betas: Ratio of impedances of medium and scattering material for each element, can be Tensor for element-wise attribution or a number for all elements
 30    :param internal_points: The internal points to use for CHIEF based BEM
 31    :param smooth_distance: amount to add to distances to avoid explosion over small values
 32    :param h: finite difference step for Burton-Miller BEM
 33    :param BM_alpha: constant alpha to use in Burton-Miller BEM
 34
 35    :return pressure: complex pressure at points
 36    '''
 37    if board is None:
 38        board = TOP_BOARD
 39
 40    if E is None:
 41        if type(scatterer) == str:
 42            scatterer = load_scatterer(scatterer)
 43        E = compute_E(scatterer,points,board,H=H, path=path,use_cache_H=use_cache_H,print_lines=print_lines,p_ref=p_ref, k=k, betas=betas, alphas=alphas, a=a, c=c, internal_points=internal_points, smooth_distance=smooth_distance, h=h, BM_alpha=BM_alpha, transducer_radius=transducer_radius)
 44    
 45    out = E@activations
 46    return out
 47
 48def propagate_BEM_pressure(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:
 49                           Tensor|None=None,E:Tensor|None=None, path:str="Media",use_cache_H:bool=True, print_lines:bool=False,p_ref=Constants.P_ref,k:float=Constants.k, betas = 0, alphas:float|Tensor=1, a=None,c=None, internal_points=None,smooth_distance=0, h=None, BM_alpha=None, transducer_radius = Constants.radius) -> Tensor:
 50    '''
 51    Propagates transducer phases to points using BEM and returns absolute value of complex pressure\n
 52    Equivalent to `torch.abs(propagate_BEM(activations,points,scatterer,board,H,E,path))` \n
 53    :param activations: Transducer hologram
 54    :param points: Points to propagate to
 55    :param scatterer: The mesh used (as a `vedo` `mesh` object)
 56    :param board: Transducers to use 
 57    :param H: Precomputed H - if None H will be computed
 58    :param E: Precomputed E - if None E will be computed 
 59    :param path: path to folder containing `BEMCache/ `
 60    :param k: wavenumber
 61    :param alphas: Absorbance of each element, can be Tensor for element-wise attribution or a number for all elements
 62    :param betas: Ratio of impedances of medium and scattering material for each element, can be Tensor for element-wise attribution or a number for all elements
 63    :param internal_points: The internal points to use for CHIEF based BEM
 64    :param smooth_distance: amount to add to distances to avoid explosion over small values
 65    :param h: finite difference step for Burton-Miller BEM
 66    :param BM_alpha: constant alpha to use in Burton-Miller BEM
 67
 68    
 69    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
 70    :param print_lines: if true prints messages detaling progress
 71    
 72    :return pressure: real pressure at points
 73    '''
 74    if board is None:
 75        board = TOP_BOARD
 76
 77    point_activations = propagate_BEM(activations,points,scatterer,board,H,E,path,use_cache_H=use_cache_H,print_lines=print_lines,p_ref=p_ref, k=k, betas=betas, alphas=alphas, a=a, c=c, internal_points=internal_points, smooth_distance=smooth_distance, h=h, BM_alpha=BM_alpha, transducer_radius=transducer_radius)
 78    pressures =  torch.abs(point_activations)
 79    return pressures
 80
 81def propagate_BEM_pressure_grad(activations: Tensor, points: Tensor,board: Tensor|None=None, scatterer:Mesh = None, 
 82                                path:str='Media', Fx=None, Fy=None, Fz=None, cat=True,p_ref=Constants.P_ref, k=Constants.k, transducer_radius = Constants.radius):
 83    '''
 84    Propagates a hologram to pressure gradient at points\n
 85    :param activations: Hologram to use
 86    :param points: Points to propagate to
 87    :param board: The Transducer array, default two 16x16 arrays
 88    :param Fx: The forward model to us for Fx, if None it is computed using `forward_model_grad`. Default:`None`
 89    :param Fy: The forward model to us for Fy, if None it is computed using `forward_model_grad`. Default:`None`
 90    :param Fz: The forward model to us for Fz, if None it is computed using `forward_model_grad`. Default:`None`
 91    :return: point velocity potential'
 92    '''
 93    
 94    if Fx is None or Fy is None or Fz is None:
 95        _Fx,_Fy,_Fz = BEM_forward_model_grad(points, scatterer ,board, p_ref=p_ref, k=k, transducer_radius=transducer_radius)
 96        if Fx is None: Fx = _Fx
 97        if Fy is None: Fy = _Fy
 98        if Fz is None: Fz = _Fz
 99    
100    Px = Fx@activations
101    Py = Fy@activations
102    Pz = Fz@activations
103
104    if cat: 
105        grad = torch.cat([Px, Py, Pz], dim=2)
106        return grad
107    return Px, Py, Pz
108
109def propagate_BEM_phase(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:Tensor|None=None,
110                  E:Tensor|None=None,path:str="Media", use_cache_H: bool=True,print_lines:bool=False,p_ref=Constants.P_ref,k:float=Constants.k, betas = 0, alphas:float|Tensor=1, a=None,c=None, internal_points=None,smooth_distance=0, h=None, BM_alpha=None, transducer_radius=Constants.radius) ->Tensor:
111    '''
112    Propagates transducer phases to phases at points using BEM\n
113    :param activations: Transducer hologram
114    :param points: Points to propagate to
115    :param scatterer: The mesh used (as a `vedo` `mesh` object)
116    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
117    :param H: Precomputed H - if None H will be computed
118    :param E: Precomputed E - if None E will be computed
119    :param path: path to folder containing `BEMCache/ `
120    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
121    :param print_lines: if true prints messages detaling progress
122    :param k: wavenumber
123    :param alphas: Absorbance of each element, can be Tensor for element-wise attribution or a number for all elements
124    :param betas: Ratio of impedances of medium and scattering material for each element, can be Tensor for element-wise attribution or a number for all elements
125    :param internal_points: The internal points to use for CHIEF based BEM
126    :param smooth_distance: amount to add to distances to avoid explosion over small values
127    :param h: finite difference step for Burton-Miller BEM
128    :param BM_alpha: constant alpha to use in Burton-Miller BEM
129
130    :return pressure: phase at points
131    '''
132    if board is None:
133        board = TOP_BOARD
134
135    if E is None:
136        if type(scatterer) == str:
137            scatterer = load_scatterer(scatterer)
138        E = compute_E(scatterer,points,board,H=H, path=path,use_cache_H=use_cache_H,print_lines=print_lines, p_ref=p_ref, k=k, betas=betas, alphas=alphas, a=a, c=c, internal_points=internal_points, smooth_distance=smooth_distance, h=h, BM_alpha=BM_alpha, transducer_radius=transducer_radius)
139    
140    out = E@activations
141    return torch.angle(out)
142
143
144def propagate_BEM_laplacian(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:Tensor|None=None,
145                  E_lap:Tensor|None=None,path:str="Media", use_cache_H: bool=True,print_lines:bool=False,p_ref=Constants.P_ref,k:float=Constants.k, internal_points=None, transducer_radius=Constants.radius):
146    '''
147    Propagate transducer hologram to the laplacian of pressure at points
148    
149    :param activations: Transducer hologram
150    :param points: Points to propagate to
151    :param scatterer: The mesh used (as a `vedo` `mesh` object)
152    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
153    :param H: Precomputed H - if None H will be computed
154    :param E: Precomputed E - if None E will be computed
155    :param path: path to folder containing `BEMCache/ `
156    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
157    :param print_lines: if true prints messages detaling progress
158    :param k: wavenumber
159    :param internal_points: The internal points to use for CHIEF based BEM
160
161    :return pressure: laplacian at points
162    '''
163    
164    if board is None:
165        board = TOP_BOARD
166
167    if E_lap is None:
168        if type(scatterer) == str:
169            scatterer = load_scatterer(scatterer)
170        E_lap = BEM_laplacian(points, scatterer ,board,H=H, path=path,use_cache_H=use_cache_H,print_lines=print_lines, p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
171
172
173    lap = E_lap @ activations
174
175    return lap
176    
177
178def propagate_BEM_laplacian_abs(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:Tensor|None=None,
179                  E_lap:Tensor|None=None,path:str="Media", use_cache_H: bool=True,print_lines:bool=False,p_ref=Constants.P_ref,k:float=Constants.k, internal_points=None, transducer_radius=Constants.radius):
180    '''
181    Propagate transducer hologram to the absolute value of the laplacian of pressure at points
182
183    :param activations: Transducer hologram
184    :param points: Points to propagate to
185    :param scatterer: The mesh used (as a `vedo` `mesh` object)
186    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
187    :param H: Precomputed H - if None H will be computed
188    :param E: Precomputed E - if None E will be computed
189    :param path: path to folder containing `BEMCache/ `
190    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
191    :param print_lines: if true prints messages detaling progress
192    :param k: wavenumber
193    :param internal_points: The internal points to use for CHIEF based BEM
194
195    :return pressure: laplacian at points
196    '''
197    
198    lap = propagate_BEM_laplacian(activations,points, scatterer ,board,H=H, path=path,use_cache_H=use_cache_H,print_lines=print_lines, p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius, E_lap=E_lap)
199
200    return torch.abs(lap)
201
202def propagate_BEM_helmholtz(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:Tensor|None=None, E:Tensor|None=None,
203                  E_lap:Tensor|None=None,path:str="Media", use_cache_H: bool=True,print_lines:bool=False,p_ref=Constants.P_ref,k:float=Constants.k, internal_points=None, transducer_radius=Constants.radius):
204    '''
205    Computes the Helmholtz eq. at points given a hologram
206
207    :param activations: Transducer hologram
208    :param points: Points to propagate to
209    :param scatterer: The mesh used (as a `vedo` `mesh` object)
210    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
211    :param H: Precomputed H - if None H will be computed
212    :param E: Precomputed E - if None E will be computed
213    :param path: path to folder containing `BEMCache/ `
214    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
215    :param print_lines: if true prints messages detaling progress
216    :param k: wavenumber
217    :param internal_points: The internal points to use for CHIEF based BEM
218
219    :return pressure: laplacian at points
220    '''
221    
222    lap = propagate_BEM_laplacian(activations,points, scatterer ,board,H=H, path=path,use_cache_H=use_cache_H,print_lines=print_lines, p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius, E_lap=E_lap)
223    pressure = propagate_BEM(activations=activations,points=points, scatterer=scatterer,board=board, H=H, E=E, use_cache_H=use_cache_H, path=path, p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius)
224
225
226    return lap + k**2 * pressure
227
228
229def propagate_BEM_helmholtz_abs(activations:Tensor,points:Tensor,scatterer:Mesh|None=None,board:Tensor|None=None,H:Tensor|None=None, E:Tensor|None=None,
230                  E_lap:Tensor|None=None,path:str="Media", use_cache_H: bool=True,print_lines:bool=False,p_ref=Constants.P_ref,k:float=Constants.k, internal_points=None, transducer_radius=Constants.radius):
231    
232    '''
233    Computes the absolute value of the Helmholtz eq. at points given a hologram
234
235    :param activations: Transducer hologram
236    :param points: Points to propagate to
237    :param scatterer: The mesh used (as a `vedo` `mesh` object)
238    :param board: Transducers to use, if `None` then uses `acoustools.Utilities.TOP_BOARD` 
239    :param H: Precomputed H - if None H will be computed
240    :param E: Precomputed E - if None E will be computed
241    :param path: path to folder containing `BEMCache/ `
242    :param use_cache_H: If True uses the cache system to load and save the H matrix. Default `True`
243    :param print_lines: if true prints messages detaling progress
244    :param k: wavenumber
245    :param internal_points: The internal points to use for CHIEF based BEM
246
247    :return pressure: laplacian at points
248    '''
249    
250    helmholtz = propagate_BEM_helmholtz(activations,points, scatterer ,board,H=H, path=path,use_cache_H=use_cache_H,print_lines=print_lines, p_ref=p_ref, k=k, internal_points=internal_points, transducer_radius=transducer_radius, E_lap=E_lap, E=E)
251    
252    return torch.abs(helmholtz)