nature biomedical engineering
https://doi.org/10.1038/s41551-023-01149-4
Artic�e
Ultrafast longitudinal imaging of
haemodynamics via single-shot volumetric
photoacoustic tomography with a single-
element detector
In the format provided by the
authors and unedited
1
Contents
Supplementary Note 1 | Fly’s eye homogenizer in PACTER.
Supplementary Note 2 | Universal calibration of PACTER.
Supplementary Note 3 | Forward model and image reconstruction.
Supplementary Note 4 | Spatial resolution analysis based on the virtual transducer array.
Supplementary Note 5 | Spatial resolution analysis based on the system matrix.
Supplementary Fig. 1 | Comparison of uniform optical absorbers for calibration.
Supplementary Fig. 2 | Fly’s eye homogenizer in PACTER.
Supplementary Fig.
3 | PACTER with different illumination areas.
Supplementary Fig.
4 | Positioning of the single
-
element ultrasonic transducer on the ER.
Supplementary Fig.
5 | Fabrication of the single
-
element ultrasonic transducer.
Supplementary Fig.
6 | Characterization of the single
-
element ultrasonic transducer.
Supplementary Fig.
7 | Comparison of the integrated and post
-
fabrication transducers.
Supplementary Fig. 8 | Object
-
dependent and
-
independent calibrations in PATER and PACTER,
respectively.
Supplementary Fig.
9 | PACTER reconstruction using the calibration data acquired at different times.
Supplementary Fig. 10 | PACTER signals from calibration using bovine blood.
Supplementary Fig. 11 | Validation of PACTER reconstruction.
Supplementary Fig. 12 | Spatial resolution analysis from the virtual
-
transducer
-
array perspective.
Supplementary Fig. 13 | Spatial resolution analysis from the system
-
matrix perspective.
Supplementary Fig. 14 | Spectrogram of the PACTER signal.
Supplementary Fig.
15 | 3D PACTER images of an object placed at different depths.
Supplementary Fig. 16 | PACTER of bovine blood flushing through a tube.
Supplementary Fig.
17 | Schematic of a PACTER system allowing handheld operation.
Supplementary Fig. 18 | Temperature stabilization in PACTER.
Supplementary Fig.
19 | Simulation of PACTER signals from ERs with different dimensions.
Supplementary Fig.
20 | Potential ways to
improve the spatial resolution of PACTER.
Supplementary Fig.
21 | Potential application of PACTER in ultrasonography.
Supplementary Fig.
22 | Effect of the iteration number in the FISTA algorithm for PACTER reconstruction.
Supplementary Fig.
23 | Effect of the image processing filters on PACTER images.
Supplementary Table 1 | Comparison of PACTER and previous techniques using ERs.
Supplementary Video 1 | Principle and implementation of PACTER (with narration).
Supplementary Video 2 | Simulation of PATER and PACTER signals.
Supplementary Video 3 | 4D PACTER image and its maximum amplitude projection of bovine blood flushing
through an S
-
shaped tube.
Supplementary Video 4 | 4D PACTER images of bovine blood flowing through a tube at different speeds.
Supplementary Video 5 | 4D PACTER image and its maximum amplitude projection of bovine blood flowing
through a tube with a speed of 272.5 mm/s, captured at 1,000 volumes per
second.
Supplementary Video 6 | 4D
in vivo
PACTER image of the abdominal vasculature of mouse 1.
Supplementary Video 7 | 4D
in vivo
PACTER image of the abdominal vasculature of mouse 2.
Supplementary Video 8 | 4D
in vivo
PACTER image of the abdominal vasculature of mouse 3.
Supplementary Video 9 | 4D
in vivo
PACTER image of the thenar vasculature of participant 1.
Supplementary Video 10 | 4D
in vivo
PACTER image of the thenar vasculature of participant 2.
Supplementary Video 11 | 4D
in vivo
PACTER image of the thenar vasculature of participant 1 in a different
region.
Supplementary Video 12 | 4D
in vivo
PACTER image of the foot vessels of participant 3.
Supplementary Video 13 | 4D
in vivo
PACTER image of the foot vessels of participant 3 in a different region.
Supplementary Video 14 | PACTER signals affected by temperature fluctuations.
Supplementary r
eferences
2
Supplementary Note
1
|
F
ly’s eye homogenizer
in
PACTER
.
We use a
fly’s eye homogenizer
to provide uniform illumination for PACTER
74
.
As shown in
Supplementary
Fig.
2
, the two microlens arrays form
multiple parallel Köhler illumination systems side
-
by
-
side. The original
beam entering the first microlens array is divided into multiple beamlets
.
Through the lenslets pairs in the
microlens arrays and the spherical lens, each beamlet of the original beam is imaged to the homogenization
plane, i.e., the top of the
ER
.
Because the
images of the beamlets are
all
superimposed on the
homogenization
plane, the intensity differences among the beamlets disappear in their superimposed
images. Therefore, the intens
ity distribution of the homogenized beam is independent of the homogeneity of
the original beam. Further,
the
square
-
type microlens arrays in
the
setup
will generate a square flat
-
top
intensity distribution in the homogenization plane, which provides a good match for the square
-
type
calibration pattern in PACTER, allowing accurate mapping between the calibration and imaging areas.
The
width
of the
homogenized beam,
푑
!
, is given by
푑
!
=
푝
"
푓
#
(
2
푓
"
−
푎
)
푓
"
$
,
(
S
1
)
where
푝
"
and
푓
"
are
the pitch
and focal length
of the
lenslets in the two identical
microlens arrays
,
푓
#
is the
focal length of the spherical lens, and
푎
is the separation between the microlens arrays. In PACTER,
푎
is set
to be identical to
푓
"
, leading to
푑
!
=
푝
"
푓
#
푓
"
.
(
S
1
)
The divergence
half
-
angle
after the homogenization plane,
휃
,
is given by
휃
=
atan
3
푑
%
+
푑
!
−
푝
"
2
푓
#
5
,
(
S
2
)
where
푑
%
is the diameter of the original beam. In PACTER, we use
the
microlens arrays with
푝
"
=
0
.
5
mm
and
푓
"
=
15
mm
, the spherical lens with
푓
#
=
250
mm
, and the original beam with a diameter of
푑
%
=
6
mm
.
Therefore, the homogenized beam has a
width
of
푑
!
≈
8
mm
, matching the size of the calibration pattern (80
by 80 steps with a step size of 0.1 mm), and the divergence half
-
angle
휃
≈
2°
. The small divergence
ensures
homogenous illumination across the whole 3D volume (
8
mm
×
8
mm
×
3
.
6
mm
) for
in vivo
imaging. Within
the 3.6 mm depth, the illumination beam merely diverges laterally by 0.
13
mm, which is much smaller than
the lateral resolution (0.5
6
mm) of PACTER. Hence, the beam divergence within the imaging volume can be
ignored.
3
Supplementary Note
2
|
Universal calibration of PACTER
.
Recently, we developed photoacoustic topography through an ergodic relay (PATER), which captured a
wide
-
field PA image with a single laser shot using a single
-
element ultrasonic transducer
38
–
40
. However,
PATER suffered from the following limitations. (1) PATER could only image a 2D projection of the object,
unable to capture
a
tomographic image in 3D. (2) The ER in PATER had to be re
-
calibrated for each different
object; because the calibration relied on
a
point
-
by
-
point scanning of a laser beam, the imagin
g procedure
was time
-
consuming. (3) PATER required that the boundary conditions
of
the object and the ER stay
unchanged throughout the experiment and thus could not be used for long
-
term imaging in unstable
environment
s
.
One major distinction between
PATER and PACTER
lies in the structures of their ERs. In PATER, the ER is
simply a prism (
Supplementary Fig.
8
a
). After a PA signal is generated from the object (
푡
&
) and detected by
the transducer (
푡
'
), it quickly propagates to and reflects from the boundary between the object and the ER
(
푡
$
), as shown in the simulation based on the
k
-
wave toolbox
75
(
Supplementary Fig.
8
b
). Once the
boundary condition between the object and the ER changes because of movements of the object, instability
of the system, or switching to a different object, the signal will also change. Therefore, if the ER is calibrated
with one object, e.g., b
ovine blood, and then used to image a different object, e.g., a black wire, the
reconstruction will fail (
Supplementary Fig.
8
c
). In other words, the PATER signal is strongly
object
-
dependent (
Supplementary Fig.
8
d
), and it requires that the boundary conditions
of
the object and the ER
remain unchanged throughout the experiment.
In comparison, the ER of PACTER consists of a prism and a fused silica rod, where the rod functions as an
acoustic delay line that temporally separates the initial and reflected PA signals
41
(
Supplementary Fig.
8
e
).
When a PA signal is generated from the object (
푡
&
) and detected by the ultrasonic transducer (
푡
'
), a part of
the acoustic signal is trapped and scrambled in the prism,
w
hereas the other part is reflected toward the
object and then reverberated to the transducer (
푡
$
), as denoted by the black dotted arrows in the simulation
(
Supplementary Fig.
8
f
). Because only the reflected part of the acoustic signal will be affected by the
boundary condition between the object and the ER, once this part is excluded from the measurement, the
acquired signal will be object
-
independent. Consequently, the ER can be
used to image any object, e.g., a
black wire, despite it has been calibrated only once using a different object, e.g., bovine blood
(
Supplementary Fig.
8
g
). Whereas almost th
e whole PATER signal is object
-
dependent (
Supplementary
Fig.
8
d
), a large segment (>100 μs) of the PACTER signal is object
-
independent (
Supplementary Fig.
8
f
),
enabling universal calibration in PACTER
despite its 3D imaging capability
(
Supplementary Video
2
).
It
should be noted that the long duration of the signals poses a limit on the maximum laser repetition rate, i.e.,
the volumetric imaging rate, the system can
support. For instance, if the duration of the signal is longer than
250 μs, the laser repetition rate should be lower than 4 kHz.
4
Supplementary Note
3
|
Forward
model and image
reconstruction
.
W
e performed
calibration
s at
pixels
on
a 2D plane
and used these
pixels
as
virtual
ultrasonic
transducers for
3D
imaging
.
If
non
-
zero initial pressure exists only on the calibration plane, the detected signal
푠
(
푡
)
at time
푡
can be expressed as
푠
(
푡
)
=
A
푝
(
푘
(
(
푡
)
)
(
*
'
,
푡
≥
0
,
(
S
3
)
w
here
푁
is the number of
calibrated virtual transducers
,
푘
(
(
푡
)
is the normalized impulse response from the
calibration at the
푛
-
th
virtual transducer
, and
푝
(
is the
root
-
mean
-
squared
PA amplitude
proportional
to the
initial pressure
at the
푛
-
th
virtual transducer
.
For
initial pressure
in a 3D volume,
we a
ssum
e
푀
source points
located at
퐫
+
,
,
푚
=
1
,
2
,
...
,
푀
, in a
n
acoustically
homogeneous
3D region
attached
to
the calibration plane.
The PA wave generated
from
the
source point
at
퐫
+
,
propagates to the
calibrated virtual transducer
퐫
(
with the speed of sound
푐
after time
푡
+
,
(
=
.
퐫
!
"
0
퐫
#
.
1
, which, through the
ER
, adds
푝
+
,
(
푘
(
K
푡
−
.
퐫
!
"
0
퐫
#
.
1
L
to the detected signal, with the PA amplitude
푝
+
,
(
quantified as
푝
+
,
(
=
2
3
4
!
,
#
5
6
%
,
!
.
퐫
!
"
0
퐫
#
.
.
Here,
휃
+
,
(
denotes the incidence angle
satisfying
cos
휃
+
,
(
=
3
퐫
!
"
0
퐫
#
5
⋅
퐧
.
퐫
!
"
0
퐫
#
.
with
퐧
being
the normal vector of the calibration plane
;
function
푤
R
휃
+
,
(
S
describes a virtual transducer’s angle
-
dependent sensitivity
;
and
푝
&
,
+
is proportional to the initial pressure at
퐫
+
,
.
We replace
푝
(
푘
(
(
푡
)
in
Eq.
(
S3
)
with
푝
+
,
(
푘
(
K
푡
−
.
퐫
!
"
0
퐫
#
.
1
L
from all the
푀
source points
and obtain
t
he detected wide
-
field PA signal
푠
(
푡
)
=
A
A
푤
R
휃
+
,
(
S
푝
&
,
+
‖
퐫
+
,
−
퐫
(
‖
푘
(
U
푡
−
‖
퐫
+
,
−
퐫
(
‖
푐
V
"
+
*
'
)
(
*
'
,
푡
≥
0
.
(
S
4
)
Here, we define
푘
(
(
푡
)
=
0
,
푛
=
1
,
2
...
,
푁
,
푡
<
0
.
For sufficiently small virtual ultrasonic transducers
, we
assume that
푤
R
휃
+
,
(
S
=
ퟏ
[
&
,
4
&
]
R
휃
+
,
(
S
cos
휃
+
,
(
.
(
S
5
)
Here
, we
use the indicator function
ퟏ
;
(
푥
)
=
[
1
,
푥
∈
퐴
0
,
푥
∉
퐴
(
S
6
)
to
rejection detections with
incidence angles greater than the critical angle
휃
<
(
quantified in
Supplementary
Note
4
)
.
Substituting Eq.
(
S5
)
into Eq. (
S4
) yields
푠
(
푡
)
=
A
A
푝
&
,
+
ퟏ
[
&
,
4
&
]
R
휃
+
,
(
S
cos
휃
+
,
(
‖
퐫
+
,
−
퐫
(
‖
푘
(
U
푡
−
‖
퐫
+
,
−
퐫
(
‖
푐
V
"
+
*
'
)
(
*
'
,
푡
≥
0
.
(
S
7
)
We let
퐿
be
the number of time points after temporal discretization. Then the
computational complexity of a
forward model based on Eq. (
S7
) is
푂
(
푀푁퐿
)
.
To accelerate the forward model in Eq. (
S7
), we split the delay term
.
퐫
!
"
0
퐫
#
.
1
from function
푘
(
(
푡
)
through
temporal convolution:
푘
(
U
푡
−
‖
퐫
+
,
−
퐫
(
‖
푐
V
=
훿
U
푡
−
‖
퐫
+
,
−
퐫
(
‖
푐
V
∗
푘
(
(
푡
)
,
푚
=
1
,
2
,
...
,
푀
,
푛
=
1
,
2
,
...
,
푁
,
푡
≥
0
.
(
S
8
)
Substituting Eq. (
S8
) into Eq. (
S7
), we obtain
푠
(
푡
)
=
A
푘
(
(
푡
)
∗
A
푝
&
,
+
ퟏ
[
&
,
4
&
]
R
휃
+
,
(
S
cos
휃
+
,
(
‖
퐫
+
,
−
퐫
(
‖
훿
U
푡
−
‖
퐫
+
,
−
퐫
(
‖
푐
V
"
+
*
'
)
(
*
'
,
푡
≥
0
.
(
S
9
)
The inner summation in Eq. (
S9
) has a complexity of
푂
(
푀푁
)
and each
temporal convolution is implemented
through three fast Fourier transforms (FFTs) with a complexity of
푂
(
퐿
log
$
퐿
)
. Thus, the forward model based
on Eq. (
S9
) has a computational complexity of
max
{
푂
(
푀푁
)
,
푂
(
푁퐿
log
$
퐿
)
}
.
We perform numerical simulations to quantify the improvement of computational efficiency brought by the
fast algorithm. Considering that
the
complexities of both the slow
(
Eq. (
S7
)
)
and fast
(
Eq. (
S9
)
)
algorithms are
linearly dependent on the number of virtual detectors
푁
, we simplify the problem to a single virtual detector
(
푁
=
1) with
퐿
= 65,536 and
푀
=
80
×
80
×
120
, and only consider the computation time of a forward
simulation. In a Windows 11 Home system with Intel® Core
™
i9
-
10900T CPU @ 1.90GHz, we perform
single
-
CPU
-
core forward simulations based on the slow and fast algorithms, respectively, 36 times in Matlab.
The average computation times of the simplified forward simulation are 20 s and
2
.
2
×
10
0
=
s (
×
9
,
100
acceleration), respectively, which correspond to 35.6 h and 14.1 s for a true forward simulation in this study
(
푁
=
80
×
80
).
5
Supplementary Note
4
|
Spatial resolution analysis based on the virtual transducer array
.
We analyze the anisotropy of the spatial resolution from two perspectives
: one with an
emphasis
on the
physical
intuition
(this Note)
, the other with an emphasi
s
on the mathematical derivation
(
Supplementary
Note
5
)
.
The two
perspective
s
clarif
y
the
differences
in value
between
the
lateral and axial
resolution
s
.
We first
analyze
the anisotropy of the spatial
resolution by
assuming that
the
PACTER signals detected by
the single
-
element transducer
are accuratel
y decoded
to the
signals
detected by
the
80
×
80
virtual
transducer
s
.
We use
the
spectr
um of the PACTER signal
detected by the integrated PMN
-
PT transducer
(30
-
MHz central frequency)
in
Supplementary Fig.
12
a
, denoted as
푠
̂
(
푓
)
,
to approximate the
detection
spectrum
of each virtual transducer element.
From
푠
̂
(
푓
)
, we obtain the central frequency
푓
&
=
∫
?
|
A
̂
(
?
)
|
'
E
?
(
)
%
∫
|
A
̂
(
?
)
|
'
E
?
(
)
%
≈
8
.
3
MHz
and
the
bandwidth
푓
F
=
n
∫
(
?
0
?
%
)
'
|
A
̂
(
?
)
|
'
E
?
(
)
%
∫
|
A
̂
(
?
)
|
'
E
?
(
)
%
≈
12
.
2
MHz
, respectively
76
. Here, we
set
t
he
upper
frequency
as
푓
G
=
125
MHz
, which is the Nyquist frequency corresponding to the 250 MHz sampling rate.
Additionally
, the speed
of sound
in
fused
silica
(
longitudinal wave speed
푐
<
=
5
.
9
mm
⋅
μ
s
0
'
)
77
is greater than
th
at in
water (
푐
=
1
.
5
mm
⋅
μ
s
0
'
)
78
,
leading to
a
critical angle
(denoted at
휃
<
)
of the ultrasonic wave’s
refraction from water into
fused silica
(
Supplementary Fig.
12
b
)
.
This leads to a limited
-
view effect in the
proposed system: waves of incidence angels greater than the critical angle
휃
<
are totally reflected and not
detected by the transducer.
The critical angle
satisfies
sin
휃
<
=
1
1
&
≈
0
.
25
, which equals the numerical
aperture
(NA)
of
the
virtual
2D array.
Here
,
the longitudinal
wave speed
is used to quantify the critical angle
due to the single
-
element transducer’s dominant sensitivity to the longitudinal waves
79
.
For
the virtual 2D
array, the
axial resolution
(along the
z
-
axis)
is quantified as
푅
;
=
0
.
88
1
?
*
≈
0
.
11
mm
,
while the lateral
resolution (in the
x
-
y
plane) is quantified as
푅
#
=
0
.
71
1
);
⋅
?
%
≈
0
.
51
mm
13
.
Both
values are close to the
measured axial (0.13 mm) and lateral (0.56 mm) resolutions shown in
Fig.
3
g
.
It is important to note
that, if a
post
-
fabrication
transducer
(10
-
MHz central frequency) is used, the central frequency
푓
&
and bandwidth
푓
2
of
the PACTER signal
will be 5.0 MHz and 4.0 MHz,
respectively (
Supplementary Fig.
7
). Consequently, the
lateral and axial resolutions of PACTER will be 0.86 mm and 0.33 mm, respectively.
Therefore, the
integrated PMN
-
PT transducer with a 30
-
MHz central frequency is more suitable for PACTER.
The
design of the integrated PMN
-
PT transducer incorporated an active aperture size of 0.4
×
0.4 mm
2
. This
decision was based on several factors. On one hand, it was essential for the transducer size to be relatively
small to ensure that the acoustic waves could be effectively scrambled in the ER before being detected by
the transducer. Conversely, if th
e transducer size is excessively small, only a minimal portion of the acoustic
waves would be captured, potentially leading to a deterioration in SNR and a corr
esponding decrease in
resolution
41
. Increasing the element size may compromise the scrambling but enhance the SNR, while
decreasing the size has the opposite effects, both of which influence image resolution. In this study, a
comprehensive analysis of the dependency of image resolution on
element size was not conducted. Instead,
through experimentation, we determined the element size of
0
.
4
×
0
.
4
m
m
$
to be a suitable balance between
sufficient scrambling and SNR, thereby ensuring high image quality.
6
Supplementary Note
5
|
Spatial resolution analysis based on the system matrix
.
Obtaining
signals detected by all the virtual transducers is an ill
-
posed problem
and we reconstruct an image
directly from the single
-
element
-
detected signals.
We iteratively invert the system matrix
퐇
(of size
퐿
×
푀
)
to
obtain the images
, which
manifest
anisotropic resolutions
(see
Methods
)
. To explain the
anisotropy
in theory,
we
analyze the
gradient
matrix
퐇
H
퐇
, whose property
eventually determines the image resolution in the
iterative method.
The
푚
-
th column of
퐇
is a temporally discretized form of the
following
transducer’s
response to the point source at
퐫
+
,
:
푠
+
(
푡
)
=
A
ퟏ
[
&
,
4
&
]
R
휃
+
,
(
S
cos
휃
+
,
(
‖
퐫
+
,
−
퐫
(
‖
푘
(
U
푡
−
‖
퐫
+
,
−
퐫
(
‖
푐
V
)
(
*
'
,
푚
=
1
,
2
,
...
,
푀
,
푡
≥
0
.
Accordingly, an element at the
푚
'
-
th
row and
푚
$
-
th
column in
퐇
H
퐇
is
expressed in temporally continuous
form as
t
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
=
t
v
A
ퟏ
[
&
,
4
&
]
R
휃
+
+
,
(
S
cos
휃
+
+
,
(
w
퐫
+
+
,
−
퐫
(
w
푘
(
U
푡
−
w
퐫
+
+
,
−
퐫
(
w
푐
V
)
(
*
'
x
v
A
ퟏ
[
&
,
4
&
]
R
휃
+
'
,
(
S
cos
휃
+
'
,
(
w
퐫
+
'
,
−
퐫
(
w
푘
(
U
푡
)
(
*
'
I
J
&
−
w
퐫
+
'
,
−
퐫
(
w
푐
V
x
d
푡
=
A
ퟏ
[
&
,
4
&
]
R
휃
+
+
,
(
+
S
cos
휃
+
+
,
(
+
w
퐫
+
+
,
−
퐫
(
+
w
ퟏ
[
&
,
4
&
]
R
휃
+
'
,
(
'
S
cos
휃
+
'
,
(
'
w
퐫
+
'
,
−
퐫
(
'
w
(
+
,
(
'
∈
{
'
,
$
,
...
,
)
}
×
t
푘
(
+
(
푡
)
푘
(
'
U
푡
+
w
퐫
+
+
,
−
퐫
(
+
w
푐
−
w
퐫
+
'
,
−
퐫
(
'
w
푐
V
d
푡
I
J
GOP
Q
&
,
.
퐫
!
'
"
0
퐫
#
'
.
1
0
.
퐫
!
+
"
0
퐫
#
+
.
1
R
=
A
훼
R
퐫
+
+
,
,
퐫
+
'
,
,
퐫
(
+
,
퐫
(
'
S
R
푘
(
+
⋆
푘
(
'
S
U
w
퐫
+
+
,
−
퐫
(
+
w
푐
−
w
퐫
+
'
,
−
퐫
(
'
w
푐
V
(
+
,
(
'
∈
{
'
,
$
,
...
,
)
}
,
푚
'
,
푚
$
∈
{
1
,
2
,
...
,
푀
}
.
Here, we define a weighting factor
훼
R
퐫
+
+
,
,
퐫
+
'
,
,
퐫
(
+
,
퐫
(
'
S
=
ퟏ
[
&
,
4
&
]
R
휃
+
+
,
(
+
S
cos
휃
+
+
,
(
+
w
퐫
+
+
,
−
퐫
(
+
w
ퟏ
[
&
,
4
&
]
R
휃
+
'
,
(
'
S
cos
휃
+
'
,
(
'
w
퐫
+
'
,
−
퐫
(
'
w
,
푚
'
,
푚
$
∈
{
1
,
2
,
...
,
푀
}
,
푛
'
,
푛
$
∈
{
1
,
2
,
...
,
푁
}
and a cross
-
correlation function
R
푘
(
+
⋆
푘
(
'
S
(
휏
)
=
t
푘
(
+
(
푡
)
푘
(
'
(
푡
+
휏
)
d
푡
I
J
GOP
{
&
,
0
S
}
,
휏
∈
ℝ
,
푛
'
,
푛
$
∈
{
1
,
2
,
...
,
푁
}
.
For any non
-
zero response
푠
+
(
푡
)
, we have
}
∫
푠
+
(
푡
)
푠
+
(
푡
)
d
푡
I
J
&
}
>
0
. A practical imaging system requires
}
∫
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
}
to be small for any
푚
'
≠
푚
$
. The image resolution is determined by the
decay
speed of
}
∫
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
}
as
(
푚
'
,
푚
$
)
moves away from the diagonal of
퐇
H
퐇
.
To explain the
anisotropy
of the resolution at
퐫
+
+
,
, we fix
푚
'
and vary
푚
$
so that
퐫
+
'
,
moves from
퐫
+
+
,
to other
locations around
퐫
+
+
,
. M
eanwhile, we estimate
the decay of
}
∫
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
}
.
We confine the movement
of
퐫
+
'
,
in a small region around
퐫
+
+
,
such that
the change of
훼
R
퐫
+
+
,
,
퐫
+
'
,
,
퐫
(
+
,
퐫
(
'
S
is negligible. Thus, we only
need to elaborate on the decay of
R
푘
(
+
⋆
푘
(
'
S
3
.
퐫
!
+
"
0
퐫
#
+
.
1
−
.
퐫
!
'
"
0
퐫
#
'
.
1
5
.
We first discuss two special cases: the maximum
-
cross
-
correlation
matrix
MCC
(
+
,
(
'
=
max
S
∈
ℝ
}
R
푘
(
+
⋆
푘
(
'
S
(
휏
)
}
,
푛
'
,
푛
$
∈
{
1
,
2
,
...
,
푁
}
and the
autocorrelation function
(
푘
(
⋆
푘
(
)
(
휏
)
,
푛
=
1
,
2
,
...
,
푁
. Both functions are determined by the bandwidth
of the transducer and the configuration of the
ER
.
Here
, we performed calibrations at
푁
=
80
×
80
pixels
with
a spacing of 0.1 mm in each dimension. We pick two
locations
퐫
(
U
+
and
퐫
(
U
'
and calculate the values of
MCC
(
U
+
,
(
,
MCC
(
U
'
,
(
,
R
푘
(
U
+
⋆
푘
(
U
+
S
(
휏
)
, and
R
푘
(
U
'
⋆
푘
(
U
'
S
(
휏
)
. The values of
MCC
(
U
+
,
(
and
MCC
(
U
'
,
(
for different
푛
form
two images shown in
Supplementary Fig.
13
a
and
b
, respectively. In each image, we draw two
perpendicular lines
퐿
'
and
퐿
$
centered at the corresponding calibration location and compare the pixel
values along the lines with the autocorrelation function values (
AC
) expressed as
R
푘
(
U
+
⋆
푘
(
U
+
S
K
V
1
L
or
R
푘
(
U
'
⋆
푘
(
U
'
S
K
V
1
L
, as shown in
Supplementary Fig.
13
c
and
d
, respectively
.
As a generalization of the full
width at half maximum (FWHM), we define
FW
W
,
X
as the full width at
훽
times the maximum of a line profile
퐹
,
7
and
FW
W
,
+
'
is equivalent to the
FWHM of
퐹
. For the first selected calibration location
퐫
(
U
+
, the decay of
R
푘
(
U
+
⋆
푘
(
U
+
S
K
V
1
L
vs.
푥
and the
decay of
MCC
(
U
+
,
(
vs.
}
퐫
(
U
+
−
퐫
(
}
are shown in the plots
AC
and
퐿
'
-
퐿
$
,
respectively, in
Supplementary Fig.
13
c
. We calculated
FW
;Y
,
X
,
FW
#
+
,
X
, and
FW
#
'
,
X
for
훽
= 0.20 and 0.40.
The values are shown in
Supplementary Fig.
13
c
. We repeat the quantification for the second selected
location
퐫
(
U
'
, as shown in
Supplementary Fig.
13
d
.
W
e observe that the increase of either
.
퐫
!
+
"
0
퐫
#
+
.
1
−
.
퐫
!
'
"
0
퐫
#
'
.
1
or
}
퐫
(
+
−
퐫
(
'
}
causes
R
푘
(
+
⋆
푘
(
'
S
3
.
퐫
!
+
"
0
퐫
#
+
.
1
−
.
퐫
!
'
"
0
퐫
#
'
.
1
5
to decay. Equivalently, avoiding the
decay
of
R
푘
(
+
⋆
푘
(
'
S
3
.
퐫
!
+
"
0
퐫
#
+
.
1
−
.
퐫
!
'
"
0
퐫
#
'
.
1
5
requires both
.
퐫
!
+
"
0
퐫
#
+
.
1
−
.
퐫
!
'
"
0
퐫
#
'
.
1
and
}
퐫
(
+
−
퐫
(
'
}
to be small.
Next, we analyze how the probability of both
.
퐫
!
+
"
0
퐫
#
+
.
1
−
.
퐫
!
'
"
0
퐫
#
'
.
1
and
}
퐫
(
+
−
퐫
(
'
}
being small decreases as
퐫
+
'
,
moves away from
퐫
+
+
,
.
We only discuss two representative cases:
(1)
퐫
+
'
,
moves on the plane crossing
퐫
+
+
,
and parallel to the calibration plane (
Supplementary Fig.
13
e
and
f
);
(2)
퐫
+
'
,
moves on the line crossing
퐫
+
+
,
and normal to the calibration plane (
Supplementary Fig.
13
g
and
h
). In the first case, we let
.
퐫
!
+
"
0
퐫
#
+
.
1
=
.
퐫
!
'
"
0
퐫
#
'
.
1
and observe how the
probability of
}
퐫
(
+
−
퐫
(
'
}
being small is violated. Because
.
퐫
!
+
"
0
퐫
#
+
.
1
=
.
퐫
!
'
"
0
퐫
#
'
.
1
and locations
퐫
+
+
,
and
퐫
+
'
,
are with the same distance to the calibration plane, locations
퐫
(
+
and
퐫
(
'
must be
on
two circles with the same radius and their centers separated by
}
퐫
+
+
,
−
퐫
+
'
,
}
. If
}
퐫
+
+
,
−
퐫
+
'
,
}
is small, for any
퐫
(
+
on one circle, there exists a
퐫
(
'
on the other circle such that
}
퐫
(
+
−
퐫
(
'
}
is small, as shown in
Supplementary Fig.
13
e
.
If
}
퐫
+
+
,
−
퐫
+
'
,
}
is large, a small value of
}
퐫
(
+
−
퐫
(
'
}
requires
both
퐫
(
+
and
퐫
(
'
to be
close to the same intersection of the two circles, as depicted in
Supplementary Fig.
13
f
.
Thus, the higher
the value of
}
퐫
+
+
,
−
퐫
+
'
,
}
, the lower the probability of
}
퐫
(
+
−
퐫
(
'
}
being small. In summary, for the first case, the
decay speed of
}
∫
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
}
as
}
퐫
+
+
,
−
퐫
+
'
,
}
increases is mainly determined by
MCC
(
+
,
(
'
’s decay
speed as
}
퐫
(
+
−
퐫
(
'
}
increases. In the second case, we let
푛
'
=
푛
$
=
푛
and observe the change of
휏
=
.
퐫
!
+
"
0
퐫
#
.
1
−
.
퐫
!
'
"
0
퐫
#
.
1
. As
퐫
+
'
,
moves away from
퐫
+
+
,
(from
Supplementary Fig.
13
g
to
h
),
푐휏
=
}
퐫
+
+
,
−
퐫
+
'
,
}
cos
휃
increases. Here
,
휃
is the angle between vectors
퐫
(
−
퐫
+
+
,
and
퐫
+
'
,
−
퐫
+
+
,
. If
퐫
(
,
퐫
+
+
,
, and
퐫
+
'
,
are on the same
line, the
calibrated virtual
transducer
at
퐫
(
is the most sensitive to the signals from
퐫
+
+
,
and
퐫
+
'
,
, and
cos
휃
≈
1
.
Thus, for the second case, the decay of
}
∫
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
}
as
}
퐫
+
+
,
−
퐫
+
'
,
}
increases is mainly
determined by the decay of
(
푘
(
⋆
푘
(
)
3
Z
퐫
!
+
"
0
퐫
!
'
"
Z
1
5
.
W
e compare the
decay
speed of
MCC
(
+
,
(
'
as
}
퐫
(
+
−
퐫
(
'
}
increases and that of
(
푘
(
⋆
푘
(
)
3
Z
퐫
!
+
"
0
퐫
!
'
"
Z
1
5
as
}
퐫
+
+
,
−
퐫
+
'
,
}
increases
by observing
t
he full
-
width values
in
Supplementary Fig.
13
c
and
d
. For
both
calibrated virtual transducers
and all choices of
훽
,
values of
FW
#
+
,
X
and
FW
#
'
,
X
are of 4 to 8 times the value of
FW
;Y
,
X
.
The
substantial
decay
-
speed difference between
FW
#
+
,
X
or
FW
#
'
,
X
with
FW
;Y
,
X
leads to the
anisotropy
of
}
∫
푠
+
+
(
푡
)
푠
+
'
(
푡
)
d
푡
I
J
&
}
’s
decay
speed as
}
퐫
+
+
,
−
퐫
+
'
,
}
increases, which
contributes to
the
anisotropy
of the
image resolution. It needs to be noted that, we
have
simplif
ied
the above analysis to identify the
cross
-
correlation
-
related and auto
correlation
-
related
limiting factors. In practice, the
anisotropy
is affected by both
the maximum
-
cross
-
correlation
matrix
MCC
(
+
,
(
'
and the autocorrelation function
(
푘
(
⋆
푘
(
)
(
휏
)
for all values of
푛
'
,
푛
$
, and
푛
with different weights. These factors are incorporated into the iterative method and lead to the
final
anisotropy
of resolution in the reconstructed image.
8
Supplementary Fig.
1
|
Comparison
of
uniform optical absorbers for calibration.
a
–
d
, PACTER signals generated by 500 consecutive laser pulses (8.6 μJ pulse energy, 1 kHz repetition rate)
using black tape (
a
), black rubber (
b
), black ink (
c
), and bovine blood (
d
) as targets. Norm., normalized.
e
,
Correlation coefficients between the PACTER signals and the average of the signals from the
first 10 laser
pulses.
9
Supplementary Fig.
2
|
Fly’s eye homogenizer in PACTER.
a
, Schematic illustrating the working principle of the fly’s eye homogenizer.
b,c
, Photographs of the original
beam from the laser (
b
) and the homogenized beam (
c
), captured using a laser alignment plate 1 cm away
from the top of the
ER
.
d
, Schematic showing the optical path of the homogenizer and the
ER
.
10
Supplementary Fig.
3
|
PACTER with different illumination areas.
a
,
Photo
graph
s showing PACTER
imaging experiments when the illumination area is smaller (left),
equivalent (middle), and larger (right) compared with the calibrated area
.
b
, Perspective views of 3D
PACTER images of
the
black wire in
a
, acquired when the illumination area is smaller (left),
equivalent
(middle), and larger (right) compared with the calibrated area.
Norm., normalized. Scale bar, 1 mm.
11
Supplementary Fig.
4
|
Positioning of the single
-
element ultrasonic transducer
on the ER
.
a
–
c
,
ER, excitation laser beam, and single
-
element ultrasonic transducer viewed at different angles. The
blue solid circle indicates the area of the excitation beam on the hypotenuse surface of the prism of the ER.
The yellow dashed circle
indicates
a clear area, where the laser beam does not interact with the ultrasonic
transducer.