Simultaneous encryption and compression of medical images based on optimized tensor compressed sensing with 3D Lorenz

Background The existing techniques for simultaneous encryption and compression of images refer lossy compression. Their reconstruction performances did not meet the accuracy of medical images because most of them have not been applicable to three-dimensional (3D) medical image volumes intrinsically represented by tensors. Methods We propose a tensor-based algorithm using tensor compressive sensing (TCS) to address these issues. Alternating least squares is further used to optimize the TCS with measurement matrices encrypted by discrete 3D Lorenz. Results The proposed method preserves the intrinsic structure of tensor-based 3D images and achieves a better balance of compression ratio, decryption accuracy, and security. Furthermore, the characteristic of the tensor product can be used as additional keys to make unauthorized decryption harder. Conclusions Numerical simulation results verify the validity and the reliability of this scheme.

Existing simultaneous compression and encryption algorithms are typically applied to ordinal images rather than medical images, because they refer to lossy compression. Alfalou et al. proposed a series of representative algorithms for the simultaneous compression and encryption of 3D images. The latest and most effective algorithm is based on spectral fusion and discrete cosine transform (DCT) [8]. However, the decryption error increases rapidly along with the increase of the number of images, indicating that it cannot handle large amounts of images simultaneously. Emerging algorithms for simultaneous encryption and compression are based mainly on compressed sensing (CS), which can activate compression during the sampling process [9]. An example being Valerio et al., who proposed a multiclass encryption by CS to withstand the common attack [10,11]. To further enhance security, CS-based encryption algorithms were constructed by combining the chaos map and some optical encryption techniques, such as double random phase encryption (DRPE), fractional Fourier transform (FrFT), and fractional Mellin transform (FrMT) [12][13][14][15][16][17][18][19].
The medical images are different from other images because of their particular properties. There are legal and strict regulations applied to medical multimedia information due to the health of a patient depending on the correctness and accuracy of this information [20]. The quality of the decrypted and compressed data must be adequate to allow for a correct diagnosis when it is reconstructed. However, the existing algorithms did not meet this requirement because most of them have not been applicable to 3D images intrinsically represented by tensors. Conventional CS theory relies on data representation in the form of one-dimensional vectors. Application of CS to higher dimensional data representation is typically performed by conversion of the data to very long vectors that must be measured using very large sampling matrices, thus destroying the intrinsic structure and imposing a huge memory burden.
Recently, Cesar et al. propose a tensor CS (TCS) based on higher order singular value decomposition (HOSVD) that introduces a direct reconstruction formula to recover a tensor from a set of multi-linear projections, which are obtained by multiplying the data tensor by a different sensing matrix in each mode [21]. This HOSVD-based TCS achieved more accurate and efficient reconstruction results when compared to other existing sparsity-based TCS methods [22][23][24]. Indeed, we believe that, if it was used to design TCS-based encryption of 3D images better performance would occur.
We further introduce alternating least squares (ALS) into HOSVD-based TCS [21], and control the measurement matrices by 3D Lorenz [25][26][27]. Such a simultaneous compression and encryption algorithm for 3D medical images has two main advantages: (1) the preservation of intrinsic structures of the tensor data for the purpose of reducing the decryption error and increasing the compression ratio; (2) the keys consist of those generated by tensor decomposition, and 3D Lorenz. Particularly, the order of the tensor product used in the TCS can be used as additional keys to make unauthorized decryption harder.
This paper is organized as follows: in "Theory" section, the related notation, definition and basic results used throughout the paper, are introduced; in "Proposed encryption" section, the encryption and decryption algorithms are proposed; in "Numerical simulation results" section, several numerical results based on 3D lung CT images are provided, to corroborate our theoretical results and evaluate the stability and robustness 1. Mode-i product: A mode-i product of a tensor A and a matrix Φ ∈ R M i ×m is denoted by A × i Φ and is of size m × (M 1 · · · M i−1 · M i+1 · · · M d ) matrix. 2. Mode-i unfolding: The mode-i unfolding A (i) of A arranges the mode-i fibers to be the columns of the resulting matrix. 3. HOSVD: The decomposition and reconstruction of A can be written as the product: where Φ i ∈ R M i ×m i , and W is a complex (m 1 × m 2 × · · · × m d )-tensor of which the subtensors obtained by corresponding singular values. 4. Tucker-TCS: in [21], a more stable, robust and accuracy tensor reconstruction of CS is proposed: where " †" stands for the MP pseudo-inverse of a matrix. We assume that the following sets of compressive multi-way measurements Z (n) are available:

