Matrix-free photoacoustic image reconstruction via sensor-token networks

1 Introduction

Photoacoustic tomography (PAT) has emerged over the past two decades as one of the most promising hybrid biomedical imaging modalities, combining the rich endogenous optical absorption contrast of hemoglobin, melanin, and lipids with the deep tissue penetration and high spatial resolution of ultrasound [1, 2]. A nanosecond pulsed laser illumination of biological tissue causes localized heating and thermoelastic expansion, generating broadband acoustic pressure waves that are detected by an array of ultrasound transducers placed around the tissue. The recovered initial pressure distribution provides a quantitative spatial map of optical absorbers and underpins emerging clinical applications, including breast cancer screening, vascular imaging, melanoma detection, and functional neuroimaging [3, 4].

The reconstruction of the initial pressure distribution from the time-resolved sensor signals is governed by the linear acoustic wave equation and, under the assumption of a homogeneous lossless medium, can be expressed compactly as a linear inverse problem y = Hx₀+n, where H is the discretized forward (system) operator, x₀ is the vectorized image (gray-scale) pixels and n is the noise. In practice, H is severely under-determined or ill-posed because (i) the sensor array geometry is sparse and frequently subtends less than 2π around the tissue, (ii) the time-sampling and grid-discretization introduce numerical artifacts, and (iii) measurement noise compounds at high spatial frequencies [5, 6]. Consequently, naive inversion (e.g., back-projection or filtered back-projection) yields reconstructions corrupted by streaking artifacts, smearing, and poor resolution of fine vascular detail [7].

To overcome ill-posedness, the PAT community has invested heavily in compressive sensing formulations that exploit the spatial sparsity of vessel-like structures. Iterative algorithms such as the Iterative Shrinkage-Thresholding Algorithm (ISTA) and the split-Bregman total-variation (SBTV) method have been shown to outperform classical back-projection in limited-view settings, recovering coarse structure with high spatial correlation [8, 9]. A systematic comparison of these reconstruction algorithms under full-view and limited-view conditions was reported in John and Barhumi [7], where iterative l₁-regularized solvers were shown to be the most robust baselines. Our own comparative analysis [7] found SBTV to be the most accurate of these classical solvers, but at the cost of a preprocessing step and a per-iteration complexity that scales with the dimensions of H, limiting its use as grid and sensor counts increase. Plug-and-play frameworks that integrate learned denoising priors with classical iterative updates further improved fidelity for limited-sample regimes [10]. More recently, deep learning has been applied to PAT reconstruction in two principal modes [11–14]. The first uses convolutional encoder-decoder networks as post-processing stages applied to an initial classical estimate, reducing artifacts but retaining the system-matrix dependence in the front-end. The second is algorithm unrolling, in which a fixed number of iterations of an iterative algorithm is unfolded into a trainable deep network with each layer corresponding to one iteration [15]. Unrolled networks preserve the physics of the forward operator while exploiting data-driven priors, and have demonstrated strong reconstruction performance under sparse-view PAT conditions [5, 11, 14]. Despite these advances, both paradigms retain repeated Hx and H^*T* multiplications at inference time, which for clinically relevant grid resolutions can dominate the per-reconstruction cost [16]. This article proposes a complementary line of attack: a direct, data-driven mapping from raw sensor measurements to the reconstructed image that completely bypasses the system matrix at inference. The proposed architecture, named the Sensor-Token Reconstruction Network (STRN), treats each sensor's full time-series as a token in a Transformer encoder [17]. Multi-head self-attention allows every sensor to interact with every other sensor, and a learned spatial decoder maps the resulting sensor representation back to the image domain. In inference, reconstruction reduces to a single forward pass through the network. The network is presented here as a direct reconstruction architecture rather than as evidence for self-attention as a reconstruction mechanism: as the ablation in Section 5.7 shows, the attention sub-layers are not the decisive component, and the learned decoder drives fidelity. Crucially, the analytical system-matrix is not discarded but repurposed: it encodes the acoustic physics that make the forward problem well-defined, and is used to synthesize physically faithful training data from which the network learns the inverse mapping. The system matrix, therefore, underpins training, while inference is freed of its per-reconstruction cost.

The principal contributions of this work are summarized as follows:

• An analytical k-space forward model is constructed and validated against the k-Wave pseudo-spectral time-domain solver under matched geometry, achieving a mean per-sensor Pearson correlation of r̄=0.919±0.049. Two targeted regularizations: k-space apodization and Gaussian temporal damping are introduced and shown by ablation to act synergistically, reducing the energy-normalized mismatch by 49%.
• The Sensor-Token Reconstruction Network (STRN) is proposed as a direct, matrix-free reconstruction architecture for PAT: a sensor-as-token Transformer encoder followed by a learned spatial decoder that maps raw measurements to the image in a single forward pass, removing the per-iteration H/H^T products required by iterative and unrolled methods and yielding an order-of-magnitude inference speed-up on the same hardware.
• A comprehensive empirical comparison is reported against three representative baselines, ISTA, SBTV, and Learned Iterative Shrinkage-Thresholding Algorithm (LISTA), on the same 46-sample held-out test set, using Structural Similarity Index Measure (SSIM), Peak Signal-to-Noise Ratio (PSNR), Normalized Mean Squared Error (NMSE), and Pearson correlation. The STRN achieves the highest mean SSIM and PSNR and the lowest mean NMSE of the four methods, with paired statistical tests confirming significance at p < 10⁻⁸ for the principal claims.
• Two analyzes establish the validity and the mechanism of the approach. An independent generalization test on k-Wave-simulated measurements quantifies a substantial but structure-preserving and rapidly closable domain gap, addressing the inverse-crime concern that arises when training and testing share a forward operator. A parameter-matched ablation then isolates the contribution of each component and shows transparently that the learned spatial decoder, rather than the self-attention sub-layers, is the dominant driver of reconstruction fidelity, a negative result reported openly to guide simpler direct-reconstruction designs. A robustness study further characterizes the trained model under additive noise and per-channel gain and timing perturbations, identifying timing jitter as its dominant sensitivity.

The remainder of the article is organized as follows: Section 2 reviews related work in PAT reconstruction. Section 3 describes the analytical H-matrix construction, its regularization, and the validation methodology. Section 4 introduces the proposed Sensor-Token Reconstruction Network. Section 5 reports the experimental setup and full numerical comparison. Section 6 discusses interpretation, limitations, and clinical implications. Section 7 concludes the article and outlines directions for future work.

2 Related works

Reconstruction algorithms for PAT can be broadly grouped into four families: analytical back-projection methods, classical iterative compressive-sensing methods, model-based deep-learning methods, and direct (end-to-end) deep-learning methods. The proposed Sensor-Token Reconstruction Network belongs to the fourth family. We briefly review each below.

2.1 Analytical and iterative compressive-sensing reconstruction

Analytical reconstruction by universal back-projection [18] provides closed-form image estimates under idealized geometry assumptions, but degrades rapidly under sparse or limited-view sensing. Compressive-sensing-based iterative methods address this limitation by exploiting the sparsity of vessel images in suitable transform domains. ISTA [19] solves an l₁-regularized least-squares problem and has been widely applied to PAT [7]. The SBTV formulation [20] enforces a piecewise-smooth structure on the reconstructed image and has been particularly effective in suppressing high-frequency artifacts arising from sparse sensor arrays. Total-variation-based reconstruction for PAT specifically was studied in [8], where the choice of regularization parameter and stopping criterion were shown to strongly influence reconstruction fidelity. A systematic comparison of full-view and limited-view iterative methods was reported in John and Barhumi [9], demonstrating that no single algorithm dominates across all sensing configurations and that the appropriate method depends on the sparsity geometry. In our earlier comparative work, ISTA-based reconstruction was examined in John and Barhumi [11] and split-Bregman total variation in John and Barhumi [7], establishing the classical iterative baselines used in the present study.

2.2 Limited-view and limited-sample PAT

