Pith. sign in

REVIEW 4 major objections 5 minor 46 references

Photo-acoustic tomographic image reconstruction from reduced data using physically inspired regularization

T0 review · 4 major / 5 minor · reviewed 2026-08-14 · deepseek-v4-flash

Pith's one-line read Joint sparsity prior lifts limited-data PAT reconstruction quality

desk verdict A promising limited-data PAT reconstruction method with a useful matrix-free implementation, but the physical prior is asserted rather than demonstrated and the paper needs an ablation. read the letter →

arxiv 1908.02313 v2 pith:YFGEZGVX submitted 2019-08-06 eess.IV

classification eess.IV
keywords photoacoustictomographylimited-datareconstructionjointsparsitynon-convexregularizationgraduatednon-convexitysecond-orderderivativesmatrix-freeforwardmodeltotalvariation
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper proposes a model-based photoacoustic tomography (PAT) reconstruction method that recovers good-quality images from far fewer transducer measurements than standard methods require. Its central idea is a regularizer built around a structural property the authors report for PAT images: high pixel intensities and high second-order derivatives are jointly sparse. Because that regularizer is non-convex, the method couples it with a graduated non-convexity (GNC) optimization scheme and with a matrix-free implementation of the PAT forward model that keeps memory use manageable. In simulated and real experiments, the method reports higher structural-similarity and figure-of-merit scores than a FISTA-based total-variation baseline in every tested configuration, with the largest gains at 16–32 transducers.

What carries the argument

The central object is a non-convex regularizer that couples pixel intensity to second-order derivatives. The first form is $$R_{h,1}(p_0,q)=\sum_r \left(\epsilon+\$\alpha$(p_0)$_r^{2}$+(1-\$\alpha$)\sum_i (D_{2,i}p_0)$_r^{2}$\right)^q$$ with $q<0.5$; the second form separates the intensity and derivative terms before raising each to the power $q$. The fractional power $q$ is what makes the cost non-convex, and the graduated non-convexity outer loop starts at the convex value $q=0.5$ and decreases $q$ toward $0.25$, warm-starting each inner solve. The inner solver is a preconditioned gradient method whose preconditioner approximates a damped Newton step and is applied with conjugate gradients. Memory is controlled by a matrix-free forward model that uses the exact time-propagator relation $p(r,t)=\mathcal{F}^{-1}\{\hat{P}_0(k)\cos(c_0\|k\|t)\}$, so the operator $H^T H$ is applied as a sequence of Fourier multiplications rather than stored as a matrix.

What would settle it

Take an unseen class of PAT images, for example vascular or organ images with diffuse backgrounds, and compute the joint histogram of pixel intensity against second-derivative magnitude; if the high-intensity, high-derivative pixels are not markedly sparser than in generic images, the proposed method should fail to beat the FISTA baseline at 16–32 transducers. A complementary check is to rerun the 16-transducer experiments with the GNC loop disabled and $q$ fixed at 0.25; if artifacts return, the convergence schedule rather than the prior is doing the work.

Watch

Extended reading notes

Core claim

The paper's central claim is that a reconstruction prior designed around the physical structure of photoacoustic images, rather than a generic image prior, can recover accurate initial-pressure distributions from a fraction of the measurements normally required. Concretely, the authors assert that in PAT images high intensity and high second-order derivatives are jointly sparse, and that a regularizer encoding this joint sparsity, minimized by their custom graduated non-convex solver, produces higher SSIM and FOM scores than the FISTA-based total-variation baseline in all tested scenarios: 16, 32, 64, and 128 transducers at 20, 30, and 40 dB input noise, plus real horsehair-phantom data. The largest reported gains occur at 16–32 transducers, where the baseline reconstructions degrade visibly while the proposed method remains close to full-array quality.

Load-bearing premise

The load-bearing premise is the asserted structural property that in photoacoustic images bright pixels and pixels with large second-order derivatives tend to occupy the same sparse locations, a claim the paper supports by observation rather than quantitative statistics; if that property is absent for a given image class, the regularizer's advantage over generic total variation would disappear.

Editorial extensions

