# Quantum Radon Transform and Its Application

Guangsheng Ma, Hongbo Li, and Jiman Zhao G. Ma, corresponding author, is with
Academy of Mathematics and Systems Science, Chinese Academy of Sciences,
Beijing, 100875, China. e-mail<EMAIL_ADDRESS>Li is with Academy of
Mathematics and Systems Science, Chinese Academy of Sciences; University of
Chinese Academy of Sciences, Beijing, 100190, China. e-mail:
<EMAIL_ADDRESS>Zhao is with School of Mathematical Sciences, Beijing
Normal University, Beijing, 100875, China. e-mail<EMAIL_ADDRESS>
###### Abstract

This paper extends the Radon transform, a classical image processing tool for
fast tomography and denoising, to the quantum computing platform. A new kind
of periodic discrete Radon transform (PDRT), called quantum Radon transform
(QRT), is proposed. The QRT has a quantum implementation that is exponentially
faster than the classical Radon transform. Based on the QRT, we design an
efficient quantum image denoising algorithm. The simulation results show that
QRT preserves the good denoising capability as in the classical PDRT. Also, a
quantum algorithm for interpolation-based discrete Radon transform (IDRT) is
proposed, which can be used for fast line detection. Both the quantum
extension of IDRT and the line detection algorithm can provide polynomial
speedups over the classical counterparts.

###### Index Terms:

Radon transform, quantum computation.

##  I Introduction

Radon transform, proposed by Johann Radon in 1917 [1], is an important image
processing tool with widespread applications in computed tomography,
geophysics, and remote sensing, etc. [2]. It changes a function $f$ defined on
the plane to a function $R_{\theta}f(\rho)$ defined on the space of lines in
the plane, whose value on the line with interception $\rho$ and slope $\theta$
equals the integral of function $f$ along the line:

$\displaystyle
R_{\theta}f(\rho):=\iint_{\mathbb{R}^{2}}f(x,y)\delta(\rho-x\sin\theta+y\cos\theta)dxdy,$
(1)

where $\theta\in(0,\pi]$, $\rho\in\mathbb{R}$, and $\delta$ is the Dirac
function. By definition, Radon transform possesses the capability of detecting
singularities along straight lines, and performs better at denoising images
with linear singularities than other image processing tools [2, 3].

To implement Radon transform, discretization is necessary. However, the
different discretization methods will result in different discrete Radon
transforms (DRTs) that have different applications. In discretization, there
are two methods to approximate the line integral in (1): the interpolation
method, and the periodic discrete grid method.

The interpolation method (IDRT) evaluates the integral along a straight line
by making interpolation among the adjacent points on the discrete image grid
of the line. The earlier DRTs are based on this method [1]. In 1987, Beylkin
discovered an exact inversion formula, and proved that if the discrete version
is based on Radon’s original formula, then the reconstruction can only be
approximate [4]. The IDRT can be used in line detection, X-ray computed
tomography [5], etc. Performing the IDRT on an $N\times N$ image often
requires at least $\Omega(N^{3})$ arithmetic operations [6].

The periodic discrete grid method (PDRT) calculates the integrals along a set
of warped lines, and does not have direct connection with the continuous Radon
transform [6, 7, 8]. Still the method possesses some very nice properties,
such as the exact reconstruction property, Fourier slice property, etc. Matus
and Flusser [9] first investigated PDRT on $\mathbb{Z}^{2}_{p}$, where $p$ is
prime. Then Hsung et al. [10] extended PDRT to $\mathbb{Z}^{2}_{p^{n}}$. The
PDRT has been used in image denoising [3], tomographic reconstruction [11],
image watermarking and encryption [12], etc. To compute the PDRT of an
$N\times N$ image, $\Omega(N^{2}\text{log}N)$ arithmetic operations are
required [13].

With the emergence of quantum computing, for many important computational
problems, it is found that quantum algorithms can provide dramatic speedup
[14], for example, exponentially fast quantum algorithms such as quantum
Fourier algorithm, Shor’s factoring algorithm [15]; polynomially fast quantum
algorithms such as Grover’s search algorithm [16], and so on [17].

To speed up image processing in the forthcoming quantum computing age, quantum
image processing (QIMP), whose topics range from quantum image representations
to image processing, has drawn a lot of attention in the last decade [18, 19,
20]. Early work on QIMP concentrated on the quantum representation of images
e.g., [21, 22, 23]. Roughly speaking, there are two typical representation
methods: (1) amplitude representation method, such as the Real Ket
representation [24] which utilizes the qubits’ amplitude and a computational
basis state to encode the grayscale and the location of a pixel respectively;
(2) basis representation method, such as the novel enhanced quantum
representation (NEQR) [25] which utilizes a register’s computational basis
state to encode both the grayscale and the location of pixels.

Another focus in QIMP is to develop image processing tools in the quantum
computation framework [26, 27], for example, Yao [28] proposed an efficient
quantum image edge detection algorithm for Real Ket images. Image
segmentation, watermarking, scrambling, and some other image processing
problems have also been investigated in QIMP [29, 30, 31].

However, there is still no extension of Radon transform to the field of
quantum image processing. In this paper, we make such extension.

We first extend PDRT to quantum computation framework. The main difficulty in
making the extension comes from the fact that there is no preference of using
unitary transforms in designing classical algorithms; conversely, quantum
algorithms mostly use unitary operators. To realize the classical ‘non-
unitary’ PDRT with unitary transformations, we design a quantum reversible
multiplication, and then use Fourier slice property of PDRT with replacing the
traditional multiplication by the quantum reversible multiplication. Finally,
we obtain a transform that is different from any existing PDRTs. We call this
new transform quantum Radon transform (QRT).

The QRT preserves many good properties as in the classical PDRT, such as
Fourier slice property. By replicating the denoising experiments designed
specially for testing PDRT [3], it is shown that the QRT is of the good
denoising capability as in the classical PDRT (cf. Fig. 8-7). The most
important advantage of QRT is that for an $N\times N$ image, the QRT can be
implemented in time $O(\log^{3}N)$, which runs exponentially faster than the
classical PDRT, which has runtime $\Omega(N^{2}\log N)$. As the application, a
QRT-based quantum denoising method is proposed, which runs exponentially
faster than the classical PDRT denoising method.

We then extend IDRT to the quantum case. The quantum extension of IDRT can
provide polynomial speedup, with the input being a Real ket quantum image and
the output being encoded with the NEQR representation. As the application, a
quantum algorithm for line detection using IDRT is given, which enables the
line detection process to be speeded up polynomially in the average case.

This paper is arranged as follows. In Section II, we introduce some background
on Radon transform and classical/quantum image processing. In Section III, we
present a reversible quantum multiplication. In Section IV, we introduce QRT
and explore some basic properties of it. In Section V, we extend
interpolation-based DRT to the quantum case. In Section VI, we present two
applications of our proposed quantum transforms.

##  II Preliminaries

###  II-A Periodic Discrete Radon Transform

Throughout this paper, we use $p$ to denote a prime number, use $n$ to denote
a positive integer, and use $N$ to denote a power of $2$. $[n]$ is the subset
of integers $\\{0,1,...,n-1\\}$. We use $I$ to denote the imaginary unit. The
$L^{2}$-norm of vector $\vec{a}=(a_{0},...,a_{n-1})$ is
$||\vec{a}||_{2}=\sqrt{\sum_{i\in[n]}|a_{i}|^{2}}$.

We begin with a specific kind of discrete ‘line’, an example of which is given
in Fig. 1 :

Discrete line $L^{n}_{l,k}$. The discrete line on lattice $\mathbb{Z}^{2}_{n}$
with interception $l$ and slope $k$ is

$L^{n}_{l,k}=\left\\{\begin{aligned} &\\{(x,y)\ \big{|}x+ky=l\mod l,\
x,y\in[n]\\},\text{if}\ k\in[n];\\\ &\\{(x,l)\
\big{|}x\in[n]\\},\quad\text{if}\quad k=n.\end{aligned}\right.$ (2)

Figure 1: Example of a discrete line $L^{11}_{1,3}=\\{(x,y)|x-3y=1\text{ mod
}11\\}$ on lattice $\mathbb{Z}^{2}_{11}$. The hollow circles indicate the
lattice points over the discrete line $L^{11}_{1,3}$.

The periodic discrete Radon transform (PDRT) – sometimes called finite Radon
transform – is defined as summations of function values at points over these
discrete lines.

######  Definition 1 (Periodic discrete Radon transform [9])

The PDRT of a function $f$ defined on lattice $\mathbb{Z}^{2}_{n}$ is

$\displaystyle r_{k}(l)=\frac{1}{\sqrt{n}}\displaystyle\sum_{(i,j)\in
L^{n}_{l,k}}f(i,j),\quad l\in[n],\ k\in[n+1].$ (3)

The following proposition gives an important property of PDRT:

######  Proposition 1 (Fourier slice property of PDRT [3])

Let $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ be the $1$-D and $2$-D discrete
Fourier transform, respectively. For a function $f$ defined on
$\mathbb{Z}^{2}_{n}$, let $r_{k}(l)$ be the PDRT of $f$, so that it is defined
on $\mathbb{Z}_{n}$. Then for any $\omega,k\in[n]$,

$\displaystyle\mathcal{F}_{1}\\{r_{k}\\}(\omega)=\mathcal{F}_{2}\\{f\\}(\omega,k\omega\emph{
mod }n);$ (4)

for $k=n$,

$\displaystyle\mathcal{F}_{1}\\{r_{n}\\}(\omega)=\mathcal{F}_{2}\\{f\\}(0,\omega).$
(5)

###### Proof:

For any $k\in[n]$, the discrete Fourier transform of function $r_{k}(l)$ in
variable $l$ is

$\displaystyle\mathcal{F}_{1}\\{r_{k}\\}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{n}\displaystyle\sum_{l,y\in[n]}f(l-ky,y)e^{-2\pi
I\frac{l\omega}{n}}$ (6) $\displaystyle=$
$\displaystyle\frac{1}{n}\displaystyle\sum_{l,y\in[n]}f(l-ky,y)e^{-2\pi
I\frac{(l-ky)\omega+yk\omega}{n}}$ $\displaystyle=$
$\displaystyle\mathcal{F}_{2}\\{f\\}(\omega,k\omega\text{ mod }n).$

For $k=n$,

$\displaystyle\mathcal{F}_{1}\\{r_{n}\\}(\omega)$
$\displaystyle=\frac{1}{n}\displaystyle\sum_{x,l\in[n]}f(x,l)e^{-2\pi
I\frac{l\omega}{n}}=\mathcal{F}_{2}\\{f\\}(0,\omega).$ (7)

∎

This Fourier slice property provides a fast implementation of PDRT. For an
$n\times n$ image $f$, each value of its PDRT can be computed in time O$(n)$,
so the whole PDRT can be obtained in time O$(n^{3})$ if directly computed by
definition. On the other hand, the Fourier slice property of PDRT allows one
to compute the PDRT in time O$(n^{2}\text{log}n)$: since the $1$-D and $2$-D
(inverse) Fourier transforms can be implemented in time O$(n\text{log}n)$ and
O$(n^{2}\text{log}n)$ [32], and the PDRT of $f$ can be obtained by performing
$1$-D inverse Fourier transform on $\mathcal{F}_{2}\\{f\\}$ according to (4),
the computational complexity of the PDRT is thus reduced to
O$(n^{2}\text{log}n)$.

Besides designing fast algorithms, another important practical issue is to
recover the original image from its Radon transform. In [9],
Mat$\acute{u}\breve{s}$ and Flusser proposed the following reconstruction
formula for PDRT on $\mathbb{Z}^{2}_{p}$:

######  Proposition 2 (Reconstruction formula for PDRT on
$\mathbb{Z}^{2}_{p}$)