The limited-view regime, in which the sensor array subtends less than the full 2π around the tissue, introduces additional ill-posedness because spatial frequencies aligned with the missing aperture directions cannot be recovered. Compressive sensing techniques tailored to this scenario have been studied in John and Barhumi [9], where iterative l₁ solvers were combined with carefully chosen sensor sub-sampling patterns to maintain reconstruction quality at reduced sensor counts. The complementary problem of limited-sample PAT, recovering the initial pressure from a small training corpus, which is the regime relevant to many clinical settings, was addressed in John and Barhumi [10], Venkatakrishnan et al. [21], Kamilov et al. [22], and Awasthi et al. [23] using a plug-and-play framework that combines a classical iterative solver with a learned denoiser. The handling of missing data in limited-view PAT, where some sensor signals are corrupted or absent, was addressed by a compressive-sensing-driven deep-learning approach in John and Barhumi [5].

2.3 Model-based deep learning and algorithm unrolling

Algorithm unrolling [15] unfolds a fixed number of iterations of an iterative reconstruction algorithm into a trainable deep network, with each layer corresponding to one iteration. The proximal operator within each layer (typically a soft-threshold for L₁-regularized problems) is replaced by a learnable mapping, while the data-fidelity step retains the forward operator H and its adjoint H^T explicitly. The LISTA [19, 24] is the canonical example. Unrolled networks for PAT have been studied in John and Barhumi [11], where the architecture was made intelligible by exposing per-layer thresholds for clinical interpretation, and in John and Barhumi [5], where unrolled networks were applied to breast-cancer-related limited-view PAT data. While these methods offer substantial improvements in reconstruction quality over classical solvers and retain physical interpretability, they share an important inference-time limitation: every layer requires an Hx and H^T multiplication (each an O(MN) matrix-vector product), so the per-image computational cost still scales as KO(MN), where K is the unrolled depth. For clinically realistic image sizes and sensor counts, this remains a bottleneck for real-time use.

2.4 Direct (end-to-end) deep learning for image reconstruction

A more recent line of work seeks to avoid the H-matrix at inference entirely by training a deep network to map sensor measurements directly to the reconstructed image. Convolutional neural networks and U-Net variants have been used in this role [2, 11, 25–27], but their treatment of the sensor array as a regular grid does not exploit the nonlocal acoustic relationships between sensors. In a related generative direction, diffusion-based models have recently been investigated for PAT reconstruction, with comparative studies benchmarking several such frameworks against classical solvers [28]. Diffusion-based priors have proven particularly effective in the sparse- and extremely-limited-view regimes: combining a diffusion model with model-based iteration [29], using multiple diffusion models with multi-scale priors for extremely limited-view reconstruction [30], and diffusion-based sparse angular recovery for isotropic-resolution imaging [31]. More recently, attention-based mechanisms have been explored in PAT reconstruction to improve feature extraction and information fusion from sparse measurements [32, 33]. However, these approaches typically apply attention within learned feature representations rather than explicitly modeling interactions among individual sensors. The proposed Sensor-Token Reconstruction Network belongs to this family and uses a Transformer-based architecture [17] in which each sensor is represented as a token. This formulation preserves the physical role of the sensor array, sensors are unordered entities communicating via the acoustic field, and enables fully parallel inference without iterative dependence on H. The extent to which the self-attention sub-layers themselves contribute to accuracy is examined directly in the ablation of Section 5.7.

2.5 Positioning of the proposed work

The present work differs from prior compressive-sensing PAT reconstruction in that the H-matrix is not invoked at inference; it differs from prior unrolled-network PAT reconstruction in that no per-layer forward/adjoint operation is required, and it differs from prior CNN-based direct reconstruction in that the sensor array is treated as a set of tokens with explicit self-attention rather than as a structured 2-D input. The validated analytical H-matrix from Section 3 is used only at the training-data generation stage to synthesize sensor measurements from phantom images, after which the network learns the inverse mapping without further reference to H.

3 Methodology

3.1 Problem formulation

In PAT, a short-pulsed laser illuminates a biological tissue sample, inducing a localized temperature rise that generates a broadband acoustic pressure wave through the thermoacoustic effect. The spatial distribution of this initial acoustic pressure, denoted x0∈ℝNxNy (where N_*x* and N_*y* are the vectorized image grid dimensions in pixels), encodes the optical absorption map of the tissue and constitutes the quantity of interest. A surrounding array of N_*s* ultrasound transducers records the time-varying pressure field over N_*t* discrete time steps, producing the measurement vector y∈ℝNsNt. Under the assumption of a lossless, homogeneous medium with constant speed of sound c, the forward problem is linear and may be written compactly as Equation 1:

y=Hx0+n(1)

where H∈ℝNsNt×NxNy is the system matrix (hereafter the H-matrix) mapping the vectorised image x₀ to all sensor measurements, and n represents additive measurement noise. The reconstruction problem estimating x₀ from y is an ill-posed linear inverse problem, the solution quality of which depends critically on the accuracy of H. The construction and validation of this matrix is the primary concern of this section.

3.2 Analytical k-space forward model construction

The H-matrix is constructed analytically using a pseudospectral k-space method, which exploits the spatial Fourier decomposition of the initial pressure field. In the k-space domain, the solution to the lossless acoustic wave equation for a two-dimensional homogeneous medium can be expressed at each sensor location r_*s* and time t_*n* as Equation 2:

p(rs,tn)=Winvdiag[cos(c|k|tn)]Wfwdx0(2)

where Wfwd∈ℂ(Nout2×NxNy) is the two-dimensional forward discrete Fourier transform (DFT) matrix that projects the image onto the k-space grid of size N_out×N_out, with m and l indexing the k-space; Winv∈ℂ(Ns×Nout2) is the inverse DFT matrix evaluated at the sensor coordinates {r_*s*}; and |k| denotes the magnitude of the two-dimensional spatial frequency vector

where m∈[-Nx2,Nx2-1] and l∈[-Ny2,Ny2-1]. The cosine term encodes the time evolution of each spectral component under free-space wave propagation. Stacking this expression across all N_*s* sensors and N_*t* time steps, the H-matrix takes the factored form given in Equation 3:

H=KWfwd(3)

where K∈ℝNsNt×Nout2 is the time-domain propagation kernel, whose row corresponding to sensor s at time t_*n* is the s-th row of W_inv scaled by the cosine propagation factor cos(c|k|t_*n*). This formulation yields the full system matrix in a single construction pass without iterative simulation, making it computationally efficient for repeated use in the inverse problem.

3.3 Simulation grid, time array, and sensor configuration

The imaging domain is discretized on a uniform Cartesian grid of N_*x*×N_*y* pixels with isotropic grid spacing d_*x* = d_*y*. The choice of d_*x* governs both the spatial resolution of the reconstructed image and the temporal sampling requirements imposed by the Courant-Friedrichs-Lewy (CFL) stability criterion. For a two-dimensional acoustic problem, the maximum stable time step is

dt,max=dxc2(4)

To ensure that all wavefronts originating within the imaging domain reach the sensor boundary before the simulation terminates, the total measurement duration T and the number of time samples N_*t* are determined jointly using Equation 5:

Nt=Tξdt,max,(5)

where ξ = 1.25 is a safety factor that provides a temporal margin to prevent signal truncation artifacts. The sensor array consists of N_*s* = 71 transducers distributed uniformly along the four edges of a 70 × 70 boundary grid that encloses the 64 × 64 imaging domain, with sensors placed at every fourth grid point; the imaging domain, the sensor boundary, and the sensors are all centered within the larger 140 × 140 k-space simulation grid (N_*out*). Because sensors are placed at grid indices 1, 5, …, 69 along each 70-point edge, the index 70 is never sampled; consequently, only the shared corner (1,1) is common to two edges, while the remaining three corners fall outside the sampling pattern. This yields 4 × 18 − 1 = 71 unique transducers. Sensor Cartesian coordinates were extracted using the k-Wave grid2cart() function and used consistently in both the analytical model and the numerical reference simulation. The complete set of simulation parameters is listed in Table 1. The system geometry, including the imaging grid and sensor array configuration, is illustrated in Figure 1 and Table 1. It is worth clarifying the sensing regime studied here, as the terms sparse-view, limited-view, and limited-sample are used with varying meanings in the PAT literature. Limited-view refers to an aperture that subtends less than the full 2π around the object, so that entire angular ranges of the acoustic field are never measured, and the associated spatial frequencies are unrecoverable. Sparse-view refers to a full surrounding aperture that is undersampled, i.e., the sensors enclose the object but are too few to satisfy spatial Nyquist sampling. Limited-sample refers to the size of the training corpus, a property of the data rather than of the measurement geometry. The configuration adopted in this work is sparse-view: the N_*s* = 71 transducers form a closed square aperture that fully surrounds the imaging domain, an approximately 2π view, but sample it sparsely, at every fourth boundary grid point. It is additionally limited-sample, since the network is trained on a small phantom corpus (Section 4.7). The present study does not address the limited-view regime, in which part of the aperture is physically absent; that setting is treated only as related work.

