SlingBAG: point cloud-based iterative algorithm for large-scale 3D photoacoustic imaging

Domain: physics
Published: 2025-12-06
Editorial outcome: major_revision

Abstract
Large-scale 3D photoacoustic imaging has become increasingly important for both clinical and pre-clinical applications. Limited by cost and system complexity, only systems with sparsely-distributed sensors can be widely implemented, which necessitates advanced reconstruction algorithms to reduce artifacts. However, the high computing memory and time consumption of traditional iterative reconstruction (IR) algorithms is practically unacceptable for large-scale 3D photoacoustic imaging. Here, we propose a point cloud-based IR algorithm that reduces memory consumption by several orders, wherein the 3D photoacoustic scene is modeled as a series of Gaussian-distributed spherical sources stored in form of point cloud. During the IR process, not only are properties of each Gaussian source, including its peak intensity (initial pressure value), standard deviation (size) and mean (position) continuously optimized, but also each Gaussian source itself adaptively undergoes destroying, splitting, and duplication along the gradient direction. This method, named SlingBAG, the sliding Gaussian ball adaptive growth algorithm, enables high-quality large-scale 3D photoacoustic reconstruction with fast iteration and extremely low memory usage. We validated the SlingBAG algorithm in both simulation study and in vivo animal experiments.

Introduction
Photoacoustic imaging (PAI), which uniquely combines ultrasound detection and optical absorption contrast, is the only non-invasive optical imaging method that can image centimeters deep in living tissue with superior spatial resolution compared to other optical imaging techniques. PAI has been widely employed in various biomedical studies, including both pre-clinical studies and clinical applications1–6.

In recent years, significant progress has been made in the development of 3D PAI through the use of 2D matrix arrays, demonstrating immense potential for in vivo applications like peripheral vessel7,8, breast9,10 and small animal imaging11. Various array configurations have been explored, ranging from spherical and planar arrays12–19 to synthetic planar arrays that rely on scanning8,20. However, limited by cost and technical difficulties, it is a great challenge to fabricate a matrix array that satisfies the spatial Nyquist sampling requirement for real-time large-scale 3D PAI, which typically requires thousands of ultrasound transducer elements, as well as an equal number of parallel data acquisition channels. Therefore, real-time 3D PAI systems only have sparsely distributed sensors in reality, while it is well known that traditional back-projection image reconstruction methods will cause severe artifacts and low signal-to-noise ratio (SNR) for extremely sparse arrays.

To address these challenges, researchers have employed iterative reconstruction (IR) methods to improve image quality with sparse sensor configurations. An early IR approach to PAI was presented by ref. 21 where the differences between observed and simulated projections were iteratively minimized to improve image quality. Subsequent IR methods have incorporated various physical imaging factors. Deán-Ben et al.22 developed a model-based 3D reconstruction algorithm utilizing a sparse matrix forward model with reduced artifacts compared to back-projection methods. Wang et al.23,24 demonstrated that iterative penalized least-squares algorithms, based on discrete-to-discrete (D-D) imaging models using expansion functions in voxel grids, with either quadratic smoothing or total variation (TV) norm penalties, could enhance the performance of small animal 3D PAI systems. Similarly, ref. 25 introduced forward and backward operators based on k-space pseudospectral methods, and recently refs. 26,27 also use k-Wave28 as their forward model in model-based IR reconstruction for planar sensor arrays. Although these IR techniques have demonstrated great performance, for large-scale 3D PAI the spatial grid and system matrix scale is substantially increased, the computational load and memory demands make these methods less feasible.

Then, subsequent research has focused on more efficient model-based IR techniques. A notable advancement involves using semi-analytical forward models29–32, often based on integral geometry, combined with on-the-fly operator calculations within iterative solvers like least squares with QR factorization (LSQR) algorithm33. This strategy avoids the explicit storage of the massive system matrix, making reconstructions of larger datasets computationally tractable. However, these methods have their own limitations. The speed of semi-analytical models often stems from approximations (e.g., assuming far-field conditions and planar surface integral instead of spherical surface integral) that can compromise reconstruction accuracy, particularly for near-field or complex geometries. Besides, their reliance on solvers like LSQR restrict the use of powerful non-smooth regularization terms, such as the L1 or TV norms, which are highly effective for noise suppression and edge preservation. Furthermore, all these methods remain a uniform voxel-grid representation, inheriting its intrinsic trade-off between spatial resolution (grid size) and computational cost and its inefficiency in modeling sparse structures like vascular networks. These persistent limitations highlight the need for a paradigm shift away from traditional grid-based representations, motivating the exploration of alternative frameworks.

In addition to traditional IR, deep learning approaches for PAI reconstruction34–36 offer improvements in computational efficiency and image quality. However, they heavily rely on large, specialized training datasets, which can be difficult to obtain, and lack generalizability across different imaging systems or under varying conditions. This limits their broader applicability, particularly in diverse clinical settings where imaging conditions are not consistent.

Similar challenges of reconstructing 3D images from sparse-view also exist in the computer vision field. To tackle this challenge, differentiable rendering technologies such as NeRF37 and 3D Gaussian splatting (3DGS)38 were developed. Differentiable rendering technology39–43 enables self-supervised reconstruction of 3D scenes directly from scene observations, usually multi-view images. The photometric difference between the rendered image generated by the differentiable rendering pipeline and the actual captured image is utilized as a loss function, guiding the reconstruction of the entire scene through gradient back-propagation. NeRF, for instance, uses multi-layer perceptrons (MLPs) to map 3D spatial positions to colors and densities, while following methods like Instant NGP use hash grids to store features at different spatial scales, allowing for faster and more memory-efficient reconstructions44. Those pioneering works inspire more IR methods for PAI. A recent work reported NeRF-based paradigms to perform 2D PAI image reconstruction45. However, that method relies on dense sampling, which makes it computationally prohibitive for large-scale 3D PAI.

