of 66
1
D
ata-driven system matrix manipulation enabling fast functional imaging and
intra-image nonrigid motion correction in tomography
Pen
g Hu
1
, X
in Tong
1
, L
i Lin
1,
2
, an
d Lihong V. Wang
1,
*
1
C
altech Optical Imaging Laboratory, Andrew and Peggy Cherng Department of Medical
Engineering, Department of Electrical Engineering, California Institute of Technology,
Pasadena, CA 91125, USA
2
P
resent address: College of Biomedical Engineering and Instrument Science, Zhejiang
University, Hangzhou 310027, China
*Correspondence should be addressed to L.V.W. (LVW@caltech.edu).
A
bstract
T
omographic imaging modalities are described by large system matrices. Sparse sampling and
tissue motion degrade system matrix and image quality. Various existing techniques improve the
image quality without correcting the system matrices. Here, we compress the system matrices to
improve computational efficiency (e.g., 42 times) using singular value decomposition and fast
Fourier transform. Enabled by the efficiency, we propose (1) fast sparsely sampling functional
imaging by incorporating a densely sampled prior image into the system matrix, which maintains
the critical linearity while mitigating artifacts and (2) intra-image nonrigid motion correction by
incorporating the motion as subdomain translations into the system matrix and reconstructing the
translations together with the image iteratively. We demonstrate the methods in 3D photoacoustic
computed tomography with significantly improved image qualities and clarify their applicability
to X-ray CT and MRI or other types of imperfections due to the similarities in system matrices.
In
troduction
T
omographic imaging modalities X-ray computed tomography (CT), magnetic resonance imaging
(MRI), and photoacoustic computed tomography (PACT) produce cross-sectional images of tissue
b
y detection of penetrating X-rays
1
, n
uclear-magnetic-resonance-induced radio waves
2,
3
,
and
light-absorption-induced ultrasonic waves
4
,
respectively. Each modality with a certain setup is
described by a system matrix
5–1
0
. A
ccurate image reconstruction poses requirements to the system
matrix, which are often violated. For example, to achieve high temporal resolution for functional
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
2
i
maging, the spatial sampling density is often sacrificed, which introduces artifacts in the
reconstructed image
11–13
an
d may affect the functional signal extraction. Also, tissue motions such
as heart beating
1
4–16
,
breathing
17–19
, ab
dominal movement
20
,21
,
and fetal movement
2
2
, cau
se
complex geometric errors in each system matrix, which introduce artifacts in the reconstructed
image and compromise valuable image features.
Numerous methods have been proposed to compensate for system-matrix imperfections from
image-domain
12
,23–27
,
signal-domain
28–33
, an
d cross-domain
34–38
p
erspectives. However, due to the
large size of each system matrix, these methods tend not to manipulate or correct the system matrix
directly and have limitations. For sparse sampling functional imaging, traditional
methods
1
2,25,26,32,38
m
itigate artifacts in images but their performances drop sharply as the sampling
density reduces. Deep neural networks (DNNs)
2
7–29,31,34,37,39
sh
ow high performance in mitigating
artifacts but tend to generate false image features when the sampling density is low, and they
require imaging-modality- and device-dependent datasets, which are not always available.
Moreover, most of the methods introduce nonlinearity while mitigating artifacts, which disrupts
the functional signals that are often much weaker than background signals.
For intra-image nonrigid motion correction, gating- and binning-based methods
1
6,30,33,35,36
ar
e
commonly used. However, they require repeated data acquisition, which is time-consuming and
infeasible for unrepeated motions. DNNs have also been used for motion correction
23
,24
,
however,
they need specific training datasets that are not universally available, and it is challenging to reject
falsely generated features in DNNs. Two system-matrix-level methods
4
0,41
h
ave been proposed for
motion correction. In the first method
40
,42
, t
he authors approximate general nonrigid motions with
localized linear translations, identify possible motion paths from multichannel navigator data, and
estimate the motion at each pixel using localized gradient-entropy metric in the image domain.
However, quantifying localized motion from only navigator data is not robust, especially when the
motion amplitude and noise level increase. In the second method
4
1
, t
he authors express breathing-
and heartbeat-induced motions in basis functions by performing singular value decomposition
(SVD) and resolve these motions in imaging. However, for general motions, especially unrepeated
motions, the method’s performance is unknown. The high computation cost of SVD in the method
also restricts its application to 3D imaging.
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
3
H
ere, we compress the system matrices using SVD and fast Fourier transform (FFT), which
enables efficient system matrix slicing and manipulation. Then, we use two new methods for
functional imaging and motion correction, respectively, by performing data-driven manipulations
of these matrices. Both methods are applicable to CT, MRI, and PACT. For sparse sampling
functional imaging, we incorporate a prior image into the system matrix to reduce unknown
variables in image reconstruction. Special configurations in the method maintain linearity in image
reconstruction while substantially mitigating artifacts, which is critical for weak functional signal
extraction. For intra-image nonrigid motion correction, we approximate the motion with localized
linear translations
4
0
.
Starting from an initially reconstructed motion-blurred image, we first
estimate the translation of each subdomain of the object by minimizing the difference between the
simulated signals from the subdomain and the detected signals. With the estimated translations,
we update the system matrix and reconstruct the image again. We iterate the correction-
reconstruction process to obtain the final image. This method does not require repeated data
acquisition and is effective for unrepeated motions.
In this work, we use 3D PACT for demonstration due to its representatively large size in
tomography: light-absorption-induced ultrasonic wave from every voxel in an image is detected
by every transducer element and the system matrix is intrinsically a 6D tensor
10
. We
apply the
proposed methods to both numerical simulations and
in vivo
experiments: mouse brain functional
imaging with sparse sampling and intra-image motion correction in human breast imaging. We
demonstrate that both methods substantially improve the functional and structural image qualities.
R
esults
Sy
stem matrix compression based on SVD and FFT
In this study, we use the 3D PACT system reported previously by Lin et al.
4
3
,
which consists of
four 256-element arc transducer arrays (central frequency of 2.25 MHz and one-way bandwidth of
98%), and we assume a homogenous medium. The rotation of the four arc-arrays forms a virtual
hemispherical array truncated at the bottom for light delivery (marked as blue arcs in
Fig. 1a
), and
the detected ultrasonic signals are used to form a 3D image (marked as a rectangular cuboid). We
denote the number of voxels in the image and the number of virtual elements in the 2D array as
and
, respectively. The system matrix is formed by
responses of all virtual elements to the
signals from all voxels, with three responses’ paths shown as black-dotted arrows in
Fig. 1a
. It is
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
4
ex
pensive to store and manipulate these responses directly from memory and computation
perspectives, respectively. Considering that the responses’ arrival times are determined by the
voxel-element distances, we temporally shift the responses so that the arrival times are aligned
along the blue-dotted line in
Fig. 1b
. Performing grouped SVDs to the shifted responses, we
approximate these responses as linear combinations (coefficients shown in
Fig. 1b
) of three
temporal singular functions with high accuracy: an efficient compression of the system matrix. In
Supplementary Note 1
, we elaborate on the process of performing SVD to the responses grouped
for each virtual element in the element’s local coordinate system. Due to the homogeneous-
medium assumption and identity of all elements in the system, SVDs for all virtual elements share
the same set of temporal singular functions, with the first three shown as red, gray, and blue curves,
respectively, in
Fig. 1b
. In
Supplementary Fig. 1e
i
, we demonstrate that the arrival-time-
dependent temporal shifting is essential for the compression based on SVD.
Fig. 1
System matrix compression and slicing.
a
A virtual 2D array (blue curves) formed by
rotation of four arc arrays, the image domain (a rectangular cuboid with red-solid edges), and paths
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
5
o
f three responses marked by black-dotted arrows from voxels (red dots) to transducer elements
(black dots).
b
Responses temporally shifted to remove the differences in arrival times (aligned
along the blue-dotted line) from the voxels to elements, and compression of these responses based
on SVD with the first three temporal singular functions shown as red, gray, and blue curves,
respectively.
c
Signals detected by an element approximated by three convolutions with the three
temporal singular functions (zero time indicted by the blue-dotted line), respectively.
d
System
matrix slicing in voxel indices (rectangular cuboids with red-dotted edges) and virtual element
indices (blue arcs with red boundaries).
The system matrix compression not only saves memory in storing the matrix but also enables
fast and accurate forward simulation. For each combination of a voxel and a virtual element,
instead of adding the whole weighted response to the signals, we only access three coefficients:
values of the first three spatial singular functions. After obtaining all the coefficients with different
arrival times for each virtual element, we perform three convolutions (implemented by FFT) to
obtain the signals, as shown in
Fig. 1c
. We summarize the implementations of the original and
compressed system matrices in Eqs.
(1)
and
(
2
)
, respectively. Through forward simulations, we
demonstrate that the compressed system matrix is approximately 42 times as fast as the original
system matrix with negligible signal-domain relative errors (maximum value of 0.005) for this
study (
Supplementary Fig. 2
). Further, we reconstruct images using a nonnegativity-constrained
iterative method (Eq.
(
S33
)
) with the compressed system matrix from signals simulated using the
original system matrix, and observe negligible image-domain relative errors (maximum value of
0.007) for this study (
Supplementary Fig. 3
).
The implementation of the system matrix based on Eq.
(
2
)
is not only efficient but also
explicit, meaning subsets of the system matrix are directly accessible. For example, to obtain the
signals from
v
oxels (enclosed in two rectangular cuboids with red-dotted edges in
Fig. 1d
)
detected by
v
irtual elements (blue arcs with red boundaries in
Fig. 1d
), we can access the
relevant system matrix elements directly and finish the computation in
o
f the time for
voxels and
elements. It needs to be noted that this scaling property critical for this study is
invalid for wave-equation-based implementations
9
o
f the system matrix.
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
6
Spa
rse sampling functional imaging using a hybrid method
We propose a hybrid method for functional imaging (Methods,
Supplementary Fig. 6
) by
combining an iterative method with the universal back-projection (UBP) method
44
.
We first
demonstrate the performance of the hybrid method by using it to reconstruct images from signals
of a numerical phantom acquired by virtual arrays with different numbers of arcs:
4
=
7
6,40,28,20,16,12
. Here
i
s the number of rotating locations of the four-arc array (with 128
transducer elements in each arc) in a virtual array. Images of the numerical phantom reconstructed
using UBP, the regularized iterative method (Eq.
(
S34
)
), and the hybrid method (Eq.
(
5
)
with a
prior image obtained by performing a smooth modulation to the numerical phantom) are shown in
the first three columns in
Fig. 2a
. The used virtual arrays are shown as blue arcs with red
boundaries in the fourth column in
Fig. 2a
. We discuss the selection of regularization parameter
in
Supplementary Note 6
and, specifically, the best regularization parameters for different values
of
4
(
Sup
plementary Fig. 5
) that are used in the iterative method in
Fig. 2a
. Maximum
amplitude projections (MAPs) of these images along the
-axis are shown in
Supplementary Fig.
7
. We see that, as
4
d
ecreases from 76 to 12, the artifacts in the images reconstructed using
UBP become more abundant. The regularized iterative method mitigates the relatively weak
artifacts in all images but failed to suppress the strong artifacts such as in the images with
4
=
1
6,12
. In contrast, with the help of the prior image, the hybrid method significantly mitigates the
artifacts. Quantitatively, for each method, we calculate the structural similarity index measures
(SSIMs) between the images with
4
= 4
0,28,20,16,12
and that with
4
= 7
6
, and
compare the values in
Fig. 2b
. The hybrid method performs the best in mitigating artifacts and
maintaining true features.
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
7
F
ig. 2
UBP, regularized iterative method, and the proposed hybrid method for sparse-sampling
imaging, and their linearity tests.
a
Images reconstructed by the three methods (first three columns)
from signals detected at sparsely distributed elements (red-bounded blue curves in the fourth
column) for
4
= 7
6,40,28,20,16,12
. Examples of maintained features and suppressed
artifacts are indicated by white-solid and white-dotted arrows, respectively.
b
SSIMs between the
reconstructed images with
4
= 4
0,28,20,16,12
and those with
4
= 7
6
for the three
methods.
c
Linearity tests of the three methods for
4
= 1
2
. Scale bar, 5 mm.
Additionally, we test the linearity of each method by reconstructing two numerical phantoms
and their summation (shown in
Supplementary Fig. 8a
c
, respectively). MAPs of the
reconstructed images and the linearity test residuals (image of the summation subtracted by images
of the two phantoms) of the three methods with
4
= 7
6,12
are shown in
Supplementary Fig.
8e
-
f
and
Supplementary Fig. 9c
, respectively, which validate that the hybrid method maintains
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint
8
l
inearity while mitigating artifacts. We summarize the MAPs in the linearity tests with
4
= 1
2
in
Fig. 2c
. In another test of the hybrid method using a prior image with dot artifacts
(
Supplementary Fig. 9b
), the linearity is still valid although artifacts appear in the reconstructed
images (
Supplementary Fig. 9d
). In summary, UBP and the hybrid method are linear, but the
regularized iterative method is nonlinear. In the hybrid method, the prior image quality does not
affect the linearity but affects the reconstructed image quality.
Due to the requirement of a prior image, the hybrid method’s practical value mainly lies in fast
functional imaging with sparse sampling. We perform numerical simulations of functional imaging
with
4
= 1
2
in
Supplementary Note 7
. After reconstructing images from simulated signals
using the three methods, we extract functional signals from them using a method based on
regularized correlation (Eq.
(
7
)
) and form functional images. As shown in
Supplementary Fig.
11
and
Supplementary Video 1
, artifacts in the UBP-reconstructed images cause artifacts in the
functional images, the regularized iterative method mitigates artifacts in functional images but also
compromises the true functional region, and the proposed hybrid has the best performance.
Moreover, we apply UBP, the regularized iterative method, and the hybrid method to mouse
brain functional imaging
in vivo
using the four-arc system. We first obtain a prior image of a mouse
brain through dense sampling (
4
= 3
96
), then electrically stimulate its right front paw and
continuously acquire signals from the mouse brain through sparse sampling (
4
= 7
6
, 2 s per
image). We use subsets of the sparsely sampled signals (
4
= 4
0,20,12
) to demonstrate the
performance of the hybrid method. For one set of sparsely sampled signals, the images
reconstructed using UBP, the regularized iterative method, and the hybrid method for
4
=
4
0,20,12
are shown in
Fig. 3a
. We observe that the iterative method mitigates the artifacts (e.g.,
those indicated by white-dotted arrows) but compromises low-amplitude features (e.g., those
indicated by white-solid arrows for
4
= 2
0
). In contrast, the hybrid method maintains low-
amplitude features while substantially mitigating the artifacts, resulting in images more similar to
the densely sampled image. Electrical stimulation of the mouse’s right front paw occurs in five
cycles, each with 12-s stimulation on and 12-s off, as shown in
Fig. 3b
. To find the best
regularization parameter (
)
in the regularized correlation (Eq.
(
7
)
), we obtain functional images
from the images reconstructed through the three methods for
4
= 4
0,20,12
using
=
0.0
2,0.08,0.32,1.28,5.12
. For UBP and the hybrid method,
= 0.3
2
is the best choice to
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted January 8, 2024.
;
https://doi.org/10.1101/2024.01.07.574504
doi:
bioRxiv preprint