CT short=CT, long=computed tomography, FC short=FC, long=fully connected,
ACRIN short=ACRIN, long=American College of Radiology Imaging Network , GT
short=GT, long=ground truth, SD short=SD, long=standard deviation, PReLU
short=PReLU, long=parametric rectified linear unit , DUGAN short=DUG,
long=double U-Net generator, GAN short=GAN, long=generative adversarial
network, SSIM short=SSIM, long=structural similarity index, MSE short=MSE,
long=mean squared error, RMSE short=RMSE, long=root mean squared error, PET
short=PET, long=positron emission tomography, SPECT short=SPECT, long=single
photon emission computed tomography,class = in, MRI short=MRI, long=magnetic
resonance imaging,class = in, DnCNN short=DnCNN, long=denoising convolutional
neural network , ROI short=ROI, long=region of interest, long-plural-form =
regions of interest, SNR short=SNR, long=signal-to-noise ratio, MLEM
short=MLEM, long=maximum likelihood expectation-maximization, MBIR short=MBIR,
long=model-based iterative reconstruction, XCAT short=XCAT, long=extended
cardiac-torso, GPU short=GPU, long=graphical processing unit, CNN short=CNN,
long=convolutional neural network, CNR short=CNR, long=contrast-to-noise
ratio, FDG short=${}^{\text{18}}$F-FDG,
long=${}^{\text{18}}$F-fludeoxyglucose, FLT short=FLT, long=fluorothymidine ,
SUV short=SUV, long=standardized uptake values, 4D short=4-D, long=four-
dimensional, class = out, 3D short=3-D, long=three-dimensional, 2D short=2-D,
long=two-dimensional, FBP short=FBP, long=filtered-backprojection, OSEM
short=OSEM, long=ordered-subsets expectation-maximisation, mumap
short=${\mu}$-map, long=attenuation map, index=$\boldmu$-map, UBO short=UBO,
long=_Université de Bretagne Occidentale_ , HU short=HU, long=Hounsfield unit,
SR short=SR, long=Super Resolution,

# DUG-RECON: A Framework for Direct Image Reconstruction using Convolutional
Generative Networks