Table 1. Simulation parameters for the analytical forward model (d_*x* = 0.100 mm).

3.4 Regularization of the analytical model

In its raw form, the analytical model exhibits two systematic artifacts. First, the truncation of the DFT matrices at the boundaries of the k-space grid introduces Gibbs-like aliasing, manifesting as large-amplitude spikes in the computed sensor signals. Second, without any temporal envelope, the model produces non-physical oscillations at late time samples that have no counterpart in a physically bounded measurement. Two targeted regularizations are applied to suppress these artifacts as stated below:

• K-Space Apodization (Windowing): A two-dimensional window function w(k) is applied element-wise to both the forward and inverse DFT matrices prior to construction of H (Equations 6a, 6b):
• where ⊙ denotes the Hadamard (element-wise) product. The window tapers the DFT coefficients smoothly to zero at the k-space boundary, suppressing Gibbs-like aliasing without introducing significant spectral leakage. Three standard apodization functions are candidates for this role: the Hann window, which provides a smooth cosine-squared taper and is the conventional choice for k-space apodization in spectral methods; the Hamming window, which raises the minimum sidelobe level relative to Hann by adding a constant offset; and the Blackman window, which achieves lower sidelobe levels at the cost of a wider main lobe.
• Gaussian Temporal Damping: A Gaussian envelope g(t) is applied to the cosine propagation kernel to enforce realistic signal decay within the measurement window (Equation 7):
• where τ is the damping time constant (measured in seconds), set to 90% of the total measurement duration T. This ensures the analytical signal decays smoothly before the end of the measurement window, preventing the non-physical late-time oscillations that arise in the undamped model. We emphasize that this damping is a numerical regulariser, not a model of physical acoustic attenuation: the medium is assumed homogeneous and lossless throughout, and the damping instead compensates for the truncation and boundary mismatch between the free-space analytical formulation and the PML-terminated k-Wave reference, whose discrepancy accumulates in the late-time tail. The value τ = 0.9T was chosen so that attenuation becomes appreciable only within the final ~10% of the window, the late-time region where the two models diverge, while leaving the physically informative early arrivals essentially unchanged; it minimized the energy-normalized mismatch against k-Wave among the values examined (Section 5.2.2).

3.5 Numerical validation via k-Wave

The analytical H-matrix is validated by comparing the predicted sensor signals *Hx_0* against those produced by the k-Wave pseudo-spectral time-domain solver [34], which numerically integrates the full linearised acoustic wave equation. k-Wave simulations were executed with a perfectly matched layer (PML) of 40 grid points on each side of the domain and an absorption coefficient α = 2, using an identical spatial grid, time array, sensor configuration, and input pressure image x₀ as the analytical model. The degree of agreement between the two models for sensors is quantified by the Pearson correlation coefficient:

rs=Cov(ysan,yskw)σ(ysan)·σ(yskw)(8)

where ysan and yskw denote the analytical and k-Wave time-series at sensor s, respectively. A value of r_*s* close to unity indicates that the analytical model faithfully reproduces the temporal waveform structure predicted by the numerical simulation. The mean correlation across all N_*s* sensors, r̄=1Ns∑rs, serves as the primary figure of merit for model accuracy.

3.6 Grid spacing study

The imaging grid spacing d_*x* is a critical design parameter because it sets a direct trade-off between spatial resolution and forward-model fidelity. A finer grid resolves finer image detail but through the CFL condition of Equation 4, requires a proportionally larger number of time samples N_*t*, which extends the temporal window over which the analytical and numerical models must remain in agreement. Since the analytical model assumes free-space propagation whereas k-Wave enforces absorbing boundary conditions through the PML, the mismatch between the two accumulates in the late-time signal tail, and the extent of that tail grows with N_*t*. A larger d_*x* therefore reduces model mismatch at the cost of resolution, while a smaller d_*x* improves resolution but degrades agreement with the reference.

To select an appropriate value, three candidate grid spacings, d_*x* = 0.075, 0.100, and 0.125 mm, were evaluated against the k-Wave reference under otherwise identical conditions. For each spacing the full H-matrix and the corresponding k-Wave signals were generated, and the per-sensor Pearson correlation was computed as the figure of merit, summarized by its mean, its minimum across the sensor array, and the fraction of sensors exceeding the acceptance threshold of r = 0.80. The spacing yielding the best agreement on these measures was carried forward to all subsequent experiments.

4 Proposed method

The classical iterative and learning-based reconstruction approaches retain a dependence on the H-matrix at inference time, either through repeated matrix-vector products or embedded forward operators. A complementary data-driven approach is proposed in this article, in which a neural network learns the mapping from sensor measurements y directly to the reconstructed image x^*. The proposed architecture, referred to as the Sensor-Token Reconstruction Network (STRN), treats each sensor's time series as a token in a Transformer encoder, allowing every sensor to attend to all others and capture the non-local acoustic relationships encoded in the PAT forward model. The steps that are adopted for the mapping of sensor data y, to the reconstructed image x^* are as listed below:

4.1 Input projection

The sensor measurement y is reshaped into a matrix y∈ℝNt×Ns, so that each column ys∈ℝNt is the time series recorded by sensor s, for s = 1, …, N_*s*. Each time series is mapped to a d_model-dimensional token embedding by a shared linear projection (Equation 9):

zs=Wprojys+bproj(9)

where W_proj is a learnable weight matrix of size d_model×N_*t* and b_proj is a learnable bias vector of size d_model. z_*s* is the resulting token embedding for sensor s. All token embeddings are assembled column-wise into the matrix Z = [z₁, …, z_{N_*s*}].

4.2 Positional encoding

Because the Transformer encoder treats its input as an unordered set, sinusoidal positional encodings [17] are added to each token to convey sensor ordering. For sensor index s and embedding dimension index k (k = 0, 1, 2, …, d_model−1), the encoding is defined in Equations 10a, 10b:

PE(s,k)=sin(s10,000k/dmodel), whenkis even(10a)

PE(s,k)=cos(s10,000(k-1)/dmodel), whenkis odd(10b)

The denominator 10,000k/dmodel grows from 1 (at k = 0) to 10,000 (at k = d_model), so that different dimensions oscillate at different frequencies: low-index dimensions change rapidly with s and encode fine position differences, while high-index dimensions change slowly and encode coarse spatial layout. Each token is updated as z_*s*←z_*s*+PE(s, :), where PE(s, :) denotes the full d_model-dimensional encoding vector for sensor s.

4.3 Multi-head self-attention encoder

The position-encoded token matrix Z is then processed by L stacked Transformer encoder blocks. Each block applies a multi-head self-attention (MHSA) sub-layer followed by a position-wise feed-forward sub-layer, with a residual connection and layer normalization after each sub-layer.

In the MHSA sub-layer, A attention heads are computed in parallel. For each head h∈{1, …, A}, the token matrix Z is projected into a query matrix Q_*h*, a key matrix K_*h*, and a value matrix V_*h*, each of per-head dimension d_head = d_model/A, via

Qh=W(Qh)Z, Kh=W(Kh)Z, Vh=W(Vh)Z,(11)

where W(Qh), W(Kh), W(Vh)∈ℝdhead×dmodel are learnable projection matrices. The query Q_*h* encodes what each sensor is searching for, the key K_*h* encodes what each sensor offers, and the value V_*h* holds the content to be aggregated. The attention output for head h is

Ah=softmax(QhTKhdhead), headh=VhAhT(12)

With the token matrix Z∈ℝdmodel×Ns (one column per sensor), the projections satisfy Qh,Kh,Vh∈ℝdhead×Ns, so the score matrix QhTKh∈ℝNs×Ns has entry (i, j) measuring the relevance of sensor j to sensor i. The row-wise softmax yields the attention matrix Ah=[αij]∈ℝNs×Ns, whose rows sum to one, and the per-head output headh=VhAhT∈ℝdhead×Ns re-expresses each sensor's value as an attention-weighted combination of all sensors. Division by dhead prevents the scores from becoming excessively large before the softmax [17]. The outputs of all A heads are concatenated and projected back to d_model dimensions via a learnable output matrix W₀ (of size d_model×d_model). A residual connection is added and the result is normalized,

