{'fraction_non_alphanumeric': 0.032595706102938146, 'fraction_numerical': 0.011579511270449582, 'mean_word_length': 4.686035613870666, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural Implicit Surface Reconstruction using Imaging Sonar

We present a technique for dense 3D reconstruction of objects using an imaging sonar, also known as forward-looking sonar (FLS). Compared to previous methods that model the scene geometry as point clouds or volumetric grids, we represent the geometry as a neural implicit function. Additionally, given such a representation, we use a differentiable volumetric renderer that models the propagation of acoustic waves to synthesize imaging sonar measurements. We perform experiments on real and synthetic datasets and show that our algorithm reconstructs high-fidelity surface geometry from multi-view FLS images at much higher quality than was possible with previous techniques and without suffering from their associated memory overhead.


I. INTRODUCTION
Imaging or forward-looking sonar (FLS) is an extensively used sensor modality by Autonomous Underwater Vehicles (AUV). The key motivation for using FLS sensors is their ability to provide long-range measurements, unlike optical cameras whose range is severely limited in turbid water-a common situation in the field. Their versatility has resulted in their incorporation as a core sensor modality in applications including robotic path planning [1,2], localization [3]- [7], and the automation of tasks potentially dangerous or mundane for humans such as underwater inspection [8] and mapping [9]- [12].
FLS outputs 2D image measurements of 3D structures by emitting acoustic pulses and measuring the energy intensity of the reflected waves. While the sonar resolves azimuth and range, the elevation angle is ambiguous, and an object at a specific range and azimuth can be located anywhere along the elevation arc. Hence, the task of 3D reconstruction using FLS measurements can be equivalently framed as the task of resolving the elevation ambiguity from the image readings.
Existing algorithms for 3D reconstruction from FLS measurements can be grouped into geometry-based, physicsbased and, more recently, learning-based methods. However, most existing approaches either place restrictions on the robotic/sensor setup (elevation aperture, motion of the vehicle, etc.); rely on volumetric grids that are prohibitively expensive for large scenes or scenes with fine-grained geometry; or, specific to learning approaches, require the use of large labeled training sets that are difficult to collect in underwater environments.
To address these shortcomings, we approach the problem of underwater FLS-based 3D reconstruction through the lens of differentiable rendering and leverage the representational This work was partially supported by the Office of Naval Research grant N00014-21-1-2482 and National Science Foundation award 1730147. Ioannis Gkioulekas was supported by a Sloan Research Fellowship. 2 The code will be made available at https://github.com/ rpl-cmu/neusis power of neural networks to encode the imaged object as an implicit surface. Our overall reconstruction approach comprises the following components: • A differentiable volumetric renderer that models the propagation of acoustic spherical wavefronts. • A representation of 3D surfaces as zero-level sets of neural implicit functions. • A regularized rendering loss for 3D reconstruction using imaging sonars. To the best of our knowledge, this work is the first to introduce a physics-based volumetric renderer suitable for dense 3D acoustic reconstruction using wide-aperture imaging sonars. We evaluate our approach against different unsupervised methods on simulated and real-world datasets, and show that it outperforms previous state of the art. We will open-source our code together with different datasets 2 .

A. 3D Reconstruction Using Imaging Sonar
Different methods have been introduced to produce both sparse [5,10,11,13]- [15] and dense 3D reconstructions using FLS. The focus of this work is on dense object-level 3D reconstruction. A number of algorithms enforced assumptions or constraints on the physical system to obtain reliable 3D models. Teixeira et al. [16] successfully reconstructed a 3D map of a ship hull by leveraging probabilistic volumetric techniques to create submaps which are later aligned using Iterative Closest Point (ICP). However, the sonar aperture was set to 1 • and all detected points were assumed to lie on the zero-elevation plane which leads to reconstruction errors and prohibits extending the method to larger apertures. A line of work [17]- [20] uses two complementary sensors (imaging and profiling sonars) and performs sensor fusion to disambiguate the elevation angle. In our work, we restrict our setup to a single imaging sonar. Westman et al. [21] proposed a method to reconstruct specific points on surfaces (aka. Fermat Paths). However, it places constraints on the vehicle's motion as it needs a view ray perpendicular to the surface at each surface point and hence, requires a large number of images collected from specific views.
Another set of methods uses generative models to obtain dense 3D reconstructions. Aykin et al. [22], [23] attempt to estimate the elevation angle of each pixel by leveraging information from both object edges and shadows which restricts the object to be on the seafloor. Westman et al. [24] further extended the idea to do away with the seafloor assumption but still required estimates of object edges. Negahdaripour et al. [25] proposed an optimization-based algorithm to refine an initial 3D reconstruction obtained using space carving by encouraging consistency between the actual sonar images and the images produced by the generative model. However, generative methods generally rely on assumptions of the surface reflectivity proprieties and on 3D estimates of object edges which makes them impractical in real scenarios.
Various volumetric methods have also been proposed. Wang et al. [26] introduced an inverse sonar sensor model to update the occupancy in a voxel grid and later extended it to handle errors in pose estimates by aligning local submaps using graph optimization [27]. Although these methods, as probabilistic frameworks, can be more robust compared to space carving techniques [9,22], they consider each voxel independently and ignore inherent surface constraints. Guerneve et al. [28] frame the problem of 3D volumetric reconstruction as a blind deconvolution with a spatially varying kernel which can be reformulated as a constrained linear least squares objective. However, the method makes a linear approximation to the vertical aperture and places restrictions on the motion of the sonar limiting its practical application. Westman et al. [29] noted the equivalence of 3D sonar reconstruction to the problem of Non-Line-of-Sight (NLOS) imaging. It introduced a regularized least square objective and solved it using the alternating direction method of multipliers (ADMM). All aforementioned volumetric methods, however, share similar limitations since extracting highfidelity surfaces from volumetric grids is difficult. These approaches can also be computationally expensive for larger scenes or a fine discretization of volumes.
More recently, learning-based methods were proposed to resolve the elevation ambiguity. DeBortoli et al. [30] proposed a self-supervised training procedure to fine-tune a Convolutional Neural Network (CNN) trained on simulated data with ground truth elevation information. Wang et al. [31] use deep networks to transfer the acoustic view to a pseudo frontal view which was shown to help with estimating the elevation angle. However, these methods are limited to simple geometries or require collecting a larger dataset of real elevation data. Arnold et al. [32] propose training a CNN to predict the signed distance and direction to the nearest surface for each cell in a 3D grid. However, the method requires ground truth Truncated Signed Distance Field (TSDF) information which can be difficult to obtain.
In this work, we propose a physics-based renderer which uses raw FLS images and known sonar pose estimates to represent objects as zero-level sets of neural networks. It does not require hand-labeled data for training nor does it place restrictions on the setup or environment (voxel size, need for reflectance information, etc.)

B. Neural Implicit Representation
NeRF [33] introduced a volume rendering method to learn a density field aimed at novel view synthesis. It samples points along optical rays and predicts an output color which is then checked against that of the ground truth pixel. IDR [34] introduced a surface rendering technique that contrary to the volume rendering technique of NeRF, only considers a single point intersection on a surface. Hence, it fails to properly capture areas of abrupt changes in the scene. NeuS [35] leveraged the volume rendering technique of NeRF to perform 3D surface reconstructions and showed impressive results against state-of-the-art neural scene representation methods for scenes with severe occlusions and complex structures. NeTF [36] applied ideas from NeRF to the problem of NLOS imaging which was shown in [29] to have close similarity to FLS 3D reconstruction. All these methods focus on 3D reconstruction using optical sensors, either intensity or time-of-flight based. Our focus is on learning surfaces from acoustic sonar images.

A. Image Formation Model
An FLS 2D image I comprises pixels corresponding to discretized range and azimuth (r i , θ i ) bins. Each pixel value is proportional to the sum of acoustic energy from all reflecting points being the elevation angle (Fig. 1c). However, the elevation angle φ i is lost since each column θ i of an FLS image is the projection onto the z = 0 plane of a circular sector π i constrained to the sonar vertical aperture (φ min , φ max ) and containing the z axis (Fig. 1b). We now present our rendering equation. Imaging sonars are active sensors that emit a pulse of sound and measure the strength of the reflected acoustic energy. Let E e be the emitted acoustic energy by the sonar. Now, consider a unique infinitesimal reflecting patch P i "illuminated" by the acoustic wave and located on the arc A(φ) ∈ π i which passes through (r i , θ i , 0) (Fig. 2). The reflected acoustic energy at P i and received by the sonar can be approximated as: where 2 is the patch thickness, σ is the particle density at P i , and the factor 1 r 2 accounts for spherical spreading on both the transmit and receive paths. T is the transmittance, corresponding to exponential attenuation of a wave due to particle absorption-equivalently, the probability that the acoustic wave travels between two points unoccluded. We note that, when the sonar emitter and receiver are collocated, this probability is identical during the transmit (sonar to P i ) and return (P i to sonar) paths; thus, transmittance is accounted for only once for both paths. This is analogous to the effect of coherent backscattering in optical wave propagation with collocated emitter and receiver [38]. Now consider a surface composed of many such patches. The received energy by the sonar is simply the sum of the reflected energy by all patches {P i } ∈ A(φ) which approximate the surface. Hence, we arrive at the following image formation model: Note that although sound propagation through liquids is fundamentally different from that of light (longitudinal vs. transverse waves), different geometric acoustic modeling techniques still borrowed heavily from graphics and ray optics [39]. These methods fundamentally rely on Huygen's principle of sound travel through mediums which approximates the spherical wavefront as many energy-carrying particles travelling at the speed of sound. Hence, analogous to the concept of a light ray, we view an acoustic ray as the ray starting at the sonar acoustic center and ending at P i (Fig. 1a).

B. Rendering Procedure
Similarly to Yariv et al. [34] and Wang et al. [35], we represent the sonar-imaged surface using two multi-layer perceptrons (MLPs): a neural implicit surface, N, which maps a spatial coordinate x = (r, θ, φ) to its signed distance; and a neural renderer, M, which outputs the outgoing radiance at x.
Once the surface S is learned, we can extract it as the zero level set of N: To train our network using the rendering loss (Eq. 2), we leverage the following equation from Wang et al. [35] to estimate the value of the density σ(x) from the SDF: where Φ s (τ ) ≡ (1+e −sτ ) −1 is the sigmoid function used as a smooth approximator of the occupancy indicator function

C. Sampling Procedure
Existing work that targets optical cameras leverages ray optics where sampling points along a ray originating at some pixel is sufficient to approximate the rendering loss. On the contrary, our rendering loss in Eq. 2 requires producing point samples along the arc at p i = (r i , θ i ) as well as samples along each acoustic ray. To obtain a balanced dataset of zero and non-zero intensity samples when processing an image, we sample N P 1 random image pixels as well as N P 2 pixels with an intensity I(r i , θ i ) greater than a threshold. Let P be the set of sampled pixels. For each pixel p i ∈ P, we use stratified sampling to obtain samples along the arc. We discretize the elevation range [−φ min , φ max ] into N A equally spaced angles. Hence, the difference between two consecutive angles is ∆φ = φmax−φmin N A . We perturb these angles by adding N A randomly generated noise values ∼ Uniform(0, 1)∆φ to obtain a set of points A p = {P p = (r i , θ i , φ Pp ))} on the arc.
For each sampled point P p , we first construct the acoustic ray R Pp which starts at the acoustic center of the sonar and terminates at P p and then sample N R − 1 points along each ray. Specifically, we first sample N R − 1 range values r such that r < r and r = i r for some i > 0 ( r being the sonar range resolution). We obtain the set of points R Pp = {p = (r , θ, φ Pp )}. The points R Pp ∪ A p constitute a set of N R points along the ray (N R − 1 points along the ray + exactly 1 point on the arc). Finally, we perturb the range value of all points by adding uniformly distributed noise ∼ Uniform(0, 1) r (Fig. 3). Note that the points R Pp ∪ A p are expressed in spherical coordinates in the local sonar coordinate frame and hence need to be re-expressed in a global reference frame common to all sonar poses. We first transform the points to Cartesian coordinates: and then transform to world frame by multiplying with the sonar to world transform T sonar The resulting set of points expressed in world frame R W Pp ∪ A W p are used as inputs to the SDF neural network N. Finally, the direction of each ray is defined by the unit vector Let X be the set of all sampled points across all pixels, arcs and rays. This is the input batch to the neural network.

