Multibeam fish detection filter
The following algorithm attempts to identify fish in a multibeam echogram as a Boolean result. It generates a bitmap to filter out background signal in multibeam and imaging sonars. It is inspired by the filter described in Balk et al. 2009.
The Window size (pings) setting governs the background noise estimation. In this example we set it to 17.
Use a Mask operator in conjunction with this filter and the original data. Then use the Multibeam Target Detection operator to detect fish targets.
List of operands
Operand 1: a Sound Metrics DIDSON or ARIS multibeam variable, or other multibeam variables.
Code
""" Echoview Code Operator source file ======================================================================== Created by Echoview(R) 10.0.218 on Monday, 29 April 2019 See the Echoview help file for Code operator documentation and examples. NumPy User Guide: https://docs.scipy.org/doc/numpy/user/ NumPy Reference: https://docs.scipy.org/doc/numpy/Reference/ SciPy Reference Guide: https://docs.scipy.org/doc/scipy/Reference/ Echoview(R) is a registered trademark of Echoview Software Pty Ltd. """ # Authorship information __author__ = "Echoview Software Pty Ltd. 2019." __disclaimer__ = ( "This example code is provided AS IS, without warranty of any " "kind, express or implied, including but not limited to the " "warranties of merchantability, fitness for a particular purpose " "and noninfringement. In no event shall Echoview Software Pty Ltd " "be liable for any claim, damages or other liability, arising " "from, out of or in connection with the use of this example code." ) __version__ = "1.0" # System Imports from typing import List # Libraries from echoview import OperatorBase, MeasurementType, OperandInput, Error import numpy as np from scipy import signal, ndimage class Operator(OperatorBase): """ Fish detection filter (background suppression) for DIDSON data ==================================================================== Generates a bitmap representing potential fish detections. Operands --------------------- * Operand 1 - DIDSON Multibeam Sv """ def result_type(self, input_types): return MeasurementType.MULTIBEAM_BOOLEAN def __init__(self): self.foreground_filter = np.full((7, 3, 3), 1.) # boxcar (range, beams, pings) self.foreground_filter /= np.sum(self.foreground_filter) self.threshold = 2.2 # Filter threshold in dB def eval(self, inputs: List[OperandInput]): input = inputs[0] data = np.concatenate([ping.data[..., np.newaxis] for ping in input.window_measurements], axis=2) # Generate a mask that reflects all valid values (i.e. not # nodata), this is then used to convert output to nodata # Note: That as signal.convolve doesn't support nodata values we # mask them to -999db (0 linearly) mask = np.isfinite(input.measurement.data) if input.measurement.type.is_db: data[~mask] = -999 # 0 linearly data = to_linear(data) else: data[~mask] = 0 # smooth data using the foreground_filter, this cuts out # the highest frequency noise and is considered the foreground foreground_linear = signal.convolve(data, self.foreground_filter, mode='same') foreground = to_db(foreground_linear[..., input.window_index]) # calculate the background signal, this is taken to be mean over # the ping/time axis background = to_db(np.mean(data, axis=(2))) # calculate the filtered signal, defined as the difference # between the foreground and background signals plus some offset filtered_data = foreground - background > self.threshold # no data values go to False since Echoview doesn't support # nodata for boolean data filtered_data[~mask] = False return filtered_data # helper functions THRESHOLD_LOWER_DB = -999 THRESHOLD_UPPER_DB = 999 THRESHOLD_UPPER_LINEAR = 10 ** (THRESHOLD_UPPER_DB / 10) def to_linear(data: np.ma.MaskedArray, out=None): """Converts the specified input to linear values""" # Equivalent to ``10 ** (data / 10)`` but more efficient due to # reusing same output array if out is not None and out.shape != data.shape: raise ValueError( 'out must be a NumPy array with the same shape as the input array') below_threshold = data <= THRESHOLD_LOWER_DB above_upper_threshold = data >= THRESHOLD_UPPER_DB out = np.divide(data, 10, out=out) out = np.power(10, out, out=out) out[below_threshold] = 0. out[above_upper_threshold] = THRESHOLD_UPPER_LINEAR return out def to_db(data: np.ma.MaskedArray): """Converts the specified linear input to dB values""" thresholded_values = data <= 0 out = 10 * np.ma.log10(data) out[thresholded_values] = THRESHOLD_LOWER_DB return out