About

Leonid Kunyansky is a faculty member in the Program in Applied Mathematics at the University of Arizona. His research interests include electromagnetic and acoustic scattering, wave propagation, photonic crystals, and computerized tomography. His work focuses on these areas, contributing to the understanding and development of mathematical models and methods related to wave phenomena and imaging technologies.

Research topics

• Mathematics
• Mathematical analysis
• Algorithm
• Computer science
• Acoustics

Selected publications

• Half-time range description for the free space wave operator and the spherical means transform

  Inverse Problems · 2025-01-31

  articleOpen accessSenior author

  Abstract The forward problem arising in several hybrid imaging modalities can be modeled by the Cauchy problem for the free space wave equation. Solution to this problems describes propagation of a pressure wave, generated by a source supported inside unit sphere S . The data g represent the time-dependent values of the pressure on the observation surface S . Finding initial pressure f from the known values of g consitutes the inverse problem. The latter is also frequently formulated in terms of the spherical means of f with centers on S . Here we consider a problem of range description of the wave operator mapping f into g . Such a problem was considered before, with data g known on time interval at least <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mrow> <mml:mo stretchy="false">[</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> <mml:mn>2</mml:mn> <mml:mo stretchy="false">]</mml:mo> </mml:mrow> </mml:math> (assuming the unit speed of sound). Range conditions were also found in terms of spherical means, with radii of integration spheres lying in the range <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mrow> <mml:mo stretchy="false">[</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> <mml:mn>2</mml:mn> <mml:mo stretchy="false">]</mml:mo> </mml:mrow> </mml:math> . However, such data are redundant. We present necessary and sufficient conditions for function g to be in the range of the wave operator, for g given on a half-time interval <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mrow> <mml:mo stretchy="false">[</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> <mml:mn>1</mml:mn> <mml:mo stretchy="false">]</mml:mo> </mml:mrow> </mml:math> . This also implies range conditions on spherical means measured for the radii in the range <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mrow> <mml:mo stretchy="false">[</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> <mml:mn>1</mml:mn> <mml:mo stretchy="false">]</mml:mo> </mml:mrow> </mml:math> .

  PublisherOA PDFDOI

• Multi-Electrode Acoustoelectric Imaging of Current Source Density via Scalar Potential Reconstruction: Feasibility Study in a Cube Phantom

  2025-10-14

  preprintOpen accessSenior author

  Acoustoelectric (AE) imaging has the potential to reconstruct current sources with high spatial and temporal resolution. However, a significant challenge is low signal-to-noise ratio (SNR), which complicates source detection and ultimately degrades spatial resolution. In this study, we propose a novel multi-electrode reconstruction method to image the current source density (CSD) in biological tissue at high image frame rates. This approach reconstructs the scalar potential of the current source with only single-step inversion regardless of the number of electrodes. It effectively removes modulation effects caused by acoustic and electrical-lead fields to obtain an unmodulated representation of the CSD. The approach was validated using a simple bio-mimicking cube phantom, with two stimulus electrodes to simulate current sources and six recording electrodes to capture AE signals. A 1.5D ultrasound transducer array operating at 0.6 MHz with focused beam scanning generated 2D AE signals for CSD reconstruction. Accurately modeled acoustic and lead fields were used in reconstruction through regularized inversion with truncated singular value decomposition (TSVD). The reconstructed scalar potential of the current source demonstrated improved spatial resolution with up to 33 percent decrease in the area of the image point spread function and a 7 dB increase in SNR compared to raw AE signals for the small set of electrodes used. Moreover, the results were in strong agreement with simulated ground truth.

  PublisherOA PDFDOI

