"""Module for the radiosity simulation."""
import numpy as np
import pyfar as pf
from . import form_factor, source_energy, receiver_energy, geometry
from . import energy_exchange_order as ee_order
[docs]
class DirectionalRadiosityFast():
"""Radiosity object for directional scattering coefficients."""
_walls_points: np.ndarray
_walls_normal: np.ndarray
_walls_up_vector: np.ndarray
_patches_points: np.ndarray
_patches_normal: np.ndarray
_patch_size: float
_n_patches: int
_speed_of_sound: float
_visibility_matrix: np.ndarray
_visible_patches: np.ndarray
_form_factors: np.ndarray
_form_factors_tilde: np.ndarray
# general data for material data
_n_bins: int
_frequencies = np.ndarray
# absorption data
_absorption: np.ndarray
_absorption_index: np.ndarray
_scattering: np.ndarray
_scattering_index: np.ndarray
_sources: list[pf.Coordinates]
_receivers: list[pf.Coordinates]
_air_attenuation: np.ndarray
_patch_to_wall_ids: np.ndarray
def __init__(
self, walls_points, walls_normal, walls_up_vector,
patches_points, patches_normal, patch_size, n_patches,
patch_to_wall_ids):
"""Create a Radiosity object for directional implementation."""
self._walls_points = walls_points
self._walls_normal = walls_normal
self._walls_up_vector = walls_up_vector
self._patches_points = patches_points
self._patches_normal = patches_normal
self._patch_size = patch_size
self._n_patches = n_patches
self._patch_to_wall_ids = patch_to_wall_ids
self._n_bins = None
self._frequencies = None
self._sources = None
self._receivers = None
self._absorption = None
self._air_attenuation = None
[docs]
@classmethod
def from_polygon(
cls, polygon_list, patch_size):
"""Create a Radiosity object for directional scattering coefficients.
Parameters
----------
polygon_list : list[PatchesDirectional]
list of patches
patch_size : float
maximal patch size in meters.
"""
# save wall information
walls_points = np.array([p.pts for p in polygon_list])
walls_normal = np.array([p.normal for p in polygon_list])
walls_up_vector = np.array([p.up_vector for p in polygon_list])
# create patches
(
patches_points, patches_normal,
n_patches, patch_to_wall_ids) = geometry.process_patches(
walls_points, walls_normal, patch_size, len(polygon_list))
# create radiosity object
return cls(
walls_points, walls_normal, walls_up_vector,
patches_points, patches_normal, patch_size, n_patches,
patch_to_wall_ids)
[docs]
def bake_geometry(self, ff_method='universal', algorithm='order'):
"""Bake the geometry by calculating all the form factors.
Parameters
----------
ff_method : str, optional
defines the form factor calculation method, by default 'universal'.
No other methods are implemented yet.
algorithm : str, optional
defines the algorithm used for Radiosity, by default 'order'.
No other algorithms are implemented yet.
"""
# Check the visibility between patches.
self._visibility_matrix = geometry.check_visibility(
self.patches_center, self.patches_normal, self.patches_points)
n_combinations = np.sum(self.visibility_matrix)
visible_patches = np.empty((n_combinations, 2), dtype=np.int32)
i_counter = 0
for i_source in range(self.n_patches):
for i_receiver in range(self.n_patches):
if not self.visibility_matrix[i_source, i_receiver]:
continue
visible_patches[i_counter, 0] = i_source
visible_patches[i_counter, 1] = i_receiver
i_counter += 1
self._visible_patches = visible_patches
if ff_method == 'universal':
self._form_factors = form_factor.universal(
self.patches_points, self.patches_normal,
self.patches_area, self._visible_patches)
else:
raise NotImplementedError()
# Calculate the form factors with directivity.
if self._sources is not None:
sources_array = np.array([s.cartesian for s in self._sources])
receivers_array = np.array([s.cartesian for s in self._receivers])
scattering_index = np.array(self._scattering_index)
scattering = np.array(self._scattering)
else:
sources_array = None
receivers_array = None
scattering_index = None
scattering = None
if self._absorption is None:
absorption = None
absorption_index = None
else:
absorption = np.atleast_2d(np.array(self._absorption))
absorption_index = np.array(self._absorption_index)
n_bins = 1 if self._n_bins is None else self._n_bins
if algorithm == 'order':
self._form_factors_tilde = \
form_factor._form_factors_with_directivity_dim(
self.visibility_matrix, self.form_factors, n_bins,
self.patches_center, self.patches_area,
self._air_attenuation, absorption,
absorption_index,
self._patch_to_wall_ids, scattering,
scattering_index,
sources_array, receivers_array)
else:
raise NotImplementedError()
[docs]
def init_source_energy(
self, source_position:np.ndarray,
ff_method="universal", algorithm='order'):
"""Initializes the energy of the source at the patches.
Note that the geometry must be baked before calling this function.
The source is assumed visible to all patches.
Parameters
----------
source_position : np.ndarray
_description_
ff_method : str, optional
defines the form factor calculation method, by default 'universal'.
No other methods are implemented yet.
algorithm : str, optional
defines the algorithm used for Radiosity, by default 'order'.
No other algorithms are implemented yet.
"""
source_position = np.array(source_position)
patch_to_wall_ids = self._patch_to_wall_ids
absorption = np.atleast_2d(np.array(self._absorption))
absorption_index = self._absorption_index
sources = np.array([s.cartesian for s in self._sources])
receivers = np.array([s.cartesian for s in self._receivers])
scattering = np.array(self._scattering)
scattering_index = self._scattering_index
patches_center = self.patches_center
if ff_method == "universal":
energy_0, distance_0 = source_energy._init_energy_universal(
source_position, patches_center, self.patches_points,
self._air_attenuation, self.n_bins)
else:
raise NotImplementedError()
if algorithm == 'order':
energy_0_dir = ee_order._add_directional(
energy_0, source_position,
patches_center, self._n_bins, patch_to_wall_ids,
absorption, absorption_index,
sources, receivers,
scattering, scattering_index)
self.energy_0_dir = energy_0_dir
self.distance_0 = distance_0
else:
raise NotImplementedError()
[docs]
def calculate_energy_exchange(
self, speed_of_sound,
histogram_time_resolution,
histogram_length, algorithm='order',
max_depth=-1, recalculate=False):
"""Calculate the energy exchange between patches.
Parameters
----------
speed_of_sound : float
speed of sound in m/s.
histogram_time_resolution : float
time resolution of the histogram, e.g. 1/1000.
histogram_length : float
duration of the histogram in seconds.
algorithm : str, optional
defines the algorithm used for Radiosity, by default 'order'.
No other algorithms are implemented yet.
max_depth : int, optional
maximum reflection order the be calculated.
e.g. 0 means the first order reflection, by default -1
recalculate : bool, optional
whether or not to recalculate energy propagation from source
position, by default False.
"""
n_samples = int(histogram_length/histogram_time_resolution)
patches_center = self.patches_center
distance_0 = self.distance_0
n_patches = self.n_patches
distance_i_j = np.empty((n_patches, n_patches))
for i in range(n_patches):
for j in range(n_patches):
distance_i_j[i, j] = np.linalg.norm(
patches_center[i, :]-patches_center[j, :])
if algorithm == 'order':
energy_0_dir = self.energy_0_dir
if not hasattr(self, 'E_matrix_total') or recalculate:
if max_depth < 1:
self.E_matrix_total = ee_order.energy_exchange_init_energy(
n_samples, energy_0_dir, distance_0,
speed_of_sound, histogram_time_resolution,
)
else:
self.E_matrix_total = ee_order.energy_exchange(
n_samples, energy_0_dir, distance_0, distance_i_j,
self._form_factors_tilde,
speed_of_sound, histogram_time_resolution, max_depth,
self._visible_patches)
else:
raise NotImplementedError()
[docs]
def collect_receiver_energy(self, receiver_pos,
speed_of_sound, histogram_time_resolution,
method="universal", propagation_fx=False):
"""Collect patch histograms as detected by receiver.
The receiver is assumed visible to all patches.
Parameters
----------
receiver_pos : np.ndarray
receiver positions in cartesian coordinates of
shape (n_receivers, 3).
speed_of_sound : float
speed of sound in m/s.
histogram_time_resolution : float
time resolution of the histogram, e.g. 1/1000.
method : str, optional
defines the algorithm used for Radiosity, by default 'order'.
No other algorithms are implemented yet.
propagation_fx : bool, optional
if true the sound propagation from patch to receiver is applied
in the final histogram, by default False.
Returns
-------
histogram_out : np.ndarray
the final energy histogram of shape
(n_receivers, n_patches, n_bins, n_samples)
"""
air_attenuation = self._air_attenuation
patches_points = self._patches_points
n_patches = self.n_patches
n_bins = self.n_bins
receiver_pos = np.atleast_2d(receiver_pos)
n_receivers = receiver_pos.shape[0]
patches_center = self.patches_center
patches_receiver_distance = np.empty([n_receivers,
self.n_patches,patches_center.shape[-1]])
E_matrix = np.empty((n_patches, n_bins, self.E_matrix_total.shape[-1]))
histogram_out = np.empty(
(n_receivers, n_patches, n_bins, self.E_matrix_total.shape[-1]) )
for i in range(n_receivers):
patches_receiver_distance = patches_center - receiver_pos[i]
if method == "universal":
# geometrical weighting
patch_receiver_energy = receiver_energy._universal(
receiver_pos[i], patches_points)
else:
raise NotImplementedError()
# access histograms with correct scattering weighting
receivers_array = np.array([s.cartesian for s in self._receivers])
receiver_idx = geometry.get_scattering_data_receiver_index(
patches_center, receiver_pos[i], receivers_array,
self._patch_to_wall_ids)
assert receiver_idx.shape[0] == self.n_patches
assert len(receiver_idx.shape) == 1
for k in range(n_patches):
E_matrix[k,:]= (self.E_matrix_total[k,receiver_idx[k],:]
* patch_receiver_energy[k])
if propagation_fx:
# accumulate the patch energies towards the receiver
# with atmospheric effects (delay, atmospheric attenuation)
histogram_out[i] = ee_order._collect_receiver_energy(
E_matrix,
np.linalg.norm(patches_receiver_distance, axis=1),
speed_of_sound, histogram_time_resolution,
air_attenuation=air_attenuation)
else:
histogram_out[i] = E_matrix
return histogram_out
[docs]
def set_wall_absorption(self, wall_indexes, absorption:pf.FrequencyData):
"""Set the wall absorption.
If this is not called before baking, a default value of 0.0 is used.
Parameters
----------
wall_indexes : list[int]
list of walls for the scattering data
absorption : pyfar.FrequencyData
absorption coefficient of cshape (1, )
"""
self._check_set_frequency(absorption.frequencies)
if self._absorption is None:
self._absorption_index = np.empty((self.n_walls), dtype=np.int64)
self._absorption_index.fill(-1)
self._absorption = []
self._absorption.append(np.atleast_1d(absorption.freq.squeeze()))
self._absorption_index[wall_indexes] = len(self._absorption)-1
[docs]
def set_air_attenuation(self, air_attenuation:pf.FrequencyData):
"""Set air attenuation factor in Np/m.
If not called, no air attenuation is applied.
Parameters
----------
air_attenuation : pyfar.FrequencyData
Air attenuation factor in Np/m.
"""
self._check_set_frequency(air_attenuation.frequencies)
self._air_attenuation = np.atleast_1d(air_attenuation.freq.squeeze())
[docs]
def set_wall_scattering(
self, wall_indexes:list[int],
scattering:pf.FrequencyData, sources:pf.Coordinates,
receivers:pf.Coordinates):
"""Set the wall scattering (brdf).
See :py:mod:`~sparrowpy.brdf` for more information. Make sure the
absorption is either included in the brdf or set independently via
wall using
:py:func:`~sparrowpy.DirectionalRadiosityFast.set_wall_absorption`.
Parameters
----------
wall_indexes : list[int]
list of walls for the scattering data
scattering : pyfar.FrequencyData
scattering data of cshape (n_sources, n_receivers)
sources : pyfar.Coordinates
source coordinates
receivers : pyfar.Coordinates
receiver coordinates
"""
assert (sources.z >= 0).all(), \
"Sources must be in the positive half space"
assert (receivers.z >= 0).all(), \
"Receivers must be in the positive half space"
self._check_set_frequency(scattering.frequencies)
if self._sources is None:
self._sources = np.empty((self.n_walls), dtype=pf.Coordinates)
self._receivers = np.empty((self.n_walls), dtype=pf.Coordinates)
self._scattering_index = np.empty((self.n_walls), dtype=np.int64)
self._scattering_index.fill(-1)
self._scattering = []
for i in wall_indexes:
sources_rot, receivers_rot = _rotate_coords_to_normal(
self.walls_normal[i], self.walls_up_vector[i],
sources, receivers)
self._sources[i] = sources_rot
self._receivers[i] = receivers_rot
self._scattering.append(scattering.freq*np.pi)
self._scattering_index[wall_indexes] = len(self._scattering)-1
def _check_set_frequency(self, frequencies:np.ndarray):
"""Check if the frequency data matches the radiosity object."""
if self._n_bins is None:
self._n_bins = frequencies.size
else:
assert self._n_bins == frequencies.size, \
"Number of bins do not match"
if self._frequencies is None:
self._frequencies = frequencies
else:
assert (self._frequencies == frequencies).all(), \
"Frequencies do not match"
@property
def n_bins(self):
"""Return the number of frequency bins."""
return self._n_bins
@property
def n_walls(self):
"""Return the number of walls."""
return self._walls_points.shape[0]
@property
def n_patches(self):
"""Return the total number of patches."""
return self._n_patches
@property
def form_factors(self):
"""Return the form factor."""
return self._form_factors
@property
def visibility_matrix(self):
"""Return the visibility matrix."""
return self._visibility_matrix
@property
def walls_area(self):
"""Return the area of the walls."""
return geometry._calculate_area(self._walls_points)
@property
def walls_points(self):
"""Return the points of the walls."""
return self._walls_points
@property
def walls_normal(self):
"""Return the normal of the walls."""
return self._walls_normal
@property
def walls_center(self):
"""Return the center of the walls."""
return geometry._calculate_center(self._walls_points)
@property
def walls_up_vector(self):
"""Return the up vector of the walls."""
return self._walls_up_vector
@property
def patches_area(self):
"""Return the area of the patches."""
return geometry._calculate_area(self._patches_points)
@property
def patches_center(self):
"""Return the center of the patches."""
return geometry._calculate_center(self._patches_points)
@property
def patches_size(self):
"""Return the size of the patches."""
return geometry._calculate_size(self._patches_points)
@property
def patches_points(self):
"""Return the points of the patches."""
return self._patches_points
@property
def patches_normal(self):
"""Return the normal of the patches."""
return self._patches_normal
@property
def patch_size(self):
"""Return the size of the patches."""
return self._patch_size
@property
def speed_of_sound(self):
"""Return the speed of sound in m/s."""
return self._speed_of_sound
def _rotate_coords_to_normal(
wall_normal:np.ndarray, wall_up_vector:np.ndarray,
sources:pf.Coordinates, receivers:pf.Coordinates):
"""Rotate the coordinates to the normal vector."""
o1 = pf.Orientations.from_view_up(
wall_normal, wall_up_vector)
o2 = pf.Orientations.from_view_up([0, 0, 1], [1, 0, 0])
o_diff = o1.inv()*o2
euler = o_diff.as_euler('xyz', True).flatten()
receivers_cp = receivers.copy()
receivers_cp.rotate('xyz', euler)
receivers_cp.radius = 1
sources_cp = sources.copy()
sources_cp.rotate('xyz', euler)
sources_cp.radius = 1
return sources_cp, receivers_cp