Let $r_{k}(l)$ be the PDRT of a function $f$ defined on $\mathbb{Z}^{2}_{p}$.
Then for any $i,j\in[p]$,

$\displaystyle
f(i,j)=\frac{1}{\sqrt{p}}\sum_{\begin{subarray}{c}\\{(l,k)|(i,j)\in
L^{p}_{l,k},\\\ l\in[p],\
k\in[p+1]\\}\end{subarray}}r_{k}(l)-\frac{1}{p}\displaystyle\sum_{x,y\in[p]}f(x,y).$
(8)

###### Proof:

The following geometric properties of discrete lines are easy to verify:

  * 1.

Every discrete line $L^{p}_{l,k}$ contains $p$ lattice points, and two
parallel discrete lines have no point of intersection;

  * 2.

For any fixed slope $k$, the $p$ parallel lines $L^{p}_{l,k}$ (where
$l\in[p]$) provide a complete cover of the lattice $\mathbb{Z}^{2}_{p}$;

  * 3.

Two discrete lines $L^{p}_{l,k}$ of different slopes will interact in exactly
one point.

For all the $p+1$ lines through a fixed point $(i,j)$, by term 1, every two of
them have only one point in common, which is just $(i,j)$. Since there are
$p(p+1)$ points on the lines, and there are $p+1$ copies of point $(i,j)$ on
such lines, there are all together $p^{2}$ different points on these lines,
which are exactly the total number of points in lattice $\mathbb{Z}_{p}^{2}$:

$\displaystyle\displaystyle\sum_{x,y\in[p]}f(x,y)+$ $\displaystyle
pf(i,j)=\sqrt{p}\sum_{\begin{subarray}{c}\\{(l,k)|(i,j)\in L^{p}_{l,k},\\\
l\in[p],\ k\in[p+1]\\}\end{subarray}}r_{k}(l).$ (9)

∎

In general, the size of an image is not the square of a prime number. There is
a series of work to extend PDRT to images of more general sizes [10, 33, 9].
For one example, Kingston [34] extends PDRT to images of size $p^{n}\times
p^{n}$. The reconstruction formula for such images are much more complicated,
e.g., (11) in [34].

###  II-B Classical Image Denoising

Image denoising is to remove noise from a noisy image, so as to restore the
true image [35, 36]. Suppose we are given a real-valued noisy signal

$\displaystyle h_{i}=f_{i}+e_{i},\quad i\in[n],$ (10)

where $\vec{f}=(f_{0},f_{1},...,f_{n-1})$ is original signal, and $e_{i}$ is
the noise sampled independently from the normal distribution
$V(0,\sigma^{2})$, where $0$ is the mean and $\sigma^{2}$ is the variance.
After performing some denoising method on $\vec{h}=(h_{0},h_{1},...,h_{n-1})$
to prepare the proceed signal $\vec{h^{\prime}}$, if the following noise level
decreases (i.e.,
$\text{Risk}(\vec{h},\vec{f})>\text{Risk}(\vec{h^{\prime}},\vec{f})$):

$\displaystyle\text{Risk}(\vec{h},\vec{f})=\frac{1}{n}E(||\vec{h}-\vec{f}\
||^{2}_{2}),$ (11)

where $E(\cdot)$ is the expectation, then we say such denoising method is
effective.

Fig. 2 shows a general procedure for denoising signal using discrete wavelet
transform (DWT). Below, we explain why DWT denoising method is effective. Let
$\vec{f}=(f_{i})$ be a pure signal where $f_{i}\equiv 1$, $i\in[n]$. Let
$\vec{h}=(h_{i})$ be the noisy signal where $h_{i}=1+e_{i}$, $e_{i}\sim
V(0,\sigma^{2})$. After denoising $\vec{h}$ using Haar wavelet and threshold
$T=\infty$ (i.e, change all wavelet coefficients to $0$), by Fig. 2, each
element of denoised signal is of the form $h^{\prime}_{i}=1+e^{\prime}_{i}$,
where

$\displaystyle
e^{\prime}_{i}\sim\frac{1}{2}V(0,\sigma^{2})+\frac{1}{2}V(0,\sigma^{2})=V(0,\sigma^{2}/2).$

Now that
$\text{Risk}(\vec{h^{\prime}},\vec{f})=\frac{\sigma^{2}}{2}<\frac{\sigma^{2}}{n}=\text{Risk}(\vec{h},\vec{f})$,
the noise level decreases, and thus the Haar denoising method works.

Figure 2: Denoising using DWT. There are three steps in total: (1) Apply the
DWT to the input noisy signal $f$, i.e., multiply $f$ by the matrix form of
some wavelet (the one used here is Haar wavelet). The resulting coefficients
can be divided into two parts: the scaling coefficients $c$, composed of the
odd rows, and the wavelet coefficients $d$, composed of even rows. (2) Apply
Hard-thresholding to wavelet coefficients $d$ with threshold $T$ to obtain new
wavelet coefficients $d^{\prime}$ (3) Apply inverse DWT to thresholded
coefficients, where in this case the inverse of the Haar wavelet matrix is
itself.

However, for general denoising algorithms, it is often hard to make such a
statistical analysis as the above. Experimentally, the following so-called
signal-noise-ratio (SNR) is often taken to measure denoising performance:

$\displaystyle\text{SNR}(\vec{h},\vec{f})=10\log_{10}(\frac{||\vec{h}||_{2}^{2}}{||\vec{h}-\vec{f}||^{2}_{2}}),$
(12)

where $f$ and $h$ are some specific test signals. In comparison with the
probability value in (11), SNR is more easily accessible by numerical
experiments. A denoising method is said to be effective if
$\text{SNR}(\vec{h^{\prime}},\vec{f})>\text{SNR}(\vec{h},\vec{f})$.

Now, we consider image denoising, i.e., $2$-dimensional signal denoising. The
simplest $2$-D DWT denoising method is to apply $1$-D denoising method to each
row in the signal matrix, then use $1$-D denoising method to each column.
However, image possesses various geometric features. By the work of Do and
Vetterli [3], the following PDRT denoising method are more effective than the
general $2$-D DWT in denoising images with obvious singularities along
straight lines, as shown in Fig. 8.

  * 1)

Apply the PDRT to the noisy image to obtain its PDRT $r_{k}(l)$,

  * 2)

For each slope $k$, perform once 1-D DWT denoising method on the PDRT
$r_{k}(l)$ along the direction of interception $l$.

  * 3)

Perform inverse PDRT by (9).

In [3], the authors explain why PDRT is better: by PDRT, the typical linear
singularities of pure image are represented by a few wavelet coefficients (in
step 2) while randomly located noisy singularities are unlikely to produce
significant coefficients. This is unlike using the $2$-D DWT where both noisy
pixels and image singularities can produce significant wavelet coefficients.
Therefore one can remove the noise with less damage to the original image by
properly thresholding wavelet coefficients.

A routine computation shows that the time complexity of donising an $n\times
n$ image using $2$-D DWT is O$(n^{2})$. The time complexity of PDRT denoising
method is O$(n^{3})$, because performing inverse PDRT has a time complexity
O$(n^{3})$ by (9).

###  II-C Quantum Image Representation and Preparation

######  Definition 2 (Real ket representation [24] )

Let $f$ be an $N\times N$ image, where $N=2^{n}$, and $f(i,j)\geq 0$ is the
image intensity at point $(i,j)$, then its Real Ket representation is

$\displaystyle\frac{1}{\sqrt{\displaystyle\sum_{i,j\in[N]}|f(i,j)|^{2}}}\sum_{i,j\in[N]}f(i,j)|i\rangle|j\rangle,$
(13)

######  Definition 3 (Novel enhanced quantum representation [25] )

The NEQR representation of image $f$ is

$\displaystyle\frac{1}{N}\sum_{i,j\in[N]}|f(i,j)\rangle|i\rangle|j\rangle.$
(14)

In this paper, we use the Real Ket representation of images in the quantum
computation framework. The famous quantum Fourier transform [17] and HHL
algorithm [37] are both based on this amplitude encoding method. Unless
otherwise specified, the term ‘quantum image’ refers to an image encoded with
the Real Ket method.

To show the connection between the above two representation methods, we need
an efficient operation called _conditional rotation_ [37, 38]:

######  Proposition 3 (Conditional rotation)

Let $a$ be the $k$-bit finite precision representation of a positive number
that is smaller than 1. Then the following mapping can be performed in time O
_(_ $k^{2.5}$_)_ :

$\displaystyle|a\rangle|0\rangle\rightarrow|a\rangle(a|0\rangle+\sqrt{1-a^{2}}|1\rangle).$
(15)

###### Proof:

Let $a:=\displaystyle\sum_{i\in[k]}2^{-i-1}a_{i}$, where $a_{i}\in\\{0,1\\}$.
Then $|a_{i}\rangle$ is the $(i+1)$-st qubit of
$|a\rangle=|a_{k-1}\rangle\cdots|a_{1}\rangle|a_{0}\rangle$ counted from
right. For $s>0$, define 1-qubit quantum gate

$R_{s}=\left[\begin{array}[]{cc}\cos s,&-\sin s\\\ \sin s,&\cos
s\end{array}\right].$ (16)

Then it holds that

$\prod_{i\in[k]}R^{\alpha_{i}}_{2^{-(i+1)}}=\prod_{i\in[k]}R_{2^{-(i+1)}\alpha_{i}}=R_{\alpha}.$
(17)

By performing the following $k$ successive 1-bit conditional rotations: for
each $i\in[k]$ the corresponding rotation is $R^{\alpha_{i}}_{2^{-(i+1)}}$
where $\alpha_{j}\in\\{0,1\\}$ is the control bit, on the last qubit of
$|a\rangle|0\rangle$, one gets

$|a\rangle|0\rangle\rightarrow|a\rangle(\cos(a)|0\rangle+\sin(a)|1\rangle).$
(18)

By Lemma $48$ in [39], given $|a\rangle$ where $a\in(0,1)$, using Taylor
series approximation allows to prepare the state $|\arccos(a)\rangle$ in time
O$(k^{2.5})$, where $\arccos(a)\in[0,\frac{\pi}{2})$. So with the help of
ancilla qubits and (18), the following sequence of mappings can be implemented
in time O$(k^{2.5})$:

$\displaystyle|a\rangle|0\rangle|0\rangle$
$\displaystyle\longrightarrow|a\rangle|\arccos(a)\rangle|0\rangle$
$\displaystyle\longrightarrow|a\rangle|\arccos(a)\rangle(a|0\rangle+\sqrt{1-a^{2}}|1\rangle)$
$\displaystyle\longrightarrow|a\rangle|0\rangle(a|0\rangle+\sqrt{1-a^{2}}|1\rangle).$

∎

By the following proposition, a Real Ket image can be prepared from its NEQR
version, and this preparation procedure is efficient if the condition number
$\kappa:=\frac{\min_{i,j\in[N]}|f(i,j)|}{\max_{i,j\in[N]}|f(i,j)|}$ is large
enough.

###### Proposition 4

Let $\vec{a}$ be a real vector realized by a unitary operator
$U:|i\rangle|0\rangle\rightarrow|i\rangle|a_{i}\rangle$, $i\in[N]$ where
$a_{i}$ is the $m$-bit finite precision representation of the vector entries,
and $U$ can be performed in time O$(T_{U})$. Given $a_{max}=\max_{j}|a_{j}|$,
the state $|\vec{a}\rangle$ can be prepared in time
O$(\frac{T_{U}+\emph{poly}(m,\log N))}{\kappa^{2}})$ by the following mapping:

$\displaystyle|0\rangle\rightarrow\cos\theta|\vec{a}\rangle|0\rangle+\sin\theta|\vec{a}^{\perp}\rangle|1\rangle,$
(19)

