A hybrid color and grayscale images encryption scheme based on the quaternion Hartley transform (QHT), the twodimensional (2D) logistic map, the double random phase encoding (DRPE) in gyrator transform (GT) domain and the threestep phaseshifting interferometry (PSI) is presented. First, we propose a new color image processing tool termed as the quaternion Hartley transform, and we develop an efficient method to calculate the QHT of a quaternion matrix. In the presented encryption scheme, the original color and grayscale images are represented by quaternion algebra and processed holistically in a vector manner using QHT. To enhance the security level, a 2D logistic mapbased scrambling technique is designed to permute the complex amplitude, which is formed by the components of the QHTtransformed original images. Subsequently, the scrambled data is encoded by the GTbased DRPE system. For the convenience of storage and transmission, the resulting encrypted signal is recorded as the realvalued interferograms using threestep PSI. The parameters of the scrambling method, the GT orders and the two random phase masks form the keys for decryption of the secret images. Simulation results demonstrate that the proposed scheme has high security level and certain robustness against data loss, noise disturbance and some attacks such as chosen plaintext attack.
I. INTRODUCTION
Optical technology for security applications has received increasing interest in the last few decades
[1

3]
. In particular, the double random phase encoding (DRPE)
[1]
technique, in which the image is encoded to be a white noise pattern with two statistically independent random phase masks (RPM), has attracted much attention. The DRPE encryption technique is constructed by using two RPMs as keys: one is located at the input domain and the other is at the Fourier domain. To enlarge the key space, it has been further extended from the Fourier transform domain to the Fresnel transform domain
[4
,
5]
, the fractional Fourier transform domain
[6]
, the Gyrator transform (GT) domain
[7]
and the Hartley transform (HT) domain
[8]
, etc. In these encryption systems, the parameters (e.g. the propagation distance, the fractional order, etc) are introduced and served as extra keys. Since the HT is mathematically equivalent to the Fourier transform but is purely real, it has a computational advantage over the aforementioned transforms for not having to manage the real and imaginary parts
[9]
. Though most of the DRPEbased encryption techniques are quite robust and secure, the recovered image loses its color information, which is useful in image processing and practical applications, and in decryption since the original image is illuminated with monochromatic light
[10]
.
Using the DRPE technique, Zhang and Karim reported a singlechannel encryption algorithm for color image in the Fourier domain
[11]
. Chen
et al
. suggested a method using coherent diffractive imaging for optical colorimage encryption and synthesis in the Fresnel domain
[3]
. Joshi
et al
.
[12]
have proposed an approach for the encryption of color images via using a fractional Fourier transform. By use of the Arnold transform and an interference method, a color image encryption technique is proposed in
[13]
. A method for securing color image based on the Arnold transform and the gyrator transform is reported in
[10]
. Based on polarization encoding, Nishchal
et al
. proposed two asymmetric color image cryptosystem techniques
[14

15]
. Li and Zhao
[16]
designed a color image encryption based on the fractional HT. However, the aforementioned encryption methods for color image ignore the potential interactions between different color channels
[17
,
18]
. In these methods, a color image is separated into the three color channels, and encryption is carried out in each channel independently. Recently, quaternion theory
[20]
has been used in DRPEbased encryption. Based on the discrete quaternion Fourier transform (DQFT) and DRPE, a color image encryption method is reported
[17]
. With an iterative phase retrieval algorithm, an algorithm to encrypt double color images in quaternion gyrator domain is described
[18]
. Using quaternion representation, the color image is processed holistically as a vector field rather than as separated color components in both of the two methods. However, the method in
[17]
is vulnerable to some attacks (such as chosenplaintext attack, etc) caused by the linearity of the cryptosystem, while the algorithm in
[18]
has high computational cost due to the phase retrieval iteration process.
In this paper, based on the quaternion Hartley transform (QHT), logisticbased scrambling technique, GTbased DRPE and threestep phaseshifting interferometry (PSI), a novel hybrid color and grayscale image encryption method is presented. In the proposed approach, the QHT is defined and its calculation for a quaternion matrix is developed first. Then a color image and a grayscale image are processed holistically in a vector manner using QHT. To resist an attack such as chosen plaintext attack, the components of the QHTtransformed images are permuted by use of the designed scrambling algorithm based on logistic maps. Since the flexible configuration of the GT based system, which means fixed distances between the generalized lenses and the manipulation of the transformation angle using lens rotation, makes the setup useful for image encryption
[19]
, the GTbased DRPE technique is used to encrypt the images after QHT and permutation. The resulting encrypted signal from the cryptosystem is recorded as interferograms by the threestep PSI. Numerical simulations have been made for demonstrating the feasibility and performance of this encryption.
II. RELATED BACKGROUND
In this section, we come back to some related theories before extending the traditional HT to the quaternion domain, showing how to calculate the QHT of a quaternion matrix and developing the logistic mapsbased scrambling technique, etc.
 2.1. Quaternion Number
Quaternions can be viewed as generalizations of complex numbers. A quaternion number can be represented as follows
[18]
:
where
a, b, c,
and
d
are real numbers, and
i, j,
and
k
are three imaginary units obeying the following rules
The conjugate and modulus of a quaternion are respectively defined by
When
a
= 0,
q
is a pure quaternion.
 2.2. Hartley Transform
In mathematics, the HT, which transforms realvalued functions to realvalued functions, is an integral transform closely related to the Fourier transform. HT has two main properties: one is that it is a real transform; the other is that it and its inverse transform are identical
[21]
. Therefore it can have computational advantages over the Fourier transform
[21]
. The HT
[8]
of a real function
f
(
x
,
y
) is defined as
and the inverse Hartley transform (IHT) is defined as
where
cas
=
cos
+
sin
,
HT
[•] and
IHT
[•] denote the Hartley and inverse Hartley operators.
 2.3. Double Random Phase Encoding in the Gyrator Domain
In mathematics, the GT with respect to parameter
α
, called as a fractional order
[7]
, of a twodimensional function
f
(
x
,
y
) can be expressed as
[7
,
9]
where the
G_{α}
(
p
,
s
) and
f
(
x
_{0}
,
y
_{0}
) are the output and input of the transform, respectively.
GTα
[•] represents gyrator operator. The inverse GT (IGT) corresponds to the GT with respect to 
α
and is given by
where
IGT
[•] denotes the inverse gyrator operator. The GT can be achieved either by using an optical system
[9]
or by using a fast algorithm
[22]
.
In
Fig. 1
, the optical setup of DRPE in the gyrator domain is shown. Three planes are defined as the input plane, the transform plane, and the output plane. The corresponding coordinates of the three planes are denoted by (
x
_{0}
,
y
_{0}
), (
p
,
s
) and (
x
_{1}
,
y
_{1}
), respectively. The two random phase masks
RPM
1 and
RPM
2 are represented by exp[
j2πφ
(
x
_{0}
,
y
_{0}
)] and exp[
j2πϕ
(
p,s
)], where
φ
(
x
_{0}
,
y
_{0}
) and
ϕ
(
p,s
) are two random functions distributed uniformly in the interval [0,1]. Using a system which is composed of three generalized lenses L
_{1}
, L
_{2}
and L
_{1}
with fixed equal distances between them, the GT can be implemented optically. GT at order
α
is reached by rotation of these lenses
[9]
. In
Fig. 1
, the dashed block GT1 which contains lenses L
_{1}
, L
_{2}
and L
_{1}
indicates the first optical GT with respect to
α
_{1}
, and the other GT2 containing lenses L'
_{1}
, L'
_{2}
and L'
_{1}
indicates the second optical GT with respect to
α
_{2}
. When the security system is illuminated perpendicularly by a plane wave, the encrypted data
g
(
x
_{1}
,
y
_{1}
) is obtained at the output plane
[9]
:
Optical setup for the double random phase encoding in the gyrator domain.
By applying IGT to the encrypted data with the conjugates of the
RPM
1 and
RPM
2, the decrypted image can be obtained as follows
[9]
:
 2.4. Phaseshifting Interferometry
Since digital holography
[2]
provides a convenient form of recording the complex encrypted images after passing through the DRPE systems, phaseshifting interferometry is employed to record the complex resulting encrypted signal in the proposed scheme. A variety of PSI techniques have been developed, including threestep, fourstep, etc
[23]
.
Let
A
(
x
,
y
)exp[
jψ
(
x
,
y
)] and
A_{r}
exp(
jδ_{k}
) be the complex amplitude distributions of the object wave in the recording plane and the reference wave in that plane at the
k
th exposure, respectively. Here,
A_{r}
is a constant,
δ_{k}
is the phase shift of the reference wave between two adjacent steps and
k
=1,2,…,
N
. The
k
^{th}
interference pattern
I_{k}
(
x
,
y
) can be represented as
[24]
,
For a known set {
δ_{k}
}(
k
=1,2,...,
N
), a digital hologram
U
(
x
,
y
) can be derived as a function {
I_{k}
} and {
δ_{k}
}
[24]
. In the threeframe case,
N
is 3. When
δ
_{1}
=0,
δ
_{2}
=
π
/2 and
δ
_{3}
=
π
, a digital hologram
U
(
x
,
y
) from the three interferograms
I
_{1}
,
I
_{2}
and
I
_{3}
can be expressed as
[24
,
25]
 2.5. Logistic Map
Chaos theory is an evolutionary theory, which describes that the nonlinear dynamical systems change from ordered state to disordered state
[26]
. The dynamical systems are established based on various chaos functions such as logistic maps, which are extremely sensitive to the initial conditions. These functions generate iterative values which are completely random in nature. In this paper, the twodimensional (2D) logistic map is used to make the change of sequence of image pixels. It is defined as
[26]
The dynamic behavior of 2D logistic is controlled by the parameters
ε
,
λ
_{1}
,
λ
_{2}
and
ρ
shown in Eq. (12). When
ε
=4,
ρ
=0.1, 0.65≤
λ
_{1}
≤0.9 and 0.65≤
λ
_{2}
≤0.9, the dynamical system is in a chaotic state and slight variations of the initial value will yield a drastically different result which is a nonperiodic and nonconverging sequence over time
[26]
.
III. QUATERNION HARTLEY TRANSFORM
 3.1. Definition
Due to the noncommutative multiplication property of quaternions, there are different types of QHT that can be defined. In this work, the leftside QHT (QHT
^{L}
) and the rightside QHT (QHT
^{R}
) are defined:
• Leftside QHT:
• Rightside QHT:
where
f_{q}
(
x
,
y
) is a twodimensional quaternion function and
µ
is a pure unit quaternion which meets the constraint that
µ
^{2}
=1.
QHT^{L}
() and
QHT^{R}
() are the leftside QHT and rightside QHT operations, respectively.
Corresponding to QHT, two forms of the inverse QHT (IQHT) are defined as follows
• Leftside IQHT (IQHTL):
• Rightside IQHT (IQHT
^{R}
):
Here,
IQHT^{L}
() and
IQHT^{R}
() are the inverse leftside QHT and rightside QHT operations, respectively. QHT and IQHT are transformation pairs of each other. They ensure that a quaternion function
f_{q}
(
x
,
y
) which is transformed into the QHT domain can be reconstructed completely by the inverse process without losing any information.
 3.2. QHT Calculation
In this subsection, the method which makes full use of the existing HT algorithm to calculate the QHT of a quaternion matrix is presented. Since the leftside QHT is used in this work, only the calculation method for it is described. By using the HT algorithm, the QHT can be implemented efficiently.
Since
f_{q}
(
x
,
y
) is a quaternion function, it can be represented as
where
f
_{0}
(
x
,
y
),
f
_{1}
(
x
,
y
),
f
_{2}
(
x
,
y
) and
f
_{3}
(
x
,
y
) are real value functions.
For leftside QHT, substituting (17) into (13), we have
Considering the general unit pure quaternion
µ
=
ξi + ηj +γk
(
ξ, η
and
γ
are real numbers), substituting
µ
into (18) and using the properties of the quaternion shown in Eq. (2), we have
where
Similarly, applying leftside IQHT to Eq. (19), the reconstructed
f’_{q}
(
x
,
y
) can be obtained.
where
It can be observed from formulas (18)~(22), the leftside QHT and IQHT of a quaternion matrix can be calculated effectively by using the traditional HT and IHT algorithms. Notice that the rightside QHT can be processed in a similar way.
IV. THE 2D LOGISTIC MAPBASED IMAGE SCRAMBLING METHOD
By using a scrambling technique, He
et al
. claimed that an optical security system can resist some potential attacks such as chosen plaintext attack
[27]
. Since the dynamic response of the chaotic system is highly sensitive to the initial values and parameters of chaotic variables, and the chaotic trajectory is unpredictable, a 2D logistic mapbased image scrambling technique is proposed to permute the position of image pixels.
Let
I
(
m
,
n
) be the tobepermuted image. The scrambling method is described as follows
1) Calculate the height
M
and the width
N
of
I
(
m
,
n
).
2) Initialize
X
(1) and
Y
(1) randomly and choose an arbitrary natural number
t
first, then iteratively generate the chaotic sequences
X
(
i
) and
Y
(
i
) whose lengths both are
MN
+
t
by using Eq. (12). Here,
i
=1,2, …,
MN
+
t
.
3) Generate two integers
s
_{1}
and
s
_{2}
which are between 1 and
t
randomly first. In other words, 1≤
s
_{1}
≤
t
and 1≤
s
_{2}
≤
t
. Then truncate
NM
elements of
X
(
i
) from the
s
_{1}
^{th}
element to obtain a chaotic sequence
L
1={
X
(
i
),
i
=
s
_{1}
,
s
_{1}
+1,…,
s
_{1}
+
MN
1}. Similarly, another chaotic sequence
L
2={
Y
(
i
),
i
=
s
_{2}
,
s
_{2}
+1,…,
s
_{2}
+
MN
1} can be obtained by truncating from
Y
(
i
).
4) Compute the chaotic sequences
L
1 and
L
2 using the following equation
where round(
Z
) is the operation that rounds the elements of
Z
to the nearest integers, and
d
is an integer. After calculation, the irregularity and distribution uniformity of the sequences
L
1 and
L
2 can be enhanced, and all the elements of
L
1 and
L
2 are larger than 0.5 and less than 0.5
[28]
. That is 0.5<
L
1<0.5 and 0.5<
L
2<0.5.
5) Sort the sequences
L
_{1}
and
L
_{2}
in a certain (ascending or descending) order to obtain two new sequence
L
’1 and
L
’2 and their corresponding permutation indices
IX
1 and
IY
2. In
IX
1 and
IY
2, there are
MN
elements respectively. The relations between
L
1 and L’1, and between
L
2 and
L
’2 are
L
’1=
L
1(
IX
1) and
L
’2=
L
2(
IY
2), respectively. For example, the
m
^{th}
element in
L
’1 corresponds to the
IX
1(
m
)
^{th}
element in
L
1.
6) Map
I
(
m
,
n
) into a onedimensional (1D) array
I
’ using the zigzag algorithm
[29]
. So the length of
I
’ is
MN
7) Use the permutation indices
IX1
and
IY2
to permute
I
’ and the scrambled vector
I
''' can be obtained as follows
8) Finally, the scrambled image
SI
can be achieved by applying inverse zigzag scan
[29]
process to
I
'''.
The inverse image scrambling process is similar to the image scrambling process. In decryption (inverse scrambling process), obtain the permutation indices
IX
1 and
IY
2 as described in steps (1)~(5) with the same parameters first. Then map the scrambled image
SI
into a 1D vector
SI
’ using the zigzag algorithm. Subsequently, permute
SI
’ back to their original position using the following equations
Finally, the decrypted image
DI
can be retrieved by applying the inverse zigzag algorithm to
SI
'''. Since the proposed logisticbased scrambling method is highly sensitive to the parameters, including the initial values
X
(1) and
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
, and
s
_{2}
, it offers a huge level of security.
V. THE PROPOSED HYBRID COLOR AND GRAYSCALE IMAGE ENCRYPTION AND DECRYPTION
Let
f_{RGB}
(
x
,
y
) be a color image with size
M×N
in the RGB color space. In principle, each image pixel of
f_{RGB}
(
x
,
y
) can be treated as a pure quaternion number with real part equal to zero
[17]
:
where
f_{R}
(
x
,
y
),
f_{G}
(
x
,
y
), and
f_{B}
(
x
,
y
) are the red, green and blue channels of
f_{RGB}
(
x
,
y
), respectively. By using this representation, each color triple can be treated as a whole unit and be applied in encryption and watermarking
[17
,
18
,
30]
. Hence, a color image which is represented by a pure quaternion can be processed holistically in a vector manner using QHT. Furthermore, it should be noticed that a quaternion contains four parts, which results that more than a simple color image can be processed in encryption. For example, the real part of the quaternion matrix can be used to represent a grayscale image while the three imaginary parts represent a color image.
Therefore, based on the QHT, the scrambling method and the GTbased DRPE, a hybrid color and grayscale images encryption scheme is proposed. In the proposed method, the three step PSI is used to record the encrypted data. The optoelectronic setup of the proposed encryption process is shown in
Fig. 2
. Supposing
f_{gs}
(
x
,
y
) denotes a grayscale image with size
M×N
, the proposed method is described as follows
Optoelectronic hybrid system implementing the proposed encryption method.
1) Normalize
f_{RGB}
(
x
,
y
) and
f_{gs}
(
x
,
y
) first. To treat the color image and grayscale image in a holistic manner,
f_{RGB}
(
x
,
y
) and
f_{gs}
(
x
,
y
) are represented by quaternion: 1) Normalize
f_{RGB}
(
x
,
y
) and
f_{gs}
(
x
,
y
) first. To treat the color image and grayscale image in a holistic manner,
f_{RGB}
(
x
,
y
) and
f_{gs}
(
x
,
y
) are represented by quaternion:
f_{q}
(
x
,
y
) =
f_{gs}
(
x
,
y
) +
if_{R}
(
x
,
y
) +
if_{G}
(
x
,
y
) +
kf_{B}
(
x
,
y
) .
2) Apply QHTL to fq(x,y):
HF
(
u
,
v
)=
QHT^{L}
[
f_{q}
(
x
,
y
)]=
A
_{0}
(
u
,
v
) +
iA
_{1}
(
u
,
v
) +
jA
_{2}
(
u
,
v
) +
jA
_{2}
(
u
,
v
) +
kA
_{3}
(
u
,
v
) .
Here,
QHT^{L}
() is the leftside QHT operation.
A
_{0}
,
A
_{1}
,
A
_{2}
and
A
_{3}
can be computed by using Eq. (20).
3) With
A
_{0}
,
A
_{1}
,
A
_{2}
and
A
_{3}
, form a new matrix using the following formula.
Since the length and width of
f_{RGB}
(
x
,
y
) and
f
_{gs}
(
x
,
y
) are both
M×N
, the size of
AI
is
2M×2N
.
4) Permute
AI
to obtain
AIS
by the proposed chaosbased scrambling method mentioned in Section 4 with the parameters
X
(1),
Y
(1),
λ
1,
λ
2,
s
_{1}
and
s
_{2}
. That is
AIS
=
scramble
(
AI, X
(1),
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
,
s
_{2}
). Here scramble() is the scrambling operation.
5) Let
EA
=
AIS
_{1}
+
jAIS
_{2}
, where
AIS
_{1}
=
AIS
[1:
M
, 1:2
N
] is the first half part of matrix
AIS
and
AIS
_{2}
=
AIS
[
M
+1:2
M
, 1:2
N
] is the second half part of
AIS
, respectively.
EA
can be regarded as the complex amplitude that will be encrypted by the DRPE technique in the gyrator domain. The size of
EA
is
M×2N
.
6) The complex amplitude
EA
is multiplied by
RPM
1, mathematically represented by function
U
(
x
_{0}
,
y
_{0}
) =
EA
(
x
_{0}
,
y
_{0}
) exp[
j2πφ
(
x
_{0}
,
y
_{0}
)].
7) As shown in
Fig. 2
, import
U
(
x
_{0}
,
y
_{0}
) into the first spatial light modulator (SLM) SLM1 which is located in the input plane under the control of the computer, and then optically transformed by first GT at order
α
_{1}
. The resulting complex distribution
U
_{1}
(
p,s
) can be represented by
8) Subsequently, the complex function
U_{1}
(
p,s
) is multiplied by
RPM
2 displayed on the SLM2 which is located in GT plane, and then optically transformed by second GT at order
α
_{2}
. So the object wave becomes
where
A_{c}
(
x
_{1}
,
y
_{1}
) and
ψ
(
x
_{1}
,
y
_{1}
) are the amplitude and the phase of
U_{c}
(
x
_{1}
,
y
_{1}
), respectively.
9) For convenience of storage and transmission, the object wave and reference wave overlap to produce realvalued interferograms. As shown in
Fig. 2
, the resulting encrypted data is recorded as three interference patterns
I_{1}
,
I_{2}
and
I_{3}
in threestep PSI case mentioned in Subsection 2.4.
In the encryption process mentioned above, the steps (1)~(6) will be performed digitally. In
Fig. 2
, M1, M2 are mirrors, and BS1 and BS2 are beam splitters. GT1 and GT2 are corresponding to the gyrator transform systems GT1 and GT2 shown in
Fig. 1
, respectively. As shown in
Fig. 2
, the complex function
U
(
x
_{0}
,
y
_{0}
) obtained in step (6) and the
RPM
2 are entered into the SLM1 and SLM2 respectively under the control of computer in order to be transformed, respectively. The phase shift
δ_{k}
(
k
=1, 2, 3) is introduced by a mirror attached with a Piezoelectric Transducer (PZT) controlled by a computer at the
k
exposure
[31]
. Finally, the resulting encrypted signal is captured and recorded as three interference patterns by a charge coupled device (CCD) which is placed at the output plane of the encryption system. The flowchart of the encryption process is shown in
Fig. 3
.
Flowchart of encryption process.
The parameters of the proposed encryption method, including
X
(1),
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
,
s
_{2}
,
RPM
1,
RPM
2 and the orders of GT
α
_{1}
and
α
_{2}
, form a very large key space enhancing the security level of the encryption system. Only when all the right keys are simultaneously used for decryption can the encrypted data be reconstructed correctly.
The decryption process, which is similar to the encryption process, but in the reversed order is depicted as follows:
1) Using the three interferograms
I
_{1}
,
I
_{2}
and
I
_{3}
, an encrypted digital hologram
UE
(
x
_{1}
,
y
_{1}
) can be achieved using Eq. (11).
2) Apply a GT to
UE
(
x
_{1}
,
y
_{1}
) with order α2 and then multiply the obtained complex distribution by the conjugate of the
RPM
2.
3) Make another GT with order 
α
_{1}
and then multiply the achieved complex distribution by the conjugate of the
RPM
1. The resulting distribution function is as follows
where
EA_{re}
(
x
_{0}
,
y
_{0}
) and
EA_{im}
(
x
_{0}
,
y
_{0}
) are the real part and imaginary part of
EA’
(
x
_{0}
,
y
_{0}
), respectively.
4) With
EA_{re}
and
EA_{im}
, form a new matrix using the following formula.
5) With the parameters
X
(1),
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
and
s
_{2}
,
AIS’
is permuted by the proposed inverse scrambling process mentioned in Section 4 to obtain
AI’
.
6) First, let
A
’
_{0}
=AI’[1:
M
, 1:
N
],
A
’
_{1}
=
AI’
[1:
M
,
N
+1:2
N
],
A
’
_{2}
=
AI’
[
M
+1:2
M
, 1:
N
] and
A’
_{3}
=
AI’
[
M
+1:2
M
,
N
+1:2
N
]. Then let
HF’
(
u
,
v
)=
A’
_{0}
+
iA’
_{1}
(
u
,
v
)+
jA’
_{2}
(
u
,
v
)+
kA’
_{3}
(
u,v)
.
7) Apply IQHT to
HF’
(
u
,
v
) to achieve
f’q
(
x
,
y
).
8) With the normalized
f’_{R}
(
x
,
y
),
f’_{G}
(
x
,
y
) and
f’_{B}
(
x
,
y
), the decrypted color image
f’_{RGB}
(
x
,
y
) can be obtained and the normalized
f’_{gs}
(
x
,
y
) is the decrypted grayscale image.
Figure 4
depicts the decryption process. In
Fig. 4
, conj(
RPM
2) and conj(
RPM
1) are the conjugate of
RPM
2 and the conjugate of
RPM
1, respectively.
Flowchart of decryption process.
VI. NUMERICAL SIMULATION RESULTS
To verify the feasibility of the proposed encryption technique, numerical simulations are performed on the color image “Peppers” and the grayscale image “Lena”, which sizes are both 512×512, shown in
Fig. 5(a)
and
Fig. 4(b)
. We carried out tests on a notebook computer with Intel(R) Core(TM) i74700HQ CPU @ 2.40GHz and 8G DDRL3 and with the MATLAB R2013a. In the experiments, the system parameters are
µ
=(
i+j+k
)/3
^{1/2}
,
α
_{1}
=0.92,
α
_{2}
=0.75,
ε
=4,
ρ
=0.1,
X
(1) = 0.24,
Y
(1) = 0.67,
t
=50000,
λ
_{1}
=0.716,
λ
_{2}
=0.881,
s
_{1}
= 12345,
s
_{2}
= 36578 and
d
=4. To evaluate the performance of image reconstruction quantitatively, the normalized mean square error (NMSE)
[18]
was used to calculate the similarity between the decrypted image and the original image, which is expressed as
Results of the proposed image encryption. (a) The original color image “Peppers”; (b) the original grayscale image “Lena”; (c), (d), (e) interferograms I_{1}, I_{2} and I_{3}, respectively; (f) correctly decrypted color image; (g) correctly decrypted grayscale image.
where,
f_{o}
(
x
,
y
) and
f_{d}
(
x
,
y
) are the normalized original and the decrypted images respectively, and
M×N
is the size of the images.
 6.1. Performance of the Encryption System
Using the proposed encryption scheme, the “Peppers” and “Lena” are encrypted and three encrypted interferograms are obtained, which are shown in
Figs. 5(c)
~
5(e)
.
Figs. 5(f)
~
5(g)
display the retrieved color and grayscale images with the correct keys respectively, which are perfect, without any noise or distortion. As shown in
Table 1
, the NMSE between
Fig. 5(a)
and
Fig. 5(f)
and that between
Fig. 5(b)
and
Fig. 5(g)
are 1.7017×10
^{30}
and 1.173×10
^{30}
, respectively.
Now we investigate the sensitivity of retrieved image to small change of the parameters of 2D logistic mapbased scrambling technique on the decrypted results.
Figure 6
showed the derivation of NMSE versus the parameters
X
(1),
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
and
s
_{2}
, respectively.
Figs. 7(a)
~
7(l)
show the decrypted color and grayscale images with wrong keys
X
(1)=0.24+1.0×10
^{15}
,
Y
(1)=0.671.0×10
^{15}
,
λ
_{1}
=0.716 +1.0×10
^{16}
,
λ
_{2}
=0.8811.0×10
^{16}
,
s
_{1}
= 12344 and
s
_{2}
= 36579, respectively.
The NMSEs of the proposed method: (a), (b) the NMSEs for deviation of X(1) and Y(1); (c), (d) the NMSEs for deviation of λ_{1} and λ_{2}; (e), (f) the NMSEs for deviation of s_{1} and s_{2}.
The decrypted color “Peppers” and grayscale “Lena” with: (a), (b) incorrect X(1)=0.24+1.0×10^{15}; (c), (d) incorrect Y(1)=0.671.0×10^{15}; (e), (f) incorrect λ _{1}=0.716+1.0×10^{16}; (g), (h) incorrect λ_{2}=0.8811.0×10^{16}; (i), (j) incorrect s_{1}= 12344; (k), (l) incorrect s_{2}= 36579.
The corresponding NMSEs are listed in the upper part of
Table 1.
Please note that in the above experiments, the other keys remain correct while a key is changed in decryption. As illustrated in
Figs. 6(a)
~
6(d)
and
Figs. 7(a)
~
7(h)
, we cannot obtain any information from the decryption images visually when the absolute values of deviations of
X
(1) and
Y
(1) are up to 10
^{15}
and those of
λ
_{1}
and
λ
_{2}
are up to 10
^{16}
. In addition, we know from
Figs. 6(e)
~
6(f)
and
Figs. 7(i)
~
7(l)
that if the parameters
s
_{1}
and
s
_{2}
are less 1 or more 1 than the correct value, the decrypted images cannot afford any valid information. From the NMSE values shown in
Table 1
, it can also be concluded that the decrypted images cannot be recognized even if the keys of the scrambling system are changed slightly. So, the parameters
X
(1),
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
and
s
_{2}
are highly sensitive to the proposed method.
Comparison between the NSMEs of the decrypted “Peppers” and “Lena” by use of correct and incorrect keys
Comparison between the NSMEs of the decrypted “Peppers” and “Lena” by use of correct and incorrect keys
In the experiments, one of the two
RPM
s was removed or just one of them was shifted transversely by one pixel. The NMSE values of the decrypted images are listed in the middle part of
Table 1
.
Figs. 8(a)
~
8(b)
display the decrypted color and grayscale images, which are totally unrecognizable, with other correct keys but without
RPM
1. The recovered color and grayscale images, which were obtained from
RPM
1 shifted transversely by one pixel, are shown in
Figs. 8(c)
~
8(d)
. These two images cannot be recognized too. The lack of the
RPM
2 or small changes to
RPM
2 in the decryption step also leads to similar results which are exhibited in
Figs. 8(e)
~
8(h)
. Experimental results convincingly demonstrate that both
RPM
s are necessary in decryption and that the illegal users cannot access the security system without the
RPM
s.
The decrypted color “Peppers” and grayscale “Lena” with: (a), (b) lack of RPM1; (c), (d) RPM1 is shifted transversely by one pixel; (e), (f) lack of RPM2; (g), (h) RPM2 is shifted transversely by one pixel.
To examine the sensitivity of small change of the orders of the GT, the decryption processes are performed by fixing one order and varying the other. The relationship curves between the NMSE and the deviation of the GT order are shown in
Fig. 9
, in which the deviation ranges from 0.1 to 0.1 and the step is 0.001. As can be seen from
Fig. 9
, the NMSE value approximates to zero when
α
_{1}
or
α
_{2}
is correct while the value sharply increases when
α
_{1}
or
α
_{2}
slightly departs from the correct value, which indicates that any tiny fluctuation will result in false decryption. The NMSE values shown in the lower part of
Table 1
also indicate that if the orders of inverse GT are changed by 0.004 from their correct values, the decrypted images are noiselike images as displayed in
Fig. 10
and cannot be recognized.
NMSEs versus (a) the deviation of order α_{1} and (b) the deviation of order α_{2}.
The decrypted color “Peppers” and grayscale “Lena” with: (a), (b) incorrect α_{1}= 0.92+0.004; (c), (d) incorrect α_{2}=0.750.004.
Next we estimate the key space of the proposed encryption scheme. According to the description of the proposed method, we know that the key space of the cryptosystem consists of the
RPM
s, the orders of GT and the parameters of 2D logistic map. In addition, the chaotic permutation module is independent from the GTbased DRPE system. In the logisticbased permutation process, the key spaces of the parameters
X
(1),
Y
(1),
λ
_{1}
,
λ
_{2}
,
s
_{1}
and
s
_{2}
should be analyzed, which are denoted by
Z
_{1}
,
Z
_{2}
,
Z
_{3}
,
Z
_{4}
,
Z
_{5}
and
Z
_{6}
, respectively. In the GTbased DRPE system, the key spaces of
RPM
1,
RPM
2 and the orders
α
_{1}
and
α
,
_{2}
, should be analyzed, which are denoted by
Z
_{7}
,
Z
_{8}
,
Z
_{9}
and
Z
_{10}
respectively. From
Figs. 6(a)
~
6(d)
and
Figs. 7(a)
~
7(h)
, the parameters
X
(1),
Y
(1),
λ
_{1}
and
λ
_{2}
maintain 15, 15, 16 and 16 digits after decimal point respectively, so
Z
_{1}
×
Z
_{2}
×
Z
_{3}
×
Z
_{4}
=10
^{62}
. Since 1≤
s
_{1}
≤
t
and 1≤
s
_{2}
≤
t
,
Z
_{5}
×
Z
_{6}
=
t
^{2}
. Therefore,
Z
_{1}
×
Z
_{2}
×
Z
_{3}
×
Z
_{4}
×
Z
_{5}
×
Z
_{6}
=
t
^{2}
×10
^{62}
. In our experiments, Z
_{1}
×Z
_{2}
×
Z
_{3}
×
Z
_{4}
×
Z
_{5}
×
Z
_{6}
=2.5×10
^{71}
since
t
=50000. Suppose that 256 levels are chosen for a pixel in initial random functions
φ
(
x
_{0}
,
y
_{0}
) and
ϕ
(
p,s
) of
RPM
1 and
RPM
2, the two
RPM
s with size of
M
×2
N
need to be decoded with 2×256
^{M×2N}
attempts. In this case,
M
=
N
=512, so
Z
_{7}
×Z
_{8}
=2×256
^{512×1024}
. Considering the sensitivity to one GT order of encryption system is 0.004, the fractional orders
α
_{1}
and
α
_{2}
can be varied with steps of 0.001 in fractional orders range [0, 1]
[7]
, which results in
Z
_{9}
×
Z
_{10}
=(1/0.001)
^{2}
=10
^{6}
. Therefore, the entire key space of the cryptosystem is
Z
_{1}
×
Z
_{2}
×
Z
_{3}
×
Z
_{4}
×
Z
_{5}
×
Z
_{6}
×
Z
_{7}
×
Z
_{8}
×
Z
_{9}
×
Z
_{10}
=5×256
^{512×1024}
×10
^{77}
≈ 256
^{512×1024}
×2
^{258}
. As stated in
[32]
, to acquire a high level of security, the size of key space should at least be larger than 2
^{100}
. It is apparent that the key space of the proposed encryption system is far larger than 2
^{100}
and enormous enough to resist brute force attack.
 6.2. Robustness of the Method Against Attacks
Information loss or noise contamination may occur during data transmission. Now, the robustness of this scheme against occlusion attack which is regarded as data loss is tested.
Figure 11(a)
shows one of the interferograms occluded by 50%. The recovered images obtained with all correct keys are illustrated in
Figs. 11(b)
~
11(c)
. It’s well known that Gaussian noise
[33]
and salt & pepper noise
[33]
are frequently appearing noises during the information transmission. The robustness test is further verified against noise attacks on the encrypted results.
Figure 12(a)
is one of the interferograms distorted by Gaussian noise with mean value 0 and standard deviation 20. The corresponding retrieved images are displayed in
Figs. 12(b)
~
12(c)
.
Fig. 12(d)
exhibits one of the interferograms damaged by salt & pepper noise with density 0.02 added.
Figs. 12(e)
~
12(f)
depicts the corresponding retrieved images. Although all the results shown in
Figs. 11(b)
~
11(c)
,
Figs. 12(b)
~
12(c)
and
Figs. 12(e)
~
12(f)
are interfered with seriously by noise and the corresponding NMSEs are big, the secret images among the noise fluctuation can still be distinguished.
Robustness against occlusion. (a) one of the interferograms with 50% occlusion; (b) the recovered “Peppers” (NMSE=0.498); (c) the recovered “Lena” (NMSE=0.4165).
Robustness against noise. (a) one of the interferograms after Gaussian noise with mean value 0 and standard deviation 20; (b) the recovered “Peppers” from interferograms that have undergone Gaussian noise addition (NMSE=0.4471); (c) the recovered “Lena” from interferograms that have undergone Gaussian noise (NMSE=0.4133); (d) one of the interferograms after salt & pepper noise with density 0.02 added; (e) the recovered “Peppers” from interferograms that have undergone salt & pepper noise (NMSE =0.4247); (f) the recovered “Lena” from interferograms that have undergone salt & pepper noise (NMSE=0.4157).
Four attacks, including cipher only attack, known plaintext attack, chosen plaintext attack and chosen ciphertext attack
[32]
, are often used to attack the DRPEbased optical security system to recover the ciphered images. Among these attacks, the chosenplaintext attack is the most powerful attack
[32]
. He
et al
. claimed that permutation techniques can prevent the attacker from accessing the ciphertext obtained from the DRPE system in such attacks
[27]
. In addition, Zhang and Xiao stated that the security system should resist other attacks if it can resist chosen plaintext attack
[32]
. For the proposed method, we assume that the attacker has acquired the GT orders
α
_{1}
and
α
_{2}
of the GTbased DRPE system but does not know the parameters of the logisticbased scrambling algorithm. Using two fake color and gray images as the plaintexts, the two random phase masks of the proposed method can be produced by utilizing the chosen plaintext attack described in
[34]
. By use of these achieved keys, the attacker can crack the original ciphertext. Two keys
RPM
1’ and
RPM
2’ are generated by choosing
Figs. 13(a)

13(b)
as the fake plaintext images and then employed to decrypt the ciphertext of “Peppers” and “Lena”. As shown in
Figs. 13(c)
~
13(d)
, the retrieved images cannot reveal any information visually. Therefore, the robustness against chosenciphertext attack is greatly enhanced because of the suggested logisticbased scrambling technique, which strengthens the nonlinearity in QHT domain. According to
[32]
, the proposed encryption approach can also resist the other attacks mentioned above.
Robustness against chosen plaintext attack. (a), (b) Two fake plaintext images “Fruits” and “Elaine”; (c), (d) decrypted “Peppers” and “Lena”, the corresponding NMSEs are 0.5161 and 0.4924, respectively.
 6.3. Complexity Analysis
In the numerical simulation, the encrypted results were obtained by computer. So the computational complexity of the proposed encryption technique is calculated. The calculation depends on three main factors. First, to deal with a color image and a grayscale image holistically in a vector manner, the QHT is computed one time. According to the developed calculation method in subsection 3.2, the QHT can be implemented via calculating HT four times. For an image with size
M
×
N
, the cost to compute the HT is
O
(
MIN
log
_{2}
MN
)
[35]
. Hence, the complexity of QHT calculation is
O
(
4MIN
log
_{2}
MN
). Second, to scramble a matrix using the proposed chaosbased technique, the 2D logistic map operation, the round operation, the quicksort algorithm
[36]
, the zigzag scan operation, the inverse zigzag scan operation and the sequence permutation operation are performed one time, two times, two times, one time, one time and one time, respectively. As shown in the steps (3) and (4) of the proposed encryption method, the size of the tobescrambled matrix
AI
is 2
M
×2
N
though the sizes of the color image and grayscale image are both
M
×
N
. Therefore, the computational cost of the scrambling procedure is
O
(2(4
MN
+
t
) + 2(4
MN
+
t
) + 2(4
MN
+
t
)
_{log2}
(4
MN
+
t
)) + (4
MN
+
t
) + (4
MN
+
t
) + (4
MN
+
t
). Here,
t
is the natural number shown in step (2) of the proposed scrambling method. Since 4
MN
+
t
≈4
MN
and log
_{2}
(4
MN
+
t
)≈2+log
_{2}
MN
as
MN
tends to infinity, the complexity of the scrambling procedure is
O
(8
MN
log
_{2}
MN
+ 44
MN
). The last factor is the complexity of the virtual optical encryption process. In the proposed encryption method, an encrypted interferogram can be obtained after the complex amplitude being modulated two times by use of random phase masks
RPM
1 and
RPM
2, transformed two times by use of GT and interfered one time by use of the reference wave. For a matrix with size
M×N
, the complexities of GT operation, modulation operation and interference operation are
O
(2
MN
log
_{2}
MN
)
[22]
,
O
(
MN
)
[37]
and
O
(4
MN
)
[37]
, respectively. Because the size of complex amplitude is
M×2N
, the computational cost of the above process is
O
(3(8
MN
log
_{2}
2
MN
+ 4
MN
+ 8
MN
)) =
O
(24
MN
log
_{2}
MN
+ 60
MN
). Thus, the total computational complexity of the proposed encryption scheme is
O
(36
MN
log
_{2}
MN
+104
MN
). Since
O
(
MN
log
_{2}
MN
) >
O
(
MN
), compared with the quaternion gyrator transformbased image encryption method whose computational complexity is
O
(4000
MN
log
_{2}
MN
)
[18]
, the proposed scheme is more efficient. In the experiments, the execution time for the encrypted results is 6.7728 seconds.
VII. CONCLUSION
In the paper, a novel definition of quaternion Hartley transform and its implementation are presented first. Then a new method for hybrid color and grayscale images encryption is proposed, in which the QHT combined with the developed chaosbased pixel scrambling technique, the GTbased DRPE and the threestep PSI is used to encrypt the color and grayscale images. By use of QHT, a color image and a grayscale image are processed holistically in a vector manner without separating the color image into three channels and manipulating the grayscale image independently, so that the complexity of the security system can be reduced effectively without any reduction of its security. The components of the QHTtransformed original images are combined as the complex amplitude which will be strengthened the nonlinearity by the designed chaosbased scrambling technique and encrypted by the GTbased DRPE technique. In addition, the encrypted interferograms which are recorded by threestep PSI are realvalued, so they are convenient for storage and transmission. One needs to specify all the right keys to recover the original images correctly. Simulation results demonstrate that the proposed scheme has a quite sensitivity to the decryption keys, a enormous key space to resist brute force attack and a good robustness against Gaussian noise, salt & pepper noise and some attacks such as chosen plaintext attack.
Acknowledgements
This work is partly supported by the Natural Science Foundation of Guangdong Province (No. 2014A030310038), the Educational Commission of Guangdong Province (No. 2013KJCX0127), the 2015 Educational Commission of Guangdong Province, the 2013 quality project of Educational Commission of Guangdong Province, the Science and Technology Planning Project of Chaozhou (No. 2013G01) and the Research Fund of Hanshan Normal University (No. LT201201).
Refregier P.
,
Javidi B.
(1995)
“Optical image encryption based on input plane and Fourier plane random encoding,”
Opt. Lett.
20
767 
769
DOI : 10.1364/OL.20.000767
Tajahuerce E.
,
Matoba O.
,
Verrall S. C.
,
Javidi B.
(2000)
“Optoelectronic information encryption with phaseshifting interferometry,”
Appl. Opt.
39
2313 
2320
DOI : 10.1364/AO.39.002313
Lee I. H.
,
Cho M. J.
(2014)
“Double random phase encryption using orthogonal encoding for multipleimage transmission,”
J. Opt. Soc. Korea
18
201 
206
DOI : 10.3807/JOSK.2014.18.3.201
Zheng J. C.
,
Li X. Q.
(2015)
“Image authentication using only partial phase information from a doublerandomphaseencrypted image in the Fresnel domain,”
J. Opt. Soc. Korea
19
241 
247
DOI : 10.3807/JOSK.2015.19.3.241
Yuan S.
,
Xin Y. H.
,
Liu M. T.
,
Yao S. X.
,
Sun X. J.
(2012)
“An improved method to enhance the security of double randomphase encoding in the Fresnel domain,”
Optics & Laser Technology
44
51 
56
DOI : 10.1016/j.optlastec.2011.05.015
Nischal N. K.
,
Joseph J.
,
Singh K.
(2004)
“Securing information using fractional Fourier transform in digital holography,”
Opt. Commun.
235
253 
259
DOI : 10.1016/j.optcom.2004.02.052
Wang Q.
,
Guo Q.
,
Lei L.
(2013)
“Double image encryption based on phase–amplitude mixed encoding and multistage phase encoding in gyrator transform domains,”
Optics & Laser Technology
48
267 
279
DOI : 10.1016/j.optlastec.2012.10.037
Muniraj I.
,
Guo C. L.
,
Lee B. G.
,
Sheridan J. T.
(2015)
“Interferometry based multispectral photonlimited 2D and 3D integral image encryption employing the Hartley transform,”
Opt. Express
23
15907 
15920
DOI : 10.1364/OE.23.015907
Abuturab M. R.
(2012)
“Securing color information using Arnold transform in gyrator transform domain,”
Optics and Lasers in Engineering
50
772 
779
DOI : 10.1016/j.optlaseng.2011.12.006
Joshi M.
,
Shakher C.
,
Singh K.
(2008)
“Color image encryption and decryption for twin images in fractional Fourier domain,”
Opt. Commun.
281
5713 
5720
DOI : 10.1016/j.optcom.2008.08.024
Chen W.
,
Quan C.
,
Tay C.
(2009)
“Optical color image encryption based on arnold transform and interference method,”
Opt. Commun.
282
3680 
3685
DOI : 10.1016/j.optcom.2009.06.014
Rajput S. K.
,
Nishchal N. K.
(2012)
“Asymmetric color cryptosystem using polarization selective diffractive optical element and structured phase mask,”
Appl. Opt.
51
5377 
5386
DOI : 10.1364/AO.51.005377
Mehra I.
,
Nishchal N. K.
(2014)
“Image fusion using wavelet transform and its application to asymmetric cryptosystem and hiding,”
Opt. Express
22
5474 
5482
DOI : 10.1364/OE.22.005474
Wang X. L.
,
Zhai H. C.
,
Li Z. L.
,
Ge Q.
(2011)
“Double randomphase encryption based on discrete quaternion Fouriertransforms,”
Optik
122
1856 
1859
DOI : 10.1016/j.ijleo.2010.11.016
Shao Z. H.
,
Shu H. Z.
,
Wu J. S.
,
Dong Z. F.
,
Coatrieux G.
,
Coatrieux J. L.
(2014)
“Double color image encryption using iterative phase retrieval algorithm in quaternion gyrator domain,”
Opt. Express
22
4932 
4943
DOI : 10.1364/OE.22.004932
Abuturab M. R.
(2015)
“Generalized Arnold mapbased optical multiple colorimage encoding in gyrator transform domain,”
Opt. Commun.
343
157 
171
DOI : 10.1016/j.optcom.2014.12.085
Hamilton W. R.
2009
Elements of Quaternions
Cambridge University Press
Cambridge, UK
Chapter 2
Gupta N.
,
Saxena D. K.
(2013)
“Implementation and contrasting comparision of transform techniques of OFDM,”
International Journal of Computer Applications
72
47 
51
Liu Z. J.
,
Chen D. Z.
,
Ma J. P.
,
Wei S. Y.
,
Zhang Y. L.
,
Dai J. M.
,
Liu S. T.
(2011)
“Fast algorithm of discrete gyrator transform based on convolution operation,”
Optik
122
864 
867
DOI : 10.1016/j.ijleo.2010.06.010
Zhang S.
,
Weide D. V. D.
,
Oliver J.
(2010)
“Superfast phaseshifting method for 3D shape measurement,”
Opt. Express
18
9684 
9689
DOI : 10.1364/OE.18.009684
Cai L. Z.
,
Liu Q.
,
Yang X. L.
(2003)
“Phaseshift extraction and wavefront reconstruction in phaseshifting interferometry with arbitrary phase steps,”
Opt. Lett.
28
1808 
1810
DOI : 10.1364/OL.28.001808
Cai L. Z.
,
He M. Z.
,
Liu Q.
,
Yang X. L.
(2004)
“Digital image encryption and watermarking by phaseshifting interferometry,”
Appl. Opt.
43
3078 
3084
DOI : 10.1364/AO.43.003078
Chen D. M.
,
Liu Y. M.
,
Chen X. D.
,
Chang Y. P.
,
Wang J.
(2011)
“A novel chaotic map and DES based file encryption algorithm,”
International Journal of Advancements in Computing Technology
3
198 
205
He M. Z.
,
Tan Q. F.
,
Cao L. C.
,
He Q. S.
,
Jin G. F.
(2009)
“Security enhanced optical encryption system by random phase key and permutation key,”
Opt. Express
17
22462 
22473
DOI : 10.1364/OE.17.022462
Wang Y.
,
Zheng D. L.
,
Ju L.
(2004)
“Digital image encryption algorithm based on threedimension Lorenz chaos,”
Journal of University of Science and Technology Beijing
(in Chinese)
26
678 
682
Smila M.
,
Sankar S.
(2014)
“Novel algorithms for finding an optimal scanning path for JPEG image compression,”
IJETCSE
8
230 
236
Tsui T. K.
,
Zhang X. P.
,
Androutsos D.
(2008)
“Color image watermarking using multidimensional Fourier transforms,”
IEEE Trans. Inform. Forensics Secur.
3
16 
28
DOI : 10.1109/TIFS.2007.916275
Wang X. G.
,
Zhao D. M.
(2006)
“Image encryption based on anamorphic fractional Fourier transform and threestep phaseshifting interferometry,”
Opt. Commun.
268
240 
244
DOI : 10.1016/j.optcom.2006.07.024
Zhang Y.
,
Xiao D.
(2013)
“Double optical image encryption using discrete Chirikov standard map and chaosbased fractional random transform,”
Opt. Lasers Eng.
51
472 
480
DOI : 10.1016/j.optlaseng.2012.11.001
Gonzalez R. C.
,
Woods R. E.
2007
Digital Image Processing
3rd ed.
Prentice Hall
New Jersey, USA
Chapter 5
Peng X.
,
Wei H. Z.
,
Zhang P.
(2006)
“Chosenplaintext attack on lensless doublerandom phase encoding in the Fresnel domain,”
Opt. Lett.
31
3261 
3263
DOI : 10.1364/OL.31.003261
Li C. B.
2013
Data Structure
Tsinghua University Press
Beijing, PRC
Chapter 9
(in Chinese)
Su X.
(2014)
“Virtual optical encryption system based on coupling sawtooth spatiotemporal chaos,” Master Thesis
Nanjing University of Posts and Telecommunications
(in Chinese)