Z←LayerNorm(Z+W0[head1;head2;…;headA])(13)

where LayerNorm standardizes each token vector to zero mean and unit variance across its d_model dimensions, followed by a learnable scale and shift. The semicolon denotes row-wise concatenation. The residual addition ensures that attention learns only the correction to the existing token representation, preserving information and stabilizing gradient flow.

4.4 Position-wise feed-forward sub-layer

Following the MHSA sub-layer, each token is passed independently through a two-layer feed-forward network (FFN),

FFN(zs)=W2·ReLU(W1zs+b1)+b2(14)

where W₁ (of size d_ff×d_model) expands each token from d_model to a hidden dimension d_ff, ReLU(α) = max(0, α) is a non-linear activation applied element-wise, and W₂ (of size d_model×d_ff) projects back to d_model dimensions. The same residual-plus-LayerNorm operation as in Equation 13 follows. Equations 11–14 together define one complete encoder block; this block is applied L times in sequence.

4.5 Spatial decoder

After L encoder blocks, the token matrix Z (of size d_model×N_*s*) is flattened column-wise into a single vector z_flat of length d_modelN_*s*, which is decoded into the reconstructed image through three fully connected layers (Equations 15a–15c):

h1=ReLU(Wd1zflat+bd1),(15a)

h2=ReLU(Wd2h1+bd2),(15b)

x*=sigmoid(Wd3h2+bd3),(15c)

where W_d₁, W_d₂, and W_d₃ are learnable weight matrices; h₁ and h₂ are intermediate hidden vectors and sigmoid(α)=11+e-α maps each output to [0, 1] to match the normalized image range. The final output x^* is reshaped to the desired image size.

4.6 Training objective

All learnable parameters are jointly optimized by minimizing a vessel-weighted mean squared error loss. Since vessel pixels occupy a small fraction of the image area relative to the background, a higher penalty weight is assigned to them. Denoting the ground truth phantom as x₀ and the network prediction as x^*, the loss is defined as

L=∑iwi(xi*-x0,i)2,(16)

where w_*i* = 8 for vessel pixels (defined as those with x_{0, *i*}>0.1) and w_*i* = 1 for background pixels (vessel mask is computed per sample from the ground-truth phantom). The process included in the proposed method is summarized in Figure 2.

4.7 Dataset preparation

A total of 191 grayscale radiological images were downloaded from Radiopaedia.org [35]. The images were selected to exhibit vessel-like, branching structures representative of the vascular contrast targeted in photoacoustic imaging, and each was treated as a surrogate initial-pressure map x₀ rather than as a calibrated optical-absorption distribution. Each source image was converted to grayscale, cast to double precision, and resized to the 64 × 64 imaging grid by bicubic interpolation. Intensity was then contrast-enhanced by saturating the lowest and highest 1% of pixel values and linearly mapping the remainder to the full [0, 1] range (MATLAB imadjust); the resulting map was scaled by the unit peak initial-pressure magnitude (p₀ = 1, arbitrary units) and embedded at the center of the 140 × 140 k-space simulation grid, from which the forward sensor signal y = Hx₀ was computed. To mitigate overfitting, given the limited sample count, data augmentation was applied by translating each phantom by small integer offsets of 2 to 3 pixels in the horizontal, vertical, and diagonal directions, and recomputing the sensor signal for each translated phantom. Of the 191 source phantoms, 44 were held out for validation and 46 for testing (neither augmented); the remaining training phantoms were augmented by translation to 488 samples, giving 578 images in total.

5 Results

The experiments reported in this section are organized into eight parts. Section 5.1 validates the analytical H-matrix by comparing its predictions against the k-Wave numerical reference. Section 5.2 presents the empirical selection of the k-space apodization window, the Gaussian damping time constant, and the imaging grid spacing. Section 5.3 establishes baseline reconstruction performance using the classical iterative methods ISTA and split-Bregman total variation (SBTV). Section 5.4 reports the training behavior and validation performance of the proposed Sensor-Token Reconstruction Network. Section 5.5 provides a direct quantitative and qualitative comparison of all four approaches. Section 5.6 evaluates the trained network on independently simulated k-Wave measurements without retraining; Section 5.7 presents a parameter-matched architectural ablation isolating the source of reconstruction fidelity; and Section 5.8 characterizes robustness to additive noise and per-channel gain and timing perturbations. Throughout, reconstruction quality is evaluated using SSIM, PSNR, NMSE and Pearson correlation r̄img, computed on normalized images in the range [0, 1].

5.1 Forward model validation

The analytical H-matrix was validated by comparing the sensor signals it predicts against those produced by the k-Wave pseudospectral solver under identical imaging conditions. For each sensor s = 1, …, N_*s*, the Pearson correlation coefficient r_*s* between the normalized analytical and k-Wave time series was computed as defined in Equation 8. Figure 3 shows representative waveform overlays for five randomly selected sensors, illustrating the close agreement between the two models across the full measurement window of N_*t* = 160 time samples.

The mean correlation across all sensors was r̄=0.919±0.049, with all 71 sensors exceeding r_*s* = 0.80.¹ The full per-sensor correlation profile is shown in Figure 4, together with the correlation histogram and summary statistics. The consistently high correlation across the entire sensor array confirms that the analytical model faithfully reproduces the temporal waveform structure of the k-Wave reference, and that the H-matrix is a reliable and computationally efficient substitute for repeated numerical simulation in the training data generation stage. The small residual discrepancy 1-r̄≈0.081 is attributable to the boundary mismatch between the two models: the analytical formulation assumes free-space propagation, whereas k-Wave enforces absorbing boundary conditions via the perfectly matched layer. This discrepancy accumulates preferentially in the late-time signal tail and is effectively suppressed by the Gaussian temporal damping applied during H-matrix construction, as discussed in Section 5.2.

5.2 Regularization and grid spacing selection

This section reports the empirical selection of the three key design parameters of the analytical forward model: the k-space apodization window, the Gaussian temporal damping constant, and the imaging grid spacing. Each parameter was evaluated independently against the k-Wave reference using Pearson correlation and NMSE as figures of merit, and the optimal configuration was carried forward to all subsequent experiments.

5.2.1 K-space apodization window

Four apodization functions were compared: a rectangular window (i.e., no apodization), the Hann, Hamming, and Blackman windows. The left panel of Figure 5 shows the mean Pearson correlation achieved by each. The rectangular window produced the lowest correlation at r̄=0.8732, confirming that unwindowed DFT matrices introduce significant Gibbs-like aliasing that degrades agreement with k-Wave. All three smooth windows improved correlation substantially and performed nearly identically; Hann: r̄=0.9187, Hamming: r̄=0.9183, Blackman: r̄=0.9197, with differences of less than 0.002 between them. The Hann window was selected for all subsequent experiments as it is the conventional standard for k-space apodization in pseudospectral methods, and its performance is indistinguishable from the alternatives within the measurement noise.

5.2.2 Gaussian temporal damping (regularization ablation)

The right panel of the Figure 5 shows an ablation study examining the individual and combined contributions of k-space apodization and Gaussian temporal damping, quantified by NMSE against k-Wave. The unregularized baseline yielded NMSE = 0.496. Applying Hann apodization alone increased NMSE to 0.573; a result that, while appearing counterintuitive, is consistent with the known behavior of k-space windowing: the Hann taper reduces signal amplitude at high spatial frequencies, lowering the total signal energy and thereby inflating the energy-normalized NMSE, even though it improves the waveform shape correlation reported in Section 5.2.1. Applying Gaussian temporal damping alone (τ = 0.9T) reduced NMSE to 0.420, demonstrating that late-time oscillation suppression has a more direct effect on energy-based error than apodization does. Applying both regularizations together yielded the lowest NMSE of 0.254, a 49% reduction from the raw baseline, confirming that the two regularizations address complementary sources of model error and act synergistically. The combined configuration was retained as the final H-matrix for all subsequent experiments.

5.2.3 Source position dependence

