Feat 980 msrp phat

opensoundscape/localization.py

@@ -9,7 +9,7 @@
 from opensoundscape import audio
 from opensoundscape.utils import cast_np_to_native
 from scipy.optimize import least_squares
-
+from scipy.spatial import ConvexHull
 
 # define defaults for physical constants
 SPEED_OF_SOUND = 343  # default value in meters per second
@@ -87,6 +87,7 @@ def __init__(
 
         # Initialize attributes
         self.tdoas = None  # time delay at each receiver
+        self.cc_matrix = None  # cross correlation matrix
         self.cc_maxs = None  # max of cross correlation for each time delay
         self.location_estimate = None  # cartesian location estimate in meters
         self.distance_residuals = None  # distance residuals in meters
@@ -341,6 +342,10 @@ def _localize_after_cross_correlation(self, localization_algorithm):
 
         return self.location_estimate
 
+    @property
+    def cc_maxs(self):
+        return self._cc_maxs
+
     @classmethod
     def from_dict(cls, dictionary):
         """Recover SpatialEvent from dictionary, eg loaded from json"""
@@ -443,6 +448,108 @@ def __call__(self, file, start_time):
         )
 
 
+class SpatialGrid:
+    """
+    Class for creating a grid of points for localizing sound events with methods that use a grid search approach.
+    """
+
+    def __init__(
+        self,
+        recorder_positions,
+        sample_rate,
+        resolution=1,
+        margin=0,
+        speed_of_sound=SPEED_OF_SOUND,
+    ):
+        """
+        Initialize a SpatialGrid object
+
+        Args:
+            recorder_positions: list of [x,y] or [x,y,z] positions of each recorder in meters
+            sample_rate: sample rate of the audio in Hz.
+            resolution: resolution of the grid in meters. Default is 1.
+            margin: margin around the convex hull of the grid in meters. Will only attempt to localize events that are inside the grid + margin.
+                    A negative margin will shrink the grid. Default is 0.
+            speed_of_sound: speed of sound in meters per second. Default is 343 m/s.
+
+
+        """
+        self.recorder_positions = np.array(recorder_positions)
+        self.sample_rate = sample_rate
+        self.resolution = resolution
+        self.margin = margin
+        self.dimensions = self.recorder_positions.shape[1]
+        self.speed_of_sound = speed_of_sound
+        self.grid = self._make_grid()
+        self.tdoa_grid = self._make_tdoa_grid()
+
+    def _make_grid(self):
+        """
+        Create a grid of points for localizing sound events
+
+        Returns:
+            grid: a list of [x,y] or [x,y,z] positions of each point in the grid
+        """
+
+        # make a grid of all the points between the min and max possible coordinates of the recorder positions
+        x = np.arange(
+            np.floor(np.min(self.recorder_positions[:, 0])) - self.margin,
+            np.ceil(np.max(self.recorder_positions[:, 0])) + self.margin,
+            self.resolution,
+        )
+        y = np.arange(
+            np.floor(np.min(self.recorder_positions[:, 1])) - self.margin,
+            np.ceil(np.max(self.recorder_positions[:, 1])) + self.margin,
+            self.resolution,
+        )
+
+        if self.dimensions == 2:
+            grid = np.array(np.meshgrid(x, y)).T.reshape(-1, 2)
+        else:
+            z = np.arange(
+                np.floor(np.min(self.recorder_positions[:, 2])) - self.margin,
+                np.ceil(np.max(self.recorder_positions[:, 2])) + self.margin,
+                self.resolution,
+            )
+            grid = np.array(np.meshgrid(x, y, z)).T.reshape(-1, 3)
+
+        hull = ConvexHull(self.recorder_positions)
+        # only keep the points that are inside the convex hull of the recorder positions
+        # self.hull is the ConvexHull object of the recorder positions
+        # self.hull.equations is the equation of the hyperplane of each face of the convex hull
+        # apply the equation of each face to the grid points to check if they are inside the convex hull
+        eps = 1e-6
+        # add a small epsilon to the margin to ensure that points on the edge of the convex hull are included
+        mask = np.all(
+            hull.equations[:, :-1].dot(grid.T) + hull.equations[:, -1][:, None]
+            <= self.margin + eps,
+            axis=0,
+        )
+        grid = grid[mask]
+
+        return grid
+
+    def _make_tdoa_grid(self):
+        """
+        Create an array containing the expected TDOAs for every recorder, for every recorder position in the grid.
+
+        Returns:
+            delays_grid: a grid, where each point has the matrix of TDOAs for each recorder position
+        """
+        # calculate the time delays for each point in the grid
+        delays_grid = []
+        for point in self.grid:
+            delays = (
+                np.linalg.norm(self.recorder_positions - point, axis=1)
+                / self.speed_of_sound
+            )
+            # now make a matrix of all the relative delays between each pair of recorders
+            # this is a matrix where the i,j-th element is the delay between recorder i and recorder j. With recorder i as the reference.
+            delays = delays - delays[:, None]
+            delays_grid.append(delays)
+        return np.array(delays_grid)
+
+
 class SynchronizedRecorderArray:
     """
     Class with utilities for localizing sound events from array of recorders

tests/test_localization.py

@@ -633,3 +633,11 @@ def test_df_to_events(LOCA_2021_aru_coords, LOCA_2021_detections):
         assert (event.tdoas == recovered_events[i].tdoas).all()
         assert (event.location_estimate == recovered_events[i].location_estimate).all()
         assert (event.cc_maxs == recovered_events[i].cc_maxs).all()
+
+
+def test_SpatialGrid():
+    # Verify that the SpatialGrid class creates the correct grid
+    receivers = np.array([[0, 0], [0, 10], [10, 0], [10, 10]])
+    grid = localization.SpatialGrid(receivers, resolution=1, margin=0)
+    # verify there are the right number of grid points
+    assert len(grid.grid) == 100