In contrast to NeRF, the 3D Gaussian splatting (3DGS)38 effectively addresses these computational challenges in 3D computer graphics. Based on point cloud instead of spatial grid, 3DGS applies object-order rendering by representing scenes as Gaussian ellipsoids, which are projected and composited onto a 2D screen using elliptical weighted average (EWA)46 for alpha blending. As a differentiable rendering strategy using explicit scene representation, 3DGS enables efficient, high-quality 3D reconstruction and real-time rendering.

Inspired by 3DGS, we propose the sliding Gaussian ball adaptive growth algorithm (SlingBAG) for 3D PAI reconstruction. Distinct from traditional voxel grids, this method stores information of the PA source in the form of a point cloud that undergoes iterative optimization, and the 3D PA scene is modeled as a series of Gaussian-distributed spherical sources. In order to fulfill the differentiable condition, we constructed a “differentiable rapid radiator" model based on the fact that a PA spherical source has the analytical solution for the wave equation in an acoustically non-viscous homogeneous medium. Each point cloud PA source is composed of multiple concentric spheres, with each sphere’s radius and initial PA pressure determined by the discretized Gaussian distribution. PA signals are the superposition of waves from spatially distributed Gaussian-distributed sources with specific peak intensity (pressure value), standard deviation (size), and mean value (spatial position). By minimizing the discrepancies between the predicted and the actual PA signals, we iteratively refine the point cloud and ultimately realize the 3D PAI reconstruction by converting the reconstruction result from a point cloud into a voxel grid.

This work introduces the SlingBAG algorithm, detailing its differentiable rapid radiator, point cloud optimization, and physics-based point cloud-to-voxel shader. Validated by both simulation and in vivo experiments, the proposed SlingBAG algorithm enables high-quality large-scale 3D PAI for broad potential applications with sparse sensor distributions.

Results
3D PA image reconstruction under an extremely sparse planar array

We first used simulated hand vessels as the target tissue to demonstrate the high performance of the SlingBAG algorithm. The simulated data, depicting human peripheral hand vessels, was used as the ground truth (Fig. 1a). At a grid point spacing of 0.2 mm, the total grid size of the imaging area is 400×520×200 (80 mm×104 mm×40 mm). To simulate 3D PA imaging, the initial dense planar matrix array has 122,500 sensors (350× 350, with a pitch of 0.4 mm) uniformly placed on a plane covering 140 mm×140 mm area, and the plane of the matrix array is about 15 mm below the imaging area (Fig. 1g). The simulation experiments were performed using k-Wave, with a sampling frequency of 40 MHz, 4900 sampling points, a medium density of 1×103 kg/m3, and a sound speed of 1500 m/s.

For sparse sensor arrangements, we undersampled the original 122,500 sensors into four different sparse matrix arrays with total sensor numbers of 49, 196, 576, and 4900, respectively. The pitches between adjacent sensors in these four different configurations are 20 mm, 10 mm, 6 mm, and 2 mm, respectively. We compared the proposed SlingBAG algorithm with the universal back-projection (UBP) algorithm47, and the model-based point-detector (MB-PD) IR algorithm augmented with L2 regularization, which uses a semi-analytical forward model with integral geometry29,30.

The reconstruction results are shown in Fig. 1b–f. For each vertical group of images, the three sub-images from top to bottom respectively show the top-view maximum amplitude projection (MAP), front-view MAP, and a cross-section slice along the green line marked in the top-view MAP of the 3D reconstruction result. Within each sub-figure, the three groups of images from left to right represent the reconstruction results of UBP, MB-PD, and SlingBAG under the same number of sensors, respectively. More 3D reconstruction results using both algorithms with different array configurations can be found in Supplementary Movies 1–3.

Under extremely sparse sensor configurations, the UBP reconstruction results have serious artifacts as shown in Fig. 1b, c, e. The reconstruction results of MB-PD are better than those of UBP; however, when the number of sensors is 49, 196, or 576, noticeable artifacts remain present. Only when the number of sensors reaches 4900 is artifact suppression significantly improved (Fig. 1f), but at this point, the advantage of MB-PD over UBP becomes less pronounced. The SlingBAG algorithm demonstrates superior 3D reconstruction capability; it even reconstructs a well-recognized outline of the hand vasculature with only 49 sensors (Fig. 1b). As the sensor matrix becomes denser, the quality of SlingBAG’s reconstruction improves rapidly. With 576 sensors (Fig. 1e), SlingBAG already achieves high-quality reconstruction; in contrast, the UBP and MB-PD results remain blurry due to reconstruction artifacts, making it difficult to discern fine vascular details (marked by the white dashed box). With a total of 4900 sensors (Fig. 1f), the SlingBAG reconstruction is nearly identical to the ground truth. Notably, for UBP reconstruction, even when using all 122,500 sensors (Fig. 1d), the result still cannot match the quality of SlingBAG’s result with only 4900 sensors, which can be easily observed from the corresponding slice images. This is because the finite size of the detection matrix array has the so-called “limited view” problem. However, the SlingBAG algorithm performs exceptionally well in limited-view conditions, which is another advantage of our method.

In Fig. 1h, we plotted the PA value distribution along the yellow dashed line in the slice images for comparison. The green solid line represents the ground truth, the red dashed line represents the reconstruction result from 4900 sensors, and the blue dashed line represents the reconstruction result from 576 sensors. As the number of sensors increases from 576 to 4900, the intensity distribution along the yellow dashed line in the slice progressively approximates the true value. The differences between the SlingBAG reconstruction results and the ground truth at several peaks are also substantially diminished with 4900 sensors, indicating that the SlingBAG algorithm’s reconstruction can achieve high quantitative accuracy with very sparse sensor distributions.

