Background Removal with Robust PCA#

Imports#

%load_ext autoreload
%autoreload 2
%matplotlib inline

sns.set(color_codes=True)
sns.set_palette(sns.color_palette("muted"))
warnings.filterwarnings("ignore")

!sudo chmod -R 777 /Landmark2/pdo/aiking/archive

# 
url = "http://backgroundmodelschallenge.eu/data/real/Video_003.zip"
path = untar_data(url); path

Path('/Landmark2/pdo/aiking/data/Video_003')

img_list = (path/"img").ls()

def fswitch(t): return t[1]+ t[0]

img_list.map(lambda e : {'name':e, 'group':fswitch(e.stem.split("_"))})

(#30) [{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00452.bmp'), 'group': '00452Img'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00351.bmp'), 'group': '00351Img'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00452.bmp'), 'group': '00452Mask'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00351.bmp'), 'group': '00351Mask'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00552.bmp'), 'group': '00552Img'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00251.bmp'), 'group': '00251Img'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00552.bmp'), 'group': '00552Mask'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00251.bmp'), 'group': '00251Mask'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00350.bmp'), 'group': '00350Mask'},{'name': Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00652.bmp'), 'group': '00652Img'}...]

(path/"img").ls()

(#30) [Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00452.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00351.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00452.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00351.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00552.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00251.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00552.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00251.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00350.bmp'),Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00652.bmp')...]

(path/"img").ls().sorted

<bound method L.sorted of [Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00452.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00351.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00452.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00351.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00552.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00251.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00552.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00251.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00350.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00652.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00652.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00350.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00250.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00250.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00550.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00550.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00651.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00450.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00651.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00450.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00551.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00252.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00551.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00252.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00451.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00650.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00352.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00451.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Mask_00352.bmp'), Path('/Landmark2/pdo/aiking/data/Video_003/img/Img_00650.bmp')]>

PILImage.open((path/"img").ls()[0])
../../_images/03_BackgroundRemovalWithRobustPCA_9_0.png 
@interact(i=(0,29))
def display_img(i=0):
    return PILImage.open((path/"img").ls()[i])

display_img()
../../_images/03_BackgroundRemovalWithRobustPCA_10_1.png

Video Clip#

video =mpe.VideoFileClip(str((path/"Video_003.avi").resolve())); video

<moviepy.video.io.VideoFileClip.VideoFileClip at 0x10a54ffa2d30>

video.subclip(0,50).ipython_display(width=300)

Moviepy - Building video __temp__.mp4.
Moviepy - Writing video __temp__.mp4

                                                                                                                                     

Moviepy - Done !
Moviepy - video ready __temp__.mp4



video.duration
# video.get_frame(0).shape[0]*video.get_frame(0).shape[1]/16
video.get_frame(0).shape

(240, 320, 3)

Preprocess Video to Matrix#

# clip = video
# k = 1
# duration = clip.duration
# # duration = 1
# scale = 50
def rgb2gray(rgb): return rgb@np.array([0.299, 0.587, 0.114])
def img2col(frame, dims): return resize(rgb2gray(frame), dims, anti_aliasing=True).astype(int).flatten()

def get_dims(frame, scale):
    w,h,c = frame.shape
    dims = (int(w*scale/100), int(h*scale/100))
    return dims
    
def vid2mat(clip, k=5, scale=50): 
    dims = get_dims(clip.get_frame(0), scale)
    return np.vstack(img2col(clip.get_frame(i/float(k)), dims) 
                     for i in range(k*int(clip.duration))).T


video.get_frame?
k=
for i in range(int(video.duration)*k):
    print(i/float(k))

0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4
2.6
2.8
3.0
3.2
3.4
3.6
3.8
4.0
4.2
4.4
4.6
4.8
5.0
5.2
5.4
5.6
5.8
6.0
6.2
6.4
6.6
6.8
7.0
7.2
7.4
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
9.4
9.6
9.8
10.0
10.2
10.4
10.6
10.8
11.0
11.2
11.4
11.6
11.8
12.0
12.2
12.4
12.6
12.8
13.0
13.2
13.4
13.6
13.8
14.0
14.2
14.4
14.6
14.8
15.0
15.2
15.4
15.6
15.8
16.0
16.2
16.4
16.6
16.8
17.0
17.2
17.4
17.6
17.8
18.0
18.2
18.4
18.6
18.8
19.0
19.2
19.4
19.6
19.8
20.0
20.2
20.4
20.6
20.8
21.0
21.2
21.4
21.6
21.8
22.0
22.2
22.4
22.6
22.8
23.0
23.2
23.4
23.6
23.8
24.0
24.2
24.4
24.6
24.8
25.0
25.2
25.4
25.6
25.8
26.0
26.2
26.4
26.6
26.8
27.0
27.2
27.4
27.6
27.8
28.0
28.2
28.4
28.6
28.8
29.0
29.2
29.4
29.6
29.8
30.0
30.2
30.4
30.6
30.8
31.0
31.2
31.4
31.6
31.8
32.0
32.2
32.4
32.6
32.8
33.0
33.2
33.4
33.6
33.8
34.0
34.2
34.4
34.6
34.8
35.0
35.2
35.4
35.6
35.8
36.0
36.2
36.4
36.6
36.8
37.0
37.2
37.4
37.6
37.8
38.0
38.2
38.4
38.6
38.8
39.0
39.2
39.4
39.6
39.8
40.0
40.2
40.4
40.6
40.8
41.0
41.2
41.4
41.6
41.8
42.0
42.2
42.4
42.6
42.8
43.0
43.2
43.4
43.6
43.8
44.0
44.2
44.4
44.6
44.8
45.0
45.2
45.4
45.6
45.8
46.0
46.2
46.4
46.6
46.8
47.0
47.2
47.4
47.6
47.8
48.0
48.2
48.4
48.6
48.8
49.0
49.2
49.4
49.6
49.8
50.0
50.2
50.4
50.6
50.8
51.0
51.2
51.4
51.6
51.8
52.0
52.2
52.4
52.6
52.8
53.0
53.2
53.4
53.6
53.8
54.0
54.2
54.4
54.6
54.8
55.0
55.2
55.4
55.6
55.8
56.0
56.2
56.4
56.6
56.8
57.0
57.2
57.4
57.6
57.8
58.0
58.2
58.4
58.6
58.8
59.0
59.2
59.4
59.6
59.8
60.0
60.2
60.4
60.6
60.8
61.0
61.2
61.4
61.6
61.8
62.0
62.2
62.4
62.6
62.8
63.0
63.2
63.4
63.6
63.8
64.0
64.2
64.4
64.6
64.8
65.0
65.2
65.4
65.6
65.8
66.0
66.2
66.4
66.6
66.8
67.0
67.2
67.4
67.6
67.8
68.0
68.2
68.4
68.6
68.8
69.0
69.2
69.4
69.6
69.8
70.0
70.2
70.4
70.6
70.8
71.0
71.2
71.4
71.6
71.8
72.0
72.2
72.4
72.6
72.8
73.0
73.2
73.4
73.6
73.8
74.0
74.2
74.4
74.6
74.8
75.0
75.2
75.4
75.6
75.8
76.0
76.2
76.4
76.6
76.8
77.0
77.2
77.4
77.6
77.8
78.0
78.2
78.4
78.6
78.8
79.0
79.2
79.4
79.6
79.8
80.0
80.2
80.4
80.6
80.8
81.0
81.2
81.4
81.6
81.8
82.0
82.2
82.4
82.6
82.8
83.0
83.2
83.4
83.6
83.8
84.0
84.2
84.4
84.6
84.8
85.0
85.2
85.4
85.6
85.8
86.0
86.2
86.4
86.6
86.8
87.0
87.2
87.4
87.6
87.8
88.0
88.2
88.4
88.6
88.8
89.0
89.2
89.4
89.6
89.8
90.0
90.2
90.4
90.6
90.8
91.0
91.2
91.4
91.6
91.8
92.0
92.2
92.4
92.6
92.8
93.0
93.2
93.4
93.6
93.8
94.0
94.2
94.4
94.6
94.8
95.0
95.2
95.4
95.6
95.8
96.0
96.2
96.4
96.6
96.8
97.0
97.2
97.4
97.6
97.8
98.0
98.2
98.4
98.6
98.8
99.0
99.2
99.4
99.6
99.8
100.0
100.2
100.4
100.6
100.8
101.0
101.2
101.4
101.6
101.8
102.0
102.2
102.4
102.6
102.8
103.0
103.2
103.4
103.6
103.8
104.0
104.2
104.4
104.6
104.8
105.0
105.2
105.4
105.6
105.8
106.0
106.2
106.4
106.6
106.8
107.0
107.2
107.4
107.6
107.8
108.0
108.2
108.4
108.6
108.8
109.0
109.2
109.4
109.6
109.8
110.0
110.2
110.4
110.6
110.8
111.0
111.2
111.4
111.6
111.8
112.0
112.2
112.4
112.6
112.8

Signature: video.get_frame(t)
Docstring:
Gets a numpy array representing the RGB picture of the clip at time t
or (mono or stereo) value for a sound clip
File:      /opt/anaconda/envs/aiking/lib/python3.9/site-packages/moviepy/Clip.py
Type:      method

dims = get_dims(video.get_frame(0), scale=25); dims
# plt.imshow(img2col(video.get_frame(0), dims).reshape(dims), cmap='gray')
# img2col(video.get_frame(0), dims).reshape(dims).min()

# plt.imshow(rgb2gray(video.get_frame(0)).astype(int), cmap='gray')

img2col(video.get_frame(0), dims)

plt.imshow(img2col(video.get_frame(0), dims).reshape(dims), cmap='gray')

<matplotlib.image.AxesImage at 0x106119faffd0>
../../_images/03_BackgroundRemovalWithRobustPCA_18_1.png 
M = vid2mat(video, k=100, scale=25); M.shape

(4800, 11300)

plt.figure(figsize=(12,12))
plt.imshow(M, cmap='gray')

<matplotlib.image.AxesImage at 0x106119fd5550>
../../_images/03_BackgroundRemovalWithRobustPCA_20_1.png 
@interact(i=IntSlider(min=0,max=11300-1, continuous_update=False))
def disp_img(i=0):
    plt.imshow(M[:,i].reshape(dims), cmap='gray')

np.save(path/"low_res_surveillance_matrix.npy", M)

SVD#

u, s, v = decomposition.randomized_svd(M,1)

u.shape, s.shape, v.shape

((4800, 1), (1,), (1, 11300))

low_rank = u@np.diag(s)@v
low_rank.shape

(4800, 11300)

plt.figure(figsize=(12,12))
plt.imshow(low_rank, cmap='gray')

<matplotlib.image.AxesImage at 0x1061198f44f0>
../../_images/03_BackgroundRemovalWithRobustPCA_27_1.png 
@interact(i=IntSlider(min=0,max=11300-1, continuous_update=False))
def disp_img(i=0):
    plt.figure()
    fig, axes= plt.subplots(2,2, figsize=(18,12))
    axes[0,0].imshow(video.get_frame(i/float(100)))
    axes[0,0].axis("off")
    axes[0,1].imshow(M[:,i].reshape(dims), cmap='gray')
    axes[0,1].axis("off")
    axes[1,1].imshow(low_rank[:,i].reshape(dims), cmap='gray')
    axes[1,1].axis("off")
    axes[1,0].imshow((M[:,i] - low_rank[:,i]).reshape(dims), cmap='gray')
    axes[1,0].axis("off")
    fig.tight_layout()
    plt.show()

disp_img(0)

<Figure size 432x288 with 0 Axes>
../../_images/03_BackgroundRemovalWithRobustPCA_29_1.png

Principal Component Analysis (PCA)#

  • Leverage the fact -> data has low intrinsic dimensionality[Dealing with high dimensionality data].

  • Alleviate the curse of dimensionality

  • Problems with scale

    • Perhaps data lies in low-dimensional subspace / low dimension manifold

  • PCA eliminates dimensions

  • Classical PCA ->

    • Best \(rank_k\) estimate L of M ( min||M-L|| where L has \(rank_k\))

    • Truncated SVD

    • Traditional PCA can handle small noise, but brittle w.r.t grossly corrupted observations–even one grossly corrupt observation can significantly mess up the answer

  • Robust PCA -> M = L+S

    • L is low-rank

    • S is sparse

    • Low-rank means that the matrix has a lot of redundant info-background

    • Sparse -> S capturing corrupted entries

  • Applications

    • Video Surveillance

    • Face Recognition

    • Latent Semantic Indexing =:

      • L captures common words used in all documents

      • S captures a few keywords that best distinguish each document from others

    • Ranking and Collaborative Filtering =: Small portion of available rankings could be noisy & even tampered with

  • The L1 norm induces sparsity :: $\(||x||_1=1\)$

    • Contours of the loss function

    • Corner of a rectangle

  • Robust PCA Optimization Problem :: \(min ||L||_* + \lambda||S||_1 \text{ subject to } L+S=M\) $\(||.||_1\)$ is the L1 norm -> results i sparse values [ max. absolute column norm]

    • $\(||.||_*\)$ is the nuclear norm -> L1 norm of singular values -> sparse singular values -> low rank