D. Discretized Image Formation Model
The discrete counterpart of the image formation model in Eq. 2 is:Î where: A p is the arc located at (r, θ), r Pp is the range of the disturbed point P p on the arc, M(P p ) is the predicted intensity at P p by the neural renderer, is the discrete opacity [35] at a point p i (p i and p i+1 being consecutive samples along the ray) which was further shown to equal: Finally, is the discrete transmittance value at P p (the endpoint of the ray). This is the product of one minus the opacity values α of all points on the acoustic ray excluding the α at P p .

E. Training Loss
Our loss function is constituted of three terms: the intensity loss in addition to eikonal and 1 regularization terms. The intensity loss encourages the predicted intensity to match the intensity of the raw input sonar images. The eikonal loss [40] is an implicit geometric regularization term used to regularize the SDF encouraging the network to produce smooth reconstructions. Finally, we draw inspiration from the NLOS volumetric albedo literature [29,41], and add the 1 loss term to help produce favorable 3D reconstructions when we use sonar images from a limited set of view directions. Hence, our final training loss term is:  Fig. 4: Our neural network architecture. The neural implicit representation N takes 3D spatial coordinates x as input and outputs their signed distance to the surface as well as a learned feature vector F(x). We use PyTorch's autodiff [42] to compute ∇N(x), the gradient of the signed distance at x.

