function fig4_10

% low rank approximations using SVD

%k=[5 10 20 40 80 160 320];
k=[10 25 50 100];

% read in image and display it
%x=imread('flower.jpg');
x=imread('butter.jpeg');
I = double(x);
% I =  im2double(x);  % requires image processing toolbox
% I =  im2uint8(x);  % requires image processing toolbox
size_of_image=size(I)

% get(gcf)
set(gcf,'Position', [4 548 560 493])
clf
% x=imread('ff.jpeg');
image(x)
set(gca,'DataAspectRatio',[1 1 1])
axis image
axis off
say=['Original Image'];
text(5,-50,say,'Color','Black','FontSize',20,'FontWeight','bold')
%pause

% separate color planes
m=length(I(:,1,1));
n=length(I(1,:,1));
IRed = I(:,:,1);
IGreen = I(:,:,2);
IBlue = I(:,:,3);

% check out some of the values in red plane
%IRed(1:5,1:5)

% compose A and then find its SVD
A=[ IRed ; IGreen; IBlue ];
size_A=size(A)
[U,S,V] = svd(A);
ss=size(S);
rs=min(ss(1),ss(2));
for i=1:rs
    sigma(i)=S(i,i);
end

% compressed versions
figure(2)
for ik=1:length(k)

    % zero out singular values bigger than k
    Sa=zeros(ss(1),ss(2));
    for ir=1:k(ik)
        Sa(ir,ir)=sigma(ir);
    end
    Aa=U*Sa*V';

    clf
    % get(gcf)
    set(gcf,'Position', [8 8 560 493])

    % reassemble color planes
    PP(:,:,1)=Aa(1:m,:);
    PP(:,:,2)=Aa(m+1:2*m,:);
    PP(:,:,3)=Aa(2*m+1:3*m,:);

    % convert floats to 8 bit integers
    AUint8=uint8(PP);
    image(AUint8);
     
    set(gca,'DataAspectRatio',[1 1 1])
    axis image
    axis off
    err=sigma(k(ik)+1)/sigma(1);
    errz=round(err,3);
    say=['k = ',num2str(k(ik)),'          Rel Error = ',num2str(errz)];
    text(5,-70,say,'Color','Black','FontSize',20,'FontWeight','bold')
    pause
    % if ik==1
    %     exportgraphics(gcf,'/Users/mark/Desktop/flower1.eps')
    % elseif ik==2
    %     exportgraphics(gcf,'/Users/mark/Desktop/flower2.eps')
    % elseif ik==3
    %     exportgraphics(gcf,'/Users/mark/Desktop/flower3.eps')
    % elseif ik==4
    %     exportgraphics(gcf,'/Users/mark/Desktop/flower4.eps')
    % end
end

% compute and then plot relative error
figure(1)
for i=1:ss(2)-1
    id(i)=i;
    err(i)=sigma(i+1)/sigma(1);
end
clf
% get(gcf)
set(gcf,'Position', [4 1115 658 230])
semilogy(id,err,'r','Linewidth',2)
hold on
set(gca,'ytick',[1e-5 1e-4 1e-3 1e-2 1e-1 1])
xlabel('k-axis')
ylabel(' Relative Error')
grid on
set(gca,'YMinorGrid','off')
box on
set(gca,'FontSize',16,'FontWeight','bold')