V.S.S. Kandarpa∗ Alexandre Bousse  Didier Benoit and Dimitris Visvikis All the
authors are affiliated to LaTIM, INSERM, UMR 1101, _Université de Bretagne
Occidentale_.
∗ (email<EMAIL_ADDRESS>
###### Abstract

This paper explores convolutional generative networks as an alternative to
iterative reconstruction algorithms in medical image reconstruction. The task
of medical image reconstruction involves mapping of projection domain data
collected from the detector to the image domain. This mapping is done
typically through iterative reconstruction algorithms which are time consuming
and computationally expensive. Trained deep learning networks provide faster
outputs as proven in various tasks across computer vision. In this work we
propose a direct reconstruction framework exclusively with deep learning
architectures. The proposed framework consists of three segments, namely
denoising, reconstruction and super resolution. The denoising and the super
resolution segments act as processing steps. The reconstruction segment
consists of a novel DUGAN which learns the sinogram-to-image transformation.
This entire network was trained on PET and CT images. The reconstruction
framework approximates 2D mapping from projection domain to image domain. The
architecture proposed in this proof-of-concept work is a novel approach to
direct image reconstruction; further improvement is required to implement it
in a clinical setting.

###### Index Terms:

Medical Image Reconstruction, Deep Learning, Generative Adversarial Networks

##  I Introduction

The use of deep learning in medical imaging has been on the rise over the last
few years. It has widely been used in various tasks across medical imaging
such as image segmentation [1, 2, 3, 4, 5], image denoising [6, 7, 8, 9],
image analysis [10, 11, 12]. The utilization of deep learning for image
reconstruction is a more challenging task. Image reconstruction using deep
learning corresponds to the task of mapping raw projection data retrieved from
the detector to image domain data. One can broadly identify three different
categories of approaches for the implementation of deep learning within the
framework of medical image reconstruction:

  * (i)

methods that use deep learning as an image processing step that improves the
quality of the raw data and/or the reconstructed image [13, 14];

  * (ii)

methods that embed deep-learning image processing techniques in the iterative
reconstruction framework to accelerate convergence or to improve image quality
[15, 16, 17];

  * (iii)

direct reconstruction with deep learning alone without any use of traditional
reconstruction methods [18, 19, 20].

The use of deep learning for the development of either data corrections or
post-reconstruction image based approaches (i) has shown potential to improve
the quality of reconstructed images. An example of data corrections by
improving the raw data through scatter correction is proposed in [14]. In this
work a modified U-Net is used to estimate scatter and correct the raw data in
order to improve CT images. Denoising the reconstructed PET images with a deep
convolutional network was done in [13]. The authors used perceptual loss along
with MSE to preserve qualitative and quantitative accuracy of the
reconstructed images. The network was initially trained on simulated data and
then fine-tuned on real patient data. Despite resulting in an improvement of
the reconstructed output, the above mentioned methods do not directly
intervene with the reconstruction process. This can be done using the two
distinct frameworks (ii) and (iii).

The first one involves the incorporation of a deep neural network into an
unrolled iterative algorithm where a trained neural network accelerates the
convergence by improving the intermediate estimates in the iterations [17, 15,
16]. The paper by Gong et al. used a modified U-Net to represent images within
the iterative reconstruction framework for PET images. The deep learning
architecture was trained on low-dose reconstructed images as input and high-
dose reconstructed images as the output. The work by Xie et al. further
extended this work by replacing the U-Net with a GAN for image representation
within the iterative framework. Kim et al incorporated a trained DnCNN along
with a novel local linear fitting function into the iterative algorithm. The
DnCNN which is trained on data with multiple noise levels improves the image
estimate at each iteration. They used simulated and real patient data in their
study. The second framework, referred to as _direct reconstruction_ , is based
on performing the whole reconstruction process (replacing the classical
framework) with a deep learning algorithm, using raw data as input and
reconstructed images as output [18, 19, 20]. In contrast to (i) and (ii) which
have been extensively investigated, direct reconstruction (iii) using deep
learning has been much less explored.

There have been to date three particular approaches that are relevant to this
strategy. The deep learning architecture proposed by Zhu et al. [19] called
AUTOMAP uses FC layers (which encode the raw data information) followed with
convolutional layers. The first three layers in this architecture are FC
layers with dimensions $2n^{2}$,$n^{2}$ and $n^{2}$ respectively where
$n\times{}n$ is the dimension of the input image. The AUTOMAP requires the
estimation of a huge number of parameters which makes it computationally
intensive. Although initially developed for MRI, AUTOMAP has been claimed to
work on other imaging modalities too. Brain images encoded into sensor-domain
sampling strategies with varying levels of additive white Gaussian noise were
reconstructed with AUTOMAP. Within the same concept of using FC layers’
architectures a three stage image reconstruction pipeline called DirectPET has
been proposed to reduce associated computational issues [18]. The first stage
down-samples the sinogram data, following which a unique Radon transform layer
encodes the transformation from sinogram to image space. Finally the estimated
image is improved using a super resolution block. This work was applied to
full body PET images and remains the only approach that can reconstruct
multiple slices simultaneously (up to 16 images). DeepPET is another approach
implemented on simulated images using encoder-decoder architecture based on
the neural network proposed by the visual geometric group [20]. Using
realistic simulated data, they demonstrated a network that could reconstruct
images faster, and with an image quality (in terms of root mean squared error)
comparable to that of conventional iterative reconstruction techniques.

In our work we explore the use of U-Net based deep learning architectures [1]
to perform a direct reconstruction from the sinogram to the image domain using
real patient datasets. Our aim is to reduce the number of trainable parameters
along with exploring a novel strategy for direct image reconstruction using
generative networks. More specifically our approach consists of a three-stage
deep-learning pipeline consisting of denoising, image reconstruction and super
resolution segments. Our experiments included training the deep learning
pipeline on PET and CT sinogram-image pairs. A single pass through the trained
network transforms the noisy sinograms to reconstructed images. The
reconstruction of both PET and CT datasets was considered and presented in the
following sections. At this stage the work presented is a proof of concept and
needs further improvement before being applied in a clinical setting.

Figure 1: Deep Learning in Medical Image Reconstruction

##  II Materials and Methods

###  II-A Image Reconstruction Model

In medical imaging, image reconstruction corresponds to the task of
reconstructing an image $\bm{x}\in\mathbb{R}^{m}$, from a scanner measurement
$\bm{y}\in\mathbb{R}^{n}$, where $m$ is the number of voxels defining the
image and $n$ is the number of detectors in the scanner. In CT, the image
$\bm{x}=\boldmu$ corresponds to X-ray attenuation, measured by the proportion
of X-rays scattered or absorbed as they pass through the object. In PET,
$\bm{x}=\bm{\lambda}$ is the distribution of a radiotracer delivered to the
patient by injection, and is measured through the detection of pairs of
$\gamma$-rays emitted in opposite directions (indirectly from the positron-
emitting radiotracer).

The measurement $\bm{y}$ is a random vector modeling the number of detection
(photon counting) at each of the $n$ detector bins, and follows a Poisson
distribution with independent entries:

$\bm{y}\sim\mathrm{Poisson}(\bm{\bar{y}}(\bm{x}))$ (1)

where $\bm{\bar{y}}(\bm{x})\in\mathbb{R}^{n}$ is the expected number of counts
(noiseless), which is a function of the image $\bm{x}$. In a simplified
setting, the expected number of counts in CT is

$\bm{\bar{y}}(\boldmu)=\exp(-\bm{L}\boldmu)$ (2)

where $\bm{L}\in\mathbb{R}^{n\times m}$ is a system matrix such that each
entry $[\bm{L}]_{i,j}$ represents the contribution of the $j$-th image voxel
to the $i$-th detector. In PET, the expected number of counts is (also in a
simplified setting) is

$\bm{\bar{y}}(\bm{\lambda})=\bm{P}\bm{\lambda}$ (3)

where $\bm{P}\in\mathbb{R}^{n\times m}$ is a system matrix such that each
entry $[\bm{P}]_{i,j}$ represents the probability that a photon pair emitted
from voxel $j$. For simplification we assume $\bm{P}=\bm{L}$ which is a
reasonable assumption in a non time-of-flight setting. Image reconstruction is
achieved by finding a suitable image $\bm{\hat{x}}=\bm{\hat{\mu}}$ or
$\bm{\hat{\lambda}}$ that approximately solves

$\bm{y}=\bm{\bar{y}}(\bm{x})\,.$ (4)

FBP techniques (see [21] for a review) can efficiently solve (4) but they are
vulnerable to noise as the system matrix $\bm{P}$ is ill-conditioned. Since
the 80’s, MBIR techniques [22, 23] became the standard approach. They consist
in iteratively approximating a solution $\bm{\hat{x}}$ such that
$\bm{\bar{y}}(\bm{\hat{x}})$ maximizes the likelihood of the measurement
$\bm{y}$. As they model the stochasticity of the system, they are more robust
to noise as compared with FBP, and can be completed with a penalty term for
additional control over the noise [24]. However, they are computationally
expensive as each iteration requires to compute a projection and a
backprojection.

Image reconstruction with deep learning however is a data driven approach
wherein there is a training and a prediction phase. Given a set of training
data which is a subset of the raw data ($\bm{y}$) and its corresponding images
($\bm{x}$), a deep learning architecture learns the mapping from raw data to
the image and improves this mapping through the training process. During the
prediction phase a subset of raw data different from the training data serves
as the input to the trained deep learning architecture. The output is a
reconstructed image which is obtained on a single forward pass through the
network. Hence making the reconstruction process through deep learning
instantaneous as opposed to an iterative process. This makes direct
reconstruction with deep learning faster and less computationally expensive
than iterative algorithms.

###  II-B Deep Learning Architectures

As shown in Figure 2, we propose a three-stage deep learning pipeline for the
task of tomographic reconstruction. In the first step the raw data (projection
space) are denoised. Next the denoised sinograms are transformed to the image
domain in the image reconstruction segment. The third and final segment
operates in the image domain to improve the image produced after domain
transformation. The following sections discuss these segments in detail.

Figure 2: Proposed Deep Learning pipeline for Direct Image Reconstruction

####  II-B1 Denoising

We used a modified U-Net architecture to denoise the Poisson sampled
sinograms, based on the work previously carried out for ultrasound denoising
[25]. The U-Net is an encoder-decoder network which was initially implemented
for segmentation but over the years its applications have broadened. As shown
in Figure 3 there are increasing number of convolutions along with max pooling
to arrive at an encoding of the input and then with convolutions followed by
upsampling, arriving at the output with an identical dimension as the input.
The important modification in the architecture mentioned in [25] with respect
to U-Net was the residual connection from the input to the final output.
Perdios et al trained the denoising architecture on simulated ultrasound
images so as to enhance ultrafast ultrasound imaging. This denoising
architecture corresponds to the first segment in our proposed framework. It
was trained on raw data pairs, i.e., low-count and high-count sinograms,
considering multiple noise levels. The detailed architecture is represented in
Figure 3. We defined the loss function between a true sinogram
$\bm{y}^{\star}=[y_{1}^{\star},\dots,y_{n}^{\star}]^{\top}\in\mathbb{R}^{n}$
and a prediction
$\bm{\hat{y}}=[\hat{y}_{1},\dots,\hat{y}_{n}]^{\top}\in\mathbb{R}^{n}$ as the
MSE:

$\mathrm{MSE}(\bm{y}^{\star},\bm{\hat{y}})=\frac{1}{n}\sum_{i=1}^{n}(y^{\star}_{i}-\hat{y}_{i})^{2}\,,$
(5)

where, $n$ is the number pixels on the sinogram, corresponding to the number
of detectors in the scanner.

####  II-B2 Image Reconstruction

The novelty in this work is the proposed U-Net based network in contrast to
previous works in direct image reconstruction using the FC layer
architectures. This design of the network draws its inspiration from
conditional GAN for image to image translation called Pix2Pix [26]. The
proposed network namely DUGAN consists of two cascaded U-Nets. The first U-Net
transforms the raw data to image while the second U-Net takes as input the
generated image and transforms it back to the raw data. The second U-Net
assesses the reconstructed image output, reiterating the relation between the
sinogram and the image. This architecture differs from the Pix2Pix, which
consists of a generator (U-Net like network) and a discriminator
(classification convolutional network). While the generator in both
architectures serves the purpose of transforming images from one domain to the
other, the discriminator with regards to Pix2Pix classifies inputs as
real/fake. The objective function for this architecture can be written as:

$\mathcal{L}_{\mathrm{total}}=\mathcal{L}_{G_{1}}+\mathcal{L}_{G_{2}}+\mathcal{L}_{G_{1}+G_{2}}$
(6)

where,

$\mathcal{L}_{G_{1}}=\frac{1}{n}\sum_{i=1}^{n}|x^{\star}_{i}-\hat{x}_{i}|$ (7)
$\mathcal{L}_{G_{2}}=\frac{1}{n}\sum_{i=1}^{n}|y^{\star}_{i}-\hat{y}_{i}|$ (8)

$G_{1},G_{2}$ are Generator 1 and Generator 2 which predict image and sinogram
respectively; $\bm{x}^{\star}$, $\bm{\hat{x}}$ the true and predicted image,
$\bm{y}^{\star}$, $\bm{\hat{y}}$ the true and predicted sinogram respectively.
$\mathcal{L}_{G_{1}+G_{2}}$ is defined similar to $\mathcal{L}_{G_{2}}$ with
the combined architecture of $G_{1}+G_{2}$, keeping the weights of $G2$ fixed.

The architecture is represented in detail in Figure 4. The training for this
architecture is summarized in Algorithm 1. A comparison of the trainable
parameters for various segments that are used to perform the domain mapping
from sinogram to image is provided in Table I, considering the AUTOMAP from
[19], the Radon inversion layer from [18] and the proposed architecture DUGAN
along with denoising and super resolution segments.

$M$ = number of epochs ;

$N$ = total training data (images/sinograms) ;

for  _i = 1,2,…, M_ do

for  _j = 1,2,…,N_ do

Train $G_{1}$: minimizing $\mathcal{L}_{G_{1}}$ ;

end for

for  _j = 1,2,…,N_ do

Train $G_{2}$: minimizing $\mathcal{L}_{G_{2}}$ ;

end for

for  _j = 1,2,…,N_ do

Train combined architecture, freezing the weights of $G_{2}$: minimizing
$\mathcal{L}_{G_{1}+G_{2}}$ ;

end for

end for

Algorithm 1 Training the DUG TABLE I: Trainable Parameters comparison Architecture | Input Size | Output Size | Trainable
---|---|---|---
|  |  | Parameters
AUTOMAP | $200\times{}168$ | $200\times{}200$ | 6,545,920,000
Radon Inversion Layer | $200\times{}168$ | $200\times{}200$ | 382,259,200
($40\times{}40$ Patch size) |  |  |
DeepPET | $128\times{}128$ | $128\times{}128$ | 62,821,473
DUG-RECON | $128\times{}128$ | $128\times{}128$ | 17,444,140
Figure 3: Representation of the denoising network. The inputs to the network
were 2D grayscale slices with resolution $128\times{}128$ and the outputs were
denoised sinograms. Figure 4: Representation of the DUGAN, the image
reconstruction block. This network was trained on denoised sinograms which
were the outputs of the previous segment. Figure 5: Representation of the
super resolution block. It consists of 8 residual blocks with Convolution,
Batch normalization and PReLu.

####  II-B3 Super Resolution

The function of the SR is to improve the estimate produced by the image
reconstruction network. Several works already exist concerning single image
super resolution [27, 28]. In this work we employed a basic super residual
network architecture to improve the reconstruction. It consists of
convolutional blocks followed by batch normalization with PReLU activation.
There were a total of 8 residual blocks in the network as represented in
Figure 5. The loss function used in this architecture was perceptual loss:

$\mathrm{PerceptualLoss}=|\mathrm{VGG}_{16}(\bm{x}^{\star})-\mathrm{VGG}_{16}(\bm{\hat{x}})|$
(9)

$\mathrm{VGG}_{16}(\bm{x}^{\star})$ and $\mathrm{VGG}_{16}(\bm{\hat{x}})$ are
the extracted features with $\mathrm{VGG}_{16}$ convolutional neural network
[29] for the true and predicted image.

The features are extracted from the 10th layer of the VGG architecture i.e.,
only the first three convolutional blocks are considered. We observed that
extracting deeper features led to the network hallucinating features in the
reconstructed images.

###  II-C Dataset Description

We applied our methodology on fluorothymidine (FLT) PET/CT images from the
American College of Radiology Imaging Network (ACRIN) FLT Breast PET/CT
database [30]. The details of the dataset are given in Table II. The sinograms
were initially generated by projecting 2-D PET and CT images slices with the
Python SKLEARN Radon transform, following the models (2) and (3) for CT and
PET respectively, with Poisson noise added. The methodology represented in
Figure 6 was used for data preparation for the PET and CT modalities
respectively. Sample pairs from the PET and CT datasets are shown in Figures 7
and 8. The CT images were downsized from $512\times{}512$, and the
reconstruction was implemented for 2-D $128\times{}128$ images.

TABLE II: Dataset Description Dataset Statistics |
---|---
Modalities | CT, PET
Number of Patients | 83
Number of PET 2D Image slices | 76,000
Number of CT 2D Image slices | 21,104
PET Matrix size | 128
CT Matrix size | 512
Scanner | GE Discovery ST
Figure 6: Data preparation Figure 7: Example PET sinogram-image pairs from the
dataset Figure 8: Example CT sinogram-image pairs from the dataset

###  II-D Training

TensorFlow [31] and Keras [32] were used for the realization of the
architectures described in the section above. These architectures were
implemented on a single Nvidia GeForce GTX 2080Ti GPU. A collection of images
$\\{\bm{x}^{\star}_{k}\\}_{k=1}^{N}$ was used to generate a corresponding
collection of noiseless sinograms $\\{\bm{y}^{\star}_{k}\\}_{k=1}^{N}$
following models (2) and (3), low-counts and high-counts sinograms,
$\\{\bm{y}^{\mathrm{LC}}_{k}\\}_{k=1}^{N}$ and
$\\{\bm{y}^{\mathrm{HC}}_{k}\\}_{k=1}^{N}$ sinograms were generated by adding
Poisson noise, where the expected number of counts was adjusted by rescaling
the intensity. The denoising segment was trained with the collection
$\\{(\bm{y}^{\mathrm{LC}}_{k},(\bm{y}^{\mathrm{HC}}_{k})\\}_{k=1}^{N}$, using
the high-count sinograms as ground truth, with total number of training
samples $N$ = 120000 and number of epochs $M$ = 100. For the training of the
second segment, we used a collection of denoised sinograms, namely
$\\{\bm{\hat{y}}_{k}\\}_{k=1}^{N}$, and their corresponding ground true images
$\\{\bm{x}^{\star}_{k}\\}_{k=1}^{N}$. The image reconstruction segment was
trained according to the Algorithm 1, that is to say by alternating between
training $G_{1}$ with $\\{(\bm{\hat{y}}_{k},\bm{x}^{\star}_{k})\\}_{k=1}^{N}$
and training $G_{2}$ with
$\\{(\bm{\hat{x}}_{k},\bm{y}^{\star}_{k})\\}_{k=1}^{N}$, where
$\bm{\hat{x}}_{k}$ is a prediction from $G_{1}$ with $\bm{\hat{x}}_{k}$ as an
input. This segment was trained with $N=40000$ and $M=50$. The CT data were
augmented by rotating the data by 90 degrees to generate the required training
data. Owing to the larger PET dataset it was not necessary to perform data
augmentation. Finally the SR segment was trained on the images predicted by
the DUGAN and the GT images
$\\{(\bm{\hat{x}}_{k},\bm{x}^{\star}_{k})\\}_{k=1}^{N}$ . The SR block was
trained with $N=20000$ for 100 epochs. For the testing, a set of 2000
sinogram-image pairs were used.

###  II-E Quantitative Analysis

Testing for the aforementioned architectures was done on samples that were not
a part of the training data. The metrics used for this analysis are RMSE and
SSIM Index. They are defined below:

$\mathrm{RMSE}(\bm{x}^{\star},\bm{\hat{x}})=\sqrt{\frac{1}{n}\sum_{j=1}^{m}(x^{\star}_{j}-\hat{x}_{j})^{2}}$
(10)

where $n$ is the number of pixels. $\hat{x}$ is the GT $x^{\star}$ the
predicted output.

$\mathrm{SSIM}(\bm{x}^{\star},\bm{x})=\frac{(2\mu_{\bm{x}^{\star}}\mu_{\bm{x}}+c_{1})(2\sigma_{\bm{x}^{\star}\bm{x}}+c_{2})}{(\mu_{\bm{x}^{\star}}^{2}+\mu_{\bm{x}}^{2}+c_{1})(\sigma_{\bm{x}^{\star}}^{2}+\sigma_{\bm{x}}^{2}+c_{2})}$
(11)

where $\mu_{\bm{x}^{\star}}$ $\mu_{\bm{x}}$ are the averages of
$\bm{x}^{\star}$ and $\bm{x}$ respectively, $\sigma_{\bm{x}^{\star}}^{2}$ and
$\sigma_{\bm{x}}^{2}$ are the variances of $\bm{x}^{\star}$ and $\bm{x}$,
$\sigma_{\bm{x}^{\star}\bm{x}}$ is the covariance between $\bm{x}^{\star}$ and
$\bm{x}$ , $c_{1}=(k_{1}L)^{2}$ and $c_{2}=(k_{2}L)^{2}$ where $k_{1}=0.01$
and $k_{2}=0.03$ by default.

###  II-F Region of Interest Analysis

The SNR and CNR were studied for four regions of interest identified within
the patient body. The SNR and CNR were evaluated by treating a region as
foreground and the other three regions as background.

$\mathrm{SNR}=\frac{\mu_{\mathrm{r}}-\mu_{\mathrm{b}}}{\sigma_{\mathrm{b}}}$
(12)

$\mathrm{CNR}=\frac{|{\mu_{\mathrm{r}}-\mu_{\mathrm{b}}}|}{\sqrt{\sigma_{\mathrm{r}}^{2}+\sigma_{\mathrm{b}}^{2}}}.$
(13)

where $\mu_{\mathrm{r}}$ and $\mu_{\mathrm{b}}$, $\sigma_{\mathrm{r}}$ and
$\sigma_{\mathrm{b}}$ correspond to the mean and standard deviation in the ROI
and the background respectively. In this study we compared the initial
reconstructed output of the DUGAN, the final reconstruction along with SR and
the original GT which was reconstructed with GE discovery ST using an OSEM
algorithm..

###  II-G Comparison with DeepPET

We implemented the architecture DeepPET [20] and compared the predictions with
our proposed approach for the reconstruction of PET images. DeepPET was
trained on $\\{(\bm{\hat{y}}_{k},\bm{x}^{\star}_{k})\\}_{k=1}^{N}$, notation
similar to the training section from above, with $N$ = 120000, exclusively on
PET data. It is worth noting that the input and output dimensions in our study
are identical while it was not originally for DeepPET. The architecture of
DeepPET is summarized in Figure 5. This architecture was trained for 100
epochs with an Adam optimiser.

Figure 9: Representation of DeepPET. The number of filters in each
convolutional layer is labeled on top of each block.

##  III Results

The predictions from the architectures along with the GT and the sinogram are
shown in Figure 10 for PET images. The results are displayed for four test
image slices across the columns. Each column shows the predicted output by the
proposed DUG-RECON architecture and the DeepPET architecture, as well as the
GT. With regards to the proposed architecture it is observed that the initial
reconstructed image i.e., the output of the DUGAN looks blurred while final
reconstructed output from the super resolution block has noticeably improved
details. The predictions by DeepPET are also visibly blurred compared to the
final reconstructed output of the proposed architecture and the GT. These
observations are further ascertained in Table III where the quantitative
metrics are tabulated. Figure 11 provides a comparison of the intensity
profiles for the predictions by DUGAN, DUG+SR and DeepPET w.r.t. to the GT for
PET images. These intensity values are observed along the line marked in
yellow in these figures. As this figure shows, the intensity profile of the
final reconstructed image of the proposed architecture is closest to the GT.
The predictions by DUGAN and DeepPET are smoother compared to the predictions
by DUGAN+SR and the GT. The ROI analysis is tabulated in Table V for the four
regions marked in Figure 12. This analysis was carried out for final
predictions by the proposed architecture and MLEM. Looking closely at Table V
we notice that the mean values of the deep learning predicted image and the
MLEM reconstructed image are comparable. The results for CT images are
displayed in Figure 13. This Figure provides a comparison between
reconstructed image predictions with the proposed architecture, DeepPET with
respect to the GT. The high-resolution nature of the CT images and a smaller
dataset presented challenges during the training of the proposed architecture.
The predictions by DeepPET appear to be better resolved than those by the
proposed architecture. However, the tissue and the bone structures are not
clearly seen in the predictions by the deep learning architectures, thereby
requiring further work to improve the reconstruction. The intensity plots are
compared for two different images in the Figure 14. The ROI analysis was
carried out for 4 regions marked in the images reconstructed with deep
learning and FBP. The image reconstructed with FBP has better SNR and CNR
compared to the image reconstructed with the proposed architecture.

TABLE III: The SSIM and RMSE for the PET images are evaluated for 4 different 2D slices. These metrics are compared amongst predictions by DUGAN, DUG+SR and DeepPET Image | Architecture |  RMSE |  SSIM
---|---|---|---
1 | DUG | 0.059 | 0.74
| DUG+SR | 0.038 | 0.84
| DeepPET | 0.047 | 0.80
2 | DUG | 0.043 | 0.76
| DUG+SR | 0.046 | 0.86
| DeepPET | 0.054 | 0.85
3 | DUG | 0.050 | 0.76
| DUG+SR | 0.038 | 0.85
| DeepPET | 0.043 | 0.83
4 | DUG | 0.061 | 0.70
| DUG+SR | 0.045 | 0.82
| DeepPET | 0.048 | 0.79
TABLE IV: The SSIM and RMSE for the CT images are evaluated for 4 different 2D slices. Here the architecture indicates the prediction by DUGAN and that of DUGAN along with SR segment Image | Architecture |  RMSE |  SSIM
---|---|---|---
1 | DUG | 0.0083 | 0.90
| DUG+SR | 0.0015 | 0.98
| DeepPET | 0.0012 | 0.99
2 | DUG | 0.0081 | 0.90
| DUG+SR | 0.0015 | 0.99
| DeepPET | 0.0014 | 0.99
3 | DUG | 0.0015 | 0.91
| DUG+SR | 0.0018 | 0.98
| DeepPET | 0.0013 | 0.99
Figure 10: Image predictions by DUGAN+SR, DeepPET and GT for four PET Images from different parts of the patient volume Figure 11: Intensity Profile across the image (highlighted by a yellow line) for a PET image prediction by DUGAN, DUGAN+SR and DeepPET compared with the GT Figure 12: SNR and CNR comparison amongst DUG+SR and OSEM for PET image along 4 regions of interest Figure 13: Image predictions by DUGAN+ SR, DeepPET and GT are displayed for 3 CT Images along different parts of the patient volume. Figure 14: Intensity Profile for two CT images (highlighted by a yellow line) predicted by DUGAN and SR compared with the GT Figure 15: SNR and CNR comparison amongst DUG+SR and OSEM for CT image along 4 regions of interest TABLE V: ROI Analysis: The mean, SD and the SNR for the 4 regions of interest marked in Figure 12 Region | Image | Mean |  SD |  SNR |  CNR
---|---|---|---|---|---
1 | DUG+SR | 0.706 | 0.024 | 7.15 | 5.71
| MLEM | 0.676 | 0.035 | 5.72 | 4.55
2 | DUG+SR | 0.713 | 0.091 | 4.81 | 3.42
| MLEM | 0.648 | 0.11 | 4.38 | 3.26
3 | DUG+SR | 0.744 | 0.071 | 2.73 | 1.65
| MLEM | 0.547 | 0.154 | 3.23 | 1.22
4 | DUG+SR | 0.117 | 0.008 | 14.96 | 10.64
| MLEM | 0.057 | 0.010 | 8.01 | 4.8
TABLE VI: ROI Analysis: The mean, SD and the SNR for the 4 regions of interest marked in Figure 15 Region | Image | Mean |  SD |  SNR |  CNR
---|---|---|---|---|---
1 | DUG+SR | 0.011 | 2.91e-4 | 7.61 | 6.66
| FBP | 0.011 | 3.45e-4 | 15.86 | 10.69
2 | DUG+SR | 0.004 | 1.53e-4 | 9.34 | 9.02
| FBP | 0.004 | 1.84e-4 | 15.36 | 13.74
3 | DUG+SR | 0.005 | 5.90e-4 | 9.40 | 5.49
| FBP | 0.005 | 4.88e-4 | 16.19 | 7.78
4 | DUG+SR | 0.012 | 8.319e-4 | 15.47 | 5.92
| FBP | 0.012 | 2.736e-4 | 14.74 | 11.47

##  IV Discussion

Deep learning has been applied to different fields of medical imaging. The
vast majority of developments concern primarily image processing and
analysis/classification tasks. Few works devoted in the field of image
reconstruction have been largely concentrated in the use of deep learning
within classical tomographic reconstruction algorithms. The main objectives of
these works have been an improvement in the speed of convergence and the
quality of the successive image estimation within the iterative reconstruction
process. The alternative approach involving direct image reconstruction
through the use of deep learning approaches to estimate images directly from
the use of raw data (sinograms or projections) has been much less explored
both for PET and CT.

Most implementations in direct image reconstruction concern the use of fully
connected layers which encode the raw data followed by convolutional layers.
In most of the proposed implementations a large number of parameters need to
be optimised which reduces the computational burden and overall robustness. In
this work we have proposed an original direct image reconstruction deep
learning framework based on an architecture inspired by convolutional
generative adversarial networks used in image to image translation. The
implementation is based on the use of a double U-Net generator (DUG)
consisting of two cascaded U-Nets. While the first network transforms the raw
data to an image the second one assesses the reconstructed image output of the
first network by reiterating the relationship between the reconstructed image
and the raw data. Two additional blocks were added; namely a network denoising
the raw data prior to their input in the DUG network and a super-resolution
block operating on the DUG output image in order to improve it’s overall
quality. The proposed network was directly trained on clinical datasets for
both PET and CT image reconstruction and its performance was assessed
qualitatively and quantitatively.

Deep neural networks usually result in blurred output. This fact is clear in
the predictions made by the DUGAN network. Both the qualitative analysis using
the profiles through the reconstructed images and the quantitative metrics
SSIM and the MSE, demonstrate the improvement of the reconstructed images
resulting from the incorporation of the SR block. The qualitative analysis
also clearly demonstrates the superiority of the proposed algorithm for direct
PET image reconstruction in comparison to alternative approaches such as
DeepPET. Finally in the ROI analysis we observed that the SNR and CNR are
higher with the deep learning approach for the PET images while they are lower
than the traditional methods for CT images. This is consistent with the
observations in the qualitative analysis, where the proposed approach was not
able to sufficiently resolve different tissue classes in the resulting
reconstructed CT images in comparison with the ground truth.

One of the potential reasons of the worse performance of DUG-RECON for CT
reconstruction relative to the superior performance observed for PET image
reconstruction may be the lower number of available CT images in the training
process. This limitation will be addressed as part of future work. Despite the
lower performance of the proposed architecture for CT images it still presents
comparable predictions and opens up avenues for deep learning architectures in
tomographic reconstruction. In general, the limitations of a deep learning
based reconstruction is the adaptability to new data which is very different
from the training sample space. Once a practical methodology is identified,
one could have a deep learning pipeline with an ensemble of networks trained
on different datasets to perform the reconstruction task.

##  V Conclusion

We have demonstrated the use of generative convolutional networks for the
tomographic image reconstruction task. More specifically we have proposed a
new architecture for direct reconstruction that approximates the 2D
reconstruction process. Also we have significantly reduced the parameters
required for the domain transform task in image reconstruction. The three-step
training pipeline based exclusively on deep learning decentralises the various
tasks involved in image reconstruction into denoising, domain transform and
super resolution. Various super resolution strategies are currently being
explored to improve the reconstructed image. Our proposed strategy for
tomographic reconstruction will eventually lead to a network based
reconstruction as we continue to improve the framework. Currently it does not
perform better than traditional methods in terms of utility metrics but still
has the advantage of instantaneous reconstruction and an effective denoising
strategy. We plan to extend the work on realistic detector data generated
through Monte Carlo simulations in addition to sinograms obtained through
Radon transform. We are also working on adapting the architecture to raw
detector data. Another important aspect of the data based deep learning
approach is that the predictions are limited by the quality of the dataset. It
becomes essential to have realistic datasets without compromising on the image
quality to improve the training of the neural networks.

## References

  * [1] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [2] Z. Guo, X. Li, H. Huang, N. Guo, and Q. Li, “Deep learning-based image segmentation on multimodal medical imaging,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 3, no. 2, pp. 162–169, 2019.
  * [3] A. Sinha and J. Dolz, “Multi-scale guided attention for medical image segmentation,” _arXiv preprint arXiv:1906.02849_ , 2019.
  * [4] J. Dolz, K. Gopinath, J. Yuan, H. Lombaert, C. Desrosiers, and I. B. Ayed, “Hyperdense-net: a hyper-densely connected CNN for multi-modal image segmentation,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 38, no. 5, pp. 1116–1126, 2018.
  * [5] M. Hatt, B. Laurent, A. Ouahabi, H. Fayad, S. Tan, L. Li, W. Lu, V. Jaouen, C. Tauber, J. Czakon _et al._ , “The first MICCAI challenge on PET tumor segmentation,” _Medical image analysis_ , vol. 44, pp. 177–195, 2018\.
  * [6] V. S. Kadimesetty, S. Gutta, S. Ganapathy, and P. K. Yalavarthy, “Convolutional neural network-based robust denoising of low-dose computed tomography perfusion maps,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 3, no. 2, pp. 137–152, 2018.
  * [7] M. Li, W. Hsu, X. Xie, J. Cong, and W. Gao, “SACNN: Self-attention convolutional neural network for low-dose CT denoising with self-supervised perceptual loss network,” _IEEE Transactions on Medical Imaging_ , 2020.
  * [8] H. Chen, Y. Zhang, W. Zhang, P. Liao, K. Li, J. Zhou, and G. Wang, “Low-dose CT denoising with convolutional neural network,” in _2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017)_. IEEE, 2017, pp. 143–146.
  * [9] Q. Yang, P. Yan, Y. Zhang, H. Yu, Y. Shi, X. Mou, M. K. Kalra, Y. Zhang, L. Sun, and G. Wang, “Low-dose CT image denoising using a generative adversarial network with wasserstein distance and perceptual loss,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 37, no. 6, pp. 1348–1357, 2018.
  * [10] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” _Medical image analysis_ , vol. 42, pp. 60–88, 2017.
  * [11] A. Amyar, S. Ruan, I. Gardin, C. Chatelain, P. Decazes, and R. Modzelewski, “3-D RPET-NET: development of a 3-D PET imaging convolutional neural network for radiomics analysis and outcome prediction,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 3, no. 2, pp. 225–231, 2019.
  * [12] S. Cui, Y. Luo, H.-H. Tseng, R. K. Ten Haken, and I. El Naqa, “Artificial neural network with composite architectures for prediction of local control in radiotherapy,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 3, no. 2, pp. 242–249, 2018.
  * [13] K. Gong, J. Guan, C.-C. Liu, and J. Qi, “PET image denoising using a deep neural network through fine tuning,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 3, no. 2, pp. 153–161, 2018.
  * [14] J. Maier, S. Sawall, M. Knaup, and M. Kachelrieß, “Deep scatter estimation (DSE): Accurate real-time scatter estimation for X-ray CT using a deep convolutional neural network,” _Journal of Nondestructive Evaluation_ , vol. 37, no. 3, p. 57, 2018.
  * [15] Z. Xie, R. Baikejiang, K. Gong, X. Zhang, and J. Qi, “Generative adversarial networks based regularized image reconstruction for PET,” in _15th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine_ , vol. 11072. International Society for Optics and Photonics, 2019, p. 110720P.
  * [16] K. Kim, D. Wu, K. Gong, J. Dutta, J. H. Kim, Y. D. Son, H. K. Kim, G. El Fakhri, and Q. Li, “Penalized PET reconstruction using deep learning prior and local linear fitting,” _IEEE Transactions on Medical Imaging_ , vol. 37, no. 6, pp. 1478–1487, 2018.
  * [17] K. Gong, J. Guan, K. Kim, X. Zhang, J. Yang, Y. Seo, G. El Fakhri, J. Qi, and Q. Li, “Iterative PET image reconstruction using convolutional neural network representation,” _IEEE Transactions on Medical Imaging_ , vol. 38, no. 3, pp. 675–685, 2019.
  * [18] W. Whiteley and J. Gregor, “Direct image reconstruction from raw measurement data using an encoding transform refinement-and-scaling neural network,” in _15th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine_ , vol. 11072. International Society for Optics and Photonics, 2019, p. 1107225.
  * [19] B. Zhu, J. Z. Liu, S. F. Cauley, B. R. Rosen, and M. S. Rosen, “Image reconstruction by domain-transform manifold learning,” _Nature_ , vol. 555, no. 7697, p. 487, 2018.
  * [20] I. Haeggstroem, C. R. Schmidtlein, G. Campanella, and T. J. Fuchs, “Deeprec: A deep encoder-decoder network for directly solving the PET reconstruction inverse problem,” _arXiv preprint arXiv:1804.07851_ , 2018.
  * [21] F. Natterer, _The Mathematics of Computerized Tomography_ , ser. Classics in Applied Mathematics. SIAM, 2001.
  * [22] L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” _IEEE Transactions on Medical Imaging_ , vol. 1, no. 2, pp. 113–122, 1982.
  * [23] J. A. Fessler, M. Sonka, and J. M. Fitzpatrick, “Statistical image reconstruction methods for transmission tomography,” _Handbook of medical imaging_ , vol. 2, pp. 1–70, 2000.
  * [24] A. R. De Pierro, “A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography,” _IEEE Transactions on Medical Imaging_ , vol. 14, no. 1, pp. 132–137, 1995.
  * [25] D. Perdios, M. Vonlanthen, A. Besson, F. Martinez, M. Arditi, and J.-P. Thiran, “Deep convolutional neural network for ultrasound image enhancement,” in _2018 IEEE International Ultrasonics Symposium (IUS)_. IEEE, 2018, pp. 1–4.
  * [26] P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros, “Image-to-image translation with conditional adversarial networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 1125–1134.
  * [27] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang _et al._ , “Photo-realistic single image super-resolution using a generative adversarial network,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 4681–4690.
  * [28] B. Lim, S. Son, H. Kim, S. Nah, and K. Mu Lee, “Enhanced deep residual networks for single image super-resolution,” in _Proceedings of the IEEE conference on computer vision and pattern recognition workshops_ , 2017, pp. 136–144.
  * [29] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [30] L. Kostakoglu, F. Duan, M. O. Idowu, P. R. Jolles, H. D. Bear, M. Muzi, J. Cormack, J. P. Muzi, D. A. Pryma, J. M. Specht _et al._ , “A phase ii study of 3-deoxy-3-18f-fluorothymidine PET in the assessment of early response of breast cancer to neoadjuvant chemotherapy: results from ACRIN 6688,” _Journal of Nuclear Medicine_ , vol. 56, no. 11, pp. 1681–1689, 2015\.
  * [31] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “Tensorflow: A system for large-scale machine learning,” in _12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16)_, 2016, pp. 265–283.
  * [32] F. Chollet _et al._ , “Keras,” https://github.com/fchollet/keras, 2015\.