IV. NETWORK ARCHITECTURE
We model N and M as two MLPs each with 4 hidden layers of size 64 (Fig 4). We additionally apply positional encoding to the input spatial coordinates and use weight normalization similar to IDR. While existing works that use optical cameras typically rely on larger networks to successfully learn high-frequency color and texture information, we found the proposed architecture to have sufficient capacity to learn different shapes from FLS images. Decreasing the size of the network was especially important to handle GPU memory overhead during training caused by the added sampling dimension (arcs).

V. EVALUATION
As our comparison metric, we use the mean and root mean square (RMS) Hausdorff distance defined as: M 1 and M 2 being respectively the ground truth (GT) and reconstructed meshes. We evaluate our method against Back-Projection (BP) and ADMM, two state-of-the-art optimization-based methods for unsupervised object-centric 3D reconstruction using imaging sonar. 1 BP is similar to the occupancy grid mapping method (OGM) as it uses the inverse sensor model to update the voxel occupancy while, however, ignoring the correlation between grid cells. We note that both ADMM and BP generate a density field F(σ). Hence, for each possible density σ (i.e., σ ∈ [0, 1]), we extract a surface using Marching Cubes and report the metrics based on the best σ value. The mesh quality generated by ADMM also depends on the regularization weight terms which we empirically tuned for each object. With our approach, extracting the zero-level set of N(x) directly generates a high-quality mesh. However, for the purpose of metric generation, we also try different level-sets near zero: We run our experiments on a system with an NVIDIA RTX Our network training time until convergence is ∼ 6 hours.

