042814Aaron

From Computational Statistics Course Wiki
Jump to navigation Jump to search
My code for the exercise
%--------------------------------------------------------------------------
%  Original Author: Isaac J. Lee
%  E-mail: ijlee2@ices.utexas.edu
%  Adapted by: John Hawkins
%  
%  This routine demonstrates the use of SVD in image compression.
%  
%  You will need to add/edit/delete some code to run this program. These
%  parts are indicated by a comment beginning with the word "PROBLEM."
%  
%  To run this program, type in:
%      svd_exercise()
%--------------------------------------------------------------------------
function svd_exercise()
    clear all
    clc; clf;
    
    % A is a h x w x 3 matrix and stores the intensities in the red, green,
    % and blue channel
    % cacti.jpg needs to be in the same directory as this file
    %A = imread('cacti.jpg');
    A = imread('Ryb_twist.jpg');
    %A = imread('Ryb.jpg');
    h = size(A, 1); w = size(A, 2);
    % Extract the red, green, and blue channel matrices
    AR = double(A(:, :, 1));
    AG = double(A(:, :, 2));
    AB = double(A(:, :, 3));
    
    % ------------------------------------------------------------
    %  Perform SVD for each color channel matrix.
    % ------------------------------------------------------------
    [UR, SR, VR] = svd(AR);
    % PROBLEM #1: Find the SVD for the green and blue channels. (Please keep
    % the variable names consistent by using UG, SG, VG, UB, SB, and VB.)
    [UG, SG, VG] = svd(AG);
    [UB, SB, VB] = svd(AB);
    figure(1)
    semilogy(diag(SB),'b-')
    hold on
    semilogy(diag(SG),'g-')
    hold on
    semilogy(diag(SR),'r-')
    B=rand(400,600);
    [U Sig Vt]=svd(B);
    semilogy(diag(Sig),'k')
    hold off
    legend('blue','green','red','rand')
    % END PROBLEM

    % ------------------------------------------------------------
    %  Find the rank of AR, AG, and AB by counting the number of
    %  positive singular values. We will use the variables r_AR,
    %  r_AG, and r_AB later to compute the relative error.
    % ------------------------------------------------------------
    r_AR = size(find(SR > 0), 1)
    % PROBLEM #2: Repeat for the green and blue channels. (Please use r_AG
    % and r_AB.)
    r_AG = size(find(SG > 0), 1)
    r_AB = size(find(SB > 0), 1)
    % END PROBLEM
    
    % ------------------------------------------------------------
    % PROBLEM #3: 
    %  We wish to find the number of singular values which contain
    %  most of the information. Plot a semilogy plot of the three 
    %  sets of singular values. Where does the added information 
    %  seem to level off?
    % ------------------------------------------------------------

    %figure

    % END PROBLEM
    
    % Display the original cacti image
    figure
    subplot(2, 3, 1);
    image(A);
    axis off; axis image;
    title('Original', 'FontSize', 13);
    set(gca, 'FontSize', 13);

    % PROBLEM #4: Choose 5 k values which you think will demonstrate the tapering off in quality of
    % the image and store them in the vector 'ranks' below, and uncomment the line.
    
    % ranks = []';
    % END PROBLEM 
ranks=[3 4 5 6 7 ]';
    for i = 1:size(ranks, 1)
        k = ranks(i);
        % ------------------------------------------------------------
        %  Compute the best rank-k approximation for each channel.
        % ------------------------------------------------------------
        AR_k = UR(:, 1:k) * SR(1:k, 1:k) * VR(:, 1:k)';
        % PROBLEM #5: Compute the best rank-k approximation for the 
        % green and blue channels.
        AB_k = UB(:, 1:k) * SB(1:k, 1:k) * VB(:, 1:k)';
        AG_k = UG(:, 1:k) * SG(1:k, 1:k) * VG(:, 1:k)';
        % END PROBLEM
        % Combine the three channels to produce the colored image
        A_k = cat(3, AR_k, AG_k, AB_k)/255;
        
        % ------------------------------------------------------------
        %  Evaluate the relative error ||A - A_{k}||_{F} / ||A||_{F}.
        %  We do this using the frobenius norm, which can be 
        %  calculated as the sum of the squares of the singular values
        % ------------------------------------------------------------
        
        errAR_k=sqrt(sum(sum(SR(k+1:end,k+1:end).^2))/sum(sum(SR.^2)));
        errAG_k=sqrt(sum(sum(SG(k+1:end,k+1:end).^2))/sum(sum(SG.^2)));
        errAB_k=sqrt(sum(sum(SB(k+1:end,k+1:end).^2))/sum(sum(SB.^2)));
        
        fprintf('k = %d relative errors in red, green, and blue channels: %1.4g, %1.4g, %1.4g.\n', k, errAR_k, errAG_k, errAB_k);
        
        % Display the rank-k approximation images
        subplot(2, 3, i + 1);
        fileName = ['rank', num2str(k), '.jpg'];
        imwrite(A_k, fileName, 'jpg');
        image(imread(fileName));
        axis off; axis image;    
        delete(fileName);
        title(['k = ', num2str(k)], 'FontSize', 13);
        set(gca, 'FontSize', 13);
    end
    maxk=(h*w)/(h+1+w)
end



r_AR =

   560


r_AG =

   560


r_AB =

   560

k = 3 relative errors in red, green, and blue channels: 0.2813, 0.2777, 0.2941.
k = 4 relative errors in red, green, and blue channels: 0.2692, 0.2648, 0.2809.
k = 5 relative errors in red, green, and blue channels: 0.2574, 0.2531, 0.2682.
k = 6 relative errors in red, green, and blue channels: 0.2454, 0.2411, 0.2556.
k = 7 relative errors in red, green, and blue channels: 0.2338, 0.2292, 0.2431.

maxk =

   280

AaronPCA.jpg AaronImageComp.jpg