Integration of first-order ODEs by Jacobi fields A. J. Pan-Collantes
Department of Mathematics, IES San Juan de Dios, Medina Sidonia, Spain; email:
<EMAIL_ADDRESS>J. A. Alvarez-Garcia Department of Mathematics, IES Jorge
Juan, San Fernando, Spain; email<EMAIL_ADDRESS>
###### Abstract

A new class of vector fields enabling the integration of first-order ordinary
differential equations (ODEs) is introduced. These vector fields are not, in
general, Lie point symmetries. The results are based on a relation between
2-dimensional Riemannian manifolds and the integrability of first-order ODEs,
which was established in a previous work of the authors. An integration
procedure is provided, together with several examples to illustrate it. A
connection between integrating factors of first-order ODEs and Schrödinger-
type equations is highlighted.

##  1 Introduction

The classical Lie method for integrating first-order ordinary differential
equations (ODEs) of the form $u^{\prime}(x)=\phi(x,u)$ relies on identifying a
vector field $X$ which is the infinitesimal generator of a non-trivial one-
parameter group of symmetries [1, 2, 3]. These vector fields are characterized
by the following condition:

$[X,A]=\rho A,$ (1)

where $[,]$ denotes the Lie bracket and $A$ is the vector field associated to
the ODE. Once this infinitesimal generator
$X=\xi\partial_{x}+\eta\partial_{u}$ is known, the integration of the ODE is
straightforward, since $(\eta-\xi\phi)^{-1}$ is an integrating factor [4, 5].