The validation in Section 5.1 reported mean statistics across all sensors for a single centrally placed source. To assess whether forward model accuracy is maintained for sources at different spatial positions within the imaging domain, three source locations were evaluated: center (32, 32), off-center (16, 16), and near-edge (8, 8). As shown in the Figure 6 (left), the mean correlation remained high across all three positions, r̄=0.9907, 0.9756, and 0.9395 respectively, all well above the acceptance threshold of r_*s* = 0.80. The modest decrease for near-edge sources is expected, as sources close to the domain boundary generate wavefronts that interact more strongly with the PML in k-Wave, producing a slightly larger free-space mismatch. Nonetheless, the consistently high correlation across all positions confirms that the H-matrix provides a reliable forward model across the full imaging domain.

5.2.4 Grid spacing selection

The imaging grid spacing was selected using the procedure described in Section 3.6, evaluating three candidate spacings, d_*x* = 0.075, 0.100, and 0.125 mm, against the k-Wave reference. Agreement was assessed using three complementary measures: the mean per-sensor Pearson correlation, the minimum correlation across the array (a worst-case indicator), and the fraction of sensors exceeding the acceptance threshold of r = 0.80.

In this grid-spacing study the number of grid points (N_*x* = N_*y* = 64, on a 70 × 70 sensor boundary) was held fixed while d_*x* was varied, so that the physical field of view scales with d_*x*; for each spacing the source positions and sensor coordinates were derived from the same grid via the k-Wave grid2cart mapping and were therefore identical between the analytical model and the k-Wave reference, so each configuration is a self-consistent geometry and the comparison isolates discretization fidelity rather than a change of aperture. The spacing d_*x* = 0.100 mm gave the best agreement on all three measures, achieving a mean correlation of r = 0.919 ± 0.049 with all 71 sensors above the r = 0.80 threshold. The finer 0.075 mm grid did not improve agreement despite its higher spatial resolution, because a smaller d_*x* reduces the CFL-stable time step and therefore increases the number of time samples N_*t* required to capture wave propagation across the domain, it increases the number of late-time samples over which the free-space analytical model and the PML-bounded k-Wave reference must agree, and these additional late-time samples are precisely where the boundary mismatch accumulates. The coarser 0.125 mm grid, conversely, shortened this window but sacrificed spatial resolution without a compensating gain in fidelity. The intermediate spacing d_*x* = 0.100 mm thus represents the best balance between resolution and forward-model accuracy, and was adopted for all subsequent experiments.

5.2.5 Singular value compaction

Figure 6 (right), shows the cumulative energy as a function of the number of singular values for both the raw H-matrix (H_raw) and the regularized matrix (H_final, Hann + damping). The number of singular values required to capture 99% of the total energy, denoted k₉₉, decreased from 3,744 for H_raw to 1,139 for H_final, representing a 3.3 × reduction in k₉₉. We emphasize, however, that energy compaction is distinct from numerical conditioning. Computed directly from the singular values, the two-norm condition number is κ(H_raw)≈15 and κ(Hfinal)≈2.3×1010: the k-space apodization and Gaussian temporal damping deliberately attenuate high-spatial-frequency and late-time components, introducing small trailing singular values that increase κ. The regularization, therefore, does not improve the condition number; rather, it concentrates the operator's energy into a low-dimensional dominant subspace. This compaction, not conditioning, is the property relevant to compressive-sensing recovery, where the sparse vessel structure resides in the dominant modes while the attenuated near-null modes are precisely those the ℓ₁/total-variation regularizers suppress during the solve.

5.3 Baseline reconstruction methods

Three baseline reconstruction methods were evaluated on the held-out test set of 46 samples to establish reference performance levels against which the proposed STRN is compared. The first two methods, ISTA and SBTV, are classical iterative compressive sensing algorithms that exploit the sparsity of x₀ and rely on repeated multiplication by H and H^*T* at every iteration. The third, LISTA, is a learned iterative method that unrolls a fixed number of ISTA iterations into a trainable network, replacing hand-tuned thresholds with learned parameters. All three methods use the validated H-matrix described in Section 3.2 as the forward operator. Reconstruction quality is reported in terms of NMSE, PSNR, SSIM, and Pearson correlation r̄img, and the metric distributions across the test set are shown in Figure 7.

5.3.1 ISTA

ISTA was run with regularization parameter λ = 0.009 and a maximum of 500 iterations, at which point the cost function had converged for all test samples (final cost 2.432 × 10¹). Across the 46 test samples, ISTA achieved a mean SSIM of 0.2116 ± 0.0888, PSNR of 10.49 ± 2.28 dB, NMSE of 3.3130 ± 4.6940, and r̄img of 0.9358 ± 0.0228. The high mean Pearson correlation indicates that the gross spatial structure of the image is recovered faithfully, but the low SSIM reflects the presence of streak artifacts and spatial smearing that degrade structural fidelity at the fine scale. The large standard deviation of the NMSE (σ = 4.6940) reveals considerable variation in reconstruction quality across test samples, with a small number of samples producing substantially higher error than the majority, as visible in the heavy right tail of the NMSE histogram. Original test samples are shown in Figure 8 and representative visual results shown in Figure 9a confirm that ISTA recovers the coarse extent of the vessel distribution but introduces ringing artifacts along high-contrast edges and fails to resolve fine vessel detail.

5.3.2 SBTV

SBTV was run with the optimal parameters identified during the parameter search (λ = 0.1000, μ = 1.0; 76 outer iterations). Across the 46 test samples, SBTV achieved a mean SSIM of 0.1733 ± 0.0750, PSNR of 9.36 ± 1.65 dB, NMSE of 4.1700 ± 6.0986, and r̄img of 0.8557 ± 0.0686. SBTV performed below ISTA on all four metrics. The lower Pearson correlation relative to ISTA (0.856 vs. 0.936) suggests that the total variation regularization, while effective at suppressing noise in standard imaging problems, introduces staircase artifacts in the vessel images used here, distorting the spatial waveform structure. The extremely large NMSE standard deviation (σ = 6.099) indicates that performance is highly image-dependent, with several samples producing NMSE values exceeding 15. The test images used are shown in Figure 8 and visual inspection of Figure 9b confirms that SBTV reconstructions exhibit over-smoothing and loss of vessel boundary sharpness relative to ISTA.

5.3.3 LISTA

LISTA, as a learned iterative method, achieved substantially higher reconstruction quality than both classical baselines. Across the 46 test samples, LISTA achieved a mean SSIM of 0.5020 ± 0.1803, PSNR of 14.26 ± 3.25 dB, NMSE of 0.8053 ± 0.1343, and r̄img of 0.4249 ± 0.1472. The SSIM improvement over ISTA is 0.290 in absolute terms, representing a 2.37 × relative gain. The lower Pearson correlation of LISTA relative to ISTA and SBTV (0.425 vs. 0.936 and 0.856) is noteworthy and reflects a fundamental difference in reconstruction character: LISTA produces spatially compact, sparse reconstructions that recover the vessel locations with good structural fidelity but do not reproduce the full spatial extent of the ground truth signal in the way that the classical methods do. The test images used are shown in Figure 8, and this trade-off is visible in Figure 9c, where LISTA reconstructions show sharper, more localized vessel structures compared to the diffuse outputs of ISTA and SBTV, though at the cost of incomplete coverage for complex multi-vessel phantoms. The NMSE standard deviation for LISTA (σ = 0.1343) is substantially lower than for ISTA and SBTV, indicating more consistent reconstruction quality across test samples. The results of all three baseline methods are summarized alongside the STRN in Table 2, presented in Section 5.5.

Table 2. Reconstruction performance of all four methods on the 46-sample test set.

Mean ± standard deviation. The best value in each column is shown in bold; ties or non-significant differences are not bolded.

5.4 Sensor-token reconstruction network: training and validation

The STRN was trained end-to-end on the 488 augmented training samples and monitored on the 44-sample validation set. The complete network and training configuration is summarized in Table 3.

Table 3. Sensor-token reconstruction network configuration.

All weights were initialized using He initialization [36]. The vessel-weighted MSE loss of Equation 16 was minimized with Adam, with the learning rate halved every 80 epochs. The checkpoint with the lowest validation loss across the 300-epoch schedule was retained for evaluation. All models were implemented in MATLAB R2025b using the Deep Learning Toolbox; the analytical H-matrix and the k-Wave reference simulations used the k-Wave toolbox [34] (version 1.4). Network training and inference were performed on a single commodity laptop CPU (Intel Core i7-1255U, 16 GB RAM); no GPU acceleration was used, and all reconstruction algorithms were evaluated on the same hardware.