A. Simulation
We use HoloOcean [43,44], an underwater simulator to collect datasets of different objects of various shapes and sizes. We use the simulator's default noise parameters; namely a multiplicative noise w sm ∼ N (0, 0.15) and additive noise w sa ∼ R(0.2) (where R is the Rayleigh distribution and parameters are in units of normalized pixel intensity in the range [0, 1]). We also enable the simulation of multipath effects. The maximum range of the sonar was set to 8m. Before feeding the raw data to the three algorithms, we perform minimal filtering of speckle noise in the images by zeroing out pixels whose intensities are less than a threshold. After generating the meshes, we align them to the GT using ICP and report in Table I the mean and RMS Hausdorff distance to the GT for different objects and two sonar vertical apertures (14 • and 28 • ). Figure 5 shows example 3D reconstructions obtained using each algorithm. We see that our method produces more accurate reconstructions compared to the baselines in terms of 3D reconstruction accuracy and mesh coverage. The neural network implicit regularization combined with the eikonal loss favors learning smooth surfaces while avoiding bad local minimas when the input images potentially do not contain enough information to completely constrain and resolve the elevation of every 3D point in space.
For large objects (specified by an asterisk in the table), we decreased the grid voxel resolution of the baseline methods by one-half (increased the voxel size from the default value of 0.025m to 0.05m) to prevent the system from running out of memory (OOM): the ADMM and BP baselines do not leverage stochastic updates and hence, need to construct the optimization objective by processing all images in one go. This leads to memory overhead for larger objects, objects that require a fine discretization of the volume, or in the presence of a large number of non-zero pixel intensities 2 . In contrast, we train our renderer on a different subset of images in every iteration and use stochastic updates (the Adam optimizer) to optimize the function which significantly reduces memory requirements.  We evaluate our method on real-world datasets of a test structure (Fig. 6a) submerged in a test tank (Fig. 6b)  IMU and a Doppler Velocity Log (DVL) to provide accurate vehicle pose information (i.e., minimal drift for the duration of data capture). Fig. 8 shows the RMS and mean Hausdorff distance error of the three methods. The quality of the mesh generated by ADMM and BP depends on the selected marching cubes threshold σ. Hence, we plot the metrics generated using different σs and report the best value. With our method, we can extract the zero-level set of N directly alleviating the need for a postprocessing step for surface generation. Since the structure is submerged and lying at the bottom of the test tank (and hence, no sonar image captures the backside of the object), we limit the matching distance of the Hausdorff metric to 0.15m, 0.2m, and 0.25m for the 1 • , 14 • , and 28 • apertures respectively. We see that our method generates higher quality reconstructions especially when using larger apertures: With 14 • , our method achieves an (RMS, Mean)=(0.058m, 0.040m) while BP and ADMM are respectively at (0.077m, 0.063m) and (0.069m, 0.052m). Similarly for a 28 • aperture, our method achieves a lower (RMS, Mean) = (0.072m, 0.055m) compared to BP (0.104m, 0.079m) and ADMM (0.091m, 0.070m). Fig. 7 shows the resulting meshes for each method. While all three methods perform well with a 1 • aperture, the difference in reconstruction quality becomes visually more apparent as the aperture angle increases. With a 14 • aperture, we begin to lose the main feature of the object with ADMM and BP: the hole, short piling and crossbar are not easily discernible. When the aperture is increased to 28 • , both baseline methods perform poorly: the hole, crossbar, and short piling are lost. On the other hand, our proposed method successfully captures the major components of the structure for all three different apertures (a base, two vertical pilings, and a crossbar).