• Digital twins enable full-reference quality assessment of photoacoustic image reconstructions

  The Journal of the Acoustical Society of America · 2025-07-01 · 1 citations

  articleOpen access

  Quantitative comparison of the quality of photoacoustic image reconstruction algorithms remains a major challenge. No-reference image quality measures are often inadequate, but full-reference measures require access to an ideal reference image. While the ground truth is known in simulations, it is unknown in vivo or in phantom studies, as the reference depends on both the phantom properties and the imaging system. This paper tackles this problem by using numerical digital twins of tissue-mimicking phantoms and the imaging system to perform a quantitative calibration to reduce the simulation gap. The contributions of this paper are twofold: First, this digital-twin framework is used to compare multiple state-of-the-art reconstruction algorithms. Second, among these is a Fourier transform-based reconstruction algorithm for circular detection geometries, which is tested on experimental data for the first time. The results demonstrate the usefulness of digital phantom twins by enabling assessment of the accuracy of the numerical forward model and enabling comparison of image reconstruction schemes with full-reference image quality assessment. This paper shows that the Fourier transform-based algorithm yields results comparable to those of iterative time reversal, but at a lower computational cost.

  PublisherDOI

• Digital twins enable full-reference quality assessment of photoacoustic image reconstructions

  ArXiv.org · 2025-05-30

  preprintOpen access

  Quantitative comparison of the quality of photoacoustic image reconstruction algorithms remains a major challenge. No-reference image quality measures are often inadequate, but full-reference measures require access to an ideal reference image. While the ground truth is known in simulations, it is unknown in vivo, or in phantom studies, as the reference depends on both the phantom properties and the imaging system. We tackle this problem by using numerical digital twins of tissue-mimicking phantoms and the imaging system to perform a quantitative calibration to reduce the simulation gap. The contributions of this paper are two-fold: First, we use this digital-twin framework to compare multiple state-of-the-art reconstruction algorithms. Second, among these is a Fourier transform-based reconstruction algorithm for circular detection geometries, which we test on experimental data for the first time. Our results demonstrate the usefulness of digital phantom twins by enabling assessment of the accuracy of the numerical forward model and enabling comparison of image reconstruction schemes with full-reference image quality assessment. We show that the Fourier transform-based algorithm yields results comparable to those of iterative time reversal, but at a lower computational cost. All data and code are publicly available on Zenodo: https://doi.org/10.5281/zenodo.15388429.

  PublisherOA PDFDOI

• Half-time Range description for the free space wave operator and the spherical means transform

  arXiv (Cornell University) · 2024-10-19

  preprintOpen accessSenior author

  The forward problem arising in several hybrid imaging modalities can be modeled by the Cauchy problem for the free space wave equation. Solution to this problems describes propagation of a pressure wave, generated by a source supported inside unit sphere $S$. The data $g$ represent the time-dependent values of the pressure on the observation surface $S$. Finding initial pressure $f$ from the known values of $g$ consitutes the inverse problem. The latter is also frequently formulated in terms of the spherical means of $f$ with centers on~$S$. Here we consider a problem of range description of the wave operator mapping $f$ into $g$. Such a problem was considered before, with data $g$ known on time interval at least $[0,2]$ (assuming the unit speed of sound). Range conditions were also found in terms of spherical means, with radii of integration spheres lying in the range $[0,2]$. However, such data are redundant. We present necessary and sufficient conditions for function $g$ to be in the range of the wave operator, for $g$ given on a half-time interval $[0,1]$. This also implies range conditions on spherical means measured for the radii in the range $[0,1]$.

  PublisherOA PDFDOI

