% from ChatGPT, this works pretty well

clear
clc

A = double(imread('peppers.bmp'));
[m,n] = size(A);

k = 150;  % <-- MUCH BETTER for 256×256

% Rank-k SVD
% [U,S,V] = svds(A, k);

% Leverage scores
% row_lev = sum(U.^2,2);
% col_lev = sum(U.^2,2);
% 
% row_prob = row_lev / sum(row_lev);
% col_prob = col_lev / sum(col_lev);

% row_idx = randsample(m, k, true, row_prob);
% col_idx = randsample(n, k, true, col_prob);

% row_idx = randsample(m, k, true, row_prob);
% col_idx = randsample(n, k, true, col_prob);
[~,~,col_idx] = qr(A,'vector');
col_idx = col_idx(1:k);
C = A(:, col_idx);

[~,~,row_idx] = qr(A','vector');
row_idx = row_idx(1:k);
R = A(row_idx, :);

% C = A(:, col_idx);
% R = A(row_idx, :);
W = A(row_idx, col_idx);

U_cur = pinv(W);
A_cur = C * U_cur * R;

imshow(uint8(A_cur));

myPSNR = psnr(A, A_cur, 255);
mySSIM = ssim(A, A_cur);



W = A(row_idx, col_idx);
U = pinv(W);
A_cur = C * U * R;