Skip to content

Commit

Permalink
Finalizes the bugfix of the angle_function function (#332)
Browse files Browse the repository at this point in the history 
* Fixes a bug in the angle_function for azimuth/colatitude computation (initial PR #329)

---------

Co-authored-by: Fabio Di Marco <[email protected]>
  • Loading branch information
@fakufaku @fabiodimarco
fakufaku and fabiodimarco authored Nov 25, 2023
1 parent 52b4eee commit 58368d9
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 114 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,8 @@ Bugfix
@hrosseel
- Fixes a bug when setting the air absorption coefficients to custom values (#191),
adds a test, and more details in the doc
- Fixes a bug in the utilities.angle_function in the calculation of the colatitude (#329)
by @fabiodimarco


`0.7.3`_ - 2022-12-05
Expand Down
139 changes: 48 additions & 91 deletions pyroomacoustics/tests/test_angle_function.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,17 +10,28 @@
a1 = np.array([[0, 0, 1], [0, 0, 1], [0, 1, 1]])
a2 = np.array([[0], [0], [1]])
a3 = np.array([0, 0, 1]).T
a4 = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]])

b1 = np.array([[0], [0], [0]])
b2 = np.array([1, -1, 1]).T
b3 = np.array([1, 0, 0]).T


a1_b1 = np.array([[0, 0, pi / 4], [0, 0, pi / 4]])
a1_b2 = np.array([[3 * pi / 4, 3 * pi / 4, pi / 2], [3 * pi / 4, pi / 2, pi / 2]])
a1_b1 = np.array([[0, 0, pi / 4], [0, 0, np.arctan(np.sqrt(2.0))]])
a1_b2 = np.array(
[
[3 * pi / 4, 3 * pi / 4, pi / 2],
[pi / 2 + np.arctan(1.0 / np.sqrt(2)), pi / 2, pi / 2],
]
)
a2_b1 = np.array([[0], [0]])
a2_b2 = np.array([[3 * pi / 4], [pi / 2]])
a3_b1 = np.array([[0], [0]])
a3_b2 = np.array([[3 * pi / 4], [pi / 2]])
a4_b1 = np.array([[0, pi / 2, 0], [0, pi / 2, pi / 2]])
a4_b3 = np.array([[pi, 3 * pi / 4, 0], [pi / 4, pi / 2, 0]])
a2_b3 = np.array([[pi], [pi / 4]])
a3_b3 = np.array([[pi], [pi / 4]])


# for 2-D coordinates
Expand All @@ -40,92 +51,38 @@
c3_d2 = np.array([[0], [pi / 2]])


class TestAngleFunction(TestCase):
def test_set_3d(self):
self.assertTrue(angle_function(a1, b1).all() == a1_b1.all())
self.assertTrue(angle_function(a1, b2).all() == a1_b2.all())

def test_point_3d_1(self):
self.assertTrue(angle_function(a2, b1).all() == a2_b1.all())
self.assertTrue(angle_function(a2, b2).all() == a2_b2.all())

def test_point_3d_2(self):
self.assertTrue(angle_function(a3, b1).all() == a3_b1.all())
self.assertTrue(angle_function(a3, b2).all() == a3_b2.all())

def test_set_2d(self):
self.assertTrue(angle_function(c1, d1).all() == c1_d1.all())
self.assertTrue(angle_function(c1, d2).all() == c1_d2.all())

def test_point_2d_1(self):
self.assertTrue(angle_function(c2, d1).all() == c2_d1.all())
self.assertTrue(angle_function(c2, d2).all() == c2_d2.all())

def test_point_2d_2(self):
self.assertTrue(angle_function(c3, d1).all() == c3_d1.all())
self.assertTrue(angle_function(c3, d2).all() == c3_d2.all())


def find_error(type_coordinates):
if type_coordinates == "3-D":
print("-" * 40)
print("type_coordinates = 3-D")
print("-" * 40)

a_range = [a1, a2, a3]
b_range = [b1, b2]
a_b_range = [a1_b1, a1_b2, a2_b1, a2_b2, a3_b1, a3_b2]

a_b_index = 0
for a in a_range:
for b in b_range:
error_azimuth = (angle_function(a, b) - a_b_range[a_b_index])[0]
error_colatitude = (angle_function(a, b) - a_b_range[a_b_index])[1]

print("for points :")
print(a)
print(b)
print(
"error in azimuth calculation: {}".format(np.average(error_azimuth))
)
print(
"error in colatitude calculation: {}".format(
np.average(error_colatitude)
)
)
print()
a_b_index += 1

elif type_coordinates == "2-D":
print("-" * 40)
print("type_coordinates = 2-D")
print("-" * 40)

c_range = [c1, c2, c3]
d_range = [d1, d2]
c_d_range = [c1_d1, c1_d2, c2_d1, c2_d2, c3_d1, c3_d2]

c_d_index = 0
for c in c_range:
for d in d_range:
error_azimuth = (angle_function(c, d) - c_d_range[c_d_index])[0]
error_colatitude = (angle_function(c, d) - c_d_range[c_d_index])[1]

print("for points :")
print(c)
print(d)
print(
"error in azimuth calculation: {}".format(np.average(error_azimuth))
)
print(
"error in colatitude calculation: {}".format(
np.average(error_colatitude)
)
)
print()
c_d_index += 1


if __name__ == "__main__":
find_error("3-D")
find_error("2-D")
def test_set_3d():
assert np.allclose(angle_function(a1, b1), a1_b1)
assert np.allclose(angle_function(a1, b2), a1_b2)


def test_point_3d_1():
assert np.allclose(angle_function(a2, b1), a2_b1)
assert np.allclose(angle_function(a2, b2), a2_b2)
assert np.allclose(angle_function(a2, b3), a2_b3)


def test_point_3d_2():
assert np.allclose(angle_function(a3, b1), a3_b1)
assert np.allclose(angle_function(a3, b2), a3_b2)
assert np.allclose(angle_function(a3, b3), a3_b3)


def test_point_3d_3():
assert np.allclose(angle_function(a4, b1), a4_b1)
assert np.allclose(angle_function(a4, b3), a4_b3)


def test_set_2d():
assert np.allclose(angle_function(c1, d1), c1_d1)
assert np.allclose(angle_function(c1, d2), c1_d2)


def test_point_2d_1():
assert np.allclose(angle_function(c2, d1), c2_d1)
assert np.allclose(angle_function(c2, d2), c2_d2)


def test_point_2d_2():
assert np.allclose(angle_function(c3, d1), c3_d1)
assert np.allclose(angle_function(c3, d2), c3_d2)
42 changes: 19 additions & 23 deletions pyroomacoustics/utilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
from scipy import signal
from scipy.io import wavfile

from .doa import cart2spher
from .parameters import constants, eps
from .sync import correlate

Expand Down Expand Up @@ -774,7 +775,8 @@ def all_combinations(lst1, lst2):

def angle_function(s1, v2):
"""
Compute azimuth and colatitude angles in radians for a given set of points `s1` and a singular point `v2`.
Compute azimuth and colatitude angles in radians for a given set of points `s1`
with respect to a reference point `v2`.
Parameters
-----------
Expand All @@ -786,8 +788,9 @@ def angle_function(s1, v2):
Returns
-----------
numpy array
2×N numpy array with azimuth and colatitude angles in radians.
2×N numpy array with azimuth and colatitude angles in radians in the
first and second row, respectively.
If the input vectors are 2-D, the colatitude is always fixed to pi/2.
"""

if len(s1.shape) == 1:
Expand All @@ -797,26 +800,19 @@ def angle_function(s1, v2):

assert s1.shape[0] == v2.shape[0]

x_vals = s1[0]
y_vals = s1[1]
x2 = v2[0]
y2 = v2[1]

# colatitude calculation for 3-D coordinates
if s1.shape[0] == 3 and v2.shape[0] == 3:
z2 = v2[2]
z_vals = s1[2]

colatitude = np.arctan2(
((x_vals - x2) ** 2 + (y_vals - y2) ** 2) ** 1 / 2, (z_vals - z2)
)
ndim = s1.shape[0]
if ndim == 2:
s1 = np.concatenate((s1, np.zeros((1, s1.shape[1]))), axis=0)
v2 = np.concatenate((v2, np.zeros((1, v2.shape[1]))), axis=0)

# colatitude calculation for 2-D coordinates
elif s1.shape[0] == 2 and v2.shape[0] == 2:
num_points = s1.shape[1]
colatitude = np.ones(num_points) * np.pi / 2
# this is slightly wasteful for 2d points, but is safer as we are relying
# on tried and tested code
az, co, r = cart2spher(s1 - v2)

# azimuth calculation (same for 2-D and 3-D)
azimuth = np.arctan2((y_vals - y2), (x_vals - x2))
if ndim == 2:
# this is only necessary to handle correctly the case s1 - v2 = 0
# in this case cart2spher returns zero, but we would like to have
# colatitude = pi/2 for consistency
co[:] = np.pi / 2

return np.vstack((azimuth, colatitude))
return np.vstack((az, co))

0 comments on commit 58368d9

Please sign in to comment.