Image Encryption Scheme Based on Multiscale Block Compressed Sensing and Markov Model

Many image encryption schemes based on compressed sensing have the problem of poor quality of decrypted images. To deal with this problem, this paper develops an image encryption scheme by multiscale block compressed sensing. The image is decomposed by a three-level wavelet transform, and the sampling rates of coefficient matrices at all levels are calculated according to multiscale block compressed sensing theory and the given compression ratio. The first round of permutation is performed on the internal elements of the coefficient matrices at all levels. Then the coefficient matrix is compressed and combined. The second round of permutation is performed on the combined matrix based on the state transition matrix. Independent diffusion and forward-backward diffusion between pixels are used to obtain the final cipher image. Different sampling rates are set by considering the difference of information between an image’s low- and high-frequency parts. Therefore, the reconstruction quality of the decrypted image is better than that of other schemes, which set one sampling rate on an entire image. The proposed scheme takes full advantage of the randomness of the Markov model and shows an excellent encryption effect to resist various attacks.


Introduction
With the rapid development of internet technology, digital images are widely used in all walks of life. As images can intuitively display information, they are widely used in banking, the military, medicine, finance, and other fields. Once images containing important information spread on the internet, they are likely to be attacked and cracked by hackers. To make meaningful images meaningless through encryption can effectively protect the security of private information in transmission and storage.
The combination of compressed sensing (CS) theory and image encryption algorithms, by which an image is simultaneously compressed and encrypted, has seen much research. The most special advantage of compressed sensing is that it can reduce data redundancy. Before data transmission, if the cipher image can be further compressed, it will not only improve the transmission efficiency but also be more suitable for transmission channels with limited bandwidth. Moreover, because the sampling process of compressed sensing is a linear random projection process, in the presence of noise, the target sparse signal can be reconstructed with high probability through a greedy algorithm or an optimized block compressed sensing to the design of an image encryption scheme. The low-frequency coefficient matrix of an image is fully sampled, and different sampling rates are set for the remaining coefficient matrices according to the amount of information they carry. There is a significant improvement in the reconstructed image after decryption. Different from the traditional encryption algorithms like DNA coding [15,16], cellular automata [9,10], and substitution box [12,13] that follow generally known rules to carry out the subsequent work, this paper introduces the Markov model in machine learning, and the scrambling process is carried out according to the information of the plain image and chaotic sequences. Experimental results show that the proposed scheme achieves a good encryption effect.
Our contributions are as follows.
1. An encryption architecture of permutation, compression, secondary scrambling, and diffusion is designed, which shows good compression performance and guarantees security; 2.
A transition probability matrix in a Markov model is introduced to scramble the image and define the state space according to the characteristics of image pixel values in the encryption process. The state transition probability matrix is constructed based on the distribution of pixel values. The process achieves good randomness, so it is difficult to predict; 3.
Information about plain images and chaotic sequences is used in the encryption process, giving the scheme high plain sensitivity to resist known-plaintext attacks (KPAs) and chosen-plaintext attacks (CPAs); 4.
Multiscale block compressed sensing theory is introduced, sampling rates of images are set by a more reasonable approach, and the reconstruction quality of decrypted images is greatly improved.
The rest of this paper is organized as follows. Preliminaries are discussed in Section 2. The steps of the proposed scheme are described in Section 3. Simulation results are shown in Section 4, and a sensitivity analysis of the scheme is discussed in Section 5. Conclusions are drawn in Section 6.

Multiscale Block Compressed Sensing
Fowler et al. [45] proposed a multiscale block compressed sensing algorithm based on a wavelet domain. The MS-BCS-SPL algorithm divides the image into blocks based on the BCS algorithm and samples each image sub-block with a different matrix. The original image is decomposed by a multilayer wavelet transform, and the wavelet coefficients of each layer are divided into blocks whose size varies with the number of layers. A measurement matrix, determined by the sampling rate of each layer, is used for measurement. Since the different levels of wavelet decomposition have different importance to the final image reconstruction quality, each layer corresponds to a different sampling rate. The smooth projection Landweber method is used to reconstruct each image block; thus, the complete reconstructed image can be obtained.
As the main information of an image is concentrated in the low-frequency coefficient after wavelet decomposition, its details are concentrated in high-frequency coefficients. To improve the reconstruction quality of an image requires one to keep the low-frequency part of the image and abandon the high-frequency part as much as possible. The measurement matrix A of the original image is decomposed into multiscale transformation matrix Ω and multiscale measurement matrix Φ , and A is represented by A = Φ Ω, so the compression process can be expressed as The multiscale transformation matrix Ω is decomposed by an L-level wavelet transform to form L measurement operators, which constitute the multiscale block measurement matrix Φ . The wavelet decomposition process of the original image x can be expressed as The s-th sub-band of x is cut into sub-blocks of size B l × B l , each sampled by a sampling matrix of corresponding size. For example, the process of compressing the j-th sub block can be expressed as After the original image is decomposed, the influence of the decomposition block of each layer on the image reconstruction is different. To improve the reconstruction quality of the image, it is necessary to set the corresponding sampling rate for each wavelet decomposition layer. Let the baseband sampling rate of wavelet decomposition be 1, so S 0 = 1, and the wavelet decomposition sub-rate of each layer can be expressed as where W l is the weight of each layer after wavelet decomposition in the whole image, i.e., For the whole image, the sampling rate can be expressed as If the sampling rate S of the whole image and the weight W l of each wavelet decomposition layer are known, then the total sampling rate S can be calculated by Equation (6). The sampling rate S l of the l-th level wavelet decomposition coefficient can be obtained by substituting S in Equation (4). When solving for the sampling rate of each layer after wavelet decomposition, there will be multiple solutions greater than 1. In practice, the sampling rate of each layer should not be greater than 1, so such solutions should be set to 1. In this case, the sampling rate is expressed as We recalculate S l from Equation (7). The above process is repeated for the wavelet decomposition of each layer l = 2, 3, . . . , L, and we check for S l > 1 to ensure that S l ≤ 1 in each layer.