If this is right

  • PAT systems could use 16 or 32 transducers instead of 128 and still recover images close to full-array quality, which would shorten scan times and lower equipment cost.
  • The advantage persists from 20 dB to 40 dB input noise, and the largest margins occur at low SNR, so the method is useful where measurements are noisy or laser fluence is limited.
  • The real horsehair-phantom test gives roughly a 3 dB improvement in the peak-to-noise figure of merit over the FISTA baseline, indicating the benefit survives experimental transducer-response and setup effects.
  • Because the second regularizer form performs better on some phantoms and the first on others, the paper gives users a practical choice rather than a single fixed prior.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Extending the joint-sparsity idea beyond PAT, the same regularizer could transfer to other high-contrast tomographic modalities such as fluorescence microscopy or diffuse optical tomography, where bright structures and sharp boundaries are also sparse.
  • A consequence the paper does not develop: the matrix-free Fourier-cosine propagator assumes a homogeneous, lossless acoustic medium and transducer positions on the imaging grid, so adapting it to heterogeneous sound speed or arbitrary detector layouts would need a different fast propagator and would likely reduce the memory savings.
  • The reported parameter settings ($q=0.25$, $n_s=10$, $\alpha=0.5$, $\lambda_p=10\lambda$) were chosen as adequate across the tested cases, so a broader tuning study on unseen anatomies and noise levels would show whether the SSIM margins, especially the large gap at 16 transducers, persist or shrink.
  • The SSIM curves suggest an adaptive acquisition strategy: add transducer positions only until the reconstruction quality plateaus, rather than always using a fixed large array.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper proposes a model-based photoacoustic tomography (PAT) reconstruction method for reduced-size datasets. The key novelty is a non-convex regularizer that couples image intensity with second-order derivatives, motivated by an asserted 'joint sparsity' property of PAT images (Section 2.1, Eqs. 13-14). The authors develop a preconditioned gradient algorithm with graduated non-convexity (GNC) for the resulting non-convex problem, and derive a matrix-free implementation of the PAT forward model based on a Fourier-domain time propagator (Section 2.4, Eq. 34). The method is evaluated on three simulated phantoms (blood vessel, Derenzo, PAT) at 16-128 transducers and 20-40 dB noise, plus one real horsehair phantom, comparing against a FISTA-based total-variation method by Huang et al. Reported SSIM and FOM gains favor the proposed method, especially at very limited transducer counts.

Significance. If the reported gains hold, the paper contributes a practical reconstruction method for limited-data PAT, with the matrix-free forward-model implementation being a genuinely useful engineering contribution that removes a memory bottleneck for large images. The comparison across transducer counts and noise levels is reasonably extensive, and the inclusion of real measured data strengthens the claims. However, the conceptual novelty—the joint-sparsity prior—is asserted rather than demonstrated. Because the regularizer changes derivative order, convexity, and intensity-derivative coupling simultaneously, the current evidence does not establish that the physical prior is the cause of the improvements. The paper would be significantly strengthened by a controlled ablation and by specifying the parameter-selection protocol, particularly for the regularization weight, to ensure a fair comparison. The central idea is defensible, but the evidence as presented is not yet conclusive on the mechanism.