5.4.1 Test-set performance

The trained network was evaluated on the 46-sample held-out test set, which was not used at any stage of training or model selection. Across the test set, the STRN achieved a mean SSIM of 0.522 ± 0.225, a mean PSNR of 22.09 ± 7.42 dB, a mean NMSE of 0.233 ± 0.172, a mean r̄img of 0.836 ± 0.103, and a mean coefficient of determination R² of 0.654 ± 0.237. Of the 46 test samples, 21 (46%) achieved SSIM >0.5 and 10 (22%) achieved SSIM >0.7, with the best sample reaching SSIM = 0.965. The standard deviation of 0.225 in SSIM reflects substantial sample-to-sample variation: the network performs strongly on phantoms whose vessel layout is well-represented in the training distribution and less well on outlier topologies, as visible in the lower tail of the distribution in Figure 10.

5.4.2 Qualitative results

Figure 11 shows reconstructions for six representative test samples (top row: ground truth; bottom row: STRN reconstruction). The network recovers the position, orientation, and branching topology of the vessel structures faithfully across all six cases, with sharp boundaries and well-localized reconstructions. Minor amplitude attenuation is observed at the thinnest vessel segments—a consequence of the limited spatial-frequency content of the training distribution and the smoothing effect of the sigmoid output activation—but the geometric fidelity remains high. Beyond reconstruction quality, the proposed network offers a substantial inference-cost advantage. Each reconstruction is produced in a single forward pass through the trained network, with no H-matrix multiplication required, and completes in a single forward pass on commodity CPU hardware. By contrast, the iterative baselines reported in Section 5.3 require either 500 (ISTA), 76 outer (SBTV), or 16 (LISTA) matrix-vector products of cost O(MN)≈4.65 × 10⁷ floating-point operations each, leading to inference times that are at least one order of magnitude longer.

5.5 STRN vs. baselines

This section presents a direct comparison of the four reconstruction methods evaluated in this work—ISTA, SBTV, LISTA, and the proposed STRN—on the same 46-sample held-out test set, using the validated H-matrix described in Section 3.2 as the forward operator. All metrics are computed on images normalized to the range [0, 1]. Aggregate results are summarized in Table 2.

The STRN achieves the highest mean SSIM and PSNR, and the lowest mean NMSE, of the four methods. To establish whether these differences are statistically significant given the modest test-set size (n = 46), each pairwise comparison was assessed using both a paired t-test on the per-sample metric values and a Wilcoxon signed-rank test, the latter providing a non-parametric check that does not assume normally distributed differences. Results are summarized in Table 4.

Table 4. Paired statistical comparison of the STRN against each baseline on the test set (n = 46).

The mean difference is computed as (STRN - baseline) for SSIM, PSNR and Pearson r̄img (where higher is better) and as (baseline - STRN) for NMSE (where lower is better), so that a positive value always favors the STRN. Significance codes: ^***p < 0.001, ^**p < 0.01, ^*p < 0.05, not significant (n.s.) p≥0.05.

Three observations follow from this analysis. First, the STRN dominates the classical iterative baselines (ISTA and SBTV) on every fidelity metric. The improvements in SSIM (+0.31 over ISTA, +0.35 over SBTV), PSNR (+11.6 dB and +12.7 dB), and NMSE (3.1 and 3.9 lower respectively) are all highly statistically significant (p < 10⁻⁴) under both the parametric and non-parametric tests. The single exception is Pearson correlation, where ISTA retains a small but significant lead (Δr̄img=-0.10,p<10-8) and SBTV is statistically indistinguishable from the STRN [p = 0.20; not significant (n.s.)]. The high Pearson correlation of ISTA reflects its well-known behavior on sparse signals: the algorithm recovers the gross spatial layout of the vessels with strong linear similarity to the ground truth while producing diffuse, low-amplitude reconstructions whose pixel-level structural fidelity (SSIM = 0.21) is poor. The STRN, by contrast, achieves only a modestly lower Pearson r̄img yet a substantially higher SSIM and PSNR, indicating reconstructions that are structurally faithful as well as numerically accurate. Second, the comparison against LISTA, the strongest baseline, establishes the proposed network's principal contribution. On SSIM, the two methods are statistically indistinguishable (Δ = +0.020, p = 0.55, n.s.), reflecting that both produce structurally compact reconstructions of comparable quality at the single-scale SSIM level. On every other metric, however, the STRN is decisively superior: PSNR is 7.83 dB higher (p < 10⁻⁸), NMSE is 3.5 × lower (p < 10⁻²¹), and Pearson r̄img is 0.41 higher (p < 10⁻²⁰). The lower Pearson r̄img of LISTA is consistent with a known limitation of unrolled networks trained with limited data: the learned proximal operator applies too strong a sparsifying threshold at each iteration, leaving incomplete coverage on multi-vessel phantoms that the STRN reconstructs more completely. Third, the qualitative comparison is shown in Figure 12, which places the ground truth and all four reconstructions (ISTA, SBTV, LISTA, and STRN) side by side on identical test samples and a common intensity scale, together with matched absolute-error maps. ISTA and SBTV recover the coarse vessel distribution but exhibit smearing and ringing along high-contrast edges, with SBTV additionally producing piecewise-constant “staircase” artifacts characteristic of total-variation regularization. LISTA produces sharper, more spatially compact reconstructions than the classical methods but tends to miss vessel segments far from the dominant signal centroid, leaving incomplete coverage on multi-vessel phantoms, visible quantitatively as the substantially lower Pearson correlation despite high SSIM. The STRN reconstructions in Figure 11 match the ground-truth vessel layout closely on most samples, with sharp boundaries and faithful branching topology. Failure modes are restricted to phantoms whose vessel patterns lie outside the augmented training distribution, corresponding to the lower tail of the SSIM histogram in Figure 10. In summary, the proposed STRN produces reconstructions of significantly higher quality than the classical iterative methods (ISTA, SBTV) on every fidelity metric, and matches or exceeds the strongest learned baseline (LISTA) on every metric except SSIM, where the two are statistically indistinguishable. This is achieved while bypassing the H-matrix entirely at inference time, providing an order-of-magnitude reduction in reconstruction time relative to all three iterative baselines.

5.6 Generalization to independently simulated k-wave measurements

A recognized risk in learning-based reconstruction is the inverse crime: when a network is trained and tested on data generated by the same discretized operator, its reported fidelity may reflect learning to invert that specific operator rather than reconstructing the underlying physics. Because the STRN is trained on analytical H-matrix measurements, this concern is addressed directly by evaluating the trained network, without modification, on measurements produced by the k-Wave pseudo-spectral solver, an independent forward model that shares the acquisition geometry but none of the analytical operator's discretization (the analytical model having been validated against k-Wave in Section 5.1). The evaluation harness was validated by confirming that, on analytical inputs, it reproduces the published test-set means of Section 5.4.

Table 5 and Figure 13a report the zero-shot outcome. Transferred without any retraining, the STRN's mean SSIM falls from 0.522 on in-distribution analytical data to 0.192 on the k-Wave data, while the Pearson correlation is largely preserved (0.836 → 0.786, a 94% retention). The gross spatial mapping therefore transfers across forward models, whereas fine amplitude calibration does not, a real, quantified domain gap (all differences p ≤ 10⁻⁴) rather than the near-perfect transfer that an inverse crime would produce.

Table 5. Zero-shot generalization of the trained STRN from analytical training data to independently k-Wave-simulated measurements (same 46 held-out phantoms, matched geometry, no retraining).

The gap is, however, rapidly closable. Fine-tuning the trained network on a small number of k-Wave samples (learning rate 10⁻⁴, initialized from the pre-trained weights) recovers fidelity quickly: as shown in Table 6 and Figure 13b, approximately 20 samples nearly double the zero-shot SSIM, and approximately 160 samples recover roughly 90% of the in-distribution SSIM. The phantoms used for fine-tuning are drawn from the training split and are disjoint from the 46 held-out test phantoms at the original-image level; adaptation fidelity is always evaluated on the k-Wave test set, so no test phantom is seen during fine-tuning. Taken together, these results support interpreting the analytical H-matrix as an efficient training surrogate whose gap to full-wave physics is substantial in SSIM yet structure-preserving and can be closed with a small number of adaptation samples, rather than as evidence that the reported accuracy is an artifact of a shared operator.

