% CUR first try
clear
clc

% read the image into memory and crop it to 250x250
I1 = imread('peppers.bmp');
I2 = double(I1(1:250,1:250));

% subplot(1,2,1), imshow(I1,[0 255]), title('Original Peppers image')
% subplot(1,2,2), imshow(I2,[0 255]), title('Cropped image to 250x250 pixels')

% number of selected columns and rows
X = 100;

% take random columns and rows, later this will be done through leverage
rndCOLS = randperm(250); rndCOLS = rndCOLS(1:X);
rndROWS = randperm(250); rndROWS = rndROWS(1:X);

% determine C (selected columns) and R (selected rows)
% show random rows and columns
C = zeros(250,X);
R = zeros(X,250);
selectedCOLS = zeros(250);
selectedROWS = zeros(250);
for i = 1 : X
    C(:,i) = I2(:,rndCOLS(i));
    R(i,:) = I2(rndROWS(i),:);
    selectedCOLS(:,rndCOLS(i)) = I2(:,rndCOLS(i));
    subplot(1,3,1), imshow(selectedCOLS,[0 255]), title('Selected columns'), xlabel(i), drawnow
    selectedROWS(rndROWS(i),:) = I2(rndROWS(i),:);
    subplot(1,3,2), imshow(selectedROWS,[0 255]), title('Selected rows'), xlabel(i), drawnow
end

% construct matrix U as an intersection of rows and columns
W = zeros(X,X);
for i = 1 : X
    for j = 1 : X
        W(i,j) = [I2(rndROWS(i),rndCOLS(j))];
    end
end
U=pinv(W);

I2_cur = C * W * R;

figure,
subplot(1,2,1), imshow(I2,[0 255]), title('Original image')
subplot(1,2,2), imshow(I2_cur,[]), title('CUR approximation')