major comments (4)
  1. [Section 2.1; Eqs. (13)-(14)] The paper's central claim is that PAT images exhibit joint sparsity of high intensities and high second-order derivatives, and that this prior is responsible for the reconstruction gains. However, this property is only asserted ('we observe' in Section 2.1) without any quantitative evidence. No histograms, scatter plots, or statistical tests are provided to show that PAT images have this structure beyond what generic first-order or second-order smoothness priors would capture. Since the regularizer is specifically designed around this property, the reader cannot verify the premise. Please add quantitative evidence of the joint-sparsity property for the phantoms and real PAT images used, and discuss whether it holds for other PAT image classes (e.g., vasculature versus diffuse absorption).
  2. [Section 3.1; Table 1] There is no controlled ablation that isolates the proposed joint-sparsity mechanism. The proposed cost (Eq. 13) differs from the FISTA/TV-1 baseline along three axes simultaneously: second-order versus first-order derivatives, non-convexity with q=0.25 solved by GNC, and intensity-derivative coupling through α. The reported SSIM gains (Table 1) could therefore be caused entirely by higher-order smoothness, by non-convex regularization, or by the specific optimization procedure, rather than by the physical prior. Please include at least the following ablations: (i) α=0 and α=1 in Eq. (13), (ii) a convex control with q=0.5 in the same GNC framework, and (iii) a second-order TV or TGV baseline. Without these, the abstract's and conclusions' attribution of improvements to the joint-sparsity prior is not supported.
  3. [Section 3.1; parameter selection] The selection of the regularization parameter λ is described only as 'determined using the model itself' (Section 3.1), which is not a reproducible or falsifiable protocol. In contrast, the FISTA baseline is explicitly tuned for best SSIM ('λ chosen for best SSIM score' in Figures 3-7). If the proposed method's λ is chosen by a different criterion or by visual inspection, the comparison may be biased in its favor. Please specify the exact λ-selection procedure for both methods, and ideally use an identical criterion (e.g., best SSIM on a validation set, or a principled discrepancy principle) for both.
  4. [Section 2.2; Eq. (17)-(18)] The use of a fixed positivity penalty with λp = 10λ is justified only by 'a series of reconstruction trials' (Section 2.2). This is an ad hoc choice that may not transfer to other phantoms, noise levels, or transducer geometries. Since the positivity constraint is part of the model (Eq. 15), the sensitivity of the results to λp should be reported, or the penalty should be systematically set (e.g., via continuation). Otherwise, the reader cannot assess whether the reported gains are contingent on a manually tuned parameter.
minor comments (5)
  1. [Section 3.1; runtime] The paper reports that the proposed method takes about 38 minutes per reconstruction versus 30 minutes for FISTA on the same machine, but provides no iteration counts, convergence curves, or memory usage measurements. A brief table of computational cost (iterations, CG calls, memory footprint with the matrix-free implementation) would make the engineering contribution of Section 2.4 concrete.
  2. [Section 3.1; Figures 6-7] Figures 6 and 7 appear to show the same experiment (16 transducers, three noise levels) for different phantoms, but the captions do not specify which phantom is in each figure. Please make the captions explicit (e.g., 'PAT phantom' and 'Derenzo phantom').
  3. [Section 3.2; Eq. (35)] The FOM defined in Eq. (35) is computed as 20*log10(S/n), where n is the standard deviation of the 'intensity' of the reconstructed image. It is unclear whether n is computed over the whole image or over a background region. Since the FOM comparison is the only quantitative metric for the real-data experiment, please define the region of interest and the noise estimation method explicitly.
  4. [Section 2.3; Algorithm 3] The GNC schedule sets q_m = 0.5 - m(0.5-q)/ns, but the text does not describe how the choice of q=0.25 and ns=10 interacts with the line-search tolerances. A sentence explaining the robustness of the GNC schedule to the inner tolerances (as claimed) would help, since Algorithm 3 is central to the non-convex optimization.
  5. [Section 2.4; Eq. (34)] The matrix-free formula in Eq. (34) is a useful contribution, but the derivation assumes transducers located on image grid points. For the real-data experiment, the transducer is rotated continuously; please clarify how the discrete transducer positions are mapped to grid points, or how the formula is adapted for off-grid positions.

Circularity Check

0 steps flagged · score 0.0 of 10

No circular derivation: the joint-sparsity regularizer is proposed and tested, not fitted or relabeled as a prediction; the minor self-citation is not load-bearing.

full rationale

The paper's derivation chain consists of (i) asserting, not deriving, that PAT images have jointly sparse intensities and second-order derivatives; (ii) constructing the nonconvex regularizer in Eqs. (13)-(14) as a modification of Ref. [34]; (iii) deriving the gradient, preconditioned CG, and GNC optimization; and (iv) evaluating on simulated and real data. No equation in this chain is equivalent by construction to a measured quantity that is then called a prediction. The regularization parameters (alpha=0.5, q=0.25, lambda via the model) are user-selected hyperparameters, not fitted values renamed as outputs. The only self-citation, Ref. [34], supplies the original log-based fluorescence regularizer and the observation that the same sparsity pattern occurs in fluorescence images; the present paper explicitly modifies the log to a fractional power and tests the modified regularizer on PAT data, so the central claim does not reduce to that citation. The simulated-data validation uses the same Fourier-cosine propagator (Eq. (4)) to generate and reconstruct data, which is an inverse-crime-style limitation, but this is a benchmarking weakness, not a circular derivation; the real measured-data experiment provides independent evidence. The main load-bearing scientific weakness is that the asserted joint-sparsity property is supported only by visual assertion and no ablation isolates the proposed prior from second-order smoothness or nonconvexity; that is an evidence gap, not a self-referential reduction. Accordingly, no circular step is identified.

