1.   引言

1982年物理学家理冯曼首次提出了量子计算的概念^[1]，这一概念为科技发展提供了新思路，然而因其对硬件的要求过高而落入低谷。直到20世纪90年代，大数质因子分解和无序数据库搜索的量子算法的相继提出^[2,3]，验证了量子并行计算的强大能力，该领域理论及应用研究才得以深入展开。

1997年俄罗斯学者研究了量子图像识别问题^[4]，随后，2003年文献[5]正式提出了量子图像处理的概念。量子图像处理的学科发展概貌可以参阅文献[6]。量子图像处理包括量子图像描述和量子图像处理算法^[7]。在量子图像描述方面，主要包括：量子比特阵列模型^[8]、实矢量模型^[9]、纠缠图像模型^[10]、量子图像灵活描述模型(Flexible Representation of Quantum Images, FRQI)^[11]以及新颖的增强描述(Novel Enhanced Quantum Representation, NEQR)模型^[12]。在量子图像处理算法方面，目前研究成果主要集中在：几何变换^[13]、色彩处理^[14]、图像置乱^[15]、特征提取^[16]、图像分割^[17]、信息隐藏^[18]、图像加密^[19]和均值滤波^[20]等。此外，量子计算与密码学的融合研究也是一个值得关注的方向^[21,22]。

虽然量子图像处理的相关研究逐步深入，但整体上还处于起步阶段，且相关方向的发展呈不平衡状态，在量子图像加密和量子图像水印^[23]方面，相关成果较为丰硕，而量子图像滤波方面的成果相对较少。近两年来陆续出现了一些量子图像空间滤波方法，例如，文献[24,25]提出了空间域量子图像滤波的一般框架，文献[26]采用冒泡排序计算像素中值，提出了量子图像上午中值滤波方法。

本文提出量子中值滤波方案，对于中值计算采用了不同的方法，与现有方法的区别在于该方案的量子线路使用了较少的基本模块，从而具有更低的复杂度。另外，本文验证了所提方案在图像降噪方面的应用，经典计算机上的仿真结果表明，本方案与经典方法的降噪效果差别不大，但在将来的量子计算机上，借助量子计算的并行性，本方案可以实现对经典计算的指数加速。

2.   预备知识

2.1   量子图像的NEQR描述

该描述是将一幅经典图像转换为量子图像时所采用的一种方法。根据NEQR描述，一个${2^n} \times {2^n}$的图像$I$可以表示为

其中$|{c_i}\rangle = |c_i^{q - 1} ,c_i^{q - 2} ,··· c_i^0\rangle ,\;\;c_i^k \in \{ 0,\;1\} ,\;\;k = q - 1, q - 2, ··· , \;0$, $|i\rangle = |y\rangle |x\rangle = |{y_{n - 1}},{y_{n - 2}}, ··· {y_0}\rangle |{x_{n - 1}} , {x_{n - 2}}, ··· {x_0}\rangle , \;{y_i},\;{x_i} \in \{ 0,\;1\}$。

在NEQR中，量子图像由颜色信息和位置信息共同表示。颜色由二值序列$|c_i^{q - 1},c_i^{q - 2} , ··· c_i^1c_i^0\rangle$来表示，灰度图像颜色值的范围是0～255。位置则由$|x\rangle$和$|y\rangle$表示，分别由$n$个量子比特组成。

图1是一个NEQR模型的例子，该图像的大小是$2 \times 2$，每一个方块代表一个像素，每个方块中的数值是颜色值。它的NEQR描述如式(2)所示

图 1  一幅2×2的灰度图像

下载: 全尺寸图片 幻灯片

2.2   量子比较器

量子比较器是量子计算的基本模块，用来比较两个无符号整数$|x\rangle$和$|y\rangle$的大小，$|{e_0}\rangle$和$|{e_1}\rangle$为量子比较器的输出。若$|{e_1}{e_0}\rangle = |00\rangle$，则 $|x\rangle = |y\rangle$，若$|{e_0}\rangle = |1\rangle$，则$|x\rangle < |y\rangle$。量子比较器的具体线路如图2左所示，线路符号如图2右所示^[27]。

图 2  量子比较器线路

下载: 全尺寸图片 幻灯片

2.3   量子模加法器

量子模加法器用来执行$|a\rangle$和$|b\rangle$两个无符号整数的模加法，量子线路如图3所示。由于量子计算的幺正性，该线路是可逆的，左右逆转该线路的方向即可实现两个无符号整数的模减法。加法器和减法器的矩阵描述互为共轭转置^[28]。

图 3  量子模加法器线路