VI. CONCLUSION AND FUTURE WORK
We proposed an approach for reconstructing 3D objects from imaging sonar which represents imaged surfaces as zerolevel sets of neural networks. We performed experiments on To easily compare against the baselines, we add the constant green dashed line to report our method's metrics. Note however that our algorithm does not depend on the σ values in the x axis. simulated and real datasets with different elevation apertures and showed that our method outperforms current state-ofthe-art techniques for unsupervised 3D reconstruction using FLS in terms of reconstruction accuracy. While existing volumetric methods can suffer from memory overhead as well as require a separate step to extract meshes from volumetric grids (a process often difficult and prone to error), our method allows for easy surface extraction from implicit representations and uses stochastic updates to lessen the computational requirements.
Our algorithm has some limitations, all of which create opportunities for future work. First, we currently focus on single-object reconstruction but plan to expand our method to large-scale reconstruction of marine environments at the scale of harbors by taking inspiration from techniques such as Block-Nerf [45]. Second, our method is currently mostly suited for offline 3D reconstructions but using techniques such as Instant-NGP [46] and Relu-Fields [47] can bring it to real-time performance needed for robotic navigation applications. Finally, all of our experiments now use only sonar but underwater robots are typically equipped with other sensors such as optical cameras. Hence, another interesting direction from future work is to fuse multi-modal sensor inputs (acoustic and optical) where, for example, optical cameras are used to obtain high resolution models of specific interest areas in the scene while a sonar, with longer range, is used elsewhere.