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.

The test image. (Please download test image at left. Only two blocks are identical.)

Input data.

PCA results.

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);

%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));

  %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;
    column = column + 1;
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);
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
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.

  • $\begingroup$ I'm not sure about your screenshots. Could you add or link the full Octave code? $\endgroup$
    – leonbloy
    Commented Apr 20, 2012 at 0:49

1 Answer 1


It's to be expected that "copied" blocks are almost equal (and more so after the PCA manipulation), so in the lexicographical sort (warning: it's understood that this lexicographic order orders first the most principal component, and so on) "copied" blocks should appear adjacent or near (the reverse is not true: adjacent lexicographicly sorted elements are not necessarily copied, nor even similar)

Here I made up a very simple example myself, in Octave, with a unidimensional signal (y) of size N=200, which has a portion of it copied (here, from 20-50 to 150-180) and a little noise added. I take a small block size (b=3). I convert to PC, sort the rows in lexicographical order (I append first the original block position in an extra column), and compute the distance between adjacent rows (notice that I'm simplifiying a lot here: I'm not discarding components, nor quantizing them; and I'm considering only adjacent rows, not a neighborhood band). I then look at the histogram of those distance, and the original offset is cleary visible.

y = filter([1],[1,-0.8,0.1],rand(1,N)-0.5); % my signal, rather arbitrary
y(20+delay:50+delay) = y(20:50);  % a portion is copied
y += (rand(1,N)-0.5)*0.1; % noise added
yy=[y(1:N-2);y(2:N-1);y(3:N)];  % octave does not have  corrmtx (this is not general in b!)
[PC, Z, W, TSQ] = princomp (yy'); % PCA
Z(:,b+1)=[1:N-2]'; % append original block position, in extra row
Z1=sortrows(Z);  % sort rows lexicographycally
Z2=abs(Z1(1:N-3,b+1)-Z1(2:N-2,b+1));  % compute temporal distances between adjacent rows
histo(Z2); % histogram: should show a peak at delay
  • $\begingroup$ Thanks! I noticed you do princomp(yy'). Taking the transpose of the B matrix in my code fixed everything. $\endgroup$
    – JCD
    Commented Apr 24, 2012 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.