下载: 全尺寸图片 幻灯片

注意，对于模加法器，Adder符号中的黑色线条的位置在右侧，根据量子计算的可逆性，当Adder符号中的黑色线条位于左侧时，即为量子模减法器。

3.   量子图像的中值滤波设计

为便于描述，本文采用模块分解的方法，将量子中值滤波线路分解为若干子模块，最后再由这些子模块组合成总线路。下面首先介绍基本模块的量子实现。

3.1   复制模块

对于两个$q$比特的量子寄存器$|{c_1}\rangle ,\;|{c_2}\rangle$，其中$|{c_1}\rangle = |{c_{1,1}}\rangle,{c_{1,2}}\rangle, ··· |{c_{1,q}}\rangle$, $|{c_2}\rangle = |{c_{2,1}}\rangle,{c_{2,2}}\rangle,···|{c_{2,q}}\rangle$=$|0{\rangle ^{ \otimes q}}$，该模块的功能是将$|{c_1}\rangle$复制到$|{c_2}\rangle$中。显然该模块可用$q$个受控非门实现，量子线路如图4所示。

图 4  复制模块的量子线路

下载: 全尺寸图片 幻灯片

3.2   循环平移模块

该模块的操作对象是量子图像的位置比特，包括循环左移、右移、上移、下移。显然该模块可采用模加法器和模减法器实现，具体线路图如图5所示。

图 5  复制模块的量子线路

下载: 全尺寸图片 幻灯片

3.3   中值计算模块

为便于描述，本文以$3 \times 3$滤波窗口为例，考察中值计算的量子线路设计方法，对于其它幅度的滤波窗口，可依此类推。具体采用循环比较的方式，设计求取9个像素中值的量子线路。设目标像素为${c_1},{c_2},{c_3},{c_4},{c_5},{c_6},{c_7},{c_8},{c_9}$。具体算法如表1所示。

表 1  循环比较算法

循环比较算法的具体实现
循环($i = 2\;\;{\rm{to}}\;\;9$)
如果$i < 9$，则$k = i - 1$，否则$k = i - 4$
循环($j = 1\;\;{\rm{to}}\;\;k$)
比较${c_j},{c_i}$，若${c_j} < {c_i}$，则交换${c_j},{c_i}$

下载: 导出CSV 

| 显示表格

在该算法中，当$i < 9$时，可以实现${c_1},{c_2}, ··· ,{c_i}$的降序排序；当$i = 9$时，可以获得${c_1},{c_2}, ··· ,{c_9}$中降序排列的前5个值，此时第5个值即为中值。

各轮循环很容易采用量子比较器和量子交换门实现，当$i = 2,\;3,\; ··· \;,\;9$时，以${C_2},\;{C_3},\; ··· \;,\;{C_9}$标记各轮排序的量子线路，具体如图6所示，为简便，图中仅给出了${C_2},\;{C_3},\;{C_9}$的量子线路。

图 6  C2, C3, C9模块的量子线路

下载: 全尺寸图片 幻灯片

将上述8个模块${C_2},\;{C_3},\; ··· \;,\;{C_9}$依次作用到像素比特${c_1},{c_2},{c_3},{c_4},{c_5},{c_6},{c_7},{c_8},{c_9}$上，即可获得该9个像素的中值，具体的量子线路如图7所示。

图 7  中值计算模块的量子线路

下载: 全尺寸图片 幻灯片

3.4   总体线路设计

量子图像中值滤波总线路的设计思想是，充分利用量子计算并行性，使所有像素的滤波操作同步进行。基于这个设计思想，本文采取的总线路的设计方法简述如下。

(1) 制备8幅空图像。首先采用Hadamard门${{H}}$和恒等门${{I}}$生成8幅与原始图像$|{I_{11}}\rangle \!=\! \dfrac{1}{{{2^n}}}\displaystyle\sum\nolimits_{{y_{11}} = 0}^{{2^n} - 1} \displaystyle\sum\nolimits_{{x_{11}} = 0}^{{2^n} - 1} {|{c_{11}}\rangle |{y_{11}}\rangle |{x_{11}}\rangle }$同大小的8幅空图像$|{I_{ij}}\rangle = \dfrac{1}{{{2^n}}}\displaystyle\sum\nolimits_{{y_{11}} = 0}^{{2^n} - 1} {\displaystyle\sum\nolimits_{{x_{11}} = 0}^{{2^n} - 1} {|0{\rangle ^{ \otimes q}}|{y_{ij}}\rangle |{x_{ij}}\rangle } }$，其中$(ij) \in \{ (12), \;(13),\;(21), \;(22), \;(23),\;(31),\;(32),\;(33)\}$，如图8左侧上下两个虚线框中线路所示，Hadamard门${{H}}$和恒等门${{I}}$的矩阵描述为