Next, to quantitatively assess the reconstruction results of the UBP, MB-PD and SlingBAG algorithms, we calculated the structural similarity index (SSIM), the mean squared error (MSE) and the peak signal-to-noise ratio (PSNR) of the top-view MAP and 3D volume for all four sensor numbers (Table 1). The SlingBAG achieved optimal performance across all evaluation metrics for each sensor number. As shown in Fig. 2, the trends and clear comparisons of all metrics for the reconstruction results of the three algorithms with respect to the number of sensors can be observed. We emphasize the following points here. First, for the 3D volume reconstruction results with 4900 sensors (Fig. 2e), the MSE of SlingBAG is only 0.000084, which is much lower than that of UBP and MB-PD, directly reflecting the high fidelity of SlingBAG’s 3D reconstructions. Second, for the reconstructions with 49, 196, and 576 sensors (Fig. 2c, f), the PSNR of SlingBAG is significantly higher than that of UBP and MB-PD; for all sensor numbers, SlingBAG achieves a PSNR above 20 dB for the top-view MAP and above 30 dB for the 3D volume, with the maximum reaching 40.78 dB, which demonstrates SlingBAG’s superior capability for artifact and noise suppression in highly sparse sensor arrays. Third, the performance metrics of SlingBAG with 196 sensors even surpass those of MB-PD with 576 sensors, indicating superior reconstruction capability among IR methods.

We further visualized the iterative results of our reconstruction algorithm’s point cloud for the 196-sensor array in Fig. 3. We first initialized a relatively dense point cloud distribution randomly within the reconstruction region. At the coarse reconstruction stage, both the position update and point duplication were disabled, only the peak initial pressure and standard deviations of each Gaussian ball underwent updates. Then, we got a preliminary convergence result in the coarse reconstruction stage. Next, we enabled the model’s position update and point cloud duplication features. It can be seen that the converged shape of the three finger vessels in the coarse reconstruction gradually becomes denser, along with finer details showing up (Fig. 3a). Figure 3b shows the 3D visualization results of the point cloud throughout the entire reconstruction process, from random initialization to the coarse reconstruction stage and finally to the fine reconstruction stage. Furthermore, through our point cloud-to-voxel-grid shader, we transformed the fine reconstruction result from point cloud form into a voxel grid (Fig. 3a). More detailed visualizations of the intermediate iterative process can be found in Supplementary Fig. 5.

3D PA image reconstruction of in vivo animal studies

Then, we performed image reconstruction of real in vivo animal experimental data using all three algorithms, the UBP, the MB-PD and the SlingBAG algorithms. The in vivo animal study data, including mouse brain, rat kidney and rat liver, was acquired by Professor Kim’s lab using a hemispherical ultrasound (US) transducer array with 1024 elements and a radius of 60 mm48. Each US transducer element in the array had an average center frequency of 2.02 MHz and a bandwidth of 54%. The effective field of view (FOV) was 12.8 mm×12.8 mm×12.8 mm, and the spatial resolutions of ~380 μm were nearly isotropic in all directions when all 1024 US transducer elements were used. All in vivo experiments were conducted with a sampling frequency of 8.33 MHz and 896 sampling points. For the mouse brain, the reconstruction region was 256×256×128; for the rat kidney and rat liver, the reconstruction region was 256×256×256. The resolution of the reconstructed 3D volume was 0.1 mm/pixel. More details regarding the 3D PAI system and animal experiments can be found in literatures49–51.

First, we compared three reconstruction results for a mouse brain as shown in Fig. 4. Figure 4a presents the reconstruction results by UBP, including the MAP results in XY Plane, XZ Plane, and YZ Plane, as well as a cross-section slice from top to bottom. Figure 4b, c presents counterpart results by the MB-PD and the SlingBAG algorithms. Figure 4e is a snapshot of the 3D SlingBAG reconstruction, which is available online (Supplementary Movie 4). According to the results, the SlingBAG reconstruction exhibits minimal artifacts and reveals much more vascular details compared with UBP and MB-PD. As in the UBP result, due to artifacts induced by sparse sensors, the fine vessels are almost indiscernible in the bottom cross-section slice. Although the MB-PD reconstruction reduced artifacts to some extent, small vessels are also blurred compared to the SlingBAG result. Further analysis showed that the contrast-to-noise ratio (CNR) of the XY Plane-MAP for the UBP and the MB-PD was 31.06 and 41.05, whereas it reached 46.39 for the SlingBAG (the selection of background and signal regions for CNR calculation can be found in Supplementary Fig. 8). We also demonstrated the SlingBAG algorithm’s superior capability for reconstruction of in vivo mouse brain under both limited-view and sparse-view conditions in Supplementary Figs. 2–4.

Next, in Figs. 5, 6, we presented the reconstruction results of a rat kidney and a rat liver, respectively. In the cross-section slice of the reconstructed kidney image in Fig. 5a, the UBP results show severe arc-shaped artifacts caused by sparse sensor configuration and limited view issues, which significantly impair the ability to resolve the fine vascular structures. The artifacts in the MB-PD results were reduced to some extent; however, obvious arc-shaped artifacts remained, impairing vascular clarity, and a slight decrease in vascular resolution was also observed (Fig. 5b). In contrast, the SlingBAG algorithm effectively suppresses these arc-shaped artifacts, producing a much clearer and cleaner vascular cross-sectional image (Fig. 5c). Similarly, the SlingBAG algorithm renders the liver vasculature with much better continuity and much fewer artifacts. The 3D reconstruction results in Figs. 5e, 6e are provided in Supplementary Movies 5, 6. For a quantitative comparison, the CNR values of the XY plane-MAP obtained with UBP reconstruction for the kidney and liver are 20.18 and 30.62, respectively; for MB-PD, they are 27.85 and 37.28, respectively; whereas for the SlingBAG results, the corresponding values are 30.26 and 42.01 (the selection of background and signal regions for CNR calculation can be found in Supplementary Fig. 8.). These results illustrate the exceptionally high quality of image reconstruction achieved by the SlingBAG algorithm.

