# Principal Component Analysis Using Structural Similarity Index for Images

Despite the advances of deep learning in specific tasks using images, the principled assessment of image fidelity and similarity is still a critical ability to develop. As it has been shown that Mean Squared Error (MSE) is insufficient for this task, other measures have been developed with one of the most effective being Structural Similarity Index (SSIM). Such measures can be used for subspace learning but existing methods in machine learning, such as Principal Component Analysis (PCA), are based on Euclidean distance or MSE and thus cannot properly capture the structural features of images. In this paper, we define an image structure subspace which discriminates different types of image distortions. We propose Image Structural Component Analysis (ISCA) and also kernel ISCA by using SSIM, rather than Euclidean distance, in the formulation of PCA. This paper provides a bridge between image quality assessment and manifold learning opening a broad new area for future research.

## Authors

• 30 publications
• 27 publications
• 32 publications
• ### Locally Linear Image Structural Embedding for Image Structure Manifold Learning

Most of existing manifold learning methods rely on Mean Squared Error (M...
08/25/2019 ∙ by Benyamin Ghojogh, et al. ∙ 0

• ### Incorporating prior knowledge about structural constraints in model identification

Model identification is a crucial problem in chemical industries. In rec...
07/08/2020 ∙ by Deepak Maurya, et al. ∙ 0

• ### Fast PET Scan Tumor Segmentation using Superpixels, Principal Component Analysis and K-means Clustering

Positron Emission Tomography scan images are extensively used in radioth...
10/18/2017 ∙ by Yeman B. Hagos, et al. ∙ 0

• ### Robust LSB Watermarking Optimized for Local Structural Similarity

Growth of the Internet and networked multimedia systems has emphasized t...
03/13/2018 ∙ by Amin Banitalebi, et al. ∙ 0

• ### Quality assessment metrics for edge detection and edge-aware filtering: A tutorial review

The quality assessment of edges in an image is an important topic as it ...
01/01/2018 ∙ by Diana Sadykova, et al. ∙ 0

• ### Principal component-based image segmentation: a new approach to outline in vitro cell colonies

The in vitro clonogenic assay is a technique to study the ability of a c...
03/10/2021 ∙ by Delmon Arous, et al. ∙ 0

• ### Spherical Principal Component Analysis

Principal Component Analysis (PCA) is one of the most important methods ...
03/16/2019 ∙ by Kai Liu, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

It has been shown that Mean Squared Error (MSE) is not a promising measure for image quality, fidelity, or similarity [wang2009mean]. The distortions of an image or similarities of two images can be divided into two main categories, i.e., structural and non-structural distortions [wang2004image]. The structural distortions, such as JPEG blocking distortion, Gaussian noise, and blurring, are the ones which are easily noticeable by Human Visual System (HVS), whereas the non-structural distortions, such as luminance enhancement and contrast change, do not have large impact on the visual quality of image.

Structural similarity index (SSIM) [wang2004image, wang2006modern] has been shown to be an effective measure for image quality assessment. It encounters luminance and contrast change as non-structural distortions and other distortions as structural ones. Due to its performance, it has recently been noticed and used in optimization problems [brunet2018optimizing] for tasks such as image denoising, image restoration, contrast enhancement, image quantization, compression, etc, noticing that the distance based on SSIM is quasi-convex under certain conditions [brunet2012mathematical].

So far, the fields of manifold learning and machine learning have largely used MSE and Euclidean distance in order to develop algorithms for subspace learning. Principal Component Analysis (PCA) is an example based on Euclidean distance or norm. However, MSE is not as promising as SSIM for image structure measurement [wang2009mean, wang2006modern] making these algorithms not effective enough in terms of capturing the structural features of image. In this paper, we introduce the new concept of image structure subspace

which is a subspace capturing the intrinsic features of an image in terms of structural similarity and distortions, and can discriminate the various types of image distortions. This subspace can also be useful for parameter estimation for (or selection between) different denoising methods, but that topic will be dealt with in future work.

The outline and contributions of the paper are as follows: We begin by defining the background methods of SSIM and PCA. We then introduce ISCA using orthonomal bases and kernals by analogy to PCA, where ISCA can be seen as PCA which uses SSIM instead of the norm. We then describe an extensive set of experiments demonstrating the performance of ISCA on projection, reconstruction and out-of-sample analysis tasks compared to various kernel PCA methods. The derivations of expressions in this paper are detailed more in the supplementary-material paper which will be released in https://arXiv.org.

## 2 Structural Similarity Index

The SSIM between two reshaped image blocks and , in color intensity range , is [wang2004image, wang2006modern]:

 R∋SSIM(˘x1,˘x2):=(2μx1μx2+c1μ2x1+μ2x2+c1)(2σx1σx2+c2σ2x1+σ2x2+c2)(σx1,x2+c3σx1σx2+c3), (1)

where , , , , , and and are defined similarly for . In this work, . The , , and are for avoidance of singularity [wang2006modern] and is the dimensionality of the reshaped image patch. Note that since , we can simplify SSIM to , where and

. If the vectors

and have zero mean, i.e., , the SSIM becomes , where [otero2014unconstrained]. We denote the reshaped vectors of the two images by and , and a reshaped block in the two images by and . The (squared) distance based on SSIM, which we denote by , is [otero2014unconstrained, brunet2012mathematical, brunet2011class]:

 R∋||˘x1−˘x2||S:=1−SSIM(˘x1,˘x2)=||˘x1−˘x2||22||˘x1||22+||˘x2||22+c, (2)

where . In ISCA and PCA which inspires ISCA, the data should be centered; therefore, the fact that and should be centered is useful.

## 3 Principal Component Analysis

Since ISCA is inspired by PCA [jolliffe2011principal] we briefly review it here. Assume that the orthonormal columns of matrix are the vectors which span the PCA subspace. Then, the projected data onto PCA subspace and the reconstructed data are and , respectively. The squared length of the projected data is where and denote the trace and Frobenius norm of matrix, respectively. Presuming that the data are already centered, the is the covariance matrix; therefore: . Maximizing the squared length of projection where the projection matrix is orthogonal is:

 maximizeU tr(U⊤SU), (3) subject to U⊤U=I,

The Lagrangian [boyd2004convex] is: , where is a diagonal matrix including Lagrange multipliers. Equating derivative of to zero gives us: . Therefore, columns of

are the eigenvectors of the covariance matrix

.

PCA can be looked at with another point of view. The reconstruction error is where is the matrix of residuals. We want to minimize the reconstruction error:

 minimizeU ||X−UU⊤X||2F, (4) subject to U⊤U=I.

The objective function is . The Lagrangian [boyd2004convex] is: , where is a diagonal matrix including Lagrange multipliers. Equating the derivative of to zero gives:

, which is again the eigenvalue problem for the covariance matrix

. Therefore, PCA subspace is the best linear projection in terms of reconstruction error.

As shown above, PCA [jolliffe2011principal] is based on (or Frobenius) norm which is not a promising measure for image quality assessment [wang2009mean]. In order to have both the minimization of reconstruction error as in PCA and using a proper measure for image fidelity, we propose ISCA.

## 4 Image Structural Component Analysis (ISCA)

### 4.1 Orthonormal Bases for One Image

Our goal is to find a subspace spanned by directions for some desired . Consider an image block which is centered (its mean is removed). We want to project it onto a -dimensional subspace and then reconstruct it back, where . Assume is a matrix whose columns are the projection directions spanning the subspace. The projection and reconstruction of are and , respectively. We want to minimize the reconstruction error with orthonormal bases of the subspace; therefore:

 minimizeU∈Rq×p ||˘x−UU⊤˘x||S, (5) subject to U⊤U=I.

According to Eq. (2) and noticing the orthonormality of projection directions, , we have:

 R∋f(U):=||˘x−UU⊤˘x||S=˘x⊤(I−UU⊤)˘x˘x⊤(I+UU⊤)˘x+c. (6)

 Rq×p∋G(U):=∂f(U)∂U=−2(1+f(U))||˘x||22+||UU⊤˘x||22+c˘x˘x⊤U. (7)

We partition a -dimensional image into non-overlapping blocks each of which is a reshaped vector . The parameter is an upper bound on the desired dimensionality of the subspace of block (). This parameter should not be a very large number due to the spatial variety of image statistics, yet also not very small so as to be able to capture the image structure. Also note that is an upper bound on the rank of .

We have instances of -dimensional subspaces, one for each of the blocks. For projecting an image into the subspace and reconstructing it back, one can project and reconstruct every block of an image separately using the bases of the block subspace. The overall bases of an image can be visualized in image-form by putting the bases of blocks next to each other (see the experiments in Section 6).

Considering all the blocks in an image, the problem in Eq. (5) becomes:

 % minimizeUi∈Rq×p b∑i=1||˘xi−UiU⊤i˘xi||S, (8) subject to U⊤iUi=I,   ∀i∈{1,…,b},

where and are the -th block and the bases of its subspace, respectively. We can embed the constraint as an indicator function in the objective function [boyd2011distributed]:

 minimizeUi,Vi∈Rq×p b∑i=1(f(Ui)+h(Vi)), (9) subject to U−V=0,

where and . The denotes the indicator function which is zero if its condition is satisfied and is infinite otherwise. The and are defined as union of partitions to form an image-form array, i.e., and [otero2018alternate].

The Eq. (9) can be solved using Alternating Direction Method of Multipliers (ADMM) [boyd2011distributed, otero2018alternate]. The augmented Lagrangian for Eq. (9) is: , where is the Lagrange multiplier, is a parameter, and . Note that the term is a constant with respect to and and can be dropped. The updates of , , and are done as [boyd2011distributed, otero2018alternate]:

 U(k+1)i :=argminUi(f(Ui)+(ρ/2)||Ui−V(k)i+J(k)i||2F), (10) V(k+1)i :=argminVi(h(Vi)+(ρ/2)||U(k+1)i−Vi+J(k)i||2F), (11) J(k+1) :=J(k)+U(k+1)−V(k+1). (12)

Considering for a matrix , the gradient of the objective function in Eq. (10) with respect to is where is defined in Eq. (7). We can use the gradient decent method [boyd2004convex] for solving the Eq. (10). Our experiments showed that even one iteration of gradient decent suffices for Eq. (10) because the ADMM itself is iterative. Hence, we can replace this equation with one iteration of gradient decent.

The proximal operator is defined as [parikh2014proximal]:

 proxλ,h(v):=argminu(h(u)+(λ/2)||u−v||22), (13)

where is the proximal parameter and is the function that the proximal algorithm wants to minimize. According to Eq. (13), the Eq. (11) is equivalent to . As is indicator function, its proximal operator is projection [parikh2014proximal]. Therefore, Eq. (11) is equivalent to where denotes projection onto a set. Here, the variable of proximal operator is a matrix and not a vector. According to [parikh2014proximal], if

is a convex and orthogonally invariant function, and it works on the singular values of a matrix variable

, i.e., where the function gives the vector of singular values of , then the proximal operator is:

 (14)

The and are the matrices of left and right singular vectors of , respectively. In our constraint , the function deals with the singular values of . The reason is that we want: , where and are because and are orthogonal matrices. Therefore, we can use Eq. (14) for Eq. (11) where sets the singular values of to one. In summary, Eqs. (10), (11), and (12) can be restated as:

 U(k+1)i :=U(k)i−ηG(U(k)i)−ηρ(U(k)i−V(k)i+J(k)i), (15) V(k+1)i J(k+1) :=J(k)+U(k+1)−V(k+1),

where columns of and are the left and right singular vectors of and is the learning rate. Iteratively solving Eq. (15) until convergence gives us the for for the image blocks indexed by . The columns of are the bases for the ISCA subspace of the -th block. Unlike in PCA, the ISCA bases do not have an order of importance but as in PCA, they are orthogonal capturing different features of image structure. The -th projected block is where its dimensions are image structural components. Note that , whether it is a block in a training image or an out-of-sample image, is centered. It is noteworthy that if we consider only one block in the images, the subscript is dropped from Eq. (15).

### 4.2 Orthonormal Bases for a Set of Images

So far, if we have a set of images, we can find the subspace bases for the -th block in each of them using Eq. (15). Now, we want to find the subspace bases for the -th block in all training images of the dataset. In other words, we want to find the subspace for the best reconstruction of the -th block in all training images. For this goal, we can look at the optimization problem in Eq. (8) or (9

) as an undercomplete auto-encoder neural network

[goodfellow2016deep] with one hidden layer where the input layer, hidden layer, and output layer have , , and neurons, respectively. The and fill the role of applying the first and second weight matrices to the input, respectively. The weights are . Therefore, we will have auto-encoders, each with one hidden layer.

For training the auto-encoder, we introduce the blocks in an image as the input to this network and update the weights based on Eq. (15). Note that we do this update of weights with only ‘one’ iteration of ADMM. Then, we move to the blocks in the next image and update the weights again by an iteration of Eq. (15

). We do this for all images one by one until an epoch is completed where an epoch is defined as introducing the block in all training images of dataset to the network. After termination of an epoch, we start another epoch to tune the weights

again. The epochs are repeated until the convergence. The termination criterion can be average reconstruction error , where is a small number and is the -th block in the -th image. After training the network, we have one -dimensional subspace for every block in all training images where the columns of the weight matrix span the subspace. Note that because of ADMM, the auto-encoders are trained simultaneously and in parallel. Again, the columns of are the bases for the ISCA subspace of the -th block.

## 5 Kernel Image Structural Component Analysis

We can map the block to higher-dimensional feature space hoping to have the data fall close to a simpler-to-analyze manifold in the feature space. Suppose is a function which maps the data to the feature space. In other words, . Let denote the dimensionality of the feature space, i.e., . We usually have . The kernel of the -th block in images and , which are and , is [hofmann2008kernel]. The kernel matrix for the -th block among the images is where . After calculating the kernel matrix, we normalize it [ah2010normalized] as where denotes the -th element of the kernel matrix. Afterwards, the kernel is double-centered as where . The reason for double-centering is that Eq. (2) requires and thus the to be centered (see Eq. (16)). Therefore, in kernel ISCA, we center the kernel rather than centering .

According to representation theory [alperin1993local], the projection matrix can be expressed as a linear combination of the projected data points. Therefore, we have where every column of is the vector of coefficients for expressing a projection direction as a linear combination of projected image blocks.

As we did for ISCA, first we consider learning the subspaces for ‘one’ image, here. Considering for the -th block in the image, the objective function of Eq. (8) in feature space is where . Note that includes the mapping of the -th block in all the images while is the mapping of the -th block in the image we are considering. The constraint of Eq. (8) in the feature space is . Therefore, the Eq. (8) in the feature space is:

 %minimizeΘi∈Rn×p b∑i=1||ϕ(˘xi)−Φ(˘Xi)ΘiΘ⊤iki||S, (16) subject to Θ⊤iKiΘi=I,   ∀i∈{1,…,b}.

Noticing the constraint and using Eq. (2), we have:

 R∋f(Θi):=||ϕ(˘xi)−Φ(˘Xi)ΘiΘ⊤iki||S=ki−k⊤iΘiΘ⊤ikiki+k⊤iΘiΘ⊤iki+c, (17)

where . The gradient of the is:

 Rn×p∋G(Θi):=∂f(Θi)∂Θi=−2(1+f(Θi))ki+k⊤iΘiΘ⊤iki+ckik⊤iΘi. (18)

We can simplify the constraint . As the kernel is positive semi-definite, we can decompose it as:

 Rn×n∋KiSVD=ΨΥΨ⊤=ΨΥ(1/2)Υ(1/2)Ψ⊤=Δ⊤Δ,

where . Therefore, the constraint can be written as: . In Eq. (16), if we embed the constraint in the objective function [boyd2011distributed], we have:

 minimizeΘi,Vi∈Rn×p b∑i=1(f(Θi)+h(ΔVi)), (19) subject to Θ−V=0,

where and and . Taking , we can restate Eq. (19) as: , subject to , where . The ADMM solution to this optimization problem is [boyd2011distributed, otero2018alternate]:

 Θ(k+1)i :=argminΘi(f(Θi)+(ρ/2)||ΔΘi−W(k)i+J(k)i||2F), (20) W(k+1)i :=argminWi(h(Wi)+(ρ/2)||ΔΘ(k+1)i−Wi+J(k)i||2F), (21) J(k+1) :=J(k)+ΔΘ(k+1)−W(k+1). (22)

With the similar explanations which we had for Eq. (15), we have:

 Θ(k+1)i :=Θ(k)i−ηG(Θ(k)i)−ηρΔ⊤(ΔΘ(k)i−W(k)i+J(k)i), (23) W(k+1)i :=Qidiag(proxρ,h(σ(ΔΘ(k+1)i+J(k)i)))Ω⊤i, J(k+1) :=J(k)+ΔΘ(k+1)−W(k+1),

where columns of and are the left and right singular vectors of . Iteratively solving Eq. (23) until convergence gives us the for for the image blocks indexed by . The columns of are the bases for the kernel ISCA subspace of the -th block. The -th projected block is and its dimensions are the kernel image structural components. Note that , whether it is the kernel over a block in a training image or an out-of-sample image, is normalized and centered. Also, note that we had considered the blocks of only one image for Eq. (23). Again, with the auto-encoder approach, we can solve these equations in successive epochs in order to find the subspaces for all the training images.

## 6 Experiments

Training Dataset: We formed a dataset out of the standard Lena image. Six different types of distortions were applied on the original Lena image (see Fig. 1), each of which has images in the dataset with different MSE values. Therefore, the size of the training set is including the original image. For every type of distortion, different levels of MSE, i.e., from to with step , were generated to have images on the equal-MSE or iso-error hypersphere [wang2006modern].

Training: In our experiments for ISCA, the parameters used were and , and for kernel ISCA, we used and . We took ( blocks inspired by [otero2014unconstrained, otero2018alternate]), , and . One of the dimensions of the trained , , and for ISCA are shown in Fig. 2. The dual variable has captured the edges because edges carry much of the structure information. As expected, and are close (Lena can be seen in them by noticing scrupulously). Note that the variables in kernel ISCA are not -dimensional and thus cannot be displayed in image form.

Projections and Comparisons:

In order to evaulate the trained ISCA and kernel ISCA subspaces, we projected the training images onto these subspaces. For projecting an image, each of its blocks is projected onto the subspace of that block. After projecting all the images, we used the 1-Nearest Neighbor (1NN) classifier to recognize the distortion type of every block. The 1NN is useful to evaluate the subspace by closeness of the projected distortions. The distortion type of an image comes from a majority vote among the blocks. The linear, Radial Basis function (RBF), and sigmoid kernels were tested for kernel ISCA. The confusion matrices for distortion recognition are shown in Fig.

3. Mostly kernel ISCA performed better than ISCA because it works in feature space; although, ISCA performed better for some distortions like contrast stretching and blurring. Moreover, we compared with PCA and kernel PCA. PCA showed weakness in contrast stretching. RBF and sigmoid kernels in kernel PCA do not perform well for JPEG distortion and contrast stretching, respectively.

Out-of-sample Projections: For out-of-sample projection, we created test images with having different distortions and some having a combination of different distortions (see Fig. 4). We did the same 1NN classification for these images. Table 1 reports the top two votes of blocks for every image with the percentage of blocks voting for those distortions. ISCA did not recognize luminance enhancement well enough because, for Eq. (2), the block is centered while in kernel ISCA, the block is centered in feature space. Overall, both ISCA and kernel ISCA performed very compelling even in recognizing the combination of distortions.

Reconstruction: The images can be reconstructed after the projection onto the ISCA subspace. For reconstruction, every block is reconstructed as where the mean of block should be added to the reconstruction. Similar to kernel PCA, reconstruction cannot be done in kernel ISCA because . Figure 5 shows reconstruction of some of training and out-of-sample images. As expected, the reconstructed images, for both training and out-of-sample images, are very similar to the original images.

## 7 Conclusion and Future Direction

This paper introduces the concept of an image structure subspace which captures the structure of an image and discriminates the distortion types. We hope this will open a broad new field for research in this area and build a greatly needed bridge between the worlds of image quality assessment and manifold learning.

For image structure subspace learning, ISCA and kernel ISCA were proposed, taking inspiration from PCA. As future work, we can consider designing deeper auto-encoder [goodfellow2016deep]

with non-linear activation functions for image structure subspace learning.

## 8 Supplementary Material

This section is the supplementary material for the paper “Principal Component Analysis Using Structural Similarity Index for Images”. In this paper, the derivation of the mathematical expressions, which were not completely detailed in the main paper, are explained. We explain the derivation of Eqs. (6), (7), (15), (17), (18), and (23).

### 8.1 Derivation of Eq. (6)

In the following, we mention the derivation of Eq. (6):

 f(U):=||˘x−UU⊤˘x||S(a)=||˘x−UU⊤˘x||22||˘x||22+||UU⊤˘x||22+c,

where is because of Eq. (2). The numerator of is simplified as:

 ||˘x−UU⊤˘x||22 =(˘x−UU⊤˘x)⊤(˘x−UU⊤˘x)=(˘x⊤−˘x⊤UU⊤)(˘x−UU⊤˘x) =˘x⊤˘x−˘x⊤UU⊤˘x−˘x⊤UU⊤˘x+˘x⊤UU⊤UIU⊤˘x (a)=˘x⊤˘x−2˘x⊤UU⊤˘x+˘x⊤UU⊤˘x =˘x⊤˘x−˘x⊤UU⊤˘x=˘x⊤(I−UU⊤)˘x,

where is because of the constraint in Eq. (5).

The first term in denominator of is simplified as:

 ||˘x||22=˘x⊤˘x,

and the second term in denominator of is simplified as:

 ||UU⊤˘x||22 =(UU⊤˘x)⊤(UU⊤˘x)=(˘x⊤UU⊤)(UU⊤˘x) =˘x⊤UU⊤UIU⊤˘x(a)=˘x⊤UU⊤˘x,

where is because of the constraint in Eq. (5). Thus, the denominator of is:

 ||˘x||22+||UU⊤˘x||22+c=˘x⊤˘x+˘x⊤UU⊤˘x+c=˘x⊤(I+UU⊤)˘x+c.

Therefore, the Eq. (6) becomes:

 R∋f(U)=