Assumptions & free parameters 5 free parameters · 5 assumptions · 0 invented entities

The method relies on standard PAT assumptions (homogeneous medium, delta-pulse source, point transducers) and on the asserted joint-sparsity property. The hyperparameters lambda, q, alpha, lambda_p, and ns are user-chosen; no new physical entities are introduced.

free parameters (5)
  • lambda (regularization weight) = not stated; determined using the model itself (Sec. 3.1)
    Controls trade-off between data fidelity and regularization; selection method not specified.
  • q (sparsity index) = 0.25
    Chosen as 'adequate for all cases' (Sec. 3.1); determines degree of non-convexity.
  • alpha (relative penalization) = 0.5
    Weight between intensity and second-order derivative terms; set to 0.5 for all test cases.
  • lambda_p (positivity penalty weight) = 10 times lambda
    Chosen by trial as sufficient to avoid significant negative values without numerical instability (Sec. 2.3).
  • ns (number of GNC steps) = 10
    Number of graduated non-convexity steps; reported as adequate for all cases (Sec. 3.1).
assumptions (5)
  • domain assumption 2D homogeneous, lossless, non-dispersive acoustic medium with point detectors
    Forward model Eq. (4) and wave equation Eq. (1) assume these conditions; deviations cause model error.
  • domain assumption Photoacoustic source is approximated as H(r)delta(t) under stress and thermal confinement
    Stated in Sec. 1.1; standard in PAT.
  • domain assumption PAT images have jointly sparse high intensity and high second-order derivatives
    Observed property claimed in Sec. 2.1 without quantitative evidence; motivates the regularization.
  • domain assumption Transducers lie on image grid points and sample pressure at points
    Eq. (34) models transducer readings as multiplication by binary mask S(r); ignores finite aperture and off-grid positions.
  • ad hoc to paper Quadratic penalty with lambda_p = 10 lambda suffices to enforce non-negativity
    Sec. 2.3 states it is a heuristic found through reconstruction trials; not analytically justified.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Photo-acoustic tomographic image reconstruction from reduced data using physically inspired regularization." pith.science (2026). https://pith.science/paper/YFGEZGVX

@misc{pith2026190802313,
  author       = {Pith},
  title        = {Pith review of: Photo-acoustic tomographic image reconstruction from reduced data using physically inspired regularization},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/YFGEZGVX}},
  note         = {Machine review of arXiv:1908.02313}
}
read the original abstract

We propose a model-based image reconstruction method for photoacoustic tomography(PAT) involving a novel form of regularization and demonstrate its ability to recover good quality images from significantly reduced size datasets. The regularization is constructed to suit the physical structure of typical PAT images. We construct it by combining second-order derivatives and intensity into a non-convex form to exploit a structural property of PAT images that we observe: in PAT images, high intensities and high second-order derivatives are jointly sparse. The specific form of regularization constructed here is a modification of the form proposed for fluorescence image restoration. This regularization is combined with a data fidelity cost, and the required image is obtained as the minimizer of this cost. As this regularization is non-convex, the efficiency of the minimization method is crucial in obtaining artifact-free reconstructions. We develop a custom minimization method for efficiently handling this non-convex minimization problem. Further, as non-convex minimization requires a large number of iterations and the PAT forward model in the data-fidelity term has to be applied in the iterations, we propose a computational structure for efficient implementation of the forward model with reduced memory requirements. We evaluate the proposed method on both simulated and real measured data sets and compare them with a recent reconstruction method that is based on a well-known fast iterative shrinkage threshold algorithm (FISTA).

Figures