Computing memory consumption analysis and comparison

Besides the quality of image reconstruction, the demand for computing memory usage is also key for practical implementation. Here, we compared the computing memory consumption of the SlingBAG algorithm with that of traditional IR methods under the same experimental setup.

For the PA imaging of hand vessels using a planar matrix array with 4900 sensors, the forward simulation combined with the backward gradient propagation of the SlingBAG algorithm required only 3.56 GB of GPU memory. In comparison, the theoretical memory consumption of traditional IR methods, such as k-Wave and the fast iterative shrinkage-thresholding algorithm (FISTA)52, is estimated to be as high as 19.32 GB, which exceeds the memory capacity limits of many consumer-grade GPUs. Similarly, for the in vivo animal imaging experiment using a hemispherical array, where the sensors are relatively far away from the reconstructed region, the SlingBAG algorithm required 2.99 GB of GPU memory, while the memory usage of traditional IR methods is estimated to be 85.87 GB, making it virtually impractical for most accessible consumer-level GPU resources. More details of the theoretical analysis for memory demand by traditional IR methods are provided in the Supplementary Note 4.

Methods
Differentiable rapid radiator based on multiorder spherical discretization

Inspired by differentiable rendering techniques from computer graphics, we propose the concept of “differentiable acoustic radiation" for PA imaging. Our model explicitly represents the 3D PA scene as a collection of isotropic Gaussian-distributed PA sources, stored as a point cloud (Fig. 7a). Each Gaussian source is characterized by three attributes: its peak initial pressure p0, standard deviation a0, and 3D coordinates (xs, ys, zs).

To build an efficient and differentiable forward model, we approximate each Gaussian source as a superposition of multiorder concentric homogeneous spheres. This discretization is based on the analytical solution for PA wave propagation from a uniform spherical source47. Specifically, we use a 10th-order spherical discretization following the 3σ rule of the Gaussian distribution (Fig. 7b). The exact radii and pressure weightings for this discretization are detailed in the Supplementary Note 10.

The analytical solution involves a non-differentiable Heaviside step function u(x). To enable gradient-based optimization, we approximate it using a high-precision, differentiable error function uϵ(x) = 1/2(1 + erf(ϵx)). This yields a fully differentiable forward propagation expression p(R, t) for a single Gaussian source, which is a linear superposition of the signals from the 10 discretized spheres:

p(R, t) = Σi=1^10 pi0/2 [uϵ(D + ai0) - uϵ(D - ai0)], (1)

where ai0 and pi0 are the radius and initial pressure of the i-th order sphere (detailed in Supplementary Note 10), R is the distance from the sphere center (xs, ys, zs) to the sensor, vs is the speed of sound, and t is time.

This formulation is fully differentiable with respect to all source attributes (p0, a0, xs, ys, zs). This allows us to analytically derive the gradients for end-to-end training via back-propagation. The complete derivations of the partial derivatives are provided in the Supplementary Note 11. Both forward and backward models are implemented with GPU acceleration using Taichi53 and CUDA.

Sliding Gaussian ball adaptive growth (SlingBAG) point cloud model

The overview of the SlingBAG reconstruction algorithm is shown in Fig. 7a. We first randomly initialize the Gaussian ball point cloud within the field of view. The model then iteratively updates all attributes by minimizing the L2 loss between the simulated signals from the differentiable rapid radiator and the experimentally acquired sensor data.

The optimization proceeds in two successive stages: a coarse reconstruction stage and a fine reconstruction stage. The coarse stage rapidly converges by updating only the standard deviation a0 and peak pressure p0 of each Gaussian ball. The fine stage updates all attributes, including the spatial coordinates (xs, ys, zs). A key feature of SlingBAG is the adaptive density optimization of the point cloud (Fig. 7c). During optimization, points are destroyed (pruned) if their pressure or size is below a threshold, saving computational resources. Conversely, points are split into smaller spheres if their size (a0) is too large, and duplicated in regions with high position gradients during the fine stage. This adaptive process dynamically adjusts the density of the point cloud, allowing the model to refine details precisely where needed. This facilitates high-quality reconstructions even from sparse sensor data. The detailed pseudocode, optimization hyperparameters, and stopping criteria for the algorithm are provided in the Supplementary Alg. 1.

Physically-based point cloud-to-voxel grid shader

The iterative reconstruction yields a 3D point cloud (Fig. 8b), which is not ideal for standard visualization. We therefore developed a “point cloud-to-voxel grid shader" to convert the final optimized point cloud data into a conventional 3D voxel grid, which is widely used for 3D visualization of PA imaging results.

The shader first initializes a 3D voxel grid of size (X, Y, Z) and aligns its coordinate space with the point cloud’s coordinate space. As illustrated in Fig. 8c, the shader then iterates over every voxel (vx, vy, vz) in the grid. For each voxel, its final intensity I(vx, vy, vz) is calculated by accumulating the pressure contributions from all N Gaussian balls in the final point cloud.

The contribution from a single Gaussian ball j (with center Pj = (Px, Py, Pz)j, peak pressure pj0, and standard deviation aj0) to a voxel at position v = (vx, vy, vz) is based on the tenth-order spherical discretization shown in Fig. 8a. The intensity is calculated as follows:

I(v) = Σj=1^N Σk=1^10 wjk I(kv - Pj k2 ≤ sk * aj0), (2)