In the recent preprint [6], the authors have introduced a theory that connects
surfaces, in the sense of 2-dimensional Riemannian manifolds, with the
integration of first-order ODEs. In particular, it is demonstrated that there
exist vector fields with a clear geometric meaning, and distinct from Lie
point symmetries, which enable the integration of the given ODE. These vector
fields arise as a modification of the notion of Jacobi field, which plays a
key role in the study of Riemannian geometry, particularly in understanding
the behavior of geodesics and its relation with the Gaussian curvature.

In this paper, we introduce the concept of relative Jacobi field in Section 2,
following an overview of the background and notation. Subsequently, in Section
3, we present a procedure for integrating first-order ODEs using a relative
Jacobi field. We explore specific scenarios, such as surfaces with constant
curvature, where the knowledge of a relative Jacobi field is obtained
automatically; and surfaces with curvature dependent solely on the variable
$x$, where the integration can be done in terms of a particular solution to a
Schrödinger-type equation. Finally, in Section 4 we present several
illustrative examples.

##  2 Relative Jacobi fields

In this work, we investigate the local integration of first-order ODEs of the
form

$u^{\prime}(x)=\phi(x,u),$ (2)

where $\phi$ is a smooth function defined on an open set
$U\subseteq\mathbb{R}^{2}$ which can be shrunk as required. Recall that the
integration of (2) is equivalent to the finding of a first integral of the
Pfaffian equation