图 8  中值计算模块的量子线路

下载: 全尺寸图片 幻灯片

(2) 图像复制。采用比较模块和复制模块，将原始图像的像素信息复制到8幅空图像中，从而得到与原始图像完全相同的8幅图像，如图8中左数第2个虚线框中线路所示。具体实现方法是，首先采用比较器比较9幅图像的位置，若位置相等就将原始图像中该位置像素值复制给空图像。

(3) 图像平移。保持图像$|{I_{22}}\rangle$不动，采用循环平移模块将其它8幅图像$|{I_{11}}\rangle$, $|{I_{12}}\rangle$, $|{I_{13}}\rangle$, $|{I_{21}}\rangle$,$|{I_{23}}\rangle$, $|{I_{31}}\rangle$, $|{I_{32}}\rangle$, $|{I_{33}}\rangle$分别向8个方向(左上、上、右上、左、右、左下、下、右下)循环平移一个单位，具体实现方法如图8中左数第3个虚线框中线路所示。特别值得指出，对于平移后的9幅图像，处于相同位置的9个像素，就是原始图像中被$3 \times 3$中值滤波罩盖住的那9个像素。此时，利用量子计算的并行性，只需对9幅图像的像素值比特实施一次中值计算，即可获得原始图像中所有像素的$3 \times 3$中值滤波结果。

(4) 中值滤波。如前所述，采用比较器和中值计算模块即可获得滤波结果。具体说来，首先采用比较器比较9幅平移后图像中的像素位置，当位置相等时，采用中值计算模块对9幅图像的像素值比特实施中值计算，由于计算得出的中值存储在$|{c_{22}}\rangle$中，所以此时，量子图像$|{I_{22}}\rangle = \dfrac{1}{{{2^n}}}\displaystyle\sum\nolimits_{{y_{22}} = 0}^{{2^n} - 1} {\displaystyle\sum\nolimits_{{x_{22}} = 0}^{{2^n} - 1} {|{c_{22}}\rangle |{y_{22}}\rangle |{x_{22}}\rangle } }$即为中值滤波后的图像，具体实现如图8中最右侧虚线框中线路所示。

4.   量复杂度分析

为便于描述，以下分析均假定图像大小为${2^n} \times {2^n}$，灰度值范围为$\{ 0,\;1,\; ··· \;,\;{2^q}\}$。

(1) 量子逻辑门的复杂度。在量子图像处理中，时间复杂度主要依据通用逻辑门(Hadamard门、非门、受控非门(Controlled -NOT gate, CNOT))的数量。为便于描述，记${\varLambda _n}({{X}})$为$n + 1$比特受控非门，其中前$n$比特为控制比特，第$n + 1$比特为目标比特，${{X}} = \left( \begin{array}{l} 0\quad 1 \\ 1\quad 0 \\ \end{array} \right)$为非门，例如CNOT为${\varLambda _1}({{X}})$, Toffoli门为${\varLambda _2}({{X}})$。文献[29]给出了几个复杂量子门与通用量子门之间的构造关系。

(a) 一个Toffoli门等价于6个CNOT门。

(b) 一个${\varLambda _n}({{X}})$门($n \ge 3$)等价于$2(n - 1)$个Toffoli门和1个CONT门。

(c) 一个交换门等价于3个CONT门。

由此可知，Toffoli门的复杂度是6，${\varLambda _n}({{X}})$门的复杂度为$12n - 11$，交换门的复杂度为3。

(2) 量子比较器和量子加法器的复杂度。由图2可知，组成量子比较器的量子门共有$n$组，分别表示为${a_0},\;{a_2},\; ··· \;,\;{a_{n - 1}}$。对于第$i$组${a_i}(0 \le i \le n - 2)$，有4个${\varLambda _{2(i + 1)}}({{X}})$门，最后一组${a_{n - 1}}$有2个${\varLambda _{2n}}({{X}})$。所以，比较器的复杂度为$4\left( {\displaystyle\sum\nolimits_{i = 0}^{n - 2} {12[2(i + 1)] - 11} } \right)$$+2[12(2n) - 11] = 48{n^2} - 44n + 22$

由图3可知，量子加法器由$n$个求和模块S和$2n - 2$个进位模块C组成，且模块S和模块C的复杂度分别为2和13。因此，量子加法器的复杂度为$2n + 13(2n - 2) = 28n - 26$。