where

1)
$|\vec{a}\rangle:=\displaystyle\sum_{i\in[N]}\frac{a_{i}}{||\vec{a}||_{2}}|i\rangle$;

2) $\vec{a}^{\perp}$ denotes the vector of entries
$a^{\perp}_{i}:=\sqrt{1-|\frac{a_{i}}{a_{max}}|^{2}}$ for $i\in[N]$;

3) the coefficient
$\cos\theta=\sqrt{\frac{||\vec{a}||^{2}_{2}}{Na^{2}_{\max}}}\geq\frac{\min\limits_{i\in[N]}|a_{i}|}{\max\limits_{i\in[N]}|a_{i}|}:=\kappa$,
and $\kappa$ is called the condition number of vector $|a\rangle$.

###### Proof:

With the number $a_{max}$ and unitary operator $U$ at hand, we can implement
the following transform [37]:

$\displaystyle|i\rangle|a_{i}\rangle\rightarrow|i\rangle|\frac{a_{i}}{a_{max}}\rangle.$
(20)

We first perform the this transform on the input state
$\displaystyle\sum_{i\in[N]}\frac{1}{\sqrt{N}}|i\rangle|a_{i}\rangle$ and then
use the rotation conditioned on $|\frac{a_{i}}{a_{max}}\rangle$ (cf.
Proposition 3):

$\displaystyle|i\rangle|\frac{a_{i}}{a_{max}}\rangle|0\rangle\rightarrow|i\rangle|\frac{a_{i}}{a_{max}}\rangle(\frac{a_{i}}{a_{max}}|0\rangle+\sqrt{1-(\frac{a_{i}}{a_{max}})^{2}}|1\rangle),$
(21)

and finally undo $U$ to clear the second register. After the above operations,
we measure the last qubit.

The running time of the above procedure is O$(T_{U}+\text{poly}(m,\log N))$.
The expected result $|0\rangle$ by measurement has probability
$\frac{\sum_{i\in[N]}a^{2}_{i}}{Na^{2}_{max}}\geq\frac{a^{2}_{min}}{a^{2}_{max}}=\kappa^{2}$,
which indicates that we have successfully prepared the state
$|\vec{a}\rangle:=\displaystyle\sum_{i\in[N]}\frac{a_{i}}{||\vec{a}||_{2}}|i\rangle$.
∎

_Remark:_ It is also possible to convert a Real Ket image to its NEQR version
by applying phase estimations [17]. The close connections between different
image representations give researchers more confidence in developing quantum
image processing tools based on a particular representation, as they are
likely to become universal once quantum techniques are sufficiently developed.

Quantum image preparation, also known as quantum (initial) state preparation,
has been extensively studied over years, e.g., [40, 41, 38]. Current
techniques allow to efficiently prepare the Real Ket (or NEQR) state of an
$N\times N$ image in time $O(\text{polylog}N)$ if the prepared image has some
special structures111such as the corresponding quantum data structure [38] is
given or the quantum state to be prepared has a well condition number [42]..
The time required to prepare an arbitrary quantum image is at most no more
than the classical preparation time, up to a logarithmic factor. Specifically,
a quantum image of form (13) or (14) can be prepared by applying at most
$N^{2}$ conditional operations, each of which has a time complexity
polylog($N$).

The following proposition will be used in Section V. It states that the inner
products of quantum states can be estimated in parallel:

######  Proposition 5 (Parallel swap test [43] )

. Given $2N$ quantum states $|\vec{u}_{0}\rangle$,
$|\vec{v}_{0}\rangle...|\vec{u}_{N-1}\rangle$, $|\vec{v}_{N-1}\rangle$, and
two state preparation unitaries:
$|k\rangle|0\rangle\rightarrow|k\rangle|\vec{u}_{k}\rangle$ and
$|k\rangle|0\rangle\rightarrow|k\rangle|\vec{v}_{k}\rangle$ (where $k\in[N]$)
that can be implemented in time O$(T_{in})$. There is a quantum algorithm with
runtime O$(\frac{T_{in}}{\epsilon})$ to achieve
$|i\rangle|0\rangle\rightarrow|i\rangle|s_{i}\rangle$ for all $i\in[N]$, where
$|s_{i}-\langle\vec{u}_{i}|\vec{v}_{i}\rangle|\leq\epsilon$.

A detailed proof can be found in Theorem 1 in [43].

##  III Quantum reversible multiplication

This section presents a reversible modular multiplication in the quantum
computation framework. Fix $N=2^{n}$. For fixed $0\leq i<2^{n-1}$ and $C=2i$,
the multiplication in $\mathbb{Z}_{N}$ by $C$ is irreversible, because the
mapping: $a\rightarrow aC\mod N$ maps both $a=0$ and $a=2^{n-1}$ to $0$. On
the other hand, for $D=2i+1$, the multiplication in $\mathbb{Z}_{N}$ by $D$ is
reversible, because if $aD=bD\text{ mod }N$ for some $a,b\in[N]$ and $a\neq
b$, then $(a-b)D\mid 2^{n}$. Since $D$ is odd, it must be that $a-b$ is a
multiple of $2^{n}$; in particular, $|a-b|\geq 2^{n}$, which is impossible for
$a,b\in[N]$.

Below, we design a unitary realization of the multiplication in
$\mathbb{Z}_{N}$ by any odd number $D\in[N]$. For $1\leq k\leq n$,

  * •

let $M_{k}$ be a to-be-realized unitary operator performing

$\displaystyle|a_{k}\rangle|b_{k}\rangle\rightarrow|a_{k}\rangle|a_{k}b_{k}\
\text{mod}\ 2^{k}\rangle,\forall a_{k},b_{k}\in[2^{k}]\ \text{and odd}\
a_{k}.$ (22)

  * •

Let $A_{k}$ be the following controlled addition in $\mathbb{Z}_{2^{k}}$:

