# Image Compression Using the Singular Value Decomposition (SVD) with MATLAB

I'm new to both linear algebra and MATLAB.
I need help with this code with the objective of compressing image using the Singular Value Decomposition (SVD).
The code is currently working, if I run with e.g sigma_threshold = 10, I get a blurry picture, if I set it to 100 it's even more blurry etc.

I have searched all over the web and found a lot of similiar examples that I have used to improve this code with, but still I know I can do it a lot more elegant.

One problem is that I dont fully understand the math behind what I'm doing, and thus I don't know really know what to do, particularly with the matrices U and V in [U S V] = svd(A).

Wondered if someone could give me a hint how to improve the sections under the comments "%Find the first index smaller than sigma_threshold" and "%Compose the red, green and blue channel again until the sigma value"

function svd_exercise(imagename,sigma_threshold)

image = imread(imagename);
image = im2double(image);

%Decompose the image in rgb values
R = image(:, :, 1);
G = image(:, :, 2);
B = image(:, :, 3);

%Compute the Single Value decomposition for each channel
[U_R, S_R, V_R] = svd(R);
[U_G, S_G, V_G] = svd(G);
[U_B, S_B, V_B] = svd(B);

%Find the first index smaller than sigma_threshold
[red_j, ~] = find(S_R > 0 & S_R < sigma_threshold, 1);
[green_j, ~] = find(S_G > 0 & S_G < sigma_threshold, 1);
[blue_j, ~] = find(S_B > 0 & S_B < sigma_threshold, 1);

%Compose the red, green and blue channel again until the sigma value
R_K = U_R(:, 1:red_j) * S_R(1:red_j, 1:red_j) * V_R(:, 1:red_j)';
G_K = U_G(:, 1:green_j) * S_G(1:green_j, 1:green_j) * V_G(:, 1:green_j)';
B_K = U_B(:, 1:blue_j) * S_B(1:blue_j, 1:blue_j) * V_B(:, 1:blue_j)';

AK = zeros(size(image));
AK(:,:,1) = R_K;
AK(:,:,2) = G_K;
AK(:,:,3) = B_K;


## 2 Answers

If you look at the structure of your matrices of singular values S_R, S_G and S_B, you'll see that the non-zero entries lie on the diagonal. You can check only these entries by using diag, removing the need for your S_R>0 condition.

This also means that for find you don't need to supply the second output argument.

To get the first non-zero entries less than sigma_threshold you can hence use

red_j = find(diag(S_R) < sigma_threshold, 1);
green_j = find(diag(S_G) < sigma_threshold, 1);
blue_j = find(diag(S_B) < sigma_threshold, 1);


Bear in mind that this is the first value which should be cut off, not the last value to include, hence you might want to do

red_j = find(diag(S_R) < sigma_threshold, 1) - 1;


or

red_j = find(diag(S_R) >= sigma_threshold, 1, 'last');


This is a good first try.
But there are some points to improve:

1. This is not the right way to compress with SVD.
When you compress with SVD you should compress blocks of the image.
See my attached implementation.
2. Important step before doing SVD is to remove the DC Level (And remember bringing it back).
3. Work in a Loop on the image channels to support Gray Scale images.
4. The names of variables should be more meaningful. It is me, but I don't like using _ in variable name. Moreover, I like differentiating Vectors and Matrices in MATLAB.

This is the function:

function [ mO ] = CompressImageSvd( mI, energyThr, blockRadius )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = CompressImageSvd( mI, energyThr, blockRadius )
%   Compresses the image using SVD.
% Input:
%   - mI            -   Input Image.
%                       Structure: Image Matrix (Signle Channel or RGB).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
%   - energyThr     -   Energy Threshold.
%                       Sets the threshold for Singular Value kept energy.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
%   - blockRadius   -   Block Radius.
%                       Sets the block radius for the compression.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, ...}.
% Output:
%   - mO            -   Output Image.
%                       Structure: Image Matrix (Signle Channel or RGB).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
% References
%   1.  SVD Wikipedia - https://en.wikipedia.org/wiki/Singular_value_decomposition.
% Remarks:
%   1.  a
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     01/09/2017  Royi
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

numRows = size(mI, 1);
numCols = size(mI, 2);
numChan = size(mI, 3); %<! Num Channels

vImageDim = [numRows, numCols];

blockLength = (2 * blockRadius) + 1;
vBlockDim   = [blockLength, blockLength];

mO = zeros([numRows, numCols, numChan]);

for ii = 1:numChan

mII     = mI(:, :, ii);
dcLevel = mean(mII(:)); %<! Extracting DC Level
mII     = mII - dcLevel;

% Decomposing the image into blocks. Each block becomes a vector in the
% Columns Images.
mColImage   = im2col(mII, vBlockDim, 'distinct');

% The SVD Step
[mU, mS, mV] = svd(mColImage);

vSingularValues = diag(mS);

vSingularValueEnergy = cumsum(vSingularValues) / sum(vSingularValues);
lastIdx = find(vSingularValueEnergy >= energyThr, 1, 'first');

vSingularValues(lastIdx + 1:end) = 0;
% mS isn't necessarily square matrix. Hence only work on its main
% diagonal.
mS(1:length(vSingularValues), 1:length(vSingularValues)) = diag(vSingularValues);

% Reconstruction of the image using "Less Energy".
mColImage = mU * mS * mV.';

% Restorig the original structure and the DC Level
mO(:, :, ii) = col2im(mColImage, vBlockDim, vImageDim, 'distinct') + dcLevel;

end

end


This is the result:

The full code can be found in my StackExchange Code Review Q157459 GitHub Repository.

• Nice, but what's the purpose of defining TRUE and FALSE and so on? I think we've had Boolean constants since MATLAB 6 or so? Jan 12, 2018 at 3:16
• Well, Just my style of coding. By the way true() / false() are functions in MATLAB and not constants.
– Royi
Feb 10, 2018 at 21:10