Table 6. Domain-adaptation recovery curve: mean SSIM on the k-Wave test set after fine-tuning the pre-trained STRN on a small number of k-Wave samples.

5.7 Architectural ablation

To identify which components of the architecture are responsible for its performance and, in particular, to test whether the self-attention sub-layers are the source of the reported fidelity, a set of parameter-matched ablations was carried out under identical training conditions (150 epochs, three random seeds, the same 46-sample test set; plain mean-squared-error loss unless stated). Results are reported as mean ± standard deviation over seeds to control for initialization variance, which is non-negligible at this training budget. These matched-budget figures are lower in absolute terms than the 300-epoch reported model of Section 5.4 and are intended for within-experiment comparison only.

The central result is that removing self-attention does not degrade reconstruction quality. Replacing the multi-head self-attention sub-layer with a parameter-matched position-wise MLP (identical parameter count, no inter-sensor mixing) matches or exceeds the full attention model on every metric and every seed (Table 7 and Figure 14a: SSIM 0.479 ± 0.003 vs. 0.458 ± 0.012). The dominant parameter block is the fully connected spatial decoder (~9M of the ~18M parameters), and it, not the attention mechanism, accounts for the fixed-array reconstruction fidelity. A decoder-size sweep confirmed this: the attention-minus-MLP advantage was −0.029, −0.005, and +0.003 SSIM for large, medium, and small decoders, respectively, so attention provided a measurable benefit only when the decoder was starved to a point at which the model no longer reconstructed usefully; the best overall configuration (large decoder, attention-free, and plain MSE) reached SSIM 0.573.

Table 7. Parameter-matched ablation (150 epochs, same 46-sample test set).

SSIM is reported as mean ± standard deviation over three seeds; PSNR, NMSE, and Pearson are reported as the seed mean. The attention-free MLP matches or exceeds the full attention model; the vessel-weighted loss does not help.

Two further findings sharpen the picture. First, the vessel-weighted loss used to train the reported model does not help under matched conditions: plain MSE yields a higher SSIM (0.546 ± 0.006 vs. 0.458 ± 0.012), so the vessel weighting is not claimed as a contribution. Second, sensor identity, not attention, is what the network exploits: positional encoding is essential (removing it lowers SSIM to 0.438, and an index-based encoding outperforms a physical-coordinate encoding, 0.52 vs. 0.44), and randomly permuting the sensor ordering at test time collapses SSIM from 0.522 to 0.074. A direct convolutional sensor-to-image baseline reached SSIM 0.382, below the token-based models but above ISTA (0.212), supporting the sensor-as-token formulation over a naive two-dimensional treatment of the measurement array.

Finally, the robustness of attention to missing sensors was tested explicitly, since inter-sensor mixing might be expected to help when part of the array is unavailable. Masked attention (attending only over present sensors) was compared with the matched MLP, both trained with random 0%–60% sensor dropout, over three seeds and all four metrics (Table 8; Figures 14b,c). Attention retains a small SSIM edge at low missing fractions that vanishes by 50%–60%, but the matched MLP is better on PSNR, NMSE, and Pearson at every missing fraction, with its lead widening as more sensors are removed. Attention, therefore, confers no defensible robustness advantage under missing sensors either. An initial single-seed screen had suggested the opposite; it did not replicate under three seeds, masking, and the full metric set, and is not relied upon.

Table 8. Missing-sensor robustness: masked attention vs. matched MLP, each separately trained with random 0%–60% sensor dropout (mean over 3 seeds).

Entries are attention/MLP. The matched MLP leads on PSNR, NMSE, and Pearson at every missing fraction.

Collectively, these ablations motivate the framing adopted throughout this article: the STRN is presented as a direct, matrix-free reconstruction architecture whose accuracy is driven by the learned decoder over a sensor-token representation, and the self-attention sub-layers are reported transparently as a component that does not, in this fixed-array setting, earn its additional cost. This negative result is itself useful, as it points toward simpler and faster direct-reconstruction designs for this class of problem.

5.8 Robustness to measurement perturbations

Robustness to realistic measurement imperfections is essential for practical deployment, since experimental acquisitions inevitably depart from the idealized noise-free measurements used for training. Three perturbation families were applied to the raw sensor measurements of the 46-sample test set, and the trained STRN was evaluated without retraining (three random realizations per sample; the unperturbed level of each family reproduces the published test-set SSIM of 0.522 exactly, providing a built-in consistency check). Additive zero-mean Gaussian noise was applied at signal-to-noise ratios of 40–10 dB; per-channel multiplicative gain perturbations, modeling detector-sensitivity variation, at standard deviations of 5%–20%; and per-channel integer timing jitter, modeling trigger and clock jitter, at standard deviations of 1–3 samples. Table 9 and Figure 15 summarize the outcome.

Table 9. Robustness of the trained STRN to measurement perturbations (mean over three realizations × 46 test samples; no retraining).

The unperturbed level of each family reproduces the published test-set values.

The network is most robust to detector-gain variation: a 20% per-channel gain spread reduces SSIM by only 0.036 (0.522 → 0.486). It is moderately robust to additive noise, retaining most of its fidelity down to 20 dB (SSIM 0.465) and degrading substantially only at 10 dB (0.195). It is most sensitive to timing jitter: a one-sample per-channel jitter already lowers SSIM to 0.461, and a three-sample jitter to 0.358. This ordering is consistent with the mechanism of a direct time-of-flight mapping; the network exploits the precise arrival-time structure of the sensor traces, so per-channel temporal misalignment is more damaging than amplitude scaling or additive noise. In practice, this identifies inter-channel timing calibration as the dominant requirement for deploying the trained model on experimental data, and motivates timing-jitter augmentation during training as a route to further robustness. The perturbations here characterize the proposed model; the same harness supports the classical and unrolled baselines for a matched comparison.

6 Discussion

The principal finding of this work is that a Transformer-based sensor-as-token architecture can produce PAT reconstructions of significantly higher fidelity than the classical iterative baselines (ISTA and SBTV) and at least comparable to the strongest unrolled-network baseline (LISTA), while bypassing the H-matrix at inference. The paired statistical comparison in Section 5.5 places this claim on a rigorous footing: STRN attains improvements of +11.6 dB PSNR over ISTA, +12.7 dB over SBTV, and +7.83 dB over LISTA, with each comparison significant at p < 10⁻⁸ under both parametric (paired t-test) and non-parametric (Wilcoxon signed-rank) tests. The NMSE reduction relative to LISTA, the strongest baseline, is a factor of 3.5 × (p < 10⁻²¹), establishing the proposed network as a clear quantitative improvement over the prior state of the art at the test-set scale studied here [11].

6.1 Interpretation of the reconstruction trade-offs