3D Lorenz
Mohmad et al. [25,26] propose a 3D discrete Lorenz system, which has a high order and also low complexity when implemented in digital hardware. The discrete Lorenz attractor employed here is given by the following difference equations where U, V, W are three state variables, A, B, C are parameters, and g1, g2, g3 are gains (step size). The calculating method of (5) is finite difference.

Proposed encryption
The proposed system can be split into encryption and decryption algorithms as illustrated in Fig. 1. The compression and decompression procedures are embedded in the encryption and decryption, respectively.

Encryption
For the initial 3D image A ∈ R M 1 ×M 2 ×M 3 , the encryption process consists of the following steps: 1. Initialize randomly the three Gaussian sensing matrices To accurately decrypt A, the optimal Φ i should satisfy: Sensing matrices  15:118 As in [28], this problem can be converted to The reconstruction algorithm as Eq. (2) of HOSVD-based TCS achieved a more accurate solution than that as Eq. (1) of HOSVD. To further improve the reconstruction accuracy, we use an alternating least squares (ALS) approach to solve Eq. (7). 2. For k = 0, iterate Eqs. (8)-(11) until Φ i converges or the maximum iteration is achieved: Then the optimal Φ * 1 , Φ * 2 and Φ * 3 are obtained. 3. Compute the compressed core tensor 4. Unfold W into its n-mode W (n) . The mode n is a private key which has three possible values: 1, 2 and 3. 5. Φ 1 , Φ 2 and Φ 3 are synchronously constructed by 3D Lorenz as Eq. (5).
Hence the compression ratio is given by: The details of how to synchronize the image by 3D Lorenz system are introduced below. One data sample ϕ i,k is inserted into U, which gives . The initial conditions U 0 , V 0 , W 0 , parameters A, B, C, and g1, g2, g3 are known by both transmitter and receiver. The transmitted signal is the U state variable, and the objective is to retrieve ϕ i,k from this signal at the receiver. Feedback is used to update the state variables at the receiver to synchronize the system and allow decryption of subsequent data values.
After φ i,k is obtained, the receiver state equations can be updated, thus achieving synchronization with the transmitter.

Decryption
The decryption process consists of the following steps: 1. E Φ i are inverse transformed by 3D Lorenz: is computed during Eq. (16). 2. W (n) and the obtained D Φ i are multiplied in the correct order to recover A ′ . There are three feasible ways to achieve this: where '⊗' represents the Kronecker product. 3. Then, fold A ′ (n) into A ′ according to the private key n.
It is obvious that besides the secret keys of measurement matrices, the unfolding model n (order of tensor product) can be used as an additional key.

Numerical simulation results
Numerical simulations were conducted with Matlab2011 on a work station with an Intel Core i7 CPU and 64 GB RAM. The decryption error and compression ratio of the proposed system are introduced in "Decryption accuracy and compression ratio" section, the histograms are analyzed in "Histograms and statistical analysis" section, the secret keys are illustrated in "Rate-distortion" section, and the robustness is stated in "Secret keys" section.