Figures reproduced from arXiv: 1908.02313 by the authors.

Figure 1
Figure 1. Numerical phantoms used for evaluation of the proposed method [PITH_FULL_IMAGE:figures/full_fig_p008_1.png] view at source ↗
Figure 2
Figure 2. a) Schematic diagram of PA data acquisition geometry with ultrasound transducers (shown by dots) around [PITH_FULL_IMAGE:figures/full_fig_p008_2.png] view at source ↗
Figure 3
Figure 3. Comparison of reconstructions obtained from simulated data with 128 transducers and 20 dB measurement [PITH_FULL_IMAGE:figures/full_fig_p010_3.png] view at source ↗
Figures from the paper (6 more)
Figure 4
Figure 4. Figure 4: Scan line based intensity profiles of reconstructed images from Fig:3. [PITH_FULL_IMAGE:figures/full_fig_p010_4.png]
Figure 5
Figure 5. Figure 5: Comparison of reconstructions obtained from simulated data with 32 transducers and 20 dB measurement [PITH_FULL_IMAGE:figures/full_fig_p011_5.png]
Figure 6
Figure 6. Figure 6: F respectively. However, the reconstruction quality was not much affected by the proposed method, although [PITH_FULL_IMAGE:figures/full_fig_p011_6.png]
Figure 6
Figure 6. Figure 6: Reconstructed images from input datasets corresponding to 16 transducers for various noise levels. (A), (B), [PITH_FULL_IMAGE:figures/full_fig_p012_6.png]
Figure 7
Figure 7. Figure 7: Reconstructed images from input data sets corresponding to 16 transducers for various noise levels. (A), (B), [PITH_FULL_IMAGE:figures/full_fig_p012_7.png]
Figure 8
Figure 8. Figure 8: Reconstructed images from horse hair phantom data using 400 transducers. (A): reconstruction from the [PITH_FULL_IMAGE:figures/full_fig_p013_8.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

46 extracted references · 46 canonical work pages

  1. [1]

    Wang L V and Yao J 2016 Nature methods 13 627

  2. [2]

    Pramanik M, Ku G, Li C and Wang L V 2008 Medical physics 35 2218–2223

  3. [3]

    Beard P 2011 Interface focus 1 602–631

  4. [4]

    Zhou Y , Yao J and Wang L V 2016Journal of biomedical optics 21 061007

  5. [5]

    Upputuri P K and Pramanik M 2016 Journal of Biomedical Optics 22 041006

  6. [6]

    2017 Nature biomedical engineering 1 0071

    Li L, Zhu L, Ma C, Lin L, Yao J, Wang L, Maslov K, Zhang R, Chen W, Shi J et al. 2017 Nature biomedical engineering 1 0071

  7. [7]

    Han Y , Ding L, Ben X L D, Razansky D, Prakash J and Ntziachristos V 2017Opt. Lett. 42 979–982

  8. [8]

    Rosenthal A, Ntziachristos V and Razansky D 2013 Current medical imaging reviews 9 318–336

Show all 46 references
  1. [9]

    Finch D and Patch S K 2004 SIAM journal on mathematical analysis 35 1213–1240

  2. [10]

    Kunyansky L A 2007 Inverse problems 23 373

  3. [11]

    Xu M and Wang L V 2002 IEEE transactions on medical imaging 21 814–822

  4. [12]

    Xu M and Wang L V 2005 Physical Review E 71 016706

  5. [13]

    Kalva S K and Pramanik M 2016 Journal of Biomedical Optics 21 086011

  6. [14]

    Pramanik M 2014 JOSA A 31 621–627

  7. [15]

    Xu Y and Wang L V 2004 Physical review letters 92 033902

  8. [16]

    Burgholzer P, Matt G J, Haltmeier M and Paltauf G 2007 Physical Review E 75 046706

  9. [17]

    Treeby B E, Zhang E Z and Cox B T 2010 Inverse Problems 26 115003

  10. [18]

    Hristova Y , Kuchment P and Nguyen L 2008Inverse Problems 24 055006

  11. [19]

    Xu Y , Wang L V , Ambartsoumian G and Kuchment P 2004Medical physics 31 724–733

  12. [20]

    Arridge S, Beard P, Betcke M, Cox B, Huynh N, Lucka F, Ogunlade O and Zhang E 2016Physics in Medicine & Biology 61 8908 14 A PREPRINT - D ECEMBER 29, 2020

  13. [21]

    Buehler A, Rosenthal A, Jetzfellner T, Dima A, Razansky D and Ntziachristos V 2011 Medical physics 38 1694–1704

  14. [22]

    Paltauf G, Viator J, Prahl S and Jacques S 2002 The Journal of the Acoustical Society of America 112 1536–1544

  15. [23]

    Schwab J, Pereverzyev Jr S and Haltmeier M 2018 SIAM Journal on Numerical Analysis 56 160–184

  16. [24]

    Prakash J, Mandal S, Razansky D and Ntziachristos V 2019 IEEE Transactions on Biomedical Engineering 1–1 ISSN 0018-9294

  17. [25]

    Huang C, Wang K, Nie L, Wang L V and Anastasio M A 2013 IEEE Transactions on Medical Imaging 32 1097–1110

  18. [26]

    Arridge S R, Betcke M M, Cox B T, Lucka F and Treeby B E 2016 Inverse Problems 32 115012

  19. [27]

    Boink Y E, Lagerwerf M J, Steenbergen W, van Gils S A, Manohar S and Brune C 2018Physics in Medicine & Biology 63 045018

  20. [28]

    Saratoon T, Tarvainen T, Cox B and Arridge S 2013 Inverse Problems 29 075006

  21. [29]

    Han Y , Tzoumas S, Nunes A, Ntziachristos V and Rosenthal A 2015Medical physics 42 5444–5452

  22. [30]

    Wang K, Huang C, Kao Y J, Chou C Y , Oraevsky A A and Anastasio M A 2013Medical physics 40

  23. [31]

    Rosenthal A, Razansky D and Ntziachristos V 2010 IEEE transactions on medical imaging 29 1275–1285

  24. [32]

    Rosenthal A, Jetzfellner T, Razansky D and Ntziachristos V 2012 IEEE transactions on medical imaging 31 1346–1357

  25. [33]

    Liu H, Wang K, Peng D, Li H, Zhu Y , Zhang S, Liu M and Tian J 2016IEEE transactions on medical imaging 35 2546–2557

  26. [34]

    Arigovindan M, Fung J C, Elnatan D, Mennella V , Chan Y H M, Pollard M, Branlund E, Sedat J W and Agard D A 2013 Proceedings of the National Academy of Sciences 201315675

  27. [35]

    Cox B T, Kara S, Arridge S R and Beard P C 2007 The Journal of the Acoustical Society of America 121 3453–3464

  28. [36]

    Cox B and Beard P 2005 The Journal of the Acoustical Society of America 117 3616–3627

  29. [37]

    Xu M and Wang L V 2006 Review of scientific instruments 77 041101

  30. [38]

    Wang K, Ermilov S A, Su R, Brecht H P, Oraevsky A A and Anastasio M A 2011IEEE transactions on medical imaging 30 203–214

  31. [39]

    Wang K, Su R, Oraevsky A A and Anastasio M A 2012 Physics in Medicine & Biology 57 5399

  32. [40]

    Tao M, Yang J and He B 2009 TR0918, Department of Mathematics, Nanjing University

  33. [41]

    Wang Y , Yang J, Yin W and Zhang Y 2008SIAM Journal on Imaging Sciences 1 248–272

  34. [42]

    Skariah D G and Arigovindan M 2017 IEEE Transactions on Image Processing 26 4471–4482

  35. [43]

    Nikolova M, Ng M K and Tam C 2010 IEEE Transactions on Image Processing 19 3073–3088

  36. [44]

    Zhou Wang, Bovik A C, Sheikh H R and Simoncelli E P 2004 IEEE Transactions on Image Processing 13 600–612

  37. [45]

    Awasthi N, Kalva S K, Pramanik M and Yalavarthy P K 2018 Journal of Biomedical Optics 23 091413

  38. [46]

    Li L, Zhu L, Shen Y and Wang L V 2017 Journal of biomedical optics 22 076017 15

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.