$-\phi dx+du\equiv 0,$ (3)

i.e., a smooth function $F\in\mathcal{C}^{\infty}(U)$ such that $dF=\mu(-\phi
dx+du)$, for certain non-vanishing function $\mu\in\mathcal{C}^{\infty}(U)$
called an integrating factor for (3).

Following the ideas developed in [6], we define in $U$ a Riemannian metric $g$
expressed in matrix form as

$\begin{pmatrix}1+\phi^{2}&-\phi\\\ -\phi&1\end{pmatrix}.$ (4)

The pair $\mathcal{S}=(U,g)$ defines a 2-dimensional Riemannian manifold,
hereafter referred to as the _surface associated with_ equation (2).

The vector field defined on $U$ given by $A=\partial_{x}+\phi\partial_{u}$ is
frequently called the associated vector field to equation (2), and its
integral curves are in correspondence with solutions to (2). It can be checked
that $g(A,A)=g(\partial_{u},\partial_{u})=1$ and $g(A,\partial_{u})=0$, so the
vector fields $A$ and $\partial_{u}$ constitute an orthonormal frame for the
tangent bundle $TU$ [6, 7].

On the surface $\mathcal{S}$ we consider the corresponding Levi-Civita
connection $\nabla$, which is the unique torsion-free connection defined on
the tangent bundle compatible with $g$ [8, 9, 10]. Given two vector fields
$X,Y$ on $U$, they can be expressed in the frame $(A,\partial_{u})$ as