(3) 中值计算模块的复杂度。如前所述，一个交换门等价于3个受控非门${\varLambda _1}({{X}})$，所以，1个${\varLambda _1}({\rm{Swap}})$等价于3个${\varLambda _2}({{X}})$。由图6可知，排序模块${C_i}(2 \le i \le 8)$包含$(i - 1)q$个受控交换门${\varLambda _1}({\rm{Swap}})$，${C_9}$包含$5q$个${\varLambda _1}({\rm{Swap}})$。根据图7，中值计算共包含$\displaystyle\sum\nolimits_{i = 2}^8 {(i - 1)q} + 5q = 33q$个${\varLambda _1}({\rm{Swap}})$，等价于$99$个${\varLambda _2}({{X}})$。

(4) 中值滤波总线路的复杂度。根据文献[30]，单比特量子门的复杂度远远低于二比特量子门，所以在计算量子线路的复杂度时一般忽略单比特量子门。根据图8，线路包括32个比较器，12个平移模块，8个受控复制模块${\varLambda _4}({\rm{Copy}})$，1个受控中值计算模块${\varLambda _{32}}({\rm{Median}})$。从图5可知平移模块的复杂度等于加法模块。1个${\varLambda _4}({\rm{Copy}})$等价于$q$个${\varLambda _5}({{X}})$。一个${\varLambda _{32}}({\rm{Median}})$等价于$33q$个${\varLambda _{33}}({\rm{Swap}})$，进一步等价于$99q$个${\varLambda _{34}}({{X}})$门。因为${\varLambda _5}({{X}})$复杂度为$12 \times 5 - 11 = 49$, ${\varLambda _{34}}({{X}})$的复杂度为$12 \times 34 - 11 = 397$，所以整个中值滤波线路的复杂度为

上述结果表明，量子中值滤波方案的计算复杂度只是图像规模(大小)的2次函数。对与经典中值滤波，由于不能执行并行计算，所以滤波操作必须逐个像素执行，因此其复杂度必然包含因子${2^{2n}}$，单就这个因子而言，经典方案的计算复杂度必然为图像规模的指数函数。因此，本文所提中值滤波方案可以实现对经典计算的加速。

5.   经典计算机上的仿真

由于目前量子计算机尚未普及，因此本仿真在经典计算机上进行，其中台式机配置为Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz 4.00 GB RAM，软件环境为Win7操作系统和Matlab (2014a)，虽不能验证量子计算的并行性，但能验证方案执行后的效果。

仿真方案为，首先在原始图像中分别加入椒盐、高斯、泊松3种噪声，然后分别采用本文所提量子中值滤波和经典中值滤波考察滤波降噪的效果。进一步我们还考察了量子中值滤波方案对于量子比特翻转噪声的滤除效果。仿真采用的5幅$512 \times 512$的灰度图像如图9所示。

图 9  仿真中用到的5幅灰度图像

下载: 全尺寸图片 幻灯片

为显示本文所提方案的优势，本文采用对比两种方案滤波前后图像峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)的方法，PSNR的定义如式(5)所示

其中，$I'(i,j)$和$I(i,j)$分别表示噪声图像和原始图像。$m$和$n$分别表示图像中每行和每列的像素数。

5幅图像滤波前后的峰值信噪比对比如表2所示，其中每种方案的左边一列是降噪后的PSNR，右边一列是与降噪之前相比PSNR的提高值。

表 2  两种方案滤波前后的峰值信噪比对比(dB)

图像	椒盐噪声						高斯噪声						泊松噪声
	滤波前	经典方案		本文方案			滤波前	经典方案		本文方案			滤波前	经典方案		本文方案
(a)	14.45	32.28	17.83	34.35	19.90		18.48	28.95	10.47	29.27	10.79		27.66	35.50	7.84	35.63	7.97
(b)	14.73	32.72	17.99	34.87	20.14		18.08	28.91	10.83	29.19	11.11		27.23	35.74	8.51	35.84	8.61
(c)	14.75	30.86	16.11	32.30	17.55		17.91	28.51	10.60	28.78	10.87		25.87	33.42	7.55	33.52	7.65
(d)	14.83	30.69	15.86	31.93	17.10		17.80	28.07	10.27	28.34	10.54		25.62	32.70	7.08	32.80	7.18
(e)	14.68	31.55	16.87	33.20	18.52		17.97	28.66	10.69	28.95	10.98		26.68	34.55	7.87	34.67	7.99
平均	14.69	31.62	16.93	33.33	18.64		18.05	28.62	10.57	28.91	10.86		26.61	34.38	7.77	34.49	7.88

下载: 导出CSV 

| 显示表格

从表2可知，本文方案的降噪效果略优于经典方案，对于椒盐、高斯、泊松3种噪声，本文方案滤波后的PSNR分别比经典方案平均提高1.71 dB, 0.29 dB, 0.11 dB。另外，与3种噪声图像降噪之前的PSNR比较，两种滤波方案均分别提升17 dB, 10 dB, 7 dB以上。这说明中值滤波方法比较适合于椒盐噪声的降噪。关于本文方案略优于经典方案的原因，是因为对于图像的边缘区域，中值的计算方法略有不同，由于滤波罩不能移出图像之外，所以经典方法一般忽略边缘像素；而本文采用循环移位的方法完整地处理了包括边缘区域的所有像素，从而导致滤波后图像的PSNR相比经典方法略有提升。

对于量子图像处理，当量子图像在量子信道中传输时，有时会发生量子比特翻转，由此导致的信息失真称为量子比特翻转噪声。下面考察本文方法对于量子比特翻转噪声的滤除效果，具体方法是，首先使$q$个灰度值比特分别按某一概率阈值翻转，翻转后的图像即为混入量子翻转噪声的图像。然后采用本文方法滤波去噪，并统计滤波前后图像的PSNR。

对于图9中的5幅图像，当概率阈值分别取0.05, 0.06, 0.08和0.10时，本文方案对于量子噪声图像滤波前后的PSNR对比如表3所示，其中“滤波后”中的左边一列是滤波后的PSNR，右边一列是与“滤波前”相比PSNR的提高值。

表 3  量子噪声图像滤波前后的峰值信噪比对比(dB)

图像	概率阈值0.05				概率阈值0.06				概率阈值0.08				概率阈值0.10
	滤波前	滤波后			滤波前	滤波后			滤波前	滤波后			滤波前	滤波后
(a)	18.89	36.07	17.18		18.14	35.36	17.22		17.07	33.90	16.83		16.19	32.27	16.08
(b)	18.88	36.36	17.48		18.15	35.88	17.73		17.03	34.42	17.39		16.16	33.00	16.84
(c)	18.63	33.29	14.66		17.96	33.03	15.07		16.91	32.36	15.45		16.07	31.42	15.35
(d)	18.64	33.70	15.06		17.98	32.41	14.43		16.91	31.78	14.87		16.07	30.80	14.73
(e)	18.79	35.52	16.73		18.10	34.24	16.14		16.93	33.19	16.26		16.09	31.99	15.9
平均	18.77	34.99	16.22		18.07	34.18	16.12		16.97	33.13	16.16		16.12	31.90	15.78

下载: 导出CSV 

| 显示表格

从表3可知，随着概率阈值的提升，加入的噪声也增加，降噪效果呈下降趋势。但对于表3中给出的4种概率阈值，滤波后图像的PSNR与滤波前相比，均有16 dB以上的提高。仿真结果验证了本文方案对于量子噪声有较好的降噪能力。以图9(b)为例，在表3中的4种概率阈值下，噪声图像分别如图10(a)—图10(d)所示，降噪效果分别如图11(a)—图11(d)所示。

图 10  不同概率阈值下的量子比特翻转噪声图像

下载: 全尺寸图片 幻灯片

图 11  不同概率阈值下量子噪声图像的滤波效果

下载: 全尺寸图片 幻灯片

从图11可知，随着概率阈值的增大，提出方案的降噪效果逐渐减弱，这与表3中的数值结果是一致的。关于本文方案适合滤除量子比特翻转噪声的仿真结果，我们给出如下解释。

本文方案中，量子图像采用NEQR描述。量子比特翻转噪声，本质上是二进制像素值每个二进制位的随机翻转，这种翻转具有脉冲噪声的效果。中值滤波是一种非线性的图像平滑方法，它能够很好地滤除脉冲噪声，同时又能够保护目标图像边缘。量子中值滤波是经典中值滤波的量子版本，具有经典中值滤波的优点。因此量子图像中值滤波适合滤除量子比特翻转噪声。

6.   结束语

针对量子图像处理中的中值滤波问题，本文研究了基于量子计算的设计方案。从设计若干基本模块的量子线路入手，通过组合基本模块，逐渐得出了总体线路的设计方法。与现有文献的不同之处在于，本文采用了新的中值计算方法。复杂性分析表明，本方案可以实现对经典方案的加速。仿真结果表明，本方案适合于经典椒盐噪声和量子比特翻转噪声图像的降噪处理。目前制约量子图像处理发展的问题之一是基本数值运算(中值、均值、均方差、协方差)的量子实现，如何设计这些基本运算的量子实现方法是我们下一步重点研究的问题。