Decryption accuracy and compression ratio
Our experiments are conducted on lung CT sequences in lung image database consortium (LIDC) [31]. Each frame of one CT sequence is preprocessed to have 512 × 512, where 512 frames were chosen. The CT sequence together is represented by a 512 × 512 × 512 tensor and has 134,217,728 voxels in total. The randomly constructed Gaussian measurement matrix for each mode is now of size 512 × m i (i = 1, 2, 3). Thus, the compression ratio is given by: The quantitative measure of the decryption error is the peak signal noise ratio (PSNR), which is based on the root mean square error (RMSE) between the decrypted data and ground truth and can be represented as: To further evaluate the performance of decryption, the structural similarity index (SSIM) is used as another indicator.
In this section, we compare the proposed algorithm with state of the art algorithms to show its superiority. These algorithms are presented briefly as follows: 1. As demonstrated previously [13], an encryption based on 2D_CS in the FrMT domain (algorithm 1) stands out for its efficient, robust, and secure encryption performance. In this algorithm, the 2D CS is based on a 2D wavelet, measuring matrices with Logistic map and a 2D NSL 0 reconstruction algorithm. Notably, although the security of FrMT is better than that of FrFT, its decryption accuracy is less than the later. In order to verity the decryption accuracy and compression ratio of our tensorbased algorithm, we replaced FrMT with FrFT in algorithm 1. 2. Additionally, we chose an encryption algorithm based on HOSVD-TCS [21] with FrFT (algorithm 2). 3. Also, as previously shown [8], an encryption algorithm based on spectral fusion and DCT obtained a better PSNR when compared with previous compression-encryption implementations. Accordingly, this became algorithm 3.
A visual evaluation of the decryption results under frames 117, 138 and 159 of one CT sequence at the compression γ = 0.125 is shown in Figs. 2, 3, 4, 5, 6. As shown in Fig. 3, all the tissues within the lung volume are clear. Our clinical experts did not find distinct differences between the decrypted and the original CT images. The quantitative summaries of the above algorithms are shown in Tables 1, 2 and Fig. 7, where the advantages of the proposed algorithm are highlighted. It is evident that the advantages of the proposed algorithm over the other methods increase with the compression ratio. Case in point is in algorithm 3, where improving the compression ratio requires a large number of frames. However, the PNSR rapidly decreases with the increase of the number of frames. Thus, the algorithm cannot handle large number of frames. We also compared the computation times (the average computation time of the experiments with all compression ratio and noise level) required in each case. The comparison shows that algorithm 3 provides a much faster computation (Table 3), while the proposed algorithm requires slightly longer computation time because it contains a procedure of iteration. Fortunately, this iteration is much simpler than those of algorithms 1 and 2.  15:118 Histograms and statistical analysis

Histograms
The histograms of two CT sequences and the encrypted images of their composed parts are shown in Figs. 8 and 9, respectively. The intensity distribution of the histograms of the encrypted images is completely dissimilar from that of the histogram of the original CT, which indicates that an intruder cannot perceive any useful information based on statistical properties. The histograms of the two original CT sequences are evidently

Statistical analysis
Statistical properties of the images can also be evaluated by the computation of the correlation between two adjacent pixels. By selecting randomly P pixels of the image, the correlation coefficient is computed as: Page 11 of 20 Wang et al. BioMed Eng OnLine (2016) 15:118 where x is the value of a selected pixel and y is the value of the correspondent adjacent pixel, D(x) is the mean square error.
It is expected that an image will have a correlation coefficient close to 1 before being submitted to the encryption; it is desirable that the correlation coefficient of a ciphered image be as close to 0 as possible. In Table 4, where the simulation results for P = 125,000 are shown, we verify that the above described premise is satisfied. This indicates that the proposed encryption scheme is secure against statistical attacks.

Normalized entropy
The normalized entropy [6] of the ciphered image is defined as    where P is the number of different values that the pixels of the ciphered image can assume, N i is the amount of pixels of the ciphered image that assume value i, and N is the total amount of pixels of the ciphered image. During the experiments, we found that the ciphered pixels' intensity of all the algorithms mentioned in the paper are almost equiprobable, so their normalized entropy are all approximate to 1.

Rate-distortion
Rate-distortion (RD) is an important indicator to evaluate the performance of compression. In a previous study [32], the following mathematical model of CS RD was constructed: The proposed TCS-based algorithm falls into the CS category, so the model in Eq. (23) applies for our algorithm. Because there are some differences between traditional vectorbased CS and the proposed tensor-based CS, we recalculate some parameters as follows:  256 of the original and fake CT sequences are shown in Fig. 12a, b, respectively. The attack result using fake decryption keys with all correct parameters is shown in Fig. 12c. It can be seen that the retrieved image is noise-like signal. We first analyze the number of pixels rate: all pixels assumed to be an even distribution, because the retrieved image