$|a_{k}\rangle|b_{k}\rangle|c\rangle\rightarrow\left\\{\begin{aligned}
&|a_{k}\rangle|a_{k}+b_{k}\text{ mod }2^{k}\rangle|c\rangle,&\text{if }c=1;\\\
&|a_{k}\rangle|b_{k}\rangle|c\rangle,&\text{if }c=0;\end{aligned}\right.$

$\hfill\forall a_{k},b_{k}\in[2^{k}],c\in[2].$

The time complexity for performing each $k$-qubit addition is O$(k^{2})$ [44],
so is the time complexity of performing unitary operator $A_{k}$ (see Section
4.2 in [17]).

For any $1\leq k\leq n$, for any $a_{k}\in[2^{k}]$, let the binary
representation of integer $a_{k}$ be

$\displaystyle
B_{i_{k-1}...i_{0}}:=\displaystyle\sum_{l\in[k]}2^{l}i_{l},\quad\text{where}\
i_{l}\in[2].$

When $k=1$, multiplication operator $M_{1}$ is the identity:

$\displaystyle|1\rangle|b_{1}\rangle\rightarrow|1\rangle|b_{1}\rangle,\quad\forall\
b_{1}\in[2].$

In the following, we realize multiplication operator $M_{k+1}$ by $A_{l}$ for
$1\leq l\leq k$ recursively.

For any $a_{k+1},b_{k+1}\in[2^{k+1}]$ where $a_{k+1}$ is odd, let their binary
representations be $B_{i_{k}...i_{1}1}$ and $B_{j_{k}...j_{1}j_{0}}$
respectively. Then their modular multiplication is

$\displaystyle(2^{k}i_{k}+B_{i_{k-1}...i_{1}1})(2B_{j_{k}...j_{1}}+j_{0})\text{
mod }2^{k+1}$ (23) $\displaystyle=$
$\displaystyle(2B_{i_{k-1}...i_{1}1}B_{j_{k}...j_{1}}+B_{i_{k}...i_{1}1}\frac{1+(-1)^{j_{0}-1}}{2})\text{
mod }2^{k+1}$ $\displaystyle=$
$\displaystyle\big{\\{}2(B_{i_{k-1}...i_{1}1}B_{j_{k}...j_{1}}\text{ mod
}2^{k})+2\delta^{1}_{j_{0}}B_{i_{k}...i_{1}}+j_{0}\big{\\}}\text{mod }2^{k+1}$

where $\delta^{1}_{j_{0}}$ is the Kronecker symbol.

So the modular multiplication $M_{k+1}$ of two integers
$B_{i_{k}...i_{1}1},B_{j_{k}...j_{1}j_{0}}\in[2^{k+1}]$, each containing $k+1$
binary digits, can be decomposed into two operators: the modular
multiplication $M_{k}$ of integers
$B_{i_{k-1}...i_{1}1},B_{j_{k}...j_{1}}\in[2^{k}]$, which occurs in the first
$k$ binary digits, followed by the controlled modular addition $A_{k}$ of the
above modular multiplication result and $B_{i_{k}...i_{1}}$, which also occurs
in the first $k$ binary digits, while the control digit $j_{0}\in[2]$ remains
in the last binary digit. The quantum circuit realizing this decomposition is
shown in Fig. 3.

Figure 3: The quantum circuit realizing $M_{k+1}$ by $M_{k}$ and $A_{k}$.

In (23) the modular multiplication $M_{k}$ can be further decomposed into
$M_{k-1}$ and $A_{k-1}$, and by doing so recursively, $M_{k+1}$ is finally
decomposed into a series of controlled modular additions:
$A_{1},A_{2},\cdots,A_{k}$. Since performing $A_{j}$ has time complexity
C$j^{2}$ for fixed constant $C>0$ and varying $1\leq j\leq k$, by
$1^{2}+2^{2}+...+k^{2}=O(k^{3})$, we get that the time complexity of realizing
$M_{k+1}$ (hence $M_{k}$) by (14), is O$(k^{3})$.

##  IV Quantum Radon Transform

The classical PDRT in (3) can be viewed as a mapping
$\mathbb{R}^{n^{2}}\rightarrow\mathbb{R}^{n(n+1)}$ as follows:

$\displaystyle\Big{(}f(0,0),f(1,0),...,$ $\displaystyle
f(n-1,n-1)\Big{)}\xrightarrow{\textnormal{PDRT}}$
$\displaystyle\Big{(}r_{0}(0),r_{0}(1),...,r_{n}(n-1)\Big{)},$ (24)

whose transformation matrix is not unitary222The lines $L^{n}_{0,0}$,
$L^{n}_{0,1}$ have intersection, so the rows of transformation matrix
corresponding to $r_{0}(0)$, $r_{1}(0)$ are not orthogonal.. It is hard to
directly design a quantum algorithm for PDRT by (3). So, we consider utilizing
the related Fourier slice property in (4).

By replacing the traditional multiplication in (4) with the quantum reversible
multiplication, we successfully design an efficient quantum algorithm that can
realize a novel transform similar to PDRT. We name it the quantum Randon
transform:

######  Definition 4 (Quantum Radon transform)

Let $f$ be a function defined on $\mathbb{Z}^{2}_{n}$, let
$\lfloor\cdot\rfloor$ be the floor function, and let $\tilde{f}$ defined on
$\mathbb{Z}^{2}_{2n}$ be related to $f$ as

$\displaystyle\tilde{f}(x^{\prime},y^{\prime}):=\frac{1}{2}(-1)^{\lfloor\frac{x^{\prime}}{n}\rfloor+\lfloor\frac{y^{\prime}}{n}\rfloor}f(x^{\prime}\text{\
mod }n,y^{\prime}\text{ mod }n),x^{\prime},y^{\prime}\in[2n].$ (25)

Recall that $L^{2n}_{l,k}$ represent discrete lines on the lattice
$\mathbb{Z}^{2}_{2n}$:

$\displaystyle
L^{2n}_{l,k}:=\\{(x^{\prime},y^{\prime})|x^{\prime}+ky^{\prime}=l\
(\text{mod}\ 2n),\ x^{\prime},\ y^{\prime}\in[2n]\\}.$

Then, the quantum Radon transform of $f$ is a function defined on
$\mathbb{Z}^{2}_{2n}$ as

$\displaystyle
QR_{f}(l,k)=\frac{1}{\sqrt{2n}}\displaystyle{\sum_{(x^{\prime},y^{\prime})\in
L^{2n}_{l,k}}}\tilde{f}(x^{\prime},y^{\prime}),\quad l,k\in[2n]$ (26)

The definition of QRT is derived from the following Algorithm 1 rather than a
deliberate construction, and it may seem a little complicated. Now, let us
give a closer observation of $\tilde{f}$ in (25). For any $x,y\in[n]$,

$\displaystyle f(x,y)$
$\displaystyle=\tilde{f}(x,y)=-\tilde{f}(x+n,y)=-\tilde{f}(x,y+n)$
$\displaystyle=\tilde{f}(x+n,y+n).$ (27)

Graphically speaking, the QRT of an image $f$ can be viewed as the classical
PDRT that performs on the symmetrized double-sized original image. Fig. 4
shows the particular summation method adopted by the QRT. Intuitively, the
alternating sums seems to tend to suppress the changes in the Radon domain
along the direction of interception, resulting in a lower sensitivity of QRT
to linear singularities. Even worse, only a half of slopes (lines) could be
detected by QRT, since when slope $k$ is even, by (IV) it holds that

$\displaystyle QR_{f}(l,k)=0.$ (28)

However, by experiments (see Fig. 8-7), it has been found that QRT preserves
the good denoising capability as in the classical PDRT; it will be further
discussed in Section VI-A later. We hold the view that the QRT has a practical
value comparable to that of the classical PDRT. On the other hand, as an
advantage, QRT is invertible, which is derived from the following invertible
QRT algorithm. The most important advantage lies in that the QRT can be
performed exponentially faster than the classical PDRT, by the following
quantum Radon transform algorithm:

(a) Sum the values of $f$ at all white points.

(b) Sum the values of $f$ at white points then subtract the values at black
points.

Figure 4: A toy example to show the different summation methods adopted to
compute the PDRT (left) of $f$ and QRT (right) of $f$, where $f$ is a function
defined on $\mathbb{Z}^{2}_{11}$. Fig. 4a shows the points used to compute
$r_{-3}(1)$ as defined in (3), i.e., summation along the discrete line
$x-3y=1\mod 11$. Fig. 4b shows the points used to compute $QR_{f}(1,-3\mod
2n)$.  Algorithm 1 Quantum Radon transform

1:An $N\times N$ quantum image
$\displaystyle{\sum_{x,y\in[N]}}f(x,y)|x\rangle|y\rangle$.

2:$\displaystyle\sum_{l,j^{\prime\prime}\in[2N]}QR_{f}(l,j^{\prime\prime})|l\rangle|j^{\prime\prime}\rangle.$

3:Prepare the state
$\displaystyle{\sum_{x,y\in[N]}}f(x,y)|x\rangle|1\rangle|y\rangle|1\rangle.$

4:Apply phase shifts conditioned to $|x\rangle,|y\rangle$ to prepare
$\displaystyle\displaystyle{\sum_{x,y\in[N]}}g(x,y)|x\rangle|1\rangle|y\rangle|1\rangle,$
(29) where $g(x,y)=f(x,y)e^{-2\pi I\frac{x+y}{2N}}$.

5:Perform quantum Fourier transforms on $|x\rangle$ and $|y\rangle$
respectively to get
$\displaystyle\displaystyle\sum_{i,j\in[N]}\Big{(}\frac{1}{N}\displaystyle\sum_{x,y\in[N]}g(x,y)e^{-2\pi
I\frac{ix+jy}{N}}\Big{)}|i\rangle|1\rangle|j\rangle|1\rangle.$ (30)

6:Perform inverse quantum reversible multiplication on registers
$|i^{\prime}\rangle=|i,1\rangle$ and $|j^{\prime}\rangle=|j,1\rangle$ to
prepare
$\displaystyle\displaystyle\sum_{i^{\prime},j^{\prime\prime}\in[2N]}\Big{(}\frac{1}{2N}\displaystyle\sum_{x^{\prime},y^{\prime}\in[2N]}\tilde{f}(x^{\prime},y^{\prime})e^{-2\pi
I\frac{i^{\prime}x^{\prime}+(i^{\prime}\odot
j^{\prime\prime})y^{\prime}}{2N}}\Big{)}$ $\displaystyle\hskip
173.56198pt|i^{\prime}\rangle|j^{\prime\prime}\rangle,$ (31) where
$i^{\prime}\odot j^{\prime\prime}=j^{\prime}$, $\odot$ is the quantum
reversible multiplication such that $i^{\prime}\odot
j^{\prime\prime}=i^{\prime}j^{\prime\prime}$ for odd $i^{\prime}\in[2N]$.

7:Apply inverse quantum Fourier transform to $|i^{\prime}\rangle$ to prepare
$\displaystyle\displaystyle\frac{1}{\sqrt{2N}}\sum_{l,j^{\prime\prime}\in[2N]}\Big{(}\frac{1}{2N}\displaystyle\sum_{i^{\prime},t\in[2N]}\displaystyle\sum_{(x^{\prime},y^{\prime})\in
L^{2N}_{t,j^{\prime\prime}}}\tilde{f}(x^{\prime},y^{\prime})$
$\displaystyle\hskip 85.35826pt\times e^{-2\pi I\frac{i^{\prime}t}{2N}}e^{2\pi
I\frac{i^{\prime}l}{2N}}\Big{)}|l\rangle|j^{\prime\prime}\rangle$ (32)
$\displaystyle=\displaystyle\sum_{l,j^{\prime\prime}\in[2N]}QR_{f}(l,j^{\prime\prime})|l\rangle|j^{\prime\prime}\rangle.$
(33)

###  IV-A Analysis of Algorithm 1

The Step 1 of Algorithm 1 is adding two extra qubits initially in the state
$|1\rangle|1\rangle$. In Step 2, the conditional phase shifts can be
implemented in time logarithmic in $N$, by a method similar to the
implementation of conditional rotation (cf. Proposition 3). After applying
Step 3, the resulting state can be viewed as the quantum Fourier transform of
some $(2N+2)$-qubit state
$\displaystyle\sum_{x^{\prime},y^{\prime}\in[2N]}g^{\prime}(x^{\prime},y^{\prime})|x^{\prime}\rangle|y^{\prime}\rangle$,
i.e.,

$\displaystyle\begin{cases}&\mathcal{F}_{2}\\{g^{\prime}\\}(2i+1,2j+1)=\mathcal{F}_{2}\\{g\\}(i,j),\quad
i,j\in[N],\\\ &\mathcal{F}_{2}\\{g^{\prime}\\}(i,j)=0,\quad i,j\in[2N]\text{\
are not both odd}.\end{cases}$ (34)

The explicit expression of function $g^{\prime}$ can be deduced by inverse
Fourier transform as follows:

$\displaystyle\hskip 24.18501ptg^{\prime}(l^{\prime},t^{\prime})=$
$\displaystyle\frac{1}{2N}\frac{1}{N}\displaystyle\sum_{x,y,i,j\in[N]}g(x,y)e^{-2\pi
I\frac{ix+jy}{N}}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\times e^{2\pi
I\frac{(2i+1)l^{\prime}+(2j+1)t^{\prime}}{2N}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}g\left(l^{\prime}\text{ mod }N,t^{\prime}\text{ mod
}N\right)e^{2\pi I\frac{l^{\prime}+t^{\prime}}{2N}}$ $\displaystyle=$
$\displaystyle\tilde{f}(l^{\prime},t^{\prime}),\qquad\qquad\qquad\qquad\qquad
l^{\prime},t^{\prime}\in[2N],$

so that the resulting state in Step 3 can be rewritten as

$\displaystyle\displaystyle\sum_{i^{\prime},j^{\prime}\in[2N]}\Big{(}\frac{1}{2N}\displaystyle\sum_{x^{\prime},y^{\prime}\in[2N]}\tilde{f}(x^{\prime},y^{\prime})e^{-2\pi
I\frac{i^{\prime}x^{\prime}+j^{\prime}y^{\prime}}{2N}}\Big{)}|i^{\prime}\rangle|j^{\prime}\rangle.$
(35)

After applying inverse quantum reversible multiplication to
$|i^{\prime}\rangle|j^{\prime}\rangle$ of (35) in Step 4, by (22), the
resulting state is as in (6). Now, from the fact that for any
$j^{\prime\prime}\in[2N]$,

$\displaystyle\begin{cases}&i^{\prime}\odot
j^{\prime\prime}=i^{\prime}j^{\prime\prime},\qquad\qquad\qquad\text{if
$i^{\prime}$ is odd,}\\\
&\mathcal{F}_{2}\\{\tilde{f}\\}(i^{\prime},j^{\prime\prime})=0,\qquad\qquad\text{
if $i^{\prime}$ is even,}\end{cases}$

one gets

$\displaystyle\hskip
24.18501pt\mathcal{F}_{2}\\{\tilde{f}\\}(i^{\prime},i^{\prime}\odot
j^{\prime\prime})=\mathcal{F}_{2}\\{\tilde{f}\\}(i^{\prime},i^{\prime}j^{\prime\prime}),\qquad
i,j\in[2N].$

Then the resulting state in Step 4 can be rewritten as
$\displaystyle\sum_{i^{\prime},j^{\prime\prime}\in[2N]}\mathcal{F}_{2}\\{\tilde{f}\\}(i^{\prime},i^{\prime}j^{\prime\prime})|i^{\prime}\rangle|j^{\prime\prime}\rangle$.
The final result in Step 5 is by the following proposition:

###### Proposition 6

[Fourier slice property of QRT] In the notations of Algorithm 1, for any given
slope $j^{\prime\prime}\in[2N]$, the $1$-D Fourier transform of
$QR_{f}(l,j^{\prime\prime})$ with respect to $l$ is

$\displaystyle\mathcal{F}_{1}\\{QR_{f}(\cdot,j^{\prime\prime})\\}(i^{\prime})=\mathcal{F}_{2}\\{\tilde{f}\\}(i^{\prime},i^{\prime}j^{\prime\prime}),\quad
i^{\prime}\in[2N].$ (36)

###### Proof:

Let $\mathcal{F}_{1}^{-1}$ be the inverse transform of $\mathcal{F}_{1}$.
Then,

$\displaystyle\mathcal{F}_{1}^{-1}\\{\mathcal{F}_{2}\\{\tilde{f}\\}(\cdot,\cdot
j^{\prime\prime})\\}(l)$ (37) $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sqrt{2N}}\displaystyle\sum_{i^{\prime}\in[2N]}\Big{(}\frac{1}{2N}\displaystyle\sum_{x^{\prime},y^{\prime}\in[2N]}\tilde{f}(x^{\prime},y^{\prime})e^{-2\pi
I\frac{i^{\prime}x^{\prime}+i^{\prime}j^{\prime\prime}y^{\prime}}{2N}}\Big{)}e^{2\pi
I\frac{i^{\prime}l}{2N}}$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sqrt{2N}}\Big{(}\frac{1}{2N}\displaystyle\sum_{i^{\prime},t\in[2N]}\displaystyle\sum_{(x^{\prime},y^{\prime})\in
L^{2N}_{t,j^{\prime\prime}}}\tilde{f}(x^{\prime},y^{\prime})e^{-2\pi
I\frac{i^{\prime}t}{2N}}e^{2\pi I\frac{i^{\prime}l}{2N}}\Big{)}$
$\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sqrt{2N}}\displaystyle\sum_{i^{\prime},t\in[2N]}\Big{(}\frac{1}{2N}\displaystyle\sum_{(x^{\prime},y^{\prime})\in
L^{2N}_{t,j^{\prime\prime}}}\tilde{f}(x^{\prime},y^{\prime})\Big{)}e^{-2\pi
I\frac{i^{\prime}t}{2N}}e^{2\pi I\frac{i^{\prime}l}{2N}}$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sqrt{2N}}\displaystyle\sum_{(x^{\prime},y^{\prime})\in
L^{2N}_{l,j^{\prime\prime}}}\tilde{f}(x^{\prime},y^{\prime})=QR_{f}(l,j^{\prime\prime}).$

