Image Compression Using the SVD in MATLAB

Question

I have to compress an image using SVD. I did a few things, but all of them produced a black & white image. However, the output must be RGB. How can I compress an image without changing its original format?

Answer 1

The typical thing to do is the low-rank approximation on separate channels. Assume that $C$ is a channel of the RGB image $I$:

rank = 10;
[U,S,V] = svd(C);
L = U(:,1:rank) * S(1:rank, 1:rank) * V(:, 1:rank)';

Now, L should be the compressed image. If you do this operation and compose the channels back, you should get a compressed RGB image.

However, such a method in my opinion is only good for mathematical understanding. It is not very practical due to the fact that it depends on image dimensions too much and cannot generate high SNR results. Or better stated, as you decrease rank, you drastically reduce the quality of the reconstruction. For this reason, one has to retain around ~50-100 components. Better methods doing this blockwise, or formulating the compression as an energy functional exist. Image denoising is a huge area of its own.

You might want to check some results with varying rank in this presentation.

---

Also, a similar question has already been asked; please follow the Stack Exchange discussion here.

Answer 2

Let's do that in a hands on approach.

We'll follow this steps:

1. Apply the SVD to the Image Blocks.
   When you compress with SVD you should compress blocks of the image.
   The reason is simple, images are objects with spatial correlation. This is exactly what we take advantage of for compression.
2. Pre Prcoessing - Remove the DC Level
   Important step before doing SVD is to remove the DC Level (And remember bringing it back).
3. RGB Images
   Work in a Loop on the image channels to support Gray Scale images. Though you can work with the SVD extension to apply it on tensors (See Apply Principal Component Analysis (PCA) for RGB Images). Yet the color correlation is not worth the trouble usually. Pay attention that usually it is good to convert to YCbCr or other channel which separates color from luminosity as we can compress the color more.

Let's have a look on the function:

function [ mO ] = CompressImageSvd( mI, energyThr, blockRadius )
% Release Notes:
%   -   1.0.000     01/09/2017  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

numRows = size(mI, 1);
numCols = size(mI, 2);
numChan = size(mI, 3); %<! Number of 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.';
    
    % Restoring the original structure and the DC Level
    mO(:, :, ii) = col2im(mColImage, vBlockDim, vImageDim, 'distinct') + dcLevel;
    
end


end

1. We iterate over the image's channels.
2. Per channel, we calculate the mean (dcLevel). We better do that per block of data. But for simplicity we did it for the whole channel.
3. We decompose the image into distinct blocks.
4. We apply the decomposition per block.
5. We keep only part of the decomposition according to the energy threshold.

The full code is available on my StackExchange Signal Processing Q18673 GitHub Repository (Look at the SignalProcessing\Q18673 folder).