function MondrianSVD

global min_r min_g min_b max_r max_g max_b

im=imread('mondrian.jpg');

im_r = double(im(:,:,1));
im_g = double(im(:,:,2));
im_b = double(im(:,:,3));

[u_r d_r v_r] = svd(im_r);
[u_g d_g v_g] = svd(im_g);
[u_b d_b v_b] = svd(im_b);

min_r = min(u_r(:));max_r = max(u_r(:));
min_g = min(u_g(:));max_g = max(u_g(:));
min_b = min(u_b(:));max_b = max(u_b(:));

for i = 1 : 80
    nBasis = i;
    im_r_new = u_r(:,1:nBasis) * d_r(1:nBasis,1:nBasis) * v_r(:,1:nBasis)';
    im_g_new = u_g(:,1:nBasis) * d_g(1:nBasis,1:nBasis) * v_g(:,1:nBasis)';
    im_b_new = u_b(:,1:nBasis) * d_b(1:nBasis,1:nBasis) * v_b(:,1:nBasis)';
    
    im_new(:,:,1) = im_r_new;
    im_new(:,:,2) = im_g_new;
    im_new(:,:,3) = im_b_new;

    h = figure(5)
    clf;
    subplot(1,5,1)
    imshow(im);
    subplot(1,5,2)
    im_new(:,:,1) = im_r_new;
    im_new(:,:,2) = im_g_new;
    im_new(:,:,3) = im_b_new;
    imshow(uint8(im_new));
    xlabel(['Number of basis: ' num2str(i)],'FontSize',20,'FontWeight','bold');
    subplot(1,5,3);
    u_new = VisualizeSingularVector(u_r(:,1:nBasis), u_g(:,1:nBasis), u_b(:,1:nBasis), [size(u_r,1), 30*nBasis]);
    imshow(u_new);
    subplot(1,5,4);
    d = d_r(1:nBasis,1:nBasis);
    d = d/d(1,1)'
    d(d==0) = -0.01;
    imagesc(d);
    axis equal
    axis off
    colormap(jet);
    subplot(1,5,5);
    v_new = VisualizeSingularVector(v_r(:,1:nBasis)', v_g(:,1:nBasis)', v_b(:,1:nBasis)', [30*nBasis, size(im,1)]);
    imshow(v_new);
    print(h, sprintf('svd_%03d.png', i), '-dpng');
end

function u_new = VisualizeSingularVector(u_r, u_g, u_b, scale)
global min_r min_g min_b max_r max_g max_b

u_new(:,:,1) = 16*255*((u_r-min_r)/(max_r-min_r)).^4;
u_new(:,:,2) = 16*255*((u_g-min_g)/(max_g-min_g)).^4;
u_new(:,:,3) = 16*255*((u_b-min_b)/(max_b-min_b)).^4;
u_new = uint8(u_new);
u_new = imresize(u_new, scale, 'nearest'); 

