15apr: Subband image coding & RPCA

Activity 1: Subband image coding

We apply a one-dimensional analysis filter (magnitude distortionless FIR filter) to process each row and column of the image separately. Then the intermediate subband signal forms another image which we call it subband image. After one time partition, we get four subband images which correspond to: lowpass signal of each row and lowpass signal of each column, lowpass signal of each row and highpass signal of each column, highpass signal of each row and lowpass signal of each column, and highpass signal of each row and highpass signal of each column. In subsequent iterations, we split the low-low component furthur.

Applying the filter to rows, followed by columns simplifies the problem as we can focus on designing good 1-D filters and then extend it to 2-D by processing rows and columns independently.

Here is the original image:

Below are the subband images, reconstructed images and error for three levels of decomposition

The code is below. A two channel synthesis filter bank is implemented which creates two subbands, one low pass the other high pass. Since the spectrum is halved, we can downsample the subband streams, to preserve the original size.

function nonunif()

clear all

close all

I = imread('lena512.bmp');

im = im2double((I));

size(im)

figure;

imshow(im);

K =3; %number of levels of decomposition

[M, N] = size(im);

%zeropad if unequal size

if (M<N)

im = [im ;repmat(zeros(1,N), (N-M), 1)];

end

if (N<M)

im = [im repmat(zeros(M,1), 1, (M-N))];

end

J = zeros(M,N);

L = zeros(M,N);

tap = 24;

fsb = 0.4;

alpha = 1;

L = im;

%define FIR filter coefficients

E00 = [-0.0087 0.0094 -0.0123 0.0235 -0.0463 0.1356 0.4623 -0.097 0.054 -0.0332 0.0221 -0.0119];

E01 = [1];

E10 = [-0.0119 0.0221 -0.0332 0.0540 -0.0970 0.4623 0.1356 -0.0463 0.0235 -0.0123 0.0094 -0.0087];

E11 = [1];

temp = [E00;E10];

h0 = temp(:)';

h1 = (-1).^[0:tap-1] .* h0;

filterlen = 24;

%subband decomposition

filtered = im;

for k = 1:K

fs = size(filtered);

blksize = fs(1)/(2^(k-1)); %size of blocks to be processed.

loop = 2^(k-1);

%newimg = zeros(fs(1)+filterlen, fs(1)+filterlen);

for p=1:loop

for q=1:loop

block = filtered((p-1)*blksize+1:p*blksize , (q-1)*blksize+1:q*blksize );

if (p == 1 && q == 1) %only process the 1st block

sb = size(block);

tempblk = zeros(sb(1)+filterlen, sb(2)+filterlen);

for i = 1:blksize

x = block(i,:);

y1 = conv(h0,x);

y2 = conv(h1,x);

sz = size(downsample(y1,2));

expandedsize = (sb(2)+filterlen)/2;

tempblk(i,1:expandedsize)= downsample(y1,2);

tempblk(i,expandedsize+1:end) = downsample(y2,2);

end

for j = 1:2*expandedsize

x = tempblk(1:blksize,j);

y1 = conv(h0,x);

y2 = conv(h1,x);

tempblk(1:expandedsize,j) = downsample(y1,2)';

tempblk(expandedsize+1:end,j) = downsample(y2,2)';

end

newimg((p-1)*2*expandedsize+1:2*p*expandedsize, (q-1)*2*expandedsize+1:2*q*expandedsize) = tempblk;

else

newimg((p-1)*2*expandedsize+1:p*2*expandedsize, (q-1)*2*expandedsize+1:q*2*expandedsize) = padarray(block, [filterlen, filterlen], 'post');

end

end

end

filtered = newimg;

figure;

imshow(filtered); %subbands image

end

%reconstruction

for k = 1:K

fs = size(filtered);

blksize = fs(1)/(2^(K-(k))); %size of blocks to be processed.

loop = 2^(K-(k));

newimg = zeros(fs(1)-loop*filterlen, fs(1)-loop*filterlen);

for p=1:loop

for q=1:loop

block = filtered((p-1)*blksize+1:p*blksize , (q-1)*blksize+1:q*blksize );

sb = size(block);

tempblk = zeros(sb(1), sb(2));

shrunksize = blksize-filterlen;

if (p == 1 && q == 1)

for i = 1:blksize

x = block(:,i);

U = [x(1:sb(1)/2)' ;x((sb(1)/2) + 1: end)'];

tempblk(1:shrunksize,i) = SynthesisBank(E00,E01,E10,E11,U);

end

for j = 1:shrunksize

x = tempblk(j,:);

U = [x(1:sb(1)/2) ; x((sb(1)/2) + 1: end)];

tempblk(j, 1:shrunksize) = SynthesisBank(E00,E01,E10,E11,U);

end

newimg((p-1)*shrunksize+1:p*shrunksize,(q-1)*shrunksize+1:q*shrunksize) = tempblk(1:shrunksize, 1:shrunksize);

else

newimg((p-1)*shrunksize+1:p*shrunksize,(q-1)*shrunksize+1:q*shrunksize) = block(1:shrunksize, 1:shrunksize);

end

end

end

filtered = newimg;

figure;

imshow(filtered); %reconstructed image

end

figure;

imshow(abs(filtered - im)); %error

end

%2 channel synthesis bank

function x = SynthesisBank(E00,E01,E10,E11,U)

N = 24; % # of taps of analysis filter

[n,m] = size(U);

if n==2 && m>2

U = U';

n = m;

end

w0 = U(:,1)+U(:,2);

w1 = U(:,1)-U(:,2);

v0 = filter(E10,E11,w0);

v1 = filter(E00,E01,w1);

y0 = upsample(v0,2);

y1 = [upsample(v1,2);0];

y0_delay = [0;y0];

x = y1+y0_delay;

x = x(N:end-2)*2;

end

Activity 2: Foreground-Background separation using RPCA

In this section we will attempt to separate the foreground and the background of surveillance videos using Robust Principle Component Analysis.

Suppose we have a matrix X, which is the sum of a low rank matrix L and a sparse matrix S. We wish to split X into these two components L and S. This is what RPCA as proposed by Candes does.

Since L is low rank and S is sparse, we wish to jointly minimize the rank of L and the L0 norm of S. But that gives rise to a non-convex problem. So RPCA bypasses this by solving convex surrogates.

(L0+S0) = arg min ||L||* + a||S||1 subject to X = L+S.

The L1 norm acts as a surrogate for the L0 norm and the nuclear norm substitutes the rank of L. the constant a acts as a tuning knob.

Now let us formulate the foreground-background separation problem into an RPCA problem. Consider N frames of a video with resolution mxn. Arrange each frame as a vector of length mn and place these vectors side by side to get a matrix (called video volume) of size mn x N

Now if the video under consideration is a surveillance video, then it is mostly stationary, with only a few activities occurring now and then. Therefore the background forms the low rank part and the activities form the sparse part. Note if we have a video where all the frames are same, then the video volume would be of rank one, and the sparse component is zero.

Formulating a surveillance video in the video volume form and applying RPCA, that is solved using the alternating direction method whose code can be found here, these are the results I got

A few drawbacks of these method are:

1. This method processes multiple frames together, that is it can't handle one frame at a time

2. Foreground objects must be moving, else they are considered part of the low rank background. For example the two people circled red are considered part of the background as they were standing still for the duration of the video.