∎

The overall runtime of Algorithm 1 is the sum of $O(\text{log}N)$ time
required to perform conditional phase shifts [45], $O(\text{log}^{2}N)$ time
required to perform quantum Fourier transform and its inverse, and
$O(\text{log}^{3}N)$ time required to perform the inverse multiplication.
Observe that all quantum gates used in Algorithm 1 are unitary and invertible.
So, an efficient inverse QRT algorithm follows immediately by reversing all
circuits in Algorithm 1. In conclusion, we have the following theorem:

###### Theorem 1

Given an $N\times N$ quantum image state
$\displaystyle{\sum_{x,y\in[N]}}f(x,y)|x\rangle|y\rangle$, then its (inverse)
QRT can be performed in time O$(\log^{3}N)$.

##  V Quantum-mechanical implementation of DRT with the interpolation method

The discretization of Radon transform needs to approximate the line integral.
A simple solution is to sample and sum along the integral line, where the
interpolation method can be employed to estimate the undefined values at non-
integer lattice sample points. Several kinds of sample and interpolation
methods are available for implementing IDRT [46, 47, 6]. Below, we define a
simplest kind of IDRT (SIDRT), and give its quantum algorithm, with the aim of
showing a general quantum approach to achieving this interpolation-based kind
of DRT.

###### Definition 5

Let the set of slopes of the basically horizontal lines and basically vertical
lines be

$\displaystyle S^{\parallel}:=\\{tan\theta_{j}|\theta_{j}=\frac{\pi
j}{N},-\frac{N}{4}\leq j<\frac{N}{4}\ \text{and}\ j\in\mathbb{Z}\\},$
$\displaystyle S^{\perp}:=\\{\frac{1}{tan\theta_{j}}|\theta_{j}=\frac{\pi
j}{N},-\frac{N}{4}<j\leq\frac{N}{4}\ \text{and}\ j\in\mathbb{Z}\\},$ (38)

respectively, where agree that $\frac{1}{0}=\infty$ and $\frac{1}{\infty}=0$.
Let $\Delta_{ki}$ be the fractional part of $ki$. Then the SIDRT of a function
$f$ on $\mathbb{Z}^{2}_{N}$ is

$\displaystyle
P_{k}(l)=\frac{1}{\sqrt{N}}\displaystyle\sum_{i\in[N]}\Big{[}\sqrt{1-|\Delta_{ki}|^{2}}f\Big{(}i,l+\lfloor
ki\rfloor\text{ mod }N\Big{)}$ $\displaystyle+\Delta_{ki}f\Big{(}i,l+\lfloor
ki\rfloor+1\text{ mod }N\Big{)}\Big{]},\quad k\in S^{\parallel},l\in[N],$ (39)
$\displaystyle
P_{k}(l)=\frac{1}{\sqrt{N}}\displaystyle\sum_{i\in[N]}\Big{[}\sqrt{1-|\Delta_{\frac{i}{k}}|^{2}}f\Big{(}l+\lfloor\frac{i}{k}\rfloor\text{
mod }N,i\Big{)}$
$\displaystyle+\Delta_{\frac{i}{k}}f\Big{(}l+\lfloor\frac{i}{k}\rfloor+1\text{
mod }N,i\Big{)}\Big{]},\quad k\in S^{\perp},l\in[N].$ (40)

By definition, $P_{k}(l)$ can be viewed as an approximate discrete line
integral along the line with interception $l$ and slope $k$.

The SIDRT is proposed for easily achieved in the quantum case. Although it
adopts a relatively simple interpolation method, the SIDRT is so useful that
enables to detect lines in complicated image, as shown in Figure 10; more
details can be found in Section VI-B.

Now, we turn to the quantum implementation of the SIDRT. We first consider
approximating the integrals along basically horizontal lines, i.e., $P_{k}(l)$
as defined in (5) where $l\in[N]$ and $k\in S^{\parallel}$.

Let $k_{\theta}:=\tan(\frac{\pi\theta}{N}-\frac{\pi}{4})$. For any $i\in[N]$,
$\theta\in[\frac{N}{2}]$, we can prepare a quantum state that contains the
location information of lattice points used to compute $P_{k_{\theta}}(l)$, by
the following sequence of mappings:

$\displaystyle|i\rangle|\theta$
$\displaystyle\rangle|0\rangle|0\rangle\longrightarrow|i\rangle|\theta\rangle|ik_{\theta}\rangle|0\rangle$
$\displaystyle\longrightarrow|i\rangle|\theta\rangle\big{|}\lfloor
ik_{\theta}\rfloor\big{\rangle}|\Delta_{ik_{\theta}}\rangle\left(\sqrt{1-|\Delta_{ik_{\theta}}|^{2}}|0\rangle+\Delta_{ik_{\theta}}|1\rangle\right)$
$\displaystyle\longrightarrow|i\rangle|\theta\rangle\big{|}\lfloor
ik_{\theta}\rfloor\big{\rangle}|0\rangle\left(\sqrt{1-|\Delta_{ik_{\theta}}|^{2}}|0\rangle+\Delta_{ik_{\theta}}|1\rangle\right).$
(41)

In (V), the first step is by trigonometric function [39] and arithmetics [17]
in computation basis. The second step is to use the decimal part of
$ik_{\theta}$ to perform control rotations (cf. Proposition 3) on the last
qubit. The third step is uncomputing $|\Delta_{ik_{\theta}}\rangle$.

Combining (V) with the addition in the computational basis gives the following
lemma:

###### Lemma 2

The following ‘location state’ can be prepared in time $O(\emph{polylog}N)$:

$\displaystyle|\theta\rangle|0\rangle\stackrel{{\scriptstyle
V}}{{\longrightarrow}}$
$\displaystyle|\theta\rangle\displaystyle\sum_{i\in[N]}\frac{1}{\sqrt{N}}\Big{(}\sqrt{1-|\Delta_{ik_{\theta}}|^{2}}|i\rangle\big{|}\lfloor
ik_{\theta}\rfloor\big{\rangle}|0\rangle$
$\displaystyle+\Delta_{ik_{\theta}}|i\rangle\big{|}\lfloor
ik_{\theta}\rfloor+1\big{\rangle}|1\rangle\Big{)},\quad\theta\in[\frac{N}{2}].$
(42)

_Remark:_ The state on the right hand side of (2) is called ‘location state’,
since its first two quibts record the locations of points used to compute
$P_{k}(0)$ for $k\in S^{\parallel}$, according to (5).

(a)

(b)

(c)

Figure 5:  Left: The ‘Lena’ image. Middle: Classical Radon transform of
‘Lena’. Right: Quantum Radon transform of ’Lena’.

Now, with the unitary capable of preparing a quantum image $f$ as follows:

$\displaystyle|0\rangle{\longrightarrow}\displaystyle\sum_{i,j\in[N]}f(i,j)|i\rangle|j\rangle,$
(43)

one can approximate the SIDRT of $f$ by the following theorem:

###### Theorem 3

Given the unitary that can prepare a quantum image
$\displaystyle\sum_{i,j\in[N]}f(i,j)|i\rangle|j\rangle$ in time O$(T_{in})$,
one can approximate $P_{k}(l)$, the SIDRT of $f$, to precision O$(\epsilon)$
in time O$(\frac{T_{in}}{\epsilon}\emph{polylog}N)$. Namely, there is a
quantum algorithm runs in time O$(\frac{T_{in}}{\epsilon}\emph{polylog}N)$ to
achieve the mapping

$\displaystyle|\theta\rangle|l\rangle|0\rangle\rightarrow|\theta\rangle|l\rangle|\tilde{P}_{k_{\theta}}(l)\rangle,\hskip
28.45274pt\forall l\in[N],\theta\in[N],$ (44)

where $|\tilde{P}_{k_{\theta}}(l)-P_{k_{\theta}}(l)|\leq\epsilon$, and
$k_{\theta}:=\tan(\frac{\pi\theta}{N}-\frac{\pi}{4})$ such that $k_{\theta}\in
S^{\parallel}\cup S^{\perp}$ as defined in (5).

###### Proof:

We first consider the realization of (44) for $\theta\in[\frac{N}{2}]$, i.e.,
the case where $k_{\theta}\in S^{\parallel}$. Observe that the following
‘image state’ can be prepared in time O($T_{in}\text{polylog}N$) by the
mapping

$\displaystyle|\theta\rangle|l\rangle|0\rangle\longrightarrow\displaystyle|\theta\rangle|l\rangle\sum_{i,j\in[N]}f(i,j)|i\rangle|j\rangle(\frac{1}{\sqrt{2}}|0\rangle+\frac{1}{\sqrt{2}}|1\rangle),$
(45)

where $\theta\in[\frac{N}{2}]$ and $l\in[N]$. Lemma 2 allows us to prepare the
‘location state’ by the mapping

$\displaystyle|\theta\rangle|l\rangle|0\rangle\longrightarrow|\theta\rangle|l\rangle\displaystyle\sum_{i\in[N]}$
$\displaystyle\frac{1}{\sqrt{N}}\Big{(}\sqrt{1-|\Delta_{ik_{\theta}}|^{2}}|i\rangle\big{|}l+\lfloor
ik_{\theta}\rfloor\big{\rangle}|0\rangle$
$\displaystyle+\Delta_{ik_{\theta}}|i\rangle\big{|}l+\lfloor
ik_{\theta}\rfloor+1\big{\rangle}|1\rangle\Big{)}$ (46)

in time O$(\text{polylog}N)$, where $\theta\in[\frac{N}{2}]$ and $l\in[N]$.

Now, by parallel swap test (Proposition 5), one can approximate the
$\frac{N^{2}}{2}$ inner products of the ‘image state’ in (45) and the
‘location state’ in (V) in parallel $\theta,\ l$. Since the result of each
inner product is just $\frac{1}{\sqrt{2}}P_{k_{\theta}}(l)$, one can perform

$\displaystyle|\theta\rangle|l\rangle|0\rangle\rightarrow|\theta\rangle|l\rangle|p_{k_{\theta}}(l)\rangle,\hskip
28.45274pt\forall\theta\in[\frac{N}{2}],\ l\in[N],$ (47)

in time O$(\sqrt{2}\frac{T_{in}}{\epsilon}\text{polylog}N)$, where
$|p_{k_{\theta}}(l)-\frac{1}{\sqrt{2}}P_{k_{\theta}}(l)|\leq\frac{\epsilon}{\sqrt{2}}$.
The theorem holds by setting
$\tilde{P}_{k_{\theta}}(l)=\sqrt{2}p_{k_{\theta}}(l)$. An efficient
implementation of multiplication by $\sqrt{2}$, a known constant, in the
computation basis can be found in [48]. The realization of (44) for
$\theta\in[N]\setminus[\frac{N}{2}]$ is similar. ∎

