I'm trying to understand an algorithm that detects similar regions of an image using PCA. The algorithm essentially divides the image into overlapping square blocks and then does PCA with each block as a dimension. Next it does a lexicographical sort of the principal components to arrange similar blocks close to one another. Finally similar regions are selected based on distance in the list and a few other measurements.
The PCA part of the algorithm is described as follows:
An image is tiled with overlapping blocks of $b$ pixles ($\sqrt{b} \times \sqrt{b}$ pixels in size), each of which is assumed to be considerably smaller than the size of the duplicated regions to be detected. Let $ \vec x_i,i=1,...,N_b,$ denote these blocks in vectorized form, where $N_b=\left(\sqrt{N}-\sqrt{b}+1 \right)^2$. We now consider an alternate representation of these blocks based on a principal component analysis (PCA)[15]. Assume that the blocks $\vec x_i$ are zero-mean, and the compute the covarience matrix as: \begin{equation*} \tag{1.1} C = \sum_{i=1}^{N_b} \vec x_i \vec x_i^T \end{equation*} The orthonormal eigenvectors,$\vec e_j$,of the matrix $C$, with corresponding eigenvalues,$\lambda_j$, satisfying: \begin{equation*} \tag{1.2} C \vec e_j = \lambda_j \vec e_j, \end{equation*} define the principal components, where $j=1,\dots,b$ and $ \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_b$. The eigenvectors, $\vec e_j$, form a new linear basis for each image block,$\vec x_i$: \begin{equation*} \tag{1.3} \vec x_i = \sum_{j=1}^b a_j \vec e_j, \end{equation*} where $a_j = \vec x_i^T \vec e_j$, and $\vec a_i = \left(a_i \dots a_b \right)$ is the new representation for each image block.
The dimensionality of this representation can be reduced by simply truncating the sum in Equation (1.3) to the first $N_t$ terms. Note that the projection onto the first $N_t$ eigenvectors of the PCA basis gives the best $N_t$-dimensional approximation in the least squares sense, if the distribution of the vectors, $\vec x_i$, is a multi-dimensional Gaussian [15]. This reduced dimension representation, therefore, provides a convenient space in which to identify similar blocks in the presence of corrupting noise, as truncation of the basis will remove minor intensity variations.
I've completed translating the algorithm into MATLAB/Octave but it does not appear to be working correctly. Below is my test image with a screen shot of the data in vectorized form and the PCA results. Here $N_b = 225$ and $ b=4$. In the PCA table all the remaining rows are either 0 or -1.
(Please download test image at left. Only two blocks are identical.)
What I don't understand is how a lexicographical sort on the PCA data will show similarity. My expectation was that row 1 and 4 in the PCA output would be identical and that PCA would allow me to remove higher columns to increase matches. Similar to the way discrete cosine transform moves 'important' information to a smaller number of values of the data set. I've verified that my PCA algorithm mostly matches (signs are different) vs MATLAB's princomp() command so I must not understand what PCA is telling me. My question is, how would a sort on PCA results show self-similarity?
%This does principal component analysis
clear all;
bSide = 2; % Size of side of one block
%Read B&W image
chanData = imread('./randS.gif');
%Setup the input values
[imageWidth,imageHeight] = size(chanData);
tic();
%Calculate the derived constants
N = imageHeight * imageWidth;
b = bSide^2;
Nb = (sqrt(N) - sqrt(b) + 1)^2;
printf("ET:%d - ",toc());
printf("Finished deriving constants: N:%d b:%d Nb:%d\n",N,b,Nb);
%A set of observations of M variables
X = zeros(b,Nb); %Block representation of input image
u = zeros(b,Nb);%Temp for zero-mean calculation
coords = zeros(Nb,2); %Coordinate of the block at this location
%Create our matrix holding the dimension data
%X = zeros([Nb,b]); %b columns by Nb rows
%Populate the dimension matrix
iWidth = imageWidth - bSide + 1;
row = 1;
column = 1;
for i = 1:Nb
%Setup local variables
Xtemp = zeros(1,b);
%Each block becomes a column of bSide rows
for j = 1:bSide
lowerX = (j - 1)*bSide + 1;
Xtemp(lowerX:lowerX+(bSide-1)) = chanData(row+(j-1),column:column+(bSide-1));
end
%Add this vector as a column to the array
X(:,i) = Xtemp;
%Now create a mean vector
u(i) = mean(X(:,i));
%Save the coordinates of this block in a two column matrix
coords(i,1) = row;
coords(i,2) = column;
%Update our source matrix indicies
if column == iWidth
row = row + 1;
column = 1;
else
column = column + 1;
end
end
clear Xtemp row column j iWidth lowerX;
printf("ET:%d - ",toc());
printf("Finished matrix construction\n");
%Get the mean for this vector and subtract it from the data
B = zeros(b,Nb);
for i = 1:Nb
B(:,i) = X(:,i) - u(i);
end
printf("ET:%d - ",toc());
printf("Allocated mean-zero matrix\n");
%clear u X;
%Find the covarience matrix
C = cov(B);
printf("ET:%d - ",toc());
printf("Found covarience matrix\n");
%clear B;
%Find the eigenvalues and eigenvectors
%Evec: Nb x Nb matrix of eigenvectors, each column is a vector
%Eval: Nb x Nb matrix where only the main diagonal contains eiganvalues
[Evec,Eval] = eig(C);
Eval = diag(Eval); %Convert eigenvalue matrix to vector
%clear C;
printf("ET:%d - ",toc());
printf("Found eigenvalues\n");
%Sort by eigenvalue. We do this by appending the eigenvalue vector as the
%first row. MATLAB won't sort by column for some reason, so we must
%transpose.
sortedEvec = [transpose(Eval); Evec];
sortedEvec = transpose(sortedEvec);
sortedEvec = sortrows(sortedEvec,-1);
sortedEvec = transpose(sortedEvec);
printf("ET:%d - ",toc());
printf("Sorted eigenvectors\n");
%Remove the eigenvalues row
sortedEvec(1,:) = [];
%Now do the final multiplication of the two transposed matrices
PCA = transpose(sortedEvec) * transpose(B);
printf("ET:%d - ",toc());
printf("Finished PCA\n");
%PCAt = PCA'; %PCAt matches princomp(B) except signs are reversed for first column
%Now integer quantize the array
PCA = floor(PCA);
printf("ET:%d - ",toc());
printf("Quantized PCA\n");
This is code I used to get the PCA matrix. princomp isn't defined for Octave but I was able to access MATLAB and run it there. My results match except the signs are inverted and the matrix from princomp is the transpose of my results.