• On the exactness of the universal backprojection formula for the spherical means Radon transform

  Inverse Problems · 2023-01-13 · 4 citations

  articleOpen accessSenior author

  Abstract The spherical means Radon transform is defined by the integral of a function f in <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:msup> <mml:mrow> <mml:mi mathvariant="double-struck">R</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>n</mml:mi> </mml:mrow> </mml:msup> </mml:math> over the sphere <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mi>S</mml:mi> <mml:mo stretchy="false">(</mml:mo> <mml:mi>x</mml:mi> <mml:mo>,</mml:mo> <mml:mi>r</mml:mi> <mml:mo stretchy="false">)</mml:mo> </mml:math> of radius r centered at a x , normalized by the area of the sphere. The problem of reconstructing f from the data where x belongs to a hypersurface <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mi mathvariant="normal">Γ</mml:mi> <mml:mo>⊂</mml:mo> <mml:msup> <mml:mrow> <mml:mi mathvariant="double-struck">R</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>n</mml:mi> </mml:mrow> </mml:msup> </mml:math> and <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mi>r</mml:mi> <mml:mo>∈</mml:mo> <mml:mo stretchy="false">(</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> <mml:mi mathvariant="normal">∞</mml:mi> <mml:mo stretchy="false">)</mml:mo> </mml:math> has important applications in modern imaging modalities, such as photo- and thermo- acoustic tomography. When Γ coincides with the boundary <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mi mathvariant="normal">∂</mml:mi> <mml:mi mathvariant="normal">Ω</mml:mi> </mml:math> of a bounded (convex) domain <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mi mathvariant="normal">Ω</mml:mi> <mml:mo>⊂</mml:mo> <mml:msup> <mml:mrow> <mml:mi mathvariant="double-struck">R</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>n</mml:mi> </mml:mrow> </mml:msup> </mml:math> , a function supported within Ω can be uniquely recovered from its spherical means known on Γ. We are interested in explicit inversion formulas for such a reconstruction. If <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll"> <mml:mi mathvariant="normal">Γ</mml:mi> <mml:mo>=</mml:mo> <mml:mi mathvariant="normal">∂</mml:mi> <mml:mi mathvariant="normal">Ω</mml:mi> </mml:math> , such formulas are only known for the case when Γ is an ellipsoid (or one of its partial cases). This gives rise to a question: can explicit inversion formulas be found for other closed hypersurfaces Γ? In this article we prove, for the so-called ‘universal backprojection inversion formulas’, that their extension to non-ellipsoidal domains Ω is impossible, and therefore ellipsoids constitute the largest class of closed convex hypersurfaces for which such formulas hold.

  PublisherDOI

• Weighted Radon transforms of vector fields, with applications to magnetoacoustoelectric tomography

  Inverse Problems · 2023-04-26 · 9 citations

  articleOpen access1st authorCorresponding

  Abstract Currently, theory of ray transforms of vector and tensor fields is well developed, but the Radon transforms of such fields have not been fully analyzed. We thus consider linearly weighted and unweighted longitudinal and transversal Radon transforms of vector fields. As usual, we use the standard Helmholtz decomposition of smooth and fast decreasing vector fields over the whole space. We show that such a decomposition produces potential and solenoidal components decreasing at infinity fast enough to guarantee the existence of the unweighted longitudinal and transversal Radon transforms of these components. It is known that reconstruction of an arbitrary vector field from only longitudinal or only transversal transforms is impossible. However, for the cases when both linearly weighted and unweighted transforms of either one of the types are known, we derive explicit inversion formulas for the full reconstruction of the field. Our interest in the inversion of such transforms stems from a certain inverse problem arising in magnetoacoustoelectric tomography (MAET). The connection between the weighted Radon transforms and MAET is exhibited in the paper. Finally, we demonstrate performance and noise sensitivity of the new inversion formulas in numerical simulations.

  PublisherOA PDFDOI

• Weighted Radon transforms of vector fields, with applications to magnetoacoustoelectric tomography

  arXiv (Cornell University) · 2023-05-10

  preprintOpen access1st authorCorresponding

  Currently, theory of ray transforms of vector and tensor fields is well developed, but the Radon transforms of such fields have not been fully analyzed. We thus consider linearly weighted and unweighted longitudinal and transversal Radon transforms of vector fields. As usual, we use the standard Helmholtz decomposition of smooth and fast decreasing vector fields over the whole space. We show that such a decomposition produces potential and solenoidal components decreasing at infinity fast enough to guarantee the existence of the unweighted longitudinal and transversal Radon transforms of these components. It is known that reconstruction of an arbitrary vector field from only longitudinal or only transversal transforms is impossible. However, for the cases when both linearly weighted and unweighted transforms of either one of the types are known, we derive explicit inversion formulas for the full reconstruction of the field. Our interest in the inversion of such transforms stems from a certain inverse problem arising in magnetoacoustoelectric tomography (MAET). The connection between the weighted Radon transforms and MAET is exhibited in the paper. Finally, we demonstrate performance and noise sensitivity of the new inversion formulas in numerical simulations.

  PublisherDOI