Three observations from Section 5.5 warrant further discussion. First, the Pearson correlation coefficient retains a small advantage for ISTA (0.936 vs. STRN's 0.836, Δr̄img=-0.10,p<10-8). This is consistent with the well-known behavior of l₁-regularized iterative solvers on sparse signals [7]: such methods recover the gross spatial layout of the vessels with high linear similarity to the ground truth, but the resulting reconstructions are diffuse and low-amplitude, which depresses the single-scale SSIM (0.21 for ISTA). The STRN, by contrast, achieves only a modestly lower Pearson r̄img but a substantially higher SSIM (0.522), reflecting reconstructions that are simultaneously structurally faithful at the pixel level and numerically accurate in the energy sense. Second, on the single-scale SSIM, the STRN and LISTA are statistically indistinguishable (Δ = +0.020, p = 0.55, n.s.). This is itself an informative result: SSIM is sensitive to fine boundary structure, and LISTA, being a learned variant of ISTA, inherits its ability to produce sparse, well-localized reconstructions. The decisive separation between the two methods appears in the energy-based metrics (PSNR, NMSE) and in Pearson correlation, where STRN's coverage of the full vessel extent and amplitude is substantially better than LISTA's. This pattern is consistent with the known limitation of unrolled networks trained on limited data: the learned proximal operator can become overly aggressive at sparsifying the iterate, leaving incomplete coverage on multi-vessel phantoms. Third, the multi-scale structural similarity index (MS-SSIM) was deliberately excluded from the principal comparison table. On the same test set, ISTA attained a higher MS-SSIM (0.866) than the STRN (0.794) despite its much lower single-scale SSIM. This counterintuitive result is a consequence of MS-SSIM's weighting of coarse-scale structural similarity, which is well preserved by diffuse low-amplitude reconstructions such as those produced by ISTA. Since fine-scale boundary structure is the metric of greatest clinical interest for vessel imaging, the single-scale SSIM is the more discriminating measure for this task, and the MS-SSIM is reported only as a STRN-specific descriptor in Section 5.4.

6.2 What drives reconstruction quality

The ablation of Section 5.7 carries a direct interpretive consequence for how the present results should be read. The network's fidelity is attributable to the learned spatial decoder operating on a sensor-token representation with positional encoding, rather than to the self-attention sub-layers, which a parameter-matched attention-free variant matches or exceeds on every metric. We therefore frame the contribution as a direct, matrix-free reconstruction architecture and report the attention finding openly. This transparency is intended to be constructive: for a fixed detector array, comparably accurate reconstruction can be obtained with simpler and cheaper per-token designs, which is a useful guide for architecture selection in this setting. It also clarifies the scope of the claim, the present evidence concerns a single fixed array, and whether self-attention becomes advantageous under variable or partially populated geometries remains an open question.

6.3 Computational efficiency at inference

The computational advantage of the STRN at inference is independent of the reconstruction-quality argument and is therefore worth separating. Each iterative baseline performs at least one matrix-vector product with H per iteration, with cost O(MN)≈4.65 × 10⁷ floating-point operations for the chosen grid size of 64 × 64 pixels and N_*s*N_*t* = 11, 360 measurements. ISTA uses 500 such iterations per reconstruction; SBTV uses 76 outer iterations; LISTA uses 16 unrolled layers. The STRN, by contrast, performs a single forward pass dominated by dense matrix multiplications totalling on the order of 10⁷ operations. On the commodity CPU used for all experiments (Intel Core i7-1255U, no GPU acceleration), and evaluated on the same hardware as the iterative baselines, this single-pass cost is at least an order of magnitude below the cumulative per-iteration cost of ISTA (500 iterations), SBTV (76 outer iterations), and LISTA (16 unrolled layers), each of which performs at least one H matrix-vector product per iteration. We emphasize that this is an operation-count comparison at the fixed 64 × 64 resolution evaluated here and does not by itself establish real-time capability: the wall-clock inference latency and memory footprint at higher resolutions (for example 128 × 128 and 256 × 256) have not been characterized, and network inference latency is in any case only one component of a complete experimental acquisition-and-reconstruction pipeline. Characterizing these, and thereby establishing real-time performance in a full experimental setting, is left to future work.

6.4 Limitations

Three limitations of the present study deserve acknowledgment. First, the training corpus comprises 191 phantoms drawn from a single public source (Radiopaedia), expanded to 488 augmented samples by small spatial translations. While the held-out test set was strictly disjoint and the validation set was used only for checkpoint selection, the absolute generalization capacity of the network to phantoms outside this distribution remains unproven. Failure cases observed in the lower tail of the SSIM histogram (Figure 10) correspond predominantly to phantoms with vessel topologies poorly represented in the training distribution. Second, the analytical H-matrix is constructed under the assumption of a homogeneous, lossless acoustic medium with constant speed of sound. Real biological tissue exhibits acoustic heterogeneity, frequency-dependent attenuation, and in some cases multiple scattering, none of which are modeled in the present forward operator. Although the validation against k-Wave confirms that the analytical model faithfully reproduces idealized forward physics, the gap between the analytical training distribution and experimental measurements must be quantified before clinical translation. Third, the present study reports results for a single grid spacing (d_*x* = 0.100 mm), sensor configuration (N_*s* = 71 on a square boundary), and image size (64 × 64). The scaling of STRN performance with image resolution, sensor count, and grid spacing has not been characterized here; this is the subject of ongoing work.

6.5 Comparison to prior self-work

The proposed architecture is positioned as a complement, not a replacement, of the unrolled-network methods previously reported by the present authors [11, 14]. Unrolled networks remain valuable for problems in which the forward operator H changes between samples (e.g., different sensor geometries, patient-specific scan protocols) and physical interpretability of the per-layer iterate is required for clinical acceptance. The STRN is most appropriate when the sensor geometry is fixed across the patient population, and inference latency is the dominant constraint. Both approaches benefit from the validated analytical H-matrix introduced in Section 3, which can be used for training data generation regardless of the downstream reconstruction architecture.

7 Conclusion and future work

This article introduced the Sensor-Token Reconstruction Network (STRN), a direct, matrix-free architecture for photoacoustic tomography image reconstruction that maps raw sensor measurements to the reconstructed image in a single forward pass, without invoking the forward operator at inference. The contributions are fourfold: (i) an analytical k-space H-matrix was constructed and validated against the k-Wave pseudo-spectral solver under matched geometry, achieving a mean per-sensor Pearson correlation of r̄=0.919±0.049, with two targeted regularisations (Hann apodization and Gaussian temporal damping) reducing the energy-normalized mismatch by 49%; (ii) the STRN architecture was introduced and trained end-to-end on 488 augmented PAT phantoms; (iii) a rigorous statistical comparison against ISTA, SBTV, and LISTA on a 46-sample held-out test set established that the STRN achieves the highest mean SSIM (0.522), the highest mean PSNR (22.09 dB), and the lowest mean NMSE (0.233), with paired t-tests and Wilcoxon signed-rank tests confirming significance at p < 10⁻⁸ for the principal claims; and (iv) two validation analyses characterized the approach, an independent k-Wave generalization test quantified a substantial but structure-preserving and rapidly closable domain gap, demonstrating that the analytical operator is an efficient training surrogate rather than a source of an inverse crime, and a parameter-matched ablation showed transparently that the learned spatial decoder, rather than the self-attention sub-layers, is the dominant driver of reconstruction fidelity. By eliminating the per-iteration Hx multiplication at inference, the STRN reduces reconstruction latency by at least an order of magnitude relative to all three iterative baselines on the same hardware. Several directions for future work are identified. First, the training corpus should be expanded with more diverse vessel topologies and, where possible, with experimentally acquired sensor signals from phantom and ex vivo tissue measurements, in order to characterize the network's generalization capacity beyond the simulated regime [11]. Second, the data augmentation strategy can be made more physically faithful by applying the corresponding sensor-index permutation when rotating or reflecting the phantom, a procedure that is exact for a sensor array with discrete rotational symmetry and that avoids the analytical/numerical distribution mismatch encountered in preliminary augmented experiments. Third, the STRN architecture should be extended to three-dimensional PAT by replacing the two-dimensional spatial decoder with a 3-D volumetric decoder, and by accommodating cylindrical or hemispherical sensor geometries common in clinical scanners [37]. Moreover, a learned MLP encoding of the physical sensor coordinates can be considered as future work. Fourth, the integration of physics-informed loss terms, penalizing deviation of the reconstructed image from the analytical forward model under perturbed sensor inputs, offers a principled means of combining the data efficiency of unrolled networks with the inference speed of the direct architecture proposed here. Fifth, the present design should be evaluated on heterogeneous acoustic media and frequency-dependent attenuation, possibly with a learned correction layer to bridge the gap between the homogeneous analytical training distribution and clinical measurements. Finally, sixth, the network's transferability to related tomographic modalities, such as thermoacoustic tomography and ultrasound computed tomography, merits investigation, as the underlying sensor-as-token paradigm is not specific to PAT and may generalize to any inverse problem in which a fixed sensor array produces a sparse-view measurement vector.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Lixin Shen, Syracuse University, United States

Reviewed by: Sudeep Mondal, University of South Florida, United States

Di Kong, Tsinghua University, China

Footnotes

The acceptance threshold r_*s* = 0.80 is adopted as a conservative criterion for waveform agreement: r_*s*≥0.80 corresponds to at least 64% (=rs2) of the per-sensor signal variance being shared between the analytical and k-Wave time series, i.e. the analytical model reproduces the dominant temporal waveform used for training-data synthesis. The reported conclusions are insensitive to reasonable changes in this threshold.