Execute MATLAB/Octave Online

%L = imread('images/lena3.jpg');
L = imread('images/Toronto2_F.png');
L=rgb2gray(L);
figure(1);
imshow(L);

b=256;  % block size used to look for duplicates
b1=sqrt(b); 
Nn = 5; % how many close rows to check
Nd = 15; %Threashold
Nc=26; % Truncate PCA at this length
Nt=26;

% calculate the total size of the image
n = numel(L);

b2 = sqrt(n)-sqrt(b)+1;
% calculate Nb
Nb= power((sqrt(n)-sqrt(b)+1),2);

% the matix of Nc plus the position
M=zeros(Nb,  Nc);
% temp index array
Mi = zeros(Nb, 2);
i=1;
disp('Starting PCA');
for r = 1:b2
   for c = 1:b2
       % Extract each block
      B = L(r:r+b1-1,c:c+b1-1);

      [pc, latent, explained] = pcacov(cov(double(B)));
      %[pc, latent, explained] = princomp(double(B), 'NumComponents', Nc);

      Z = pc(1:Nc);
      Mi(i,:) = [r c];
      M(i,:) = Z;      
      i = i+1;
   end
end

disp('Sorting M -> S');
%Sort M array in lexicographic order -> S
[S, index] = sortrows(M);
P= zeros(1,3);
disp('Finding Duplicates');
for i = 1:Nb
    iv = index(i);
    xi=mod(iv,b2) + 1;
    yi=ceil(iv/b2);
    j = i+1;    
    while j < Nb && abs(i - j) < Nn
        jv=index(j);
        xj=mod(jv,b2) + 1;
        yj=ceil(jv/b2);

        z=sqrt(power(xi-xj,2) + power(yi-yj,2));
        % only process those whose size is above Nd 

        if z > Nd

            idx = find(P(:,1)== xi & P(:,2)==yi, 1, 'last');
            if isempty(idx)==1
               P = [P; [xi, yi, 1]];                 
            else
               P(idx,3) = P(idx,3) + 1;
            end

            idx = find(P(:,1)== xi & P(:,2)==yi, 1, 'last');
            if isempty(idx)==1
               P = [P; [xj, yj, 1]];
            else
               P(idx,3) = P(idx,3) + 1;                  
            end            
        end                        
       j = j + 1;
    end
end

disp('Sorting findings');
rows = size(P,1);
% sort descending order
P = sortrows(P, -3);
% Mark the found blocks
disp('Creating Image');
idx = 1;
% Create a black image
RI = zeros(sqrt(n), sqrt(n));
while idx < rows && P(idx,3) > 5
    x = P(idx,1);
    y = P(idx,2);
    RI(x,y) = 1;
    idx = idx + 1;
end

figure(2);
imshow(RI);