Given a Real Ket quantum image $f$ and its preparation unitary as in (43), by
Theorem 3, after performing the transform of (44) on input state
$\sum_{\theta,l\in[N]}\frac{1}{N}|\theta\rangle|l\rangle|0\rangle$, one can
prepare the SIDRT of $f$ in NEQR encoded form. This NEQR encoded output is
already can be used for practical application, such as line detection shown in
Section VI-B later. Moreover, if necessary, one can continue to transform the
NEQR encoded output into its Real Ket version by Proposition 4, thus keeping
the input and output encoded in the same way.

###  V-A Efficiency Analysis of Theorem 3

Below, we discuss what is a reasonable choice of the precision $\epsilon$ in
Theorem 3.

Random image. The term ‘random image’ refers to an image whose each pixel
value is sampled from the uniform distribution $\mathbb{U}[0,1]$
independently.

###### Proposition 7

Let $f$ be an $N\times N$ random image, i.e., $f(i,j)\sim\mathbb{U}[0,1]$ for
$i$, $j\in[N]$. Then the expectation of the minimal value of the SIDRT of
normalized $f$ is no less than $\frac{\sqrt{3}}{2N^{0.5}}$.

_Remark:_ The reason why we consider the normalized rather than the original
image here is that the quantum image is a normalized state.

###### Proof:

We begin with a probability inequality. Let $a$, $b$ be two discrete random
variables whose density function are $P(a=a_{i})=p_{i},P(b=b_{i})=q_{i}$,
where $i\in[N]$, respectively. Then by Cauchy-Schwarz inequality,

$\displaystyle
E(\frac{a}{b})=\sum_{i,j\in[N]}p_{i}q_{j}\frac{a_{i}}{b_{j}}=\frac{(\sum_{i,j}p_{i}q_{j}b_{j})(\sum_{i,j}p_{i}q_{j}\frac{a_{i}}{b_{j}})}{E(b)}$
$\displaystyle\geq\frac{(\sum_{i,j}q_{j}p_{i}\sqrt{a_{i}})^{2}}{E(b)}=\frac{(\sum_{i}p_{i}\sqrt{a_{i}})^{2}}{E(b)}=\frac{E^{2}(\sqrt{a})}{E(b)}.$
(48)

Let $\mathbb{U}[0,1]$ be the uniform distribution on $[0,1]$, and let random
variables $x_{k}$, $y_{k}\sim\mathbb{U}[0,1]$ for $k\in[N^{2}]$. Notice that
$f(i,j)\sim\mathbb{U}[0,1]$. Let $P_{k_{0}}(l_{0})$ be the SIDRT of $f$ as in
(5), then for any $k_{0}$, $l_{0}\in[N]$, the expectation of the SIDRT of
normalized $f$ at point $(k_{0},l_{0})$ is

$\displaystyle
E(\frac{P_{k_{0}}(l_{0})}{\sqrt{\displaystyle\sum_{i,j\in[N]}f^{2}(i,j)}})=E(\frac{\frac{1}{\sqrt{N}}\displaystyle\sum_{k\in[N]}\Delta_{k}x_{k}+\sqrt{1-|\Delta_{k}|^{2}}x_{k}}{\sqrt{\displaystyle\sum_{k\in[N]}x_{k}^{2}+\displaystyle\sum_{k\in[N^{2}-N]}y_{k}^{2}}})$
$\displaystyle\geq
E(\frac{\frac{1}{\sqrt{N}}\displaystyle\sum_{k\in[N]}x_{k}}{\sqrt{\displaystyle\sum_{k\in[N]}1+\displaystyle\sum_{k\in[N^{2}-N]}y_{k}^{2}}})\overset{(\ref{vba})}{\geq}\frac{E^{2}(\sqrt{\frac{1}{\sqrt{N}}\displaystyle\sum_{k\in[N]}x_{k}})}{E(\sqrt{N+\displaystyle\sum_{k\in[N^{2}-N]}y_{k}^{2}})}$
$\displaystyle\geq\frac{E^{2}(\sqrt{\frac{1}{\sqrt{N}}\displaystyle\sum_{k\in[N]}x_{k}})}{\sqrt{E(N+\displaystyle\sum_{k\in[N^{2}-N]}y_{k}^{2})}}\approx\frac{\sqrt{3}}{2N^{0.5}}\quad\quad(N\rightarrow\infty),$
(49)

where the last inequality follows from that $E^{2}(\nu)\leq E(\nu^{2})$ for
arbitrary random variable $\nu$, and the last approximate equality is by the
central limit theorem [49] which states that for any independent and
identically distributed (IID) random variables $\nu_{0},\nu_{1}...,\nu_{n-1}$
with mean $\mu$ and variance $\sigma^{2}$,
E$(\sqrt{|\sum_{j\in[n]}\nu_{j}|})\approx\sqrt{n\mu-\frac{\sigma^{2}}{\mu}}$
($n\rightarrow\infty$), and so
$E(\sqrt{\displaystyle\sum_{k\in[N]}x_{k}})\approx\sqrt{\frac{N}{2}-\frac{1}{6}}$.
∎

The above Proposition implies that a reasonable choice of precision in Theorem
3 is $\epsilon=\Theta(\frac{1}{N^{0.5}})$, e.g., $\frac{1}{100N^{0.5}}$, for
producing a good approximation to SIDRT. In this case, when compared with the
classical SIDRT whose running time is $\Omega(N^{3})$, quantum SIDRT achieves
a polynomial speedup by Theorem 3, because the quantum image preparation time
$T_{in}\leq O(N^{2})$ in the worst case, thus $\frac{T_{in}}{\epsilon}\leq
O(N^{2.5})<N^{3}$.

##  VI Quantum application

We present two potential applications—quantum image denoising using QRT and
quantum line detection using SIDRT.

###  VI-A Denoising using QRT

We replicate the denoising experiments in [3], which are specially designed
for testing PDRT. The experiment results shown in Fig. 8-7 suggest that our
QRT is of the good denoising capability as in the classical PDRT. Below, we
give an efficient quantum image denoising algorithm using QRT.

(a)

Figure 6: Comparison of denoising half-plane truncated Gaussian function
$f(x_{1},x_{2})=1_{\\{x_{1}>x_{2}\\}}e^{(x_{1}-128)^{2}-(x_{2}-128)^{2}}$. By
denoising using QRT, the SNR increases 5.0602. The QRT, with a comparable
performance to PDRT, is better than DWT at handling the image which is
piecewise smooth with singularities along a straight line (here is the
diagonal line).

(a)

(b)

(c)

(d)

(e)

Figure 7: In the area marked by the red circles, it can be seen that the QRT,
similar to PDRT, enables to restore the straight edges of geometric solids
more clearly than DWT, even in a lower total denoising performance.

(a) Case I

(b) Case I

(c) Case I

(d) Case I

(e) Case I

Figure 8: Comparison of denoising half-plane truncated Gaussian function
$f(x_{1},x_{2})=1_{\\{x_{1}>x_{2}\\}}e^{(x_{1}-128)^{2}-(x_{2}-128)^{2}}$. By
(12), a higher singal-noise raito (SNR) indicates a better performance of
denoising. So, QRT, with a comparable performance to PDRT, is better than DWT
at handling the image which is piecewise smooth with singularities along a
straight line (here is the diagonal line).

  * Step 1.

Apply the QRT algorithm on the input $N\times N$ quantum image $f$ to perform

$\displaystyle\displaystyle{\sum_{x,y\in[N]}}f(x,y)|x\rangle|y\rangle\longrightarrow\displaystyle\sum_{l,k\in[2N]}QR_{f}(l,k)|l\rangle|k\rangle.$
(50)

  * Step 2.

Perform Hadamard gate (namely, $2\times 2$ Haar transform as shown in Fig. 2)
on the least significant qubit (LSB) of $|l\rangle$ in (50) to produce a state

  * $\displaystyle\frac{1}{\sqrt{2}}$ $\displaystyle\displaystyle\sum_{l^{\prime}\in[N],k\in[2N]}\big{(}QR_{f}(2l^{\prime},k)+QR_{f}(2l^{\prime}+1,k)\big{)}|l^{\prime}0\rangle|k\rangle$ $\displaystyle+\big{(}QR_{f}(2l^{\prime},k)-QR_{f}(2l^{\prime}+1,k)\big{)}|l^{\prime}1\rangle|k\rangle.$ (51)

  * Step 3.

Measure the LSB of the first register. If the outcome is $|0\rangle$, then
apply the inverse Hadamard gate to the measured qubit, and then execute the
inverse QRT algorithm. The resulting state is (up to a normalization factor):

$\displaystyle QR^{-1}\Big{(}$
$\displaystyle\displaystyle\sum_{\begin{subarray}{c}l^{\prime}\in[N],s\in[2]\\\
k\in[2N]\end{subarray}}\frac{1}{2}\Big{(}QR_{f}(2l^{\prime},k)+QR_{f}(2l^{\prime}+1,k)\Big{)}|l^{\prime}s\rangle|k\rangle\Big{)}$

The above is a quantum analogue of PDRT denoising method introduced in Section
II-B. Notice that the Hadamard gate in Step 2 plays the role of Haar wavelet,
and measuring LSB to obtain $|0\rangle$ in Step 3 has the effect of making all
wavelet coefficients zero. So, this quantum QRT denoising method indeed
simulates denoising using QRT$+$Haar wavelet and threshold $\infty$.

To evaluate the efficiency of this denoising method, we consider the success
probability (namely, the probability of measuring $|0\rangle$ in Step 3).
Since $QR_{f}(l^{\prime},k)=0$ holds for any even $k\in[2N]$, the success
probability is thus reduced to:

$\displaystyle\frac{1}{2}\displaystyle\sum_{l^{\prime},k\in[N]}|QR_{f}(2l^{\prime},2k+1)+QR_{f}(2l^{\prime}+1,2k+1)|^{2}.$
(52)

By the following theoretical analysis (cf. Proposition 8) and numerical
simulation (cf. Fig. 9), this probability of success is quite high,
particularly, it has no tendency to be tiny as the image size $N$ increases.

###### Proposition 8

Let $f$ be an $N\times N$ random image, i.e., $f(i,j)\sim\mathbb{U}[0,1]$ for
$i$, $j\in[N]$. Given the noisy image $h(i,j)=f(i,j)+\epsilon e_{ij}$,
$i,j\in[N]$, where $e_{ij}$ is the noise sampled from the normal distribution
$V(0,\sigma^{2})$ independently, and $\epsilon\in\mathbb{R}$ is noise level,
then the QRT denoising method can be implemented with an average probability
of success $p>1/2$.

###### Proof:

By Definition 4, for any fixed $l,k\in[2N]$,

$\displaystyle QR_{h}(l,k)-QR_{f}(l,k)=\frac{1}{2\sqrt{2N}}\sum_{(i,j)\in
L^{2N}_{l,k}}\epsilon e_{ij}.$ (53)

Notice that for two independent random variables $x\sim
V(u_{1},\sigma_{1}^{2})$, $y\sim V(u_{2},\sigma_{2}^{2})$, it holds that $x\pm
y\sim V(u_{1}\pm u_{2},\sigma_{1}^{2}+\sigma_{2}^{2})$. So, the distribution
of the right-hand side of (53) is $\epsilon V(0,\sigma^{2}/4)$. Let
$E_{f,e}(QR_{h})$ denote the expectation of $QR_{h}$ with respect to the
random variables $f,e$. Then, for any $l,k\in[2N]$,

$\displaystyle E_{f,e}(|QR_{h}(l,k)-QR_{h}(l+1,k)|^{2})$ $\displaystyle=$
$\displaystyle E_{f}(|QR_{f}(l,k)-QR_{f}(l+1,k)|^{2})$ $\displaystyle\hskip
56.9055pt+E_{e}(|\frac{\epsilon}{2\sqrt{2N}}\sum_{s\in[2N]}e_{s}-\sum_{s\in[2N]}e_{s}^{\prime}|^{2})$
(54) $\displaystyle=$ $\displaystyle
E_{f}(|QR_{f}(l,k)-QR_{f}(l+1,k)|^{2})+\epsilon^{2}\sigma^{2}/2,$ (55)