• On the exactness of the universal backprojection formula for the spherical means Radon transform

  arXiv (Cornell University) · 2022-07-17

  preprintOpen accessSenior author

  The spherical means Radon transform $\mathcal{M}f(x,r)$ is defined by the integral of a function $f$ in $\mathbb{R}^{n}$ over the sphere $S(x,r)$ of radius $r$ centered at a $x$, normalized by the area of the sphere. The problem of reconstructing $f$ from the data $\mathcal{M}f(x,r)$ where $x$ belongs to a hypersurface $Γ\subset\mathbb{R}^{n}$ and $r \in(0,\infty)$ has important applications in modern imaging modalities, such as photo- and thermo- acoustic tomography. When $Γ$ coincides with the boundary $\partialΩ$ of a bounded (convex) domain $Ω\subset\mathbb{R}^{n}$, a function supported within $Ω$ can be uniquely recovered from its spherical means known on $Γ$. We are interested in explicit inversion formulas for such a reconstruction. If $Γ=\partialΩ$, such formulas are only known for the case when $Γ$ is an ellipsoid (or one of its partial cases). This gives rise to the natural question: can explicit inversion formulas be found for other closed hypersurfaces $Γ$? In this article we prove, for the so-called "universal backprojection inversion formulas", that their extension to non-ellipsoidal domains $Ω$ is impossible, and therefore ellipsoids constitute the largest class of closed convex hypersurfaces for which such formulas hold.

  PublisherOA PDFDOI

• Current Source Density Imaging Using Regularized Inversion of Acoustoelectric Signals

  IEEE Transactions on Medical Imaging · 2022-11-03 · 3 citations

  articleOpen access

  Acoustoelectric (AE) imaging can potentially image biological currents at high spatial (~mm) and temporal (~ms) resolution. However, it does not directly map the current field distribution due to signal modulation by the acoustic field and electric lead fields. Here we present a new method for current source density (CSD) imaging. The fundamental AE equation is inverted using truncated singular value decomposition (TSVD) combined with Tikhonov regularization, where the optimal regularization parameter is found based on a modified L-curve criterion with TSVD. After deconvolution of acoustic fields, the current field can be directly reconstructed from lead field projections and the CSD image computed from the divergence of that field. A cube phantom model with a single dipole source was used for both simulation and bench-top phantom studies, where 2D AE signals generated by a 0.6 MHz 1.5D array transducer were recorded by orthogonal leads in a 3D Cartesian coordinate system. In simulations, the CSD reconstruction had significantly improved image quality and current source localization compared to AE images, and performance further improved as the fractional bandwidth (BW) increased. Similar results were obtained in the phantom with a time-varying current injected. Finally, a feasibility study using an in vivo swine heart model showed that optimally reconstructed CSD images better localized the current source than AE images over the cardiac cycle.

  PublisherOA PDFDOI

Recent grants

• Collaborative Research: Mathematical Methods for Emerging Techniques of Biomedical Imaging

  NSF · $109k · 2009–2012

• Inverse Problems Arising in Novel Modalities of Biomedical Imaging

  NSF · $300k · 2018–2022

• Development of Efficient Algorithms for Computational Electromagnetism and Computerized Tomography

  NSF · $32k · 2003–2005

• Collaborative research: Mathematics of emerging imaging methods in medicine and homeland security

  NSF · $183k · 2012–2017

Frequent coauthors

• Peter Kuchment

  22 shared

• Mark Agranovsky

  6 shared

• Benjamin Robert Holman

  University of Arizona

  3 shared

• Philip Hoskins

  3 shared

• Linh V. Nguyen

  3 shared

• Gaik Ambartsoumian

  3 shared

• Oscar P. Bruno

  California Institute of Technology

  3 shared

• Russell S. Witte

  University of Arizona

  2 shared

Labs

• Program in Applied MathematicsPI