$\displaystyle X$ $\displaystyle=x_{1}A+x_{2}\partial_{u},$ $\displaystyle Y$
$\displaystyle=y_{1}A+y_{2}\partial_{u},$

where $x_{i},y_{i}\in\mathcal{C}^{\infty}(U)$, $i=1,2$. We have that (see [6]
and references therein)

$\nabla_{X}Y=(X(y_{1})-\phi_{u}x_{2}y_{2})A+(X(y_{2})+\phi_{u}x_{2}y_{1})\partial_{u}.$

It is well-known that an important notion in the study of Riemannian manifolds
is the concept of geodesic. Moreover, the study of geodesic vector fields,
i.e., vector fields whose integral curves are geodesics, have attracted
increasing attention in recent research [11, 12, 13]. They are characterized
as follows.

######  Definition 2.1.

A vector field $X$ will be called a geodesic vector field if

$\nabla_{X}X=0.$

Remarkably, note that in the surface $\mathcal{S}$ associated to equation (2),

$\nabla_{A}A=(A(1)-\phi_{u}\cdot 0\cdot 0)A+(A(0)+\phi_{u}\cdot 0\cdot
1)\partial_{u}=0,$

so $A$ is a geodesic vector field.

The Gaussian curvature of a surface is a pivotal notion, arising in multiple
mathematical and physical contexts. In the case of the surface $\mathcal{S}$
associated to equation (2), Gaussian curvature can be calculated with the
following simple formula [7, 6]:

$\mathcal{K}(x,u)=-\partial_{u}(A(\phi(x,u))).$ (5)

Gaussian curvature and geodesics are related by the notion of Jacobi vector
fields along geodesic curves [10, 14], which represent the infinitesimal
separation between nearby geodesics. They can be thought of as measuring how
geodesics spread out or come together as they move through the space. In [6]
it was introduced the following notion, deeply related to Jacobi vector fields
along geodesics:

######  Definition 2.2.

Given a geodesic vector field $X$, a vector field $J$ will be called a Jacobi
vector field relative to $X$ if the following equation holds

$\nabla_{X}\nabla_{X}J+\mathcal{K}\left(g(X,X)J-g(J,X)X\right)=0.$ (6)

We now introduce two linear operators: given $h\in\mathcal{C}^{\infty}(U)$ we
define:

$\displaystyle\mathfrak{S}(h):=A(h)-\phi_{u}h,$
$\displaystyle\mathfrak{T}(h):=A(h)+\phi_{u}h.$

These operators satisfy the following properties:

######  Proposition 2.1.

Consider the first-order ODE (2) with the corresponding Pfaffian equation (3).
Given a non-vanishing smooth function $h\in\mathcal{C}^{\infty}(U)$:

  1. a)

If $\mathfrak{T}(h)=0$ then $h$ is an integrating factor for (3).

  2. b)

If $\mathfrak{S}(h)=0$ then $h^{-1}$ is an integrating factor for (3).

  3. c)

$(\mathfrak{T}\circ\mathfrak{S})(h)=(A^{2}+\mathcal{K})(h)$, where
$\mathcal{K}$ is the curvature of the associated surface $\mathcal{S}$.