where $e_{s},\ e^{\prime}_{s}\sim V(0,\sigma^{2})$ are both the short hands of
IID noises $e_{ij}$ in (53). (54) is by combining (53) and the fact that the
mean of random variables $e_{s},\ e^{\prime}_{s}$ is $0$. Similarly,

$\displaystyle E_{f,e}(|QR_{h}(l,k)+QR_{h}(l+1,k)|^{2})$
$\displaystyle=E_{f}(|QR_{f}(l,k)+QR_{f}(l+1,k)|^{2})+\epsilon^{2}\sigma^{2}/2.$
(56)

Now, we denote

$\displaystyle\Delta_{+}=\displaystyle\sum_{l,k\in[N]}|QR_{h}(2l,2k+1)+QR_{h}(2l+1,2k+1)|^{2}$
$\displaystyle\Delta_{-}=\displaystyle\sum_{l,k\in[N]}|QR_{h}(2l,2k+1)-QR_{h}(2l+1,2k+1)|^{2}$
$\displaystyle\Delta^{\prime}_{+}=\displaystyle\sum_{l,k\in[N]}|QR_{f}(2l,2k+1)+QR_{f}(2l+1,2k+1)|^{2}$
$\displaystyle\Delta^{\prime}_{-}=\displaystyle\sum_{l,k\in[N]}|QR_{f}(2l,2k+1)-QR_{f}(2l+1,2k+1)|^{2}.$

From (52), our aim is to lower bound the following kind of average probability
of success:

$\displaystyle
P_{\text{success}}=\frac{E_{f,e}(\Delta_{+})}{E_{f,e}(\Delta_{+}+\Delta_{-})}$
$\displaystyle\xlongequal{(\ref{kk1}),\
(\ref{kk2})}\frac{E_{f}(\Delta^{\prime}_{+})+N^{2}\epsilon^{2}\sigma^{2}/2}{E_{f}(2\sum_{i,j\in[N]}|f(i,j)|^{2})+N^{2}\epsilon^{2}\sigma^{2}},$
(57)

where the denominator in (VI-A) is by the following relations:

$\displaystyle\sum_{l,k\in[N]}\Delta^{\prime}_{+}+\Delta^{\prime}_{-}=\sum_{l,k\in[2N]}2|QR_{f}(l,k)|^{2}=2\sum_{i,j\in[N]}|f(i,j)|^{2}$
(58)

Notice that if $k$ is even, then $QR_{f}(l,k)\equiv 0$. We first consider
estimating $E_{f}(|QR_{f}(l,k)+QR_{f}(l+1,k)|^{2})$ and
$E_{f}(|QR_{f}(l,k)-QR_{f}(l+1,k)|^{2})$ for any odd $k\in[2N]$. Denote
intervals

$\displaystyle P_{1}:=[0,N)\times[0,N)\bigcup[N,2N)\times[N,2N),$
$\displaystyle P_{2}:=[0,N)\times[N,2N)\bigcup[N,2N)\times[0,N),$ (59)

and denote the number of points in $L^{2N}_{l,k}\bigcap P_{1}$,
$L^{2N}_{l,k}\bigcap P_{2}$ by

$\displaystyle C_{1}(L^{2N}_{l,k})=\text{card}(\\{(x,y)\in
L^{2N}_{l,k}|(x,y)\in P_{1}\\}),$ $\displaystyle
C_{2}(L^{2N}_{l,k})=\text{card}(\\{(x,y)\in L^{2N}_{l,k}|(x,y)\in P_{2}\\}),$
(60)

respectively. By geometry, it can be verified that

$\displaystyle C_{1}(L^{2N}_{l,k})+C_{2}(L^{2N}_{l,k})=2N,\quad\forall
l,k\in[2N],$ (61)
$\displaystyle|C_{1}(L^{2N}_{l,k})-C_{1}(L^{2N}_{l+1,k})|=2,\quad\forall\textnormal{
odd }k\in[2N].$ (62)

By (26), $QR_{f}(l,k)$ is a sum of values $\tilde{f}(i,j)$ over line
$L^{2N}_{l,k}$ on $\mathbb{Z}^{2}_{2N}$, where each $\tilde{f}(i,j)$ takes
values from uniform distribution on $[0,1]$ for $(i,j)\in P_{1}$, and each
$\tilde{f}(i,j)$ takes values from uniform distribution on $[-1,0]$ for
$(i,j)\in P_{2}$. Combining (61) and (62) gives that for any odd $k$,

$\displaystyle E_{f}(|QR_{f}(l,k)-QR_{f}(l+1,k)|^{2})$ (63) $\displaystyle=$
$\displaystyle\frac{1}{8N}E(\big{|}\sum\limits_{C_{1}(l)+C_{2}(l+1)}\nu-\displaystyle\sum_{C_{2}(l)+C_{1}(l+1)}\nu\big{|}^{2})$
$\displaystyle=$
$\displaystyle\frac{1}{8N}E(\big{|}\displaystyle\sum_{2N+2}\nu-\displaystyle\sum_{2N-2}\nu\big{|}^{2}),$

where $C_{1}(l)$ is the short hand of $C_{1}(L^{2N}_{l,k})$, and $\sum_{k}\nu$
denotes the summation of $k$ numbers of independent random variables
$\nu\sim\mathbb{U}[0,1]$. Notice that for IID random variables
$\nu_{1},\nu_{2}\sim\mathbb{U}[0,1]$, we have
$E(|\nu_{1}|^{2})=|E(\nu_{1})|^{2}+D(\nu_{1})$, where $D$ is the variance,
$D(\nu_{1}\pm\nu_{2})=2D(\nu_{1})$, and $E(\nu_{1}\pm\nu_{2})=E(\nu_{1})\pm
E(\nu_{1})$. So,

$\displaystyle
E(\big{|}\displaystyle\sum_{2N+2}\nu-\displaystyle\sum_{2N-2}\nu$
$\displaystyle\big{|}^{2})=D(\sum_{2N+2}\nu-\displaystyle\sum_{2N-2}\nu)$
$\displaystyle+E(\sum_{2N+2}\nu-\displaystyle\sum_{2N-2}\nu)=\frac{N}{3}+2.$
(64)

Combining (VI-A) and (63) gives

$\displaystyle
E_{f}(\displaystyle\sum_{l\in[N],k\in[N]}|QR_{f}(2l,2k+1)-QR_{f}(2l+1,2k+1)|^{2})$
$\displaystyle=N^{2}\times\frac{1}{8N}\times(\frac{N}{3}+2)=\frac{N^{2}}{24}+\frac{N}{4}$
(65)

Then by (VI-A), (58), (VI-A), and
$E(\sum\limits_{i,j\in[N]}|f(i,j)|^{2})=\frac{N^{2}}{3}$, the lower bound is
achieved as follows:

$\displaystyle P_{\text{success}}$
$\displaystyle=\frac{E(2|f|^{2})-(\frac{N^{2}}{24}+\frac{N}{4})+N^{2}\frac{\epsilon^{2}\sigma^{2}}{2}}{E(2|f|^{2})+N^{2}\epsilon^{2}\sigma^{2}}$
$\displaystyle=1-\frac{\frac{1}{24}+\frac{1}{4N}+\frac{\epsilon^{2}\sigma^{2}}{2}}{\frac{2}{3}+\epsilon^{2}\sigma^{2}}>1/2.$
(66)

∎

_Remark:_ For a real-world image, due to the smoothness of images, the ratio
of wavelet coefficients to filter coefficients is likely to be very small, so
that the probability of success of quantum denoising method is likely to be
high, such as 95$\%$ shown in Figure 9.

By treating the probability of success as a constant, we conclude that the
QRT-based quantum image denoising method has a time complexity O$(\log^{3}n)$,
which is exponentially faster than the classical PDRT denoising method.

Also, there are more quantum techniques that allow to simulate denoising using
other Daubechies wavelets [26] and thresholds [38] in the quantum computation
framework.

###  VI-B Line detection using SIDRT

The SIDRT enables one to detect possible line in images: when the SIDRT
$P_{k}(l)$ of image $f$ reaches maximum for some pair $l$ and $k$, there is
likely to be a straight line with gradient $k$ and interception $l$ in the
image. For example, suppose that $f$ depicts a line composed of points with
grayscale $1$, while the grayscales of all other points in the picture are
$0$. Then the SIDRT reaches maximum at the corresponding line, with the
coordinates of the maximum in the Radon domain giving the slope and
interception of the line respectively. Fig. 10 shows a powerful line detection
capability possessed by SIDRT.

In classical case, performing SIDRT has a running time $\Omega(N^{3})$, and
finding the maximum of SIDRT has a time complexity O$(N^{2})$. So, the total
running time of classical SIDRT-based line detection algorithm is
$\Omega(N^{3})$; in contrast, there is a quantum algorithm that can run in
$O(N^{1.5})$ time to execute line detections, by the following proposition:

###### Proposition 9

Let $f$ be an $N\times N$ quantum image that can be prepared in time
O$(T_{in})$ by the unitary $U:|0\rangle\rightarrow|\vec{f}\rangle$. There
exist a quantum algorithm that can detect line in $f$ by using SIDRT in time
$O(\frac{N^{1.5}T_{in}}{\epsilon}\emph{polylog N})$, where $\epsilon$ is the
desired precision. This algorithm outputs a pair $(\theta_{0},l_{0})$ such
that
$\tilde{P}_{k_{\theta_{0}}}(l_{0})\geq(1-\frac{2}{\sqrt{3}}\epsilon)\displaystyle\max_{\theta,l}P_{k_{\theta}}(l)$.

###### Proof:

Applying Theorem 3 to the input state
$\sum_{\theta,l\in[N]}\frac{1}{N}|\theta\rangle|l\rangle|0\rangle$ with
setting precision $\frac{\epsilon}{N^{0.5}}$, one can prepare the following
state in time $O(\frac{N^{0.5}T_{in}\text{polylog N}}{\epsilon})$:

$\displaystyle\sum_{\theta,l\in[N]}\frac{1}{N}|\theta\rangle|l\rangle|\tilde{P}_{k_{\theta}}(l)\rangle,$
(67)

where
$|\tilde{P}_{k_{\theta}}(l)-P_{k_{\theta}}(l)|\leq\frac{\epsilon}{N^{0.5}}$.
Then by the quantum algorithm for finding the maximum among $N$ terms with
query complexity O$(\sqrt{N})$ [50], one can find $(l_{0},\theta_{0})$ such
that
$\tilde{P}_{k_{\theta_{0}}}(l_{0})=\max\limits_{\theta,l\in[N]}\tilde{P}_{k_{\theta}}(l)$
in time O$(\frac{N^{1.5}T_{in}\text{polylog N}}{\epsilon})$.

By Proposition 7, in the average case, it holds that
$\displaystyle\max_{\theta,l}P_{k_{\theta}}(l)\geq\frac{\sqrt{3}}{2N^{0.5}}$.
Now that

$\displaystyle|\tilde{P}_{k_{\theta}}(l)-P_{k_{\theta}}(l)|\leq\frac{\epsilon}{N^{0.5}}\leq\frac{2\epsilon\max
P_{k_{\theta}}(l)}{\sqrt{3}},$ (68)

it holds that

$\displaystyle\tilde{P}_{k_{\theta_{0}}}(l_{0})=\max\tilde{P}_{k_{\theta}}(l)\geq(1-\frac{2\epsilon}{\sqrt{3}})\max
P_{k_{\theta}}(l).$ (69)

∎