where I(·) is the indicator function, which is 1 if the condition (the voxel center is within the k-th sphere) holds, and 0 otherwise. The distance between the voxel center and the point center is dj = ∥v - Pj∥2. The parameters wjk and sk represent the pressure weights and radius scales for the ten spheres.

This process effectively renders the physical pressure distribution of each Gaussian ball onto the 3D voxel grid. By accumulating these contributions, we achieve a comprehensive and detailed 3D visualization that accurately reflects the optimized point cloud data. The detailed pseudocode for this shader is available in the Supplementary Alg. 2.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In the work of Li et al, the authors proposed an image reconstruction method using the point cloud-based representation for\n3D photoacoustic imaging. This idea was translated from the Gaussian splatting technique in the field of computer vision\nwhich has been a research hotspot in this field in recent years. It is very interesting to see its successful translation to\nphotoacoustic imaging which is novel. The authors demonstrated that the SlingBAG algorithm enables impressively low\nmemory usage and fast iteration for high-quality large-scale 3D PA imaging with sparse sensor distribution, thanks to the\nsparsity characteristics in PA sources and the developed differentiable rapid radiator model for forward simulation. Overall,\nthe paper is interesting, clearly presented and is novel in methodology and has shown some\npotentials in addressing\nimportant challenges in large-scale 3D PA reconstruction. However, I think this works need further improvements to warrant\npublishing in Nature Communications. Particularly, the authors should compare their proposed method with current state-of-\nthe-art algorithms, the model-based iterative reconstruction methods, instead of only with the UBP algorithms. After all,\nmodel-based iterative reconstruction methods have been studied and developed for 3D PA imaging in the last ten years,\nand some of them were proposed to address the same problem just like this work, while the addressed problem size may be\na bit smaller. I recommend the following works [1-3] for the authors' reference.\nSpecific comments\n1, Details on how the PA signals of the hand vessel were simulated must be provided. I guess the PA signals were obtained\nwith the proposed differentiable rapid radiator model. While the basis of this model (Eq. 1) has been widely used for UBP-\nbased PA image reconstruction, it does not simulate real wave propagation, as well as the spatial and electrical response\nfunction of the receiver. The authors should at least acknowledge this limitation. Using k-wave to simulate the signals is a\nbetter choice, the authors are encouraged to demonstrate the capability of the rapid radiator model in modeling real PA\nsignals. The title 'real PA signal' of Supplementary Fig. 1a is misleading. It gives readers the impression that the PA signal is\nfrom experimental measurement, while this is not the fact.\n2, Page 4, line 153. \"Because current traditional IR methods cannot handle such large-scale cases, we compared the\nproposed SlingBAG algorithm with the widely-used universal back-projection (UBP) algorithm [44].\" I do not totally agree\nwith the authors that a comparison with SlingBAG and traditional IR method is unfeasible. Using k-wave as the forward\nmodel is not a good choice for large-scale PA image reconstruction problem. This has been well discussed by the authors in\nSupplementary Note 4. In ref [1], the reconstruction problem of 200x200x200 voxels has been presented, with 64 s\nreconstruction time on the GPU. Their later work further reduced the computational complexity of the algorithm for the\ninversion of large-scale problem for 3D PA imaging [2]. I think the key is that they adaptively used semi-analytical PA\nmodel\n[3] instead of wave-propagation models such k-Wave, and that the model-matrix is calculated\non-the-fly instead of pre-\nstored, so that memory consumption can be maintained in a reasonable range. The authors are encouraged the compare\ntheir method with the model-based iterative reconstruction method, at least for the reconstruction problem of a smaller and\nfeasible size, to fully demonstrate the advantage of SlingBAG over the state-of-the-art method. This open-source code\n(https://github.com/berkanlafci/pyoat) might be a good reference material for the authors to implement the model-based\niterative reconstruction method.\n3, Fig. 2, the reported metrics MAP and Slice SSIM are not defined. Is the slice SSIM calculated over the whole imaging\nvolume (for all the slices)?\n4, Fig. 4c the reference system is not properly defined, where are x, y and z?\n5, Comparing the cross-section slices in Fig. 5, it looks like the artifacts relative to limited-view problem were well-mitigated\nby SlingBAG. Could the authors discuss why the point-cloud representation can better work for addressing the limited-view\nproblem, while UBP and SlingBAG basically use the same forward model?\n6, Page 12, line 239-240, duplicated definition of the acronym CNR.\n7, the CNR of the images formed by SlingBAG look exceptionally high, how was the CNR computed? Is the CNR value over\n61.55 dB realistic considering the SNR of the imaging system?\n8, Page 12, line 264-272, The image reconstruction problem size for the experimental data is larger than the simulation data\n(1132 × 1132 × 760 vs 200 × 520 × 400), but the required GPU memory is smaller (2.99 GB vs 3.56 GB). Why is that? I\nassume the memory consumption is related to the number of the Gaussian balls. For the experimental data, should more\nGaussian balls be initialized for the experimental data? Please also specify how many Gaussian balls were initialized for the\nsimulated data and each experimental data, and how many balls were eventually retained after the reconstruction is done for\nall the cases. Is there a guideline on how many Gaussian balls need to be initialized for good image reconstruction quality?\nWould the number affect the image quality? These details will help interested readers to reproduce this work and improve its\nimpact.\n9, Page 13, line 278-279. \"This work proposes SlingBAG, a point cloud-based iterative 3D PAI reconstruction algorithm,\nwhich substantially increases speed with much less memory consumption. \" What does the \"speed\" represent? Image\nreconstruction speed? Compared to what method? There is no proper data presented in the manuscript to support this\nstatement.\n10, Page 13, line 308-312. \"It takes less than 4 hours for SlingBAG to complete 3D PA image reconstruction with 196\nsensors\". Contrasting SlingBAG with the method proposed in [1-2], the image reconstruction speed of SlingBAG is not very\nimpressive.\n\"The significant improvement in both reconstruction range and speed is largely attributable to the adaptive point cloud\nrepresentation, which inherently aligns better with the sparsity characteristics in PA sources. \" In the scenarios that not only\nthe vessels, but also other light observable materials are imaged, the sparse PA sources assumption might not be valid.\nWould SlingBAG also works efficiently in these scenarios?\n11, it seems that the rapid radiator model cannot account for the characteristics of the transducer. Would the frequency\nresponse and the spatial impulse response of the transducers affect the reconstruction?\n12, Can the authors provide more clues on how to determine the coefficients for the radius and pressure shown in Table 1.\nIt\nis understood that they resemble the spherical discretization shown in Fig 8a to approximate the Gaussian function. But how\nto find the specific value for the discretized presentation of the radius and pressure is unclear to me.\n13, Algorithms 1, please provide more details on the implementation of the SlingBAG algorithm. How to define the coarse\nreconstruction stage and fine reconstruction stage? What are the convergence conditions of each stage? Why are the\ngradients of loss function over a0 and p0 used firstly in the coarse stage? Have the authors carried out some sensitivity\nanalysis, i.e. the change in loss is more sensitive to the change in a0 and p0?\n14, Page 21, line 486-498. What is the criterion to determine \"excessively low\" or \"\"excessively large\" for discarding p0 and\na0? Can the authors add more details on spatially splitting? Does the spatial splitting only change the position of the new\nballs? What about the radius and pressure of the new balls? For the duplication feature, all the balls are duplicated every 40\niteration (Supplementary Fig. 5)? Then, why the Gaussian balls count at 160th iteration is less than that at 40th iteration? Is it\nbecause many were discarded according to the value of p0 and a0?\n15, it is a pity that Fig. 8c does not well illustrate the processing of point cloud to voxel grid shader. I do not understand how\nthe alignments of coordinate spaces were carried out and how to relate the values of p0 and a0 with the gray-scale intensity\nand brush size. After the gray-scale intensity and brush size are determined, how to apply the painting in the voxel\nrepresentation of the images? All these listed points need to be better addressed both in the text and in an updated Fig. 8c.\nReferences\n[1] L. Ding, X. L. Deán-Ben and D. Razansky, \"Efficient 3-D Model-Based Reconstruction Scheme for Arbitrary Optoacoustic\nAcquisition Geometries,\" in IEEE Transactions on Medical Imaging, vol. 36, no. 9, pp. 1858-1867, Sept. 2017", "reviewer_2": "Summary:\nThe manuscript presents the Sliding Gaussian Ball Adaptive Growth (SlingBAG) algorithm, an optimization-based\nreconstruction algorithm for large-scale 3D photoacoustic (PA) imaging. The algorithm models the to-be-reconstructed object\nas a collection of isotropic Gaussian blobs whose initial pressure values, sizes, and positions are iteratively optimized\nduring reconstruction. The primary motivation is to mitigate the computational time and GPU memory constraints associated\nwith certain optimization-based reconstruction approaches. In both presented studies, the reconstructed image quality of\nSlingBAG is compared to the universal backprojection (UBP) algorithm, with the authors reporting overall superior\nperformance of SlingBAG. The computational time for the forward simulation is assessed against k-Wave (pseudo-spectral\nmethod) in the first study, and GPU memory reduction is analyzed. However, while deploying adaptive point-cloud-based\nexpansion functions presents a promising approach to reducing computational load, the study lacks rigorous validation and\nfair comparisons with other optimization-based methods. Additionally, a comprehensive\ncomparative analysis against\nregularized optimization-based reconstruction methods is missing, leaving the claim of \"fast iteration\" and \"extremely low\nmemory usage\" insufficiently supported. Given these concerns, major revisions are necessary before the manuscript can be\nconsidered for publication.\nMajor comments:\n1. The choice of the UBP method as a reference method presents a significant weakness of the paper. The presented\nrational for this was that other iterative reconstruction methods “cannot handle such large scale cases” (line 154). Such a\nclaim is not true (see Comment #2). It is well known that the UBP method will not perform well under the considered\nconditions. Although conventional iterative reconstruction methods may be slower, they can still be implemented for the\npresented studies. The accuracy and computation time of the proposed method should be compared against standard\niterative reconstruction methods in the field of photoacoustic tomography.\n2. As a follow up to the above comment: the authors assume an acoustically homogeneous medium. Under such conditions,\nthere are a variety of integral geometry-based forward models that can be employed. These models can be parallelized in\ndifferent ways on GPUs to improve computational efficiency. Additionally, the action of these forward models can be\ncomputed on-the-fly, thereby eliminating the need to store a large system matrix. The presented research is not put in the\ncontext of these and other methods.\n3. The authors claim that SlingBAG achieves \"fast iteration\" and \"extremely low memory usage.\" Its computational\nperformance, in terms of execution time and GPU memory consumption, is benchmarked against the k-space\npseudospectral algorithm in k-Wave, demonstrating reductions in both aspects. However, the reported results are\nincomplete, as only partial performance metrics are provided. A comprehensive breakdown of reconstruction times across\ndifferent experiments is needed to substantiate these claims. Additionally, the authors should report the iteration count,\ntotal reconstruction time, and stopping criteria to clarify the computational demands of SlingBAG.\n4. The proposed method is an optimization-based reconstruction algorithm, which is often deployed for noisy and sparse-view\nimaging scenarios. However, the implementation of the method does not incorporate any regularization techniques,\nwhich are typically essential for stabilizing solutions in ill-posed inverse problems. The authors should explicitly clarify\nwhether their discretization approach is compatible with standard regularization methods. If regularization is unnecessary in\ntheir framework, a thorough justification and analysis should be provided to explain how the proposed method inherently\nensures stability and robustness in the presence of noise and limited data.\n5. Quantitative accuracy in photoacoustic (PA) imaging is essential, as functional properties are derived from the\nreconstructed initial pressure distributions. Optimization-based approaches are particularly valuable in this regard, as they\nenable precise initial value estimation. In Figure 1(l), it would be preferable to use true values for the y-axis or explicitly state\nthat the curves are scaled with the same parameter rather than individually. Additionally, incorporating evaluation metrics\nsensitive to global scaling, such as mean squared error (MSE), would provide a more comprehensive assessment. This\nwould supplement the current metrics (SSIM and CNR), which remain invariant under uniform scaling and may not fully\ncapture differences in absolute pressure values.\nOther comments:\n1. It would be beneficial to include mathematical equations in the point cloud to voxel conversion section to quantitatively\ndescribe the transformation process. This would enhance clarity and provide a more rigorous foundation for understanding\nhow the reconstructed pressure distributions are mapped onto a voxel grid."}, "author_response": "Authors’ Response\nfor the submitted paper,\n“Sliding Gaussian ball adaptive growth: point cloud-based iterative\nalgorithm for large-scale 3D photoacoustic imaging”\nWe would like to thank both the editor and the reviewers for your valuable comments and\nsuggestions, which have helped us to significantly improve our manuscript. We have carefully\nrevised our manuscript accordingly, and thoroughly addressed all comments point by point in\ndetail, including providing additional experiments. The key revisions include:\n1. Achieving faster reconstruction speed: We substantially improved the iterative\nreconstruction speed of SlingBAG by implementing shared memory access techniques,\nachieving doubling in speed.\n2. Providing comparative experiments with state-of-art IR method: As suggested by\nReviewers 1 and 2, we reproduced the recommended iterative methods using\nPython+GPU and added a comparative analysis of integral-based iterative reconstruction\nwith SlingBAG algorithm in the main text.\n3. Analyzing the compatibility with Regularization: In response to Reviewer 2’s concern\nregarding SlingBAG’s compatibility with regularization, we included supplementary\nexperiments demonstrating SlingBAG’s effective integration with regularization\ntechniques.\nNotably, after optimizing the CUDA kernels, in vivo 3D reconstruction for mouse brain, rat\nliver, and kidney datasets was completed within 1–5 hours. Specifically, the mouse brain\nreconstruction finished in 1.89 hours, with the maximum amplitude projection (MAP) contrast-\nto-noise ratio (CNR) increasing from 28.49 (using universal back-projection (UBP)) to 46.39\n(using SlingBAG algorithm). We believe the SlingBAG algorithm will emerge as one of the\nmost competitive algorithms for large-scale 3D photoacoustic reconstruction. Below we\nprovide detailed responses to each reviewer's comments. Note that the comments from editor and\nreviewers are in italic font, whereas our reply is not italicized.\nThe Point-by-Point Responses to the Comments\nReviewer #1\nIn the work of Li et al, the authors proposed an image reconstruction method using the point\ncloud-based representation for 3D photoacoustic imaging. This idea was translated from the\nGaussian splatting technique in the field of computer vision which has been a research hotspot\nin this field in recent years. It is very interesting to see its successful translation to\nphotoacoustic imaging which is novel. The authors demonstrated that the SlingBAG algorithm\nenables impressively low memory usage and fast iteration for high-quality large-scale 3D PA\nimaging with sparse sensor distribution, thanks to the sparsity characteristics in PA sources\nand the developed differentiable rapid radiator model for forward simulation. Overall,\nthe paper is interesting, clearly presented and is novel in methodology and has shown some\npotentials in addressing important challenges in large-scale 3D PA reconstruction. However,\nI think this works need further improvements to warrant publishing in Nature Communications.\nParticularly, the authors should compare their proposed method with current state-of-the-art\nalgorithms, the model-based iterative reconstruction methods, instead of only with the UBP\nalgorithms. After all, model-based iterative reconstruction methods have been studied and\ndeveloped for 3D PA imaging in the last ten years, and some of them were proposed to address\nthe same problem just like this work, while the addressed problem size may be a bit smaller. I\nrecommend the following works [1-3] for the authors' reference.\nReferences\n[1] L. Ding, X. L. Deán-Ben and D. Razansky, \"Efficient 3-D Model-Based Reconstruction\nScheme for Arbitrary Optoacoustic Acquisition Geometries,\" in IEEE Transactions on Medical\nImaging, vol. 36, no. 9, pp. 1858-1867, Sept. 2017"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have done a great job in addressing most of my previous comments. I particularly appreciate their efforts in\nimplementing the MB-PD method for comparing the methods using simulation data. However, significant technical\nadvancement over the state-of-the-art must be thoroughly demonstrated to warrant publication in Nature Communications. I\nbelieve the comparison should be extended to experimental data to fully demonstrate the advantages of the proposed\nmethod over MB-PD. With the already implemented code and the presented data, it does not require too much extra efforts to\ncomplete the investigation and make the work more rigorous.\nIt is very interesting to see that the simulated PA signals well match those predicted by the differentiable rapid radiator in a\ncontrolled environment (Supplementary Fig. 1&9). Readers may also be interested to know such discrepancy in the in vivo\ndata given that the electrical impulse response of the transducer element was not considered in the forward model for PA\nsignal prediction, i.e. modeling errors (perhaps small) existed in the reconstruction. I presume the discrepancy should be\nsmall, as SlingBAG has achieved high-quality reconstruction.", "reviewer_2": "The revised manuscript presents the Sliding Gaussian Ball Adaptive Growth (SlingBAG) algorithm, an optimization-based\nreconstruction algorithm for large-scale 3D photoacoustic (PA) imaging. The authors have been highly responsive to the\nprovided critiques.\nThey have implemented alternative model-based image reconstruction methods to serve as more meaningful references\nthan the naïve UBP method.\nThey have also addressed the issue of regularization. Specifically, they have shown that their framework can incorporate\ncertain standard regularization functionals. They have also emphasized that their method of discretizing the object also\nprovides a means of implicit regularization that can improve image quality.\nThey have provided the requested additional data regarding memory requirements and reconstruction times to better support\ntheir claims.\nThey have also provided additional mathematical equations that describe the point cloud to voxel conversion.\nIn summary, the manuscript is significantly improved and I find it suitable for publication."}, "author_response": "Authors’ Response\nfor the submitted paper,\n“Sliding Gaussian ball adaptive growth: point cloud-based iterative\nalgorithm for large-scale 3D photoacoustic imaging”\nWe would like to thank both the editor and the reviewers for your valuable comments and\nsuggestions, which have helped us to significantly improve our manuscript. We have carefully\nrevised our manuscript accordingly, and thoroughly addressed all comments point by point in\ndetail, as well as providing additional experiments. The key revisions include:\n1. We extended the comparison between the MB-PD method with the proposed point-cloud\nbased SlingBAG approach to all in vivo experimental datasets, results further demonstrate\nthe superiority of the SlingBAG.\n2. We added discussion of the influence of the electrical impulse response (EIR) on the\nSlingBAG’s performance. Both simulated hand-vessel data (576 elements) and in vivo rat\nliver and kidney datasets were used to compare SlingBAG reconstructions performed with\nand without signal deconvolution with respect to the EIR. Results show the influence is\nminor.\nWe believe the SlingBAG algorithm will emerge as one of the most competitive algorithms for\nlarge-scale 3D photoacoustic reconstruction. Below we provide detailed responses to each\nreviewer's comments. Note that the comments from editor and reviewers are in italic font, whereas\nour reply is not italicized.\nThe Point-by-Point Responses to the Comments\nReviewer #1\nIn the work of Li et al, the authors proposed an image reconstruction method using the point\ncloud-based representation for 3D photoacoustic imaging. This idea was translated from the\nGaussian splatting technique in the field of computer vision which has been a research hotspot\nin this field in recent years. It is very interesting to see its successful translation to\nphotoacoustic imaging which is novel. The authors demonstrated that the SlingBAG algorithm\nenables impressively low memory usage and fast iteration for high-quality large-scale 3D PA\nimaging with sparse sensor distribution, thanks to the sparsity characteristics in PA sources\nand the developed differentiable rapid radiator model for forward simulation. Overall,\nthe paper is interesting, clearly presented and is novel in methodology and has shown some\npotentials in addressing important challenges in large-scale 3D PA reconstruction. However,\nI think this works need further improvements to warrant publishing in Nature Communications.\nParticularly, the authors should compare their proposed method with current state-of-the-art\nalgorithms, the model-based iterative reconstruction methods, instead of only with the UBP\nalgorithms. After all, model-based iterative reconstruction methods have been studied and\ndeveloped for 3D PA imaging in the last ten years, and some of them were proposed to address\nthe same problem just like this work, while the addressed problem size may be a bit smaller. I\nrecommend the following works [1-3] for the authors' reference.\nReferences\n[1] L. Ding, X. L. Deán-Ben and D. Razansky, \"Efficient 3-D Model-Based Reconstruction\nScheme for Arbitrary Optoacoustic Acquisition Geometries,\" in IEEE Transactions on Medical\nImaging, vol. 36, no. 9, pp. 1858-1867, Sept. 2017"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have well addressed all of my concerns. In this revision, the authors have extended the comparision between MB-PD and SlingBAG to the experimental data. They have also nicely demonstrated the role of EIR of the transducer elements in the reconstructed images using simulated and in vivo data. I believe that the work is ready for publishing in Nature Communications. Congratulations to the authors on this nice work!"}, "author_response": "Authors’ Response\nfor the submitted paper,\n“Sliding Gaussian ball adaptive growth: point cloud-based iterative\nalgorithm for large-scale 3D photoacoustic imaging”\nWe are writing to express our sincere gratitude for the positive feedback from both the editor and the reviewer\ns on our revised m anuscript titled \"Sliding Gaussian ball adaptive growth: point cloud-based\niterative algorithm for large-scale 3D photoacoustic imaging\".\nWe are pleased that our revisions have fully addressed all concerns and that the reviewer finds\nour work ready for publication.\nPlease find our point-by-point response to the reviewer's kind remarks below.\nThe Point-by-Point Responses to the Comments\nReviewer #1\n(Remarks to the Author):\nThe authors have well addressed all of my concerns. In this revision, the authors have extended\nthe comparison between MB-PD and SlingBAG to the experimental data. They have also nicely\ndemonstrated the role of EIR of the transducer elements in the reconstructed images using\nsimulated and in vivo data. I believe that the work is ready for publishing in Nature\nCommunications. Congratulations to the authors on this nice work!\n(Remarks on code availability):\nIt is working.\nResponse:\nWe sincerely thank Reviewer #1 for their positive and encouraging assessment. We are\ndelighted that they found all previous concerns to be well-addressed. We are particularly\ngrateful for their acknowledgment of the key additions in this revision, including the extended\ncomparison of MB-PD and SlingBAG on experimental data and the demonstration of the\ntransducer EIR's role using simulated and in vivo data.\nWe also appreciate the reviewer taking the time to verify the functionality of our code. We are\ndeeply grateful for their strong endorsement that our work is ready for publication in Nature\nCommunications and for their kind congratulations."}]

Decision letter