###### Proof.

To prove part a) observe that, by Poincaré Lemma, a smooth function $h$ is an
integrating factor for (3) if and only if

$d(-\phi hdx+hdu)=(\phi_{u}h+\phi h_{u}+h_{x})dx\wedge du=0,$

and therefore if and only if

$\phi_{u}h+\phi h_{u}+h_{x}=(A+\phi_{u})(h)=\mathfrak{T}(h)=0.$

For part b), note that

$d(-\phi h^{-1}dx+h^{-1}du)=\left(\frac{\phi_{u}}{h}-\frac{\phi
h_{u}}{h^{2}}-\frac{h_{x}}{h^{2}}\right)dx\wedge du.$

Then, $h^{-1}$ is an integrating factor if and only if

$\phi_{u}h-\phi h_{u}-h_{x}=-\mathfrak{S}(h)=0.$

Finally, for part c),

$\displaystyle(\mathfrak{T}\circ\mathfrak{S})(h)$
$\displaystyle=\mathfrak{T}(\mathfrak{S}(h))$
$\displaystyle=\mathfrak{T}(A(h)-\phi_{u}h)$
$\displaystyle=A(A(h)-\phi_{u}h)+\phi_{u}(A(h)-\phi_{u}h)$
$\displaystyle=A^{2}(h)-A(\phi_{u}h)+\phi_{u}A(h)-\phi_{u}^{2}h$
$\displaystyle=A^{2}(h)-A(\phi_{u})h-\phi_{u}A(h)+\phi_{u}A(h)-\phi_{u}^{2}h$
$\displaystyle=A^{2}(h)-(A(\phi_{u})+\phi_{u}^{2})h$
$\displaystyle=A^{2}(h)-(\phi_{ux}+\phi\phi_{uu}+\phi_{u}^{2})h$
$\displaystyle=A^{2}(h)-\partial_{u}(\phi_{x}+\phi\phi_{u})h$
$\displaystyle=(A^{2}+\mathcal{K})(h),$

where we have used equation (5). ∎

In the next section we will show the relevance of these operators in the
integration of a first-order ODE like (2), once a Jacobi field relative to $A$
is known.

##  3 Integration of first-order ODEs

Given a first-order ODE of the form (2), the associated vector field $A$ is a
geodesic vector field with $g(A,A)=1$, and therefore a Jacobi field $J$
relative to $A$ satisfies

$\nabla_{A}\nabla_{A}J+\mathcal{K}\left(J-g(J,A)A\right)=0$ (7)

Observe that the vector fields of the form $J=(ax+b)A$, where
$a,b\in\mathbb{R}$, trivially satisfy equation (7), and thus we will refer to
them as trivial Jacobi fields relative to $A$.

Given a non-trivial Jacobi field $J$ relative to $A$, we can express it in the
frame $\\{A,\partial_{u}\\}$ as

$J=\sigma A+\delta\partial_{u},$

for certain smooth functions $\sigma,\delta\in\mathcal{C}^{\infty}(U)$. By
substitution into equation (7) we obtain

$\nabla_{A}\nabla_{A}\left(\sigma
A+\delta\partial_{u}\right)+\mathcal{K}\delta\partial_{u}=0.$

Since $\nabla_{A}(A)=0$ and $\nabla_{A}(\partial_{u})=0$ we have

$A^{2}(\sigma)A+A^{2}(\delta)\partial_{u}+\mathcal{K}\delta\partial_{u}=0,$

and thus,

$\displaystyle A^{2}(\sigma)$ $\displaystyle=0,$ (8) $\displaystyle
A^{2}(\delta)+\mathcal{K}\delta$ $\displaystyle=0.$ (9)

Thus, to integrate equation (2) using a non-trivial Jacobi field $J$ relative
to $A$, one may follow the following steps:

  1. 1.

Decompose the vector field $J$ as

$J=\sigma A+\delta\partial_{u},$

for certain smooth functions $\sigma,\delta\in\mathcal{C}^{\infty}(U)$.

  2. 2.

Define $F:=A(\sigma)$, which according to equation (8) is a first integral of
the vector field $A$, and therefore a first integral of (3). If $F$ is not a
constant function, then equation (2) is solved.

  3. 3.

On the other hand, if $A(\sigma)=a$, where $a\in\mathbb{R}$, then define the
function $G:=\sigma-ax$, which constitutes a first integral of equation (3),
since

$A(G)=A(\sigma-ax)=A(\sigma)-a=0.$

If $G$ is not a constant function, the first-order ODE (2) is solved. On the
other hand, in the case where $G=\sigma-ax=b$ for some $b\in\mathbb{R}$, we
have $\sigma=ax+b$. Consequently, it is necessary that $\delta\neq 0$, as
$\delta=0$ would imply that $J$ is a trivial Jacobi field relative to $A$.

  4. 4.

Assuming $\delta\neq 0$, and considering equation (9) and part c) in
Proposition 2.1, the following relation holds

$A^{2}(\delta)+\mathcal{K}\delta=(\mathfrak{T}\circ\mathfrak{S})(\delta)=0.$
(10)

There are two possibilities:

    1. (a)

If $\mathfrak{S}(\delta)=0$, then, as indicated in part b) of Proposition 2.1,
$\delta^{-1}$ serves as an integrating factor for (3). Compute a smooth
function $I=I(x,u)$ satisfying:

$dI=\delta^{-1}(-\phi dx+du)$

via quadrature, and equation (2) is solved.

    2. (b)

If $\mathfrak{S}(\delta)\neq 0$, part a) of Proposition 2.1 permits the use of
$\mathfrak{S}(\delta)$ as an integrating factor for (3), since
$\mathfrak{T}(\mathfrak{S}(\delta))=0$ because of (10). Then, the 1-form

$\mathfrak{S}(\delta)(-\phi dx+du)$

is closed, and a primitive $I=I(x,u)$ provides the solution to equation (2).

In summary, _the knowledge of a non-trivial Jacobi field relative to $A$_
leads to the integration of the first-order ODE (2).

###  3.1 Constant curvature surfaces

If the curvature of the associated surface is constant, a non-trivial Jacobi
field can always be found, and the ODE can be integrated. Suppose
$\mathcal{K}=k$, with $k\in\mathbb{R}$. We can assume a Jacobi field relative
to $A$ of the form $J=\delta(x)\partial_{u}$. In this case, equation (7) takes
the particular form

$\nabla_{A}\nabla_{A}\delta(x)\partial_{u}+k\delta(x)\partial_{u}=0,$

and since $\nabla_{A}\partial_{u}=0$, the function $\delta$ must satisfy the
second-order ODE

$\ddot{\delta}(x)+k\delta(x)=0.$

The general solution to this ODE is given by

$\delta(x)=\begin{cases}A\cos{\sqrt{k}x}+B\sin{\sqrt{k}x},&\text{ if }k>0,\\\
A+Bx,&\text{ if }k=0,\\\ A\cosh{\sqrt{-k}x}+B\sinh{\sqrt{-k}x},&\text{ if
}k<0,\end{cases}$ (11)

with $A,B\in\mathbb{R}$, so we can choose a particular solution and proceed as
in Step 4 of the procedure described in Section 3.

###  3.2 Relation to Schrödinger equation

In the case where the curvature of the associated surface $\mathcal{K}=k(x)$
depends only on the variable $x$, a Jacobi field relative to $A$ could be
found by solving the following Schrödinger-type equation:

$\ddot{\delta}(x)+k(x)\delta(x)=0.$ (12)

Therefore, according to our procedure described in Section 3, an integrating
factor for equation (2) could be expressed in terms of a particular solution
of equation (12). Recall that some specific values of $k(x)$ lead to solutions
$\delta(x)$ expressed in terms of well-known special functions like Bessel,
Mathieu, Airy or hypergeometric functions. In these cases, finding the
integrating factor $\mathfrak{T}(\delta)$ provided by our method through a
direct approach appears to be challenging.

Observe that in the general case, when $\mathcal{K}=\mathcal{K}(x,u)$ depends
on both $x$ and $u$, our approach provides a relationship between integrating
factors of first-order ODEs and the partial differential equation for
$\delta=\delta(x,u)$:

$A^{2}(\delta)+\mathcal{K}\delta=0,$ (13)

which could be interpreted as a generalized Schrödinger equation.

##  4 Examples

We will show in this section several examples to illustrate the previous
results.

######  Example 4.1.

Consider the first-order ODE given by

$u^{\prime}(x)=\frac{1}{u},$ (14)

whose associated vector field is $A=\partial_{x}+\frac{1}{u}\partial_{u}$. The
curvature corresponding to the surface associated to equation (14) is,
according to equation (5),

$\mathcal{K}=-\partial_{u}\left(A\left(\frac{1}{u}\right)\right)=-\frac{3}{u^{4}}.$

On the other hand, it can be checked that the vector field

$J=x\partial_{x}+\frac{x(u^{2}-x+1)}{u}\partial_{u}$

is a Jacobi field relative to $A$, since it satisfies equation (7). We will
use $J$ to integrate equation (14), following the Jacobi field-based method
developed in Section 3.

First, $J$ is expressed in the frame $\\{A,\partial_{u}\\}$ as

$J=xA+\frac{x(u^{2}-x)}{u}\partial_{u}.$

Secondly, we note that the first integral presented in Step 2, $F:=A(x)=1$, is
constant, as is the first integral described in Step 3, $G:=x-x=0$.

Thus, we proceed with Step 4, where $\delta=\frac{x(u^{2}-x)}{u}\neq 0$
satisfies

$(\mathfrak{T}\circ\mathfrak{S})\left(\frac{x(u^{2}-x)}{u}\right)=0.$

Then, $\mathfrak{S}\left(\frac{x(u^{2}-x)}{u}\right)=u\neq 0$ must serve as an
integrating factor for the Pfaffian equation ${-\frac{1}{u}dx+du}\equiv 0$.
This is indeed the case, and it enables us to obtain, by quadrature, the first
integral

$I(x,u)=-x+u^{2}.$

The general solution to (14) is given by $-x+u^{2}=C$, where $C\in\mathbb{R}$.

It is interesting to observe that, according to the characterization given in
equation (1), $J$ is not a symmetry of equation (14), since:

$[J,A]=-\partial_{x}-\frac{u^{2}+1}{u}\partial_{u}.$

######  Example 4.2.

In this example, we consider the first-order ODE

$u_{1}=e^{x}+\sqrt{2ue^{x}-e^{2x}-u^{2}+1},$ (15)

with associated vector field given by

$A=\partial_{x}+\left(e^{x}+\sqrt{2ue^{x}-e^{2x}-u^{2}+1}\right)\partial_{u}.$

The curvature of the associated surface $\mathcal{S}$ is given by

$\mathcal{K}=-\partial_{u}\left(A\left(e^{x}+\sqrt{2ue^{x}-e^{2x}-u^{2}+1}\right)\right)=1.$

Since the curvature is constant and positive, according to equation (11) in
Section 3.1, a Jacobi field relative to $A$ is given by

$J=\cos(x)\partial_{u}.$

It can be checked that this vector field is not a Lie point symmetry, by
performing the corresponding commutation relation $[J,A]$.

Next, to integrate equation (15) by means of $J$ with the procedure described
in Section 3, we start at Step 4, since $\sigma=0$. In this case,
$\delta=cos(x)$ and

$\mu:=\mathfrak{S}\left(\cos(x)\right)=-\sin(x)+\frac{cos(x)(u-e^{x})}{\sqrt{2ue^{x}-e^{2x}-u^{2}+1}}\neq
0.$

Then, $\mu$ is an integrating factor for the Pfaffian equation

$-\left(e^{x}+\sqrt{2ue^{x}-e^{2x}-u^{2}+1}\right)dx+du\equiv 0.$

A first integral can therefore be obtained by quadrature,

$I(x,u)=-\arcsin(e^{x}-u)-x,$

so the general solution of equation (15) is given in explicit form by

$u(x)=e^{x}+\sin(x+C),\quad C\in\mathbb{R}.$

######  Example 4.3.

Consider the class of first-order linear ODEs

$u^{\prime}(x)=q(x)u+p(x),$ (16)

where $p,q$ are smooth real functions. The vector field associated with
equation (16) is given by

$A=\partial_{x}+(q(x)u+p(x))\partial_{u}.$

Note that the curvature of the associated surface depends only on $x$:

$\mathcal{K}=-\partial_{u}(A(\phi))=-q^{\prime}(x)-q(x)^{2}.$

Then, to find a Jacobi field relative to $A$ we consider a particular solution
of the following Schrödinger equation

$\ddot{\delta}(x)-\left(q^{\prime}(x)+q(x)^{2}\right)\delta(x)=0,$ (17)

which allows us to construct the vector field $J=\delta(x)\partial_{u}$, which
is orthogonal to $A$.

The reader can check that the solutions to equation (17) are given by

$\delta(x)=C_{1}\left(H(x)+C_{2}\right)e^{Q(x)},$

where $Q$ is a smooth function satisfying $Q^{\prime}(x)=q(x)$, $H$ is a
smooth function satisfying $H^{\prime}(x)=e^{-2Q(x)}$, and
$C_{1},C_{2}\in\mathbb{R}$.

The choice $C_{1}=1,C_{2}=0$ provides the following Jacobi field relative to
$A$:

$J=H(x)e^{Q(x)}\partial_{u},$ (18)

which can be used to integrate equation (16). We have $\sigma=0$ and
$\delta=H(x)e^{Q(x)}$, so we continue in Step 4. We have that

$\mathfrak{S}\left(H(x)e^{Q(x)}\right)=A\left(H(x)e^{Q(x)}\right)-\phi_{u}H(x)e^{Q(x)}=e^{-Q(x)}\neq
0,$

and therefore, $\mu:=e^{-Q(x)}$ is an integrating factor for the Pfaffian
equation

$-(q(x)u+p(x))dx+du\equiv 0.$

Indeed, the general solution to equation (16) can be expressed as

$-P(x)+ue^{-Q(x)}=C,$

where $C\in\mathbb{R}$ and $P$ is a smooth function satisfying
$P^{\prime}(x)=p(x)e^{-Q(x)}$.

## References

  * [1] N. H. Ibragimov.  A Practical Course in Differential Equations and Mathematical Modelling: Classical and New Methods, Nonlinear Mathematical Models, Symmetry and Invariance Principles.  World Scientific, Beijing, 2010.
  * [2] H. Stephani.  Differential Equations: Their Solutions Using Symmetry.  Cambridge University Press, New York, 1989.
  * [3] P. J. Olver.  Applications of Lie groups to differential equations, volume 107\.  Springer-Verlag, New York, 1986.
  * [4] J. Sherring and G. Prince.  Geometric aspects of reduction of order.  Trans. Amer. Math. Soc., 334(1):433–453, 1992.
  * [5] A. J. Pan-Collantes, A. Ruiz, C. Muriel, and J. L. Romero.  $\mathcal{C}^{\infty}$-structures in the integration of involutive distributions.  Phys. Scr., 98(8):085222, jul 2023.
  * [6] A. J. Pan-Collantes and J. A. Álvarez García.  Surfaces associated with first-order ODEs, 2023.  DOI: 10.48550/arXiv.2312.04489.
  * [7] Z. O. Bayrakdar and T. Bayrakdar.  Burgers’ Equations in the Riemannian Geometry Associated with First-Order Differential Equations.  Advances in Mathematical Physics, 2018:1–8, 2018.
  * [8] W. Chen, S. S. Chern, and K. S. Lam.  Lectures on differential geometry, volume 1.  World Scientific Publishing Company, 1999.
  * [9] S. Morita.  Geometry of Differential Forms.  American Mathematical Society, Rhode Island, 2001.
  * [10] M. P. Do Carmo.  Riemannian geometry, volume 6.  Springer, 1992.
  * [11] J. F. Cariñena and M. C. Muñoz-Lecanda.  Geodesic and Newtonian vector fields and symmetries of mechanical systems.  Symmetry, 15(1):181, jan 2023.
  * [12] S. Deshmukh and V. Azam Khan.  Geodesic vector fields and Eikonal equation on a Riemannian manifold.  Indagationes Mathematicae, 30(4):542–552, 2019.
  * [13] S. Deshmukh, P. Peska, and N. Bin Turki.  Geodesic vector fields on a riemannian manifold.  Mathematics, 8(1):137, 2020.
  * [14] J. M. Lee.  Introduction to Riemannian manifolds, volume 2.  Springer, 2018.