Figure 9: Denoise using QRT$+$Haar wavelet with setting threshold value
$\infty$ (i.e., make all wavelet coefficients zeros). The original and noisy
image are as shown in Fig. 8. By denoising, the SNR increases 2.674. The ratio
of ($L^{2}$-norm of) wavelet coefficients to filter coefficients is 0.0422,
which means the probability of success of the quantum denoising algorithm on
Section VI-A is higher than 95$\%$ for this example.

Figure 10: Line detection using SIDRT. The upper is the detected image, where
a distinct straight line connecting (228,53) and (97,217) is of the gradient
approximately -1.251 and $y$-interception $(\frac{53}{1.251}+228\mod
256\approx)\ 14.366$. The bottom shows the SIDRT of the upper image. The SIDRT
reaches maximum at point (13,246), which implies there may be a straight line
of the $y$-interception $13$ and gradient
$\frac{1}{\tan(\frac{256-246}{256}-\frac{1}{4})\pi}\approx-1.284$ in the upper
image. This detection result is in good agreement with the observation.

##  VII Conclusion

This paper presents a novel discrete Radon transform for efficient quantum
implementation. By theoretical analysis and numerical experiments, it is shown
that our new proposed quantum Radon transform has similar functionality to the
classical PDRT, and is more suitable for quantum implementation than PDRT. In
addition to its quantum advantage, the QRT has a classical application value
that it can be used as a ‘reversible PDRT’.

Also, a polynomially fast quantum implementation of another interpolation-
based kind of DRT is given. There are two problems that deserve further
investigation: 1. is it possible to provide exponential speedup for performing
interpolated-based DRTs or line detection? 2. Considering that the
computational basis state of the maximum amplitude is most likely to be
observed by quantum measurement, can this property be used to design a better
line detection algorithm?

## Acknowledgment

This work is supported by Chinese Postdoctoral Science Foundation Grant No.
2020M680716, China National Key Research and Development Projects
2020YFA0712300, 2018YFA0704705, and National Natural Science Foundation of
China (Grant Nos. 11471040 and 11761131002).

## References

  * [1] J. Radon, “Über die bestimmung von funktionen durch ihre integralwerte längs gewisser mannigfaltigkeiten,” 1917.
  * [2] S. R. Deans, _The Radon transform and some of its applications_. Courier Corporation, 2007.
  * [3] M. N. Do and M. Vetterli, “Image denoising using orthonormal finite ridgelet transform,” in _Wavelet applications in signal and image processing VIII_ , vol. 4119. International Society for Optics and Photonics, 2000, pp. 831–842.
  * [4] G. Beylkin, “Discrete Radon transform,” _IEEE transactions on acoustics, speech, and signal processing_ , vol. 35, no. 2, pp. 162–172, 1987\.
  * [5] P. A. Toft, “The Radon transform-theory and implementation,” 1996.
  * [6] B. T. Kelley and V. K. Madisetti, “The fast discrete Radon transform. i. theory,” _IEEE Transactions on Image Processing_ , vol. 2, no. 3, pp. 382–400, 1993.
  * [7] D. P. Lun _et al._ , “Discrete periodic radon transform based weighted nuclear norm minimization for image denoising,” in _2017 Fifth International Symposium on Computing and Networking (CANDAR)_. IEEE, 2017, pp. 395–400.
  * [8] C. Carranza, M. Pattichis, and D. Llamocca, “Fast and parallel computation of the discrete periodic radon transform on gpus, multicore cpus and fpgas,” in _2018 25th IEEE International Conference on Image Processing (ICIP)_. IEEE, 2018, pp. 4158–4162.
  * [9] F. Matus and J. Flusser, “Image representation via a finite Radon transform,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 15, no. 10, pp. 996–1006, 1993.
  * [10] T. Hsung, D. P. Lun, and W. Siu, “The discrete periodic Radon transform,” _IEEE Transactions on Signal Processing_ , vol. 44, no. 10, pp. 2651–2657, 1996.
  * [11] I. Svalbe and D. van der Spek, “Reconstruction of tomographic images using analog projections and the digital Radon transform,” _Linear Algebra and Its Applications_ , vol. 339, no. 1-3, pp. 125–145, 2001.
  * [12] A. Kingston, I. Svalbe _et al._ , “Projective transforms on periodic discrete image arrays,” _Advances in Imaging and Electron Physics_ , vol. 139, p. 76, 2006.
  * [13] T. M. Khanipov, “Computational complexity lower bounds of certain discrete Radon transform approximations,” _arXiv preprint arXiv:1801.01054_ , 2018\.
  * [14] C. Shao, Y. Li, and H. Li, “Quantum algorithm design: Techniques and applications,” _Journal of Systems Science and Complexity_ , vol. 32, no. 1, pp. 375–452, 2019.
  * [15] P. W. Shor, “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,” _SIAM Review_ , vol. 41, no. 2, pp. 303–332, 1999.
  * [16] C. Zalka, “Grover’s quantum searching algorithm is optimal,” _Physical Review A_ , vol. 60, no. 4, p. 2746, 1999.
  * [17] M. E. Nielsen, M. A. Nielsen, and I. L. Chuang, _Quantum Computation and Quantum Information_. Cambridge University Press, 2000.
  * [18] N. Abura’ed, F. S. Khan, and H. Bhaskar, “Advances in the quantum theoretical approach to image processing applications,” _ACM Computing Surveys (CSUR)_ , vol. 49, no. 4, pp. 1–49, 2017.
  * [19] G. Beach, C. Lomont, and C. Cohen, “Quantum image processing (quip),” in _32nd Applied Imagery Pattern Recognition Workshop, 2003. Proceedings._ IEEE, 2003, pp. 39–44.
  * [20] F. Yan, A. M. Iliyasu, and S. E. Venegas-Andraca, “A survey of quantum image representations,” _Quantum Information Processing_ , vol. 15, no. 1, pp. 1–35, 2016.
  * [21] P. Q. Le, F. Dong, and K. Hirota, “A flexible representation of quantum images for polynomial preparation, image compression, and processing operations,” _Quantum Information Processing_ , vol. 10, no. 1, pp. 63–84, 2011.
  * [22] H. S. Li, Q. Zhu, R. G. Zhou, L. Song, and X. J. Yang, “Multi-dimensional color image storage and retrieval for a normal arbitrary quantum superposition state,” _Quantum Information Processing_ , vol. 13, no. 4, pp. 991–1011, 2014.
  * [23] S. E. Venegas-Andraca and S. Bose, “Storing, processing, and retrieving an image using quantum mechanics,” in _Quantum Information and Computation_ , vol. 5105. International Society for Optics and Photonics, 2003, pp. 137–147.
  * [24] J. I. Latorre, “Image compression and entanglement,” _arXiv preprint quant-ph/0510031_ , 2005.
  * [25] Y. Zhang, K. Lu, Y. Gao, and M. Wang, “NEQR: a novel enhanced quantum representation of digital images,” _Quantum Information Processing_ , vol. 12, no. 8, pp. 2833–2860, 2013.
  * [26] A. Fijany and C. P. Williams, “Quantum wavelet transforms: Fast algorithms and complete circuits,” in _NASA International Conference on Quantum Computing and Quantum Communications_. Springer, 1998, pp. 10–33.
  * [27] Y. Zhang, K. Lu, and Y. Gao, “Qsobel: a novel quantum image edge extraction algorithm,” _Science China Information Sciences_ , vol. 58, no. 1, pp. 1–13, 2015.
  * [28] X.-W. Yao, H. Wang, Z. Liao, M.-C. Chen, J. Pan, J. Li, K. Zhang, X. Lin, Z. Wang, Z. Luo _et al._ , “Quantum image processing and its application to edge detection: theory and experiment,” _Physical Review X_ , vol. 7, no. 3, p. 031041, 2017.
  * [29] S. Caraiman and V. I. Manta, “Histogram-based segmentation of quantum images,” _Theoretical Computer Science_ , vol. 529, pp. 46–60, 2014.
  * [30] X. H. Song, S. Wang, S. Liu, A. A. A. El Latif, and X. M. Niu, “A dynamic watermarking scheme for quantum images using quantum wavelet transform,” _Quantum information processing_ , vol. 12, no. 12, pp. 3689–3706, 2013.
  * [31] H. S. Li, X. Chen, S. Song, Z. Liao, and J. Fang, “A block-based quantum image scrambling for gneqr,” _IEEE Access_ , vol. 7, pp. 138 233–138 243, 2019\.
  * [32] H. J. Nussbaumer, “The fast Fourier transform,” in _Fast Fourier Transform and Convolution Algorithms_. Springer, 1981, pp. 80–111.
  * [33] A. Kingston and I. Svalbe, “Generalised finite Radon transform for n$\times$ n images,” _Image and Vision Computing_ , vol. 25, no. 10, pp. 1620–1630, 2007.
  * [34] A. Kingston, “Orthogonal discrete Radon transform over $p^{n}\times p^{n}$ images,” _Signal processing_ , vol. 86, no. 8, pp. 2040–2050, 2006.
  * [35] A. Buades, B. Coll, and J.-M. Morel, “A review of image denoising algorithms, with a new one,” _Multiscale Modeling & Simulation_, vol. 4, no. 2, pp. 490–530, 2005.
  * [36] D. L. Donoho and J. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” _Biometrika_ , vol. 81, no. 3, pp. 425–455, 1994.
  * [37] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” _Physical Review Letters_ , vol. 103, no. 15, p. 150502, 2009\.
  * [38] I. Kerenidis and A. Prakash, “Quantum recommendation systems,” in _8th Innovations in Theoretical Computer Science Conference (ITCS 2017)_. Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017.
  * [39] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, “Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics,” in _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_. ACM, 2019, pp. 193–204.
  * [40] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, “Quantum machine learning,” _Nature_ , vol. 549, no. 7671, pp. 195–202, 2017.
  * [41] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, “Preconditioned quantum linear system algorithm,” _Physical Review Letters_ , vol. 110, no. 25, p. 250504, 2013.
  * [42] C. Shao, “From linear combination of quantum states to grover’s searching algorithm,” _arXiv preprint arXiv:1807.09693_ , 2018.
  * [43] ——, “A quantum model for multilayer perceptron,” _arXiv preprint arXiv:1808.10561_ , 2018.
  * [44] L. Ruiz Perez and J. C. Garcia Escartin, “Quantum arithmetic with the quantum Fourier transform,” _Quantum Information Processing_ , vol. 16, no. 6, p. 152, 2017.
  * [45] C. Wang and L. Wossnig, “A quantum algorithm for simulating non-sparse Hamiltonians,” _arXiv preprint arXiv:1803.08273_ , 2018.
  * [46] A. Averbuch, R. Coifman, D. Donoho, M. Israeli, and J. Walden, “Fast slant stack: A notion of Radon transform for data in a cartesian grid which is rapidly computible, algebraically exact, geometrically faithful and invertible,” _SIAM Scientific Computing_ , 2001.
  * [47] D. L. Donoho and A. G. Flesia, “Digital ridgelet transform based on true ridge functions,” in _Studies in Computational Mathematics_. Elsevier, 2003, vol. 10, pp. 1–30.
  * [48] I. L. Markov and M. Saeedi, “Constant-optimized quantum circuits for modular multiplication and exponentiation,” _arXiv preprint arXiv:1202.6614_ , 2012\.
  * [49] S. Catterall, “Central limit theorem for square roots of sums of i.i.d. random variables,” https://stats.stackexchange.com/questions/241504/central-limit-theorem-for-square-roots-of-sums-of-i-i-d-random-variables; accessed 19-May-2021.
  * [50] A. Ahuja and S. Kapoor, “A quantum algorithm for finding the maximum,” _arXiv preprint quant-ph/9911082_ , 1999.