Chaotic Systems
The tent-logistic-tent system (TLTS) and tent-sine-tent system (TSTS) are, respectively, obtained as [32] X n+1 = f TLT (X n , r) =    r 2 2 X n 1 − r 2 X n + r 2 X n r 14 mod1, X n < 0.5 X n+1 = f TST (X n , r) = r 4 sin(π r 2 X n ) + r 2 X n r 14 mod1, X n < 0.5 r 4 sin π r 2 (1 − X n + r 2 (1 − X n ) r 14 mod1, X n ≥ 0.5 (9) where r ∈ (1.05.4] is the control parameter of the TLT and TST chaotic systems. When r is greater than 1.05, the LE value of the chaotic system is positive, i.e., the system is in a chaotic state. The hybrid chaotic system [12] is defined as where r ∈ [1. 4,4] is the control parameter of the hybrid chaotic map. When r is greater than 1.1, the LE value of the chaotic system is positive, i.e., the system is in a chaotic state. We use the three chaotic systems to generate three measurement matrices, respectively. A new chaotic system-improved sine-exponent-logistic(ISEL) is obtained as [46] X n+1 = (sin(πX n )) a ln (4bX n (1−X n )+c) where a ∈ [0, 1], b ∈ [0, 5], c ∈ [1.5, 2.8] are the control parameters of the ISEL system. We use this system to generate chaotic sequences for diffusion.

Markov Model
The Markov model can be used to simulate random processes [47]. In our scheme, a Markov chain is used to construct the state transition matrix for scrambling the image matrix.

State Space
In the Markov chain, every variable has several possible values, and the set of all these is called the state space. To generate this, numbers are divided into four categories in our scheme: If the integer part of a number is positive and odd, then it is a positive-odd (po-od) number; If the integer part of a number is negative and odd, then it is a negative-odd (ne-od) number; If the integer part of a number is positive and even, then it is a positive-even (po-ev) number; If the integer part of a number is negative and even, then it is a negative-even (ne-ev) number.

Markov State Transition Matrix
Obviously, since there is more than one state at each time, there are several cases of transferring from the previous to the current state. All conditional probabilities form a transition probability matrix, as shown in Table 1. Here, a i , i = 1, 2, . . . is the state at the previous time, a j , j = 1, 2, . . . is the state at the current time, and p ij , i = 1, 2, . . . ; j = 1, 2, . . . is the conditional probability from a i to a j . According to the above definition, the header of the state transition matrix is constructed.
Next, we initialize a matrix f of size 4 × 4, which is used to record the frequency of state transition. For each column vector of the matrix to be measured, we first determine the type of the first element, i.e., the row coordinate of the state transition matrix. Then the type of the next element is determined, i.e., the column coordinate of the state transition matrix. The position in the matrix f corresponding to the coordinate point is incremented by 1 to obtain the updated matrix f. For example, the matrix to record the frequency of state transition is shown in Table 2. We calculate the sum of four frequencies in each row of the matrix f, respectively, which means the total frequency of transitions from one type of number to other types of number. Each element is divided by the sum of the corresponding row to get a probability, which means the transition probability from one type of number to another type of number. After all probabilities are calculated, the row-based state transition probability matrix (RSTPM) is obtained, as shown in Table 3. The generation process of the column-based state transition probability matrix (CSTPM) is similar to that of RSTPM, where we get a row vector instead of a column vector of the matrix.  For RSTPM, we find all positions of the probabilities greater than 0.25 of the matrix shown in Table 3. The row coordinates indicate whether odd-or even-column vectors are to be selected, and whether the direction of movement is up or down. If the row coordinate is odd (even), then odd (even)-column vectors are selected; if the row coordinate is positive (negative), then all selected elements move up (down). Specifically, if the row coordinate is a po-od (ne-od) number, then all odd-column vectors move up (down), and if the row coordinate is a po-ev (ne-ev) number, then all even-column vectors move up (down). The column coordinates indicate the number of cyclic shifts, whose values are generated by the chaotic map.
For CSTPM, we find all positions of the probabilities greater than 0.25 of the matrix like Table 3. The row coordinates indicate whether odd-or even-row vectors are to be selected, and whether the direction of movement is left or right. If the row coordinate is odd (even), then the odd (even)-row vectors are selected, and if the row coordinate is positive (negative), then all selected elements move left (right). Specifically, if the row coordinate is a po-od (ne-od) number, then all odd-row vectors move left (right), and if the row coordinate is a po-ev (ne-ev) number, then all even-row vectors move left (right). The column coordinates indicate the number of cyclic shifts, whose values are generated by the chaotic map.

Scheme Based on Multiscale Block Compressed Sensing
The proposed encryption scheme is shown in Figure 1. The image is decomposed by discrete wavelet transform to get the coefficient matrices, and each coefficient matrix is scrambled for the first time and measured by the corresponding measurement matrix.
After compression, all matrices of measurement values are merged, and that matrix is scrambled by state transition matrices. The final cipher image is obtained by diffusion.  The encryption process should be closely related to plain information to make it resistant to known-and chosen-plaintext attacks. The hash value K of the plain image is generated by the SHA256 function. K is converted to binary numbers, and 32 groups of these are generated by dividing every eight bits into a group, K = k 1 , k 2 , . . . , k 32 (12) We count the number of zeros, number L of ones, length of the longest continuous zero sequence, and length of the longest continuous ones sequence in the hash value K as L 0 , L 1 , l 0 , and l 1 , respectively, and      r 1 = 1.05 + mod (max(k 1 ,k 2 ,k 3 ,k 4 ,k 5 ,k 6 )×t 1 ) (min(k 1 ,k 2 ,k 3 ,k 4 ,k 5 ,k 6 )×t2) × 10 14 , 2.95 r 2 = 1.05 + mod (max(k 7 ,k 8 ,k 9 ,k 10 ,k 11 ,k 12 )×t 3 ) (min(k 7 ,k 8 ,k 9 ,k 10 ,k 11 ,k 12 )×t4) × 10 14 , 2.95 We construct and find their Kronecker product, where ⊗ represents the Kronecker product operator, 14 , 1 · represents rounding down, and Z 3 (i), i = 1, 2, 3, 4 represents the i-th element of matrix Z 3 .

Calculating Sampling Rates Based on Multiscale Block Compressed Sensing Theory
Given the target sampling rate, the sampling rate sequences of the coefficient matrices of each layer after wavelet decomposition are calculated according to Equations (4)- (7), and denoted as subrates. In this paper, the image is decomposed by three-level wavelet transform, and there are three sampling rates.

Generating the Measurement Matrix
In Section 3.1.1, all parameters and initial values of chaotic maps were generated, then r 2 and y 00 were brought into Equation (9), with t + n 1 × n 1 iterations; r 1 and x 00 were brought into Equation (8), with t + n 2 × n 2 iterations; r 0 and x 0 were brought into Equation (10), with t + n 3 × n 3 iterations. The sinusoidal value of the initial value is multiplied by a small coefficient every 2000 iterations to disturb the initial value of the next iteration, and three chaotic sequences X 1 , X 2 , X 3 are finally generated after the first t numbers are discarded. X 1 , X 2 , and X 3 are one-dimensional vectors of length n 1 × n 1 , n 2 × n 2 , and n 3 × n 3 , respectively, where n 1 , n 2 , n 3 are determined by the block sizes given in advance. To enhance the randomness of chaotic sequences, the generated chaotic sequences are further processed as The generated chaotic sequences are transformed to matrix form, and the corresponding orthogonal bases Φ 1 , Φ 2 , Φ 3 of X 1 , X 2 , X 3 , respectively, are used as redundant measurement matrices. The new row dimensions m i , i = 1, 2, 3 are obtained by multiplying the row dimensions n i , i = 1, 2, 3 of the redundant measurement matrices by the corresponding sampling rates. The first m i rows of Φ 1 , Φ 2 , Φ 3 are extracted as formal measurement matrices Φ 1 , Φ 2 , Φ 3 .

Encrypting the Plain Image
Step 1: After discrete wavelet decomposition of the original image, a low-frequency coefficient and nine high-frequency coefficients in horizontal, vertical, and diagonal directions are obtained. After three-layer wavelet decomposition, the ratio of the original image size M 0 to the block size of each layer is M 0 : M 1 : M 2 : M 3 = 8 : 4 : 2 : 1. The three-level wavelet decomposition is shown in Figure 2. Step 2: Given an array of length 3, the values of elements represent the block sizes. The third-, second-, and first-level wavelet decomposition coefficient matrices H 3 , V 3 , D 3 , H 2 , V 2 , D 2 , and H 1 , V 1 , D 1 are divided into blocks of block_size 3 × block_size 3 , block_size 2 × block_size 2 , and block_size 1 × block_size 1 , respectively. Each block matrix is expanded to a column vector after partitioning. In this way, each block matrix in the third-, second-, and third-level wavelet decomposition coefficient matrices is transformed to a column vector of length block_size 3 2 , block_size 2 2 , and block_size 1 2 , respectively, and nine new coefficient matrices NH 3 , NV 3 , ND 3 , NH 2 , NV 2 , ND 2 ,NH 1 , NV 1 , ND 1 are constructed by merging the corresponding column vectors.
Step 3: The three chaotic sequences X 1 , X 2 , X 3 generated in Section 3.1.3 are sorted in ascending order to obtain the corresponding index sequences Ind 1 , Ind 2 , Ind 3 . The coefficient matrices A 3 , NH 3 , NV 3 , ND 3 , NH 2 , NV 2 , ND 2 , and NH 1 , NV 1 , ND 1 of the third-, second-, and first-level wavelet decomposition, respectively, are scrambled with Ind 3 , Ind 2 , and Ind 1 , respectively, to obtain A 3 , H 3 , V 3 , D 3 , H 2 , V 2 , D 2 , and H 1 , V 1 , D 1 . For example, if A 3 is a matrix of size m × n, we expand it to a sequence of length m × n, i is the index of the i-th element in the sequence, and where A 3 is the scrambled sequence of A 3 .
Step 4: The coefficient matrix A 3 remains unchanged. We compress the coefficient matrices Step 5: Coefficient matrices with the frequency of the three-level wavelet decomposition are combined to form the matrix to be measured, To enhance the randomness, T is processed according to the information l 0 , l 1 , L 0 , L 1 obtained in Section 3.1.1 to obtain Step 6: According to the method proposed in Section 2.3.2, RSTPM and CSTPM are constructed according to the element information of matrix NT.
Step 7: The maximum and minimum values of nine coefficient matrices, where min and max are the minimum and maximum values, respectively, of M. The coefficient matrices H i , V i , D i , i = 1, 2, 3 are quantized to obtain corresponding matrices The low-frequency coefficient matrix A 3 is decomposed by SVD to obtain three sub-matrices u,s,v of the same dimension, whose maximum and minimum values are obtained. The sub-matrices are quantized to the interval [0,255] by Equation (24) to obtain corresponding matrices U, S, VT.
Step 8: We combine the high frequency coefficients of each layer into groups in the same direction: H 1 , H 2 , H 3 are in the first group, V 1 , V 2 , V 3 are in the second group, and D 1 , D 2 , D 3 are in the third group. The index values are obtained using the chaotic sequences X 1 , X 2 , X 3 generated in Section 3.1.3, Lind 1 = f loor mod abs(Lind 1 ) × 10 8 , 6 + 1 Lind 2 = f loor mod abs(Lind 2 ) × 10 8 , 6 + 1 Lind 3 = f loor mod abs(Lind 3 ) × 10 8 , 6 + 1 (25) That is, the index values Lind 1 , Lind 2 , and Lind 3 of the first through third groups, respectively, are determined by the last, penultimate, and antepenultimate elements of X 1 , X 2 , X 3 . The three index values are mapped to the interval [1,6], which indicates that there are six possible permutations for each group, The first through third groups after sorting are being recorded as Y 1 , Y 2 , and Y 3 , respectively, whose internal three block matrices are recorded as Step 9: The quantization matrices U, S, and VT obtained from Step 7 are inserted in the first through third groups of matrices, respectively, after the first round of combination.
The index values are obtained using the three chaotic sequences X 1 , X 2 , X 3 generated in Section 3.1.3, Hind 1 = f loor mod abs(Hind 1 ) × 10 8 , 4 + 1 Hind 2 = f loor mod abs(Hind 2 ) × 10 8 , 4 + 1 That is, the index values Hind 1 , Hind 2 , and Hind 3 of the first through third groups, respectively, are determined by the first through third elements, respectively, of X 1 , X 2 , X 3 . At the same time, the three index values are mapped to the interval [1,4], which indicates that there are four possible permutations for each group: the top, the gaps between two adjacent block matrices, and the bottom, We mark the first through third group matrices after combination as Y 1 , Y 2 , and Y 3 , respectively, and combine them as Step 10: The dimension of the matrix should be adjusted to prevent attackers from obtaining the information of the encryption scheme through the dimension of cipher images. Suppose the dimension of the merged matrix is m × n. We find two factors m 1 and n 1 of m × n such that m 1 × n 1 = m × n, minimize |m 1 − n 1 |, and adjust the dimension of the matrix to m 1 × n 1 , obtaining the matrix information Step 11: The parameters a, b, c and initial value v 0 generated in Section 3.1.1 are brought into Equation (11), iterating t + (d 12 + 1) × m 1 × n 1 times, and the sinusoidal value of the initial state is multiplied by a small coefficient every 2000 iterations to disturb the initial value of the next iteration. The chaotic sequence V is generated after the first t numbers are discarded. To enhance the randomness, we determine the chaotic sequence according to which we generate sub-sequences V 1 , The control parameters of the scrambling process are generated according to V 1 , V 2 , V 3 , where w 1 ∼ w 4 are used for subsequent shift operations, fix is a function that rounds toward zero, and m 1 and n 1 are the largest prime factors of m 1 and n 1 , respectively.
Step 12: The scrambling operations are based on RSTPM, as generated in Step 6. The shift numbers are w 1 -w 4 when column coordinates are the po-od, po-ev, ne-od, and ne-ev numbers, respectively. The scrambling operations are shown in Table 4. The up arrow means move up, the down arrow means move down. We record the scrambled matrix as T and set the initial transition flag bit matrix and all element values to zero. If the element of a position is shifted, then the flag bits of the corresponding position change from 0 to 1. Taking the state transition probability matrix generated in Table 3 as an example, the positions with probability values greater than 0.25 are in the ne-od number and po-ev number columns. That is, the probability values of the two middle columns are selected. Therefore, the flag bits of these two columns are set to 1. Step 13: The scrambling operations are based on CSTPM, as generated in Step 6. The shift numbers are w 1 , w 3 , w 2 , and w 4 , if the column coordinates are the po-od, ne-od, po-ev, and ne-ev numbers, respectively. We record the scrambled matrix as T . We set the state transition flag bits in the same way. If the element of a position is shifted, then the flag bits of the corresponding position change from 0 to 1.
Step 14: The chaotic sequence V generated in Step 11 is quantized to the interval [0,255] to get the new matrix V . The matrix information d 1 and d 2 obtained in Step 10 is used to generate the chaotic sequences, V is sampled at an interval of d 1 to obtain sequence V 1 , and V is sampled at an interval of d 2 to obtain sequence V 1 , We take the position of the last element of sequence V 2 as the start position, and take the consecutive m 1 × n 1 elements from sequence V and record them as sequence V 0 , We perform an XOR operation between matrix T and sequence V 0 , and denote T as A.
Step 15: The matrix A is changed with front addition and a modular operation to obtain the matrix B, using the sequence V 1 generated in Step 14. We perform a cyclic left shift on B to obtain B , i.e., and the definition of BitCircShift is shown in Algorithm 1.
B is changed through back addition and a modular operation to get C, using the sequence V 2 generated in Step 14. We perform a cyclic left shift on C to obtain C , i.e., At this point, the final cipher image is obtained. y2←(x mod 2 −k )×2 8+k 9: end if 10: y←y1 + y2 11: end

Decryption Process
The image decryption scheme is the reverse of image encryption. The process is controlled by the key, including the initial values of chaotic systems, state transition flag bit matrix, and maximum and minimum of matrices, as received by the sender. The process is shown in Figure 3. Step 1: The cipher image matrix is expanded to a sequence C, and then we perform a cyclic right shift operation on C to obtain D, which is changed with inverse back addition and a modular operation to get the sequence D , using the sequence V 2 generated before, We perform a cyclic right shift operation on D to obtain E, which is changed with inverse front addition and a modular operation to get the sequence E by using the sequence V 1 generated before, Step 2: We perform an XOR operation between sequences E and V 0 to obtain sequence E , and transform E to a matrix of size m 1 × n 1 .
Step 3: E is inversely scrambled according to the previously generated transition flag bit matrix, and we find the positions of its elements whose values are 1. We inversely scramble the matrix according to the rule described in Section 2.3.2. The matrix is recorded as E after inversely scrambling row vectors, and that matrix's column vectors are inversely scrambled to obtain E .
Step 4: Find the positions of the U, S, VT and coefficient matrices H i , V i , D i , i = 1, 2, 3 in matrix E according to Hind 1 ∼ Hind 3 and Lind 1 ∼ Lind 3 , then we divide the matrix E into blocks to obtain them. The sub-matrices of low-frequency coefficient matrices U,S,VT are inversely quantized to obtain corresponding matrices U , S , VT , and the coefficient matrices HH i , VV i , DD i , i = 1, 2, 3 are obtained by inverse quantization of corresponding matrices H i , V i , D i , i = 1, 2, 3, where quantization is expressed as where max and min are the maximum and minimum values, respectively, of M. The low-frequency coefficient matrix A 3 is obtained by calculating the product of the three sub-matrices, Step 5: We first generate index vectors Ind 1 ∼ Ind 3 . Coefficient matrices A 3 , HH 3 ,VV 3 , DD 3 , HH 2 , VV 2 , DD 2 , and HH 1 , VV 1 , DD 1 are inversely scrambled with Ind 3 , Ind 2 , and Ind 1 , respectively. For example, we expand A 3 of size m × n to a sequence of length m × n, whose i-th element is where A 3 is the scrambled sequence of A 3 .
Step 6: We reconstruct the image with the SPL algorithm. The column vectors of the reconstructed coefficient matrices of each layer are restored to block matrices according to the corresponding block size. Block matrices of size block_size 3 × block_size 3 , block_size 2 × block_size 2 , and block_size 1 × block_size 1 are, respectively, combined into third-, second-, and first-level wavelet decomposition coefficient matrices A 3 , H 3 , V 3 , D 3 , H 2 , V 2 , D 2 , and H 1 , V 1 , D 1 , which are combined into a new matrix according to Figure 2.
Step 7: The matrix after inverse wavelet transform is the final decrypted image.

Encryption and Decryption Results
Simulation experiments were carried out on a laptop computer with an Intel Core i5-6200U CPU at 2.3 GHz and 4 GB RAM, on the MATLAB R2015a platform. Images of size 512 × 512, including Lena, Goldhill, Cameraman, Peppers, Barbara, and Jet, were selected as test images. The external keys were t 1 = 0.11, t 2 = 0.22, t 3 = 0.33, t 4 = 0.44, t 5 = 2.723, t 6 = 0.618, t 7 = 3.141, t 8 = 4.6692, t = 600. The array block_size, which determines the block sizes of the image matrix, was set to [8,16,32], i.e., block_size 1 = 32, block_size 2 = 16, block_size 3 = 8. The target sampling rate was 0.25. The experimental results of plain, cipher, and decrypted images are shown in Figure 4. The cipher images were meaningless and unrecognizable noise-like images with little connection to the plain images. Valid information about the original images cannot be obtained from the corresponding cipher images. The encryption scheme adjusted the dimensions of cipher images to 288 × 256, which further hid the information of the plain image and the encryption method. The decrypted images were meaningful images that could be clearly identified and were quite similar to the original images. Hence, the proposed scheme had good encryption and decryption effects.

PSNR between Plain and Decrypted Images under Different Sampling Rates
The peak signal to noise ratio (PSNR) is used to objectively judge the quality of a decrypted image. A larger PSNR indicates a smaller difference between plain and decrypted images, and higher reconstruction accuracy. For gray images, where M and N are the row dimension and column dimension of the image, and X(i, j) and Y(i, j) are the pixel values of the plain image X and decrypted image Y, respectively, at position (i, j). In this experiment, 512 × 512 images, including Lena, Cameraman, Peppers, and Couple, were selected as test images, and PSNR values between plain and decrypted images were measured at different sampling rates. Table 5 compares partial experimental results with those in Gan et al. [10]. The scheme of Gan et al. [10] adopts traditional compressed sensing theory, using one sampling rate for the whole image, which reveals that the distribution differences of image information between low-and high-frequency coefficients are not fully considered, and the resulting decrypted image does not have high reconstruction quality. The image reconstruction effect of the proposed scheme is better than that of the scheme in Gan et al. [10]. Table 6 compares some experimental results with those in Luo et al. [28]. In the latter, the plain image is first decomposed into approximate and detail components by discrete wavelet transform (DWT). All pixels in the approximate component are retained, and the remaining detailed components are measured by measurement matrices. A lower sampling rate is adapted to the horizontal direction decomposition coefficient (LH) and diagonal direction decomposition coefficient (HH), and a larger sampling rate is adapted to the vertical direction decomposition coefficient (HL). The reconstruction quality of the decrypted image is improved. However, there are some limitations because the sampling rates are set artificially. From Table 6, it can be seen that the proposed algorithm can better improve the reconstruction quality of the decrypted images.

Influence of Wavelet Basis on Image Reconstruction Effect (PSNR)
In this experiment, 12 images were decomposed by DWT in different wavelet bases, including the commonly used Symlets8, Haar, and CDF9/7. The PSNR values of the decrypted images were measured, with results as shown in Table 7. The experimental results show that the CDF9/7 wavelet base performs better in most cases; hence, this was chosen as the wavelet base.

Time Complexity Analysis
An algorithm should have fast encryption and decryption speeds to meet real-time needs. In this experiment, the plain image Lena with size of 512 × 512 was used as a test image, the sampling rate was set to 0.25, and the running time of each process of the scheme was measured. The results are shown in Table 8. In addition, the running times of the whole encryption process and decryption process at different sampling rates were also measured and the results are shown in Table 9. We can learn from the experimental result that the processes of iteration of chaotic systems and image reconstruction consume most of the time in the algorithm. Additionally, as the sampling rate increases, the encryption time also increases, but the decryption time is basically the same. More concretely, in encryption algorithm illustrated in Section 3.1, iteration of chaos includes Section 3.1.3 and Step 11 in Section 3.1.4, the compression process is Step 4 in Section 3.1.4, the first round of permutation is Step 3 in Section 3.1.4, the second round of permutation is Step 12 and Step 13 in Section 3.1.4, and the diffusion includes Step 14 and Step 15. The decryption process illustrated in Section 3.2 is the inverse process of the encryption process, and the compression process is replaced by the reconstruction process.  Table 9. Runtime statistics of encryption process and decryption process at different sampling rates. In what follows, we analyzed the time complexity of our encryption algorithm in Section 3.1 in detail. Assume the size of plain image is m × n, the block sizes of the third-, second-, and first-level wavelet decomposition coefficient matrices are n 1 , n 2 , n 3 , respectively, and the target sampling rate is CR. Section 3.1.3 is to generate three chaotic sequences for scrambling, and time complexity is Θ n 1 4 + n 2 4 + n 3 4 .

Sampling
Step 11 in Section 3.1.4 is to generate the chaotic sequence for diffusion, and time complexity is Θ(CR × m × n).
Step 3 in Section 3.1.4 is to scramble the coefficient matrices with index sequences, and time complexity is Θ(m × n).
Step 6 in Section 3.1.4 is to generate the RSTPM and CSTPM, and time complexity is Θ(2 × CR × m × n).
Step 12 and Step 13 in Section 3.1.4 are to scramble the matrix of measurement values after CS with RSTPM and CSTPM, and time complexity is Θ(CR × m × n).
Step 14 and Step 15 in Section 3.1.4 are to diffuse the matrix after scrambling, and time complexity is Θ(3 × CR × m × n). For the computational cost of the proposed method determined in other steps, the time complexity is about Θ(3 × CR × m × n). From the result shown in Table 8, the total time complexity is approximately equal to Θ(5 × m × n). Comparing with the encryption algorithms in [49,50] listed in Table 10, our scheme has a smaller time complexity. However, the time complexity of our scheme is equal to that of [10]. Table 10. Comparison results on time complexity of Encryption algorithm.

Key Space
To withstand a brute-force attack, an encryption scheme should have a large key space. If the calculation accuracy of the computer is 10 −14 , the external key is t 1 ∼ t 8 , accounted for 10 14 8 = 10 112 key space. Table 11 compares the proposed scheme and other schemes. Since the ideal key space is suggested to be at least 2 100 < 10 31 for a good cryptosystem [9], the result illustrates that the key space of the proposed scheme is large enough to resist all kinds of attacks.

Histogram Analysis
A histogram is an effective index to evaluate the distribution of pixel values. With an effective image encryption scheme, the histogram of a cipher image should be evenly distributed, so as to effectively resist statistical attacks. In this experiment, the histograms of the images in Figure 4 were drawn, with results as shown in Figure 5, from which it can be seen that the pixels of cipher images are uniformly distributed, and are quite different from those of the original images, making it impossible to obtain useful information. The histograms of the reconstructed images are similar to those of the corresponding plain images, which indicates that the reconstruction effect of decrypted images is good. All in all, attackers can obtain no useful information about plain images through statistical attacks when the proposed scheme is used.

Sensitivity Analysis
Key sensitivity and plain sensitivity are important metrics to evaluate cryptosystems. Weak sensitivity enables the easy attack of a cryptosystem. We tested the key sensitivity and plain sensitivity of our scheme.

Plain Sensitivity
In a differential attack, the attacker encrypts an original and modified plain image with the same key and determines the relationship between the plain and encrypted images by comparing the corresponding cipher images. To resist such attacks, an encryption algorithm should have strong plain sensitivity, i.e., two plain images with small differences should have significant differences after encryption. The main indicators to measure the difference between two images include number of pixels change rate (NPCR) and unified average changing intensity (UACI) [51], which respectively reflect the proportion of the number of different pixels to the size of the image and the average ratio of the differences between pixels at corresponding positions and 255. These are expressed as NPCR(P 1 , where sgn(·) is the sign function, and UACI(P 1 , The critical values of NPCR and UACI of the two random images were 99.6094% and 33.4635%, respectively. For NPCR, the more the critical value is exceeded, the stronger the plain sensitivity of the encryption scheme. For UACI, the closer to the critical value, the stronger the plain image sensitivity of the encryption scheme. In the experiment, one pixel in plain image P 1 was randomly selected and the change of its pixel value was set to 1 to obtain the modified plain image P 2 . Cipher images C 1 and C 2 were obtained by encrypting plain images P 1 and P 2 , respectively, with the same key, and NPCR and UACI were calculated, as shown in Table 12. The experimental results show that NPCR exceeded the critical value for all images except Peppers, and UACI was close to the critical value for all images. Hence, the proposed encryption scheme has a strong resistance to differential attacks.

Key Sensitivity
An encryption scheme should have high sensitivity to the secret key in both the encryption and decryption processes [9], i.e., a small change of the secret key should produce a completely different cipher image. Similarly, in decryption, the plain image should not be recovered by a decryption key differing slightly from the encryption key. We tested the key sensitivity from the aspects of sensitivity in both encryption and decryption.
The initial values of the four chaotic systems used in this paper were changed slightly from K = x 0 , x 00 , y 00 , v 0 to K1 = x 0 + 10 −14 , x 00 , y 00 , v 0 , K2 = x 0 , x 00 + 10 −14 , y 00 , v 0 , K3 = x 0 , x 00 , y 00 + 10 −14 , v 0 , and K4 = x 0 , x 00 , y 00 , v 0 + 10 −14 . The Lena image was encrypted with the correct and wrong keys, with results as shown in Figure 6. From the experimental results, it can be seen that a small change in the key can cause a huge change in the cipher image but will not reveal information related to the plain image. To quantify the differences between cipher images obtained from the same plain image, NPCR and UACI were calculated using the correct and wrong keys to encrypt images, with results as shown in Table 13. It can be seen that NPCI and UACI are all close to the theoretical values, and more than 99.5% of pixels were modified, which implies that the encryption process is highly sensitive to the secret key. Figure 6. Experimental results of key sensitivity in encryption process: (a) cipher image using correct key K; (b,d,f,h) cipher images using wrong keys K1-K4, respectively; (c,e,g,i) differential images between (a,b), (a,d), (a,f), (a,h), respectively. Next, a cipher image was decrypted with both the correct and wrong keys, with results as shown in Figure 7, from which it can be seen that only through the correct key can we get the correct decrypted image. Even a slight key change produces a visually meaningless image. Table 14 lists the values of NPCR from Figure 7b-f. When the key has a tiny change, more than 99.8% of pixels are altered, which indicates that the decryption process is highly sensitive to the secret key.

Correlation Coefficients
Strong correlations, close to 1, exist between the pixels of natural images. The correlation coefficient reflects the linear relationship between adjacent pixels and is an important index to evaluate image encryption schemes. For a good encryption algorithm, correlations between adjacent pixels of cipher images should be very weak, tending to zero, which means that the correlation between pixels is largely eliminated. The correlation coefficient of u and v is calculated as [52] where N is the number of adjacent pixel pairs from the image, and (u i , v i ), i = 1, 2, . . . , N is the gray value of a pair. In this experiment, 512 × 512 images Lena, Goldhill, and Peppers were selected as test images, and 5000, 6000 or 8000 pairs of adjacent pixels were randomly selected from the cipher images in the horizontal, vertical and diagonal directions for calculation. To reduce randomness, each calculation was carried out 100 times, the final result was taken as the average value, and this was compared with other schemes. The experimental results are shown in Table 15, and show that the correlations between adjacent pixels of the cipher images were small. We also determined adjacent pixel distributions for plain and cipher images, with results as shown in Figure 8, showing that adjacent pixels of the plain images were basically linearly distributed. However, there were weak correlations between adjacent pixels of the compressed cipher images, with correlation coefficients close to 0, showing that the proposed scheme has a good encryption effect.

Information Entropy
Information entropy reflects the randomness of the distribution of image pixels and is calculated as [37] where p(s i ) is the probability of the i-th pixel value s i of image p, and n is its total number of digits. For an 8-bit grayscale image, its cipher image should have information entropy near 8 bits. In this experiment, the 512 × 512 images Lena, Peppers, Cameraman, Barbara, and Jet were selected as test images. The experiment was carried out at sampling rates of 0.5 and 0.25, with results as shown in Table 16, from which it can be seen that the information entropies of cipher images were all greater than 7.99. The results show the encryption effect of the proposed scheme was better than that of comparison schemes [10,40,43], which indicates that our scheme has better randomness than other schemes based on compressed sensing. Besides, compared to the algorithm proposed by Brindha et al. [53], our scheme can achieve higher information entropy. It shows that although our lossy compression algorithm using compressed sensing cannot fully restore the original image from the cipher image, it can obtain more secure cipher images than the lossless compression algorithm. What is more, it can be seen that the cipher images generated from our encryption scheme have greater information entropies than those of schemes based on substitution box [14] and quantum key image [28] in most instances, which means our scheme is more resistant to entropy attack than some encryption schemes based on non-chaotic techniques.

Robustness Analysis
During the transmission process, cipher images can be contaminated by noises and data loss, making it hard to decrypt them and recover the corresponding plain images. In this experiment, noise and crop attacks were utilized to test the robustness of the proposed scheme.

Noise Attack
Common types of noise include Gaussian noise (GN), speckle noise (SN), and salt-andpepper noise (SPN). To evaluate the robustness of the proposed scheme to noise attacks, Lena was encrypted as the test image, different levels of three noises were added, and decrypted images were obtained, with experimental results as shown in Figure 9 from which it can be seen that the proposed scheme has a stronger resistance to speckle noise and salt-and-pepper noise than Gaussian noise.

Crop Attack
Different areas of cipher images were randomly selected after encrypting Lena, cutting image blocks of size 16 × 16, 32 × 32, 64 × 64, 288 × 16,16 × 256 to obtain cipher images with some data loss, along with corresponding decrypted images. The results are shown in Figure 10, and numerical quantization results of decrypted images are shown in Table 17. From the results, although a crop attack will worsen the visual effect of decrypted images, the proposed scheme still retains most information of the image, i.e., the scheme has a certain robustness to crop attacks.

Conclusions
An image encryption scheme based on multiscale block compressed sensing theory was proposed. In the scheme, considering that different information is carried by the lowand high-frequency coefficients of an image, different sampling rates were set for the lowand high-frequency coefficients of the image, so as to improve the reconstruction quality of decrypted images. Our scheme was experimentally compared with a scheme using the traditional compressed sensing theory and setting a sampling rate in each encryption process. Under the same compression ratio, the PSNR values between natural images and corresponding decrypted images of our scheme were better by more than 8 dB and were better by about 3-5 dB over a scheme that sets different sampling rates. Therefore, the proposed encryption scheme is suitable for natural image transmission with complex structures and large amounts of information. The safety of images is ensured, image information can be better preserved, and better visual effects can be obtained. With the combination of chaotic systems and a Markov model, the image is scrambled inside each coefficient matrix and is then scrambled among the coefficient matrices, and the encryption is completed by the strategies of independent and global diffusion. Experimental results show that the proposed scheme has a large key space, high plain sensitivity, and high key sensitivity in both the encryption and decryption processes. In particular, compared with an encryption scheme designed for problems of low entropy, the experimental scheme in this paper has a certain increase in information entropy values of most natural images, which shows the Markov probability model has certain advantages over traditional scrambling methods in simulating random processes. Therefore, the scheme has a certain guiding role for our future research work. At the same time, the scheme can effectively resist brute-force, differential, and statistical attacks. It also has a certain robustness to noise and crop attacks. However, high-intensity attacks cause large information loss, resulting in a poor visual effect. Hence, the design of a more robust encryption scheme warrants further study. Due to the use of multiscale block compressed sensing theory that better suits images, this scheme is suitable for mass image data transmission and can effectively save transmission space. To be honest, the times of the encryption and decryption processes in our scheme are long. Therefore, it remains to study the replacement of existing methods by more efficient chaotic systems and reconstruction methods while maintaining good encryption and decryption effects. Data Availability Statement: Data are available upon reasonable request. Additional data (beyond those included in the main text) are available from the corresponding author upon request. The data are not publicly available due to [the data also forms part of an ongoing study].

Conflicts of Interest:
The authors declare no conflict of interest.