DECOLOR: Moving Object Detection by Detecting Contiguous Outliers
in the Low-Rank Representation
Xiaowei Zhou et al.
Abstract—Object detection is a fundamental step for automated video analysis in many vision applications. Object detection in a video
is usually performed by object detectors or background subtraction techniques. Often, an object detector requires manually labeled
examples to train a binary classifier, while background subtraction needs a training sequence that contains no objects to build a background model. To automate the analysis, object detection without a separate training phase becomes a critical task. People have tried to tackle this task by using motion information. But existing motion-based methods are usually limited when coping with complex scenarios such as nonrigid motion and dynamic background. In this paper, we show that the above challenges can be addressed in a unified framework named DEtecting Contiguous Outliers in the LOw-rank Representation (DECOLOR). This formulation integrates object detection and background learning into a single process of optimization, which can be solved by an alternating algorithm efficiently. We explain the relations between DECOLOR and other sparsity-based methods. Experiments on both simulated data and real sequences demonstrate that DECOLOR outperforms the state-of-the-art approaches and it can work effectively on a wide range of complex scenarios.
本文的三个贡献点，分别是：
１．提出了一个在低秩框架下的离群点检测方法，即：DECOLOR．
２．引入了连续先验，用ＭＲＦ来对相邻像素点之间的联系进行建模．
３．结合了参数运动模型，来弥补相机的运动．使得DECOLOR在动态背景下，仍可以取得不错的效果．
3 C ONTIGUOUS O UTLIER D ETECTION IN THE L OW -R ANK R EPRESENTATION
３.２. Formulation
本文将输入的视频记为Ｄ，那么要做的工作就是将Ｄ分解为前景F和背景B,其模型的构建分别为：
背景模型：在不考虑动态背景的情况下，视频的背景在没有光照的变化和周期性的纹理的影响时，背景的强度应该是不变的．所以，背景图像可以看作是现行相关的，构成了一个低秩矩阵B,于是就有了：
rank(B) <= K;
前景模型：
所谓前景，也就是相对于背景而言运动的物体，也可以说是 any object that moves differently from the background. 前景目标的运动导致了强度的变化，这就使得低秩的背景无法模拟此类像素点．所以，在低秩表示的时候，这些像素点就成了离群点．另外，就是文章引入了一个先验，即：前景目标应该是连续的区域并且相对较小．可以用马尔科夫随机场来对前景像素点进行建模．马尔科夫随机场是一种图模型(Graph Model)，谈到图模型，自然少不了边和点．这里，就是将单个像素点比作是图的一个顶点，而两个像素点之间的联系，就可看作是图的边．根据Ising model 可以得到如下的关于前景S的能量函数：
其中，Uij 代表前景支持Sij的一元势能， Lambda ij,kl 控制着像素点Sij 和Skl之间相关性的强度．一元势能Uij的定义为：
从这个定义看出，ij位置的像素点为前景时，其一元能量为 Lambda ij, 当其为背景时，则是０．当最小化公式（４）时，会使得 Uij 最小，也就使得前景像素点尽量稀疏了．( sparse foreground ). 从公式（４）的第二项可以看出，当ij 和kl位置处的像素点相同时，即 Sij = Skl=0或者＝１时，第二项为０;否则就是１．这样最小化公式（４）也使得前景区域和背景区域尽可能的平滑，
信号模型：
信号模型描述了Ｄ的构成．在背景区域，Sij = 0, 我们认为 Dij = Bij + εij, 其中εij代表独立同分布的高斯噪声．在前景区域，Sij = 1,背景区域被前景区域所遮挡，所以，Dij等于前景的强度．
将上述三个模型统一起来，就是提出的如下model:
为了使得上述最小化能量可解，作了相关的松弛，即： 用核范数取代了 rank(B)<=K, 将公式(6)写成其对偶的形式，并且引入矩阵操作符，得到最终的能量函数：
其中，A 是图的邻结矩阵．
3.3. Algorithm
公式（７）的目标函数是非凸的，包括连续的和离散的变量．联合的优化是非常困难的，所以，在这里采用了迭代的方式，分别进行前景S和背景B的求解．B-step是一个凸优化问题，S-step是联合优化问题．
3.3.1 Estimation of the Low-Rank Matrix B
给定前景支持S的预测，那么最小化公式（７）就变成了如下的矩阵完全问题：
这是从部分观察来学习一个低秩的矩阵．该问题可以用 soft-impute 算法来求解，并且用到了如下的引理：
重写公式(8),我们有：
利用引理１，通过迭代的利用下列式子，可以得到公式（８）的最优解．
SOFT-IMPUTE 算法对应的ｃｏｄｅ如下：
1 %% function for soft-impute
2 function [Z,Znorm,alpha] = softImpute(X,Z,Omega,alpha0,maxRank)
3 %
4 % This program implements the soft-impute algorithm followed by
5 % postprocessing in the Matrix completion paper Mazumder'10 IJML
6 % min || Z - X ||_Omega + \alpha || Z ||_Nulear
7 % \alpha is decrease from alpha0 to the minima value that makes rank(Z) <= maxRank
8
9 % X is the incomplete matrix
10 % maxRank is the desired rank in the constraint
11 % Omega is the mask with value 1 for data and 0 for missing part
12 if isempty(Z)
13 Z = X;
14 end
15 if isempty(Omega)
16 Omega = true(size(X));
17 end
18 if isempty(alpha0)
19 [U,D] = svd(X,'econ');
20 alpha0 = D(2,2);
21 end
22 if isempty(maxRank)
23 maxRank = -1;
24 end
25 % parameters
26 eta = 0.707;
27 epsilon = 1e-4;
28 maxInnerIts = 20;
29 %% trivial
30 % no rank constraint
31 if maxRank >= min(size(X))
32 Z = X;
33 [U,D] = svd(Z,'econ');
34 Znorm = sum(diag(D));
35 alpha = 0;
36 return;
37 end
38 % no observation
39 if sum(Omega(:)) == 0
40 % no data
41 Z = zeros(size(X));
42 Znorm = 0;
43 alpha = alpha0;
44 return;
45 end
46 %% soft-impute
47 % 1. initialize
48 outIts = 0;
49 alpha = alpha0;
50 % 2. Do for alpha = alpha0 > alpha_1 > alpha_2 > ... > alpha_maxRank
51 disp('begin soft-impute iterations');
52 while 1
53 outIts = outIts + 1;
54 energy = inf;
55 for innerIts = 1:maxInnerIts
56 % (a)i
57 C = X.*Omega + Z.*(1-Omega);
58 [U,D,V] = svd(C,'econ');
59 VT = V';
60 % soft impute
61 d = diag(D);
62 idx = find(d > alpha);
63 Z = U(:,idx) * diag( d(idx) - alpha ) * VT(idx,:);
64 % (a)ii
65 Znorm = sum(d(idx)-alpha);
66 energy_old = energy;
67 energy = alpha*Znorm + norm(Z(Omega(:))-X(Omega(:)),'fro')/2;
68 if abs(energy - energy_old) / energy_old < epsilon
69 break
70 end
71 end
72 % check termination condition of alpha
73 k = length(idx); % rank of Z
74 disp(['alpha = ' num2str(alpha) '; rank = ' num2str(k) '; number of iteration: ' num2str(innerIts)]);
75 if k <= maxRank && alpha > 1e-3
76 alpha = alpha*eta;
77 else
78 break;
79 end
80 end
81 end
3.3.2 Estimation of the Outlier Support S
在假设给定低秩矩阵B的情况下，如何求前景支持S ? 我们可以得到如下的公式：
从文章中得知，当给定B的时候，C 是一个常数，那么可以忽略，Sij前面的那一大堆项，也是定量，于是就成了 Sij 和 ||A vec(S)||1两家独大，这两项即是马尔科夫随机场的两项．
于是，通过上述两种方式，就可以迭代的求解出前景S和背景B.
对应的code 为：
1 for i = 1:size(Dtau,2) 2 GCO_SetDataCost( hMRF, (amplify/gamma)*[ 0.5*(E(:,i)).^2, ~OmegaOut(:,i)*beta + OmegaOut(:,i)*0.5*max(E(:)).^2]' ); 3 GCO_Expansion(hMRF); 4 Omega(:,i) = ( GCO_GetLabeling(hMRF) == 1 )'; 5 energy_cut = energy_cut + double( GCO_ComputeEnergy(hMRF) ); 6 end
额，说累了，歇会．．．出去走走．．．
4 EXTENSION TO MOVING ACKGROUND
Here, we use the 2D parametric transforms [60] to model the translation, rotation, and planar deformation of the background.
Dj ○ τj 代表被向量 τj ∈IRp转换后的第ｊ帧图像，其中ｐ表示运动模型的参数的编号，如：p=6代表仿射运动，p=8代表投影运动．所以，本文所提出的分解变成了 D○τ = B+Ｅ＋ε，其中, D○τ = [D1○τ1, ... , Dn○τn], τ 是一个向量代表所有 τj 的集合．
接下来，我们用D○τ代替公式（７）中的Ｄ，并且用B,S来预测 τ,迭代的最小化下列公式：
现在开始讨论如何通过 τ 来最小化公式（２０）：
此处，我们利用增量优化的方法来解决这个参数运动估计问题：在每一次迭代，我们以小的增幅 △τ 更新 τ， 将 D○τ线性化为D○τ + J△τ ，
其中，J为雅克比矩阵．所以，τ 可以用下列方式进行更新：
通过△τ 的最小化过程时一个最小二乘问题，有闭合解．
实际上，τ１，．．．，τn 的更新可以分别独立的进行，为了加速DECOLOR的收敛，我们初始化τ，粗略的将每一帧Dj和中间帧Dn/2进行对比．这个预对齐的过程是用robust multiresolution method来进行的．
总的算法，归结为：
－－－－－－－－－－－－－－ 算法部分结束 －－－－－－－－－－－－－－－
代码解读：
1.RUN_REAL_MOVING.m
1 clear all
2 close all
3
4 addpath('internal');
5 addpath(genpath('gco-v3.0'));
6
7 %% data
8 % static background
9 dataList = {'people1','people2','cars6','cars7'};
10
11 for dataID = 1:4
12
13 dataName = dataList{dataID};
14 load(['data\' dataName],'ImData');
15
16 %% run DECOLOR
17 opt.flagAlign = 1;
18 opt.tol = 1e-3;
19 [LowRank,Mask,tau,info] = ObjDetection_DECOLOR(ImData,opt);
20
21 % warp masks to match the original images
22 for i = 1:size(ImData,3)
23 % use [Iwarp,Omega] = warpImg(I,tau,mode,extrapval)
24 Mask(:,:,i) = warpImg(double(Mask(:,:,i)),tau(:,i),1,0)>0.5;
25 cropRatio = 0.01;
26 Mask([1:round(cropRatio*size(Mask,1)),round((1-cropRatio)*size(Mask,1)):end],:,i) = false;
27 Mask(:,[1:round(cropRatio*size(Mask,2)),round((1-cropRatio)*size(Mask,2)):end],i) = false;
28 end
29 save(['result\' dataName '_DECOLOR.mat'],'dataName','Mask','LowRank','tau','info');
30
31
32 %% displaying
33 load(['result\' dataName '_DECOLOR.mat'],'dataName','Mask','LowRank','tau');
34 moviename = ['result\' dataName,'_DECOLOR.avi']; fps = 12;
35 mov = avifile(moviename,'fps',fps,'compression','none');
36 for i = 1:size(ImData,3)
37 figure(1); clf;
38 subplot(2,2,1);
39 imshow(ImData(:,:,i)), axis off, colormap gray; axis off;
40 title('Original image','fontsize',12);
41 subplot(2,2,2);
42 imshow(LowRank(:,:,i)), axis off,colormap gray; axis off;
43 title('Low Rank','fontsize',12);
44 subplot(2,2,3);
45 imshow(ImData(:,:,i)), axis off,colormap gray; axis off;
46 hold on; contour(Mask(:,:,i),[0 0],'y','linewidth',5);
47 title('Segmentation','fontsize',12);
48 subplot(2,2,4);
49 imshow(ImData(:,:,i).*uint8(Mask(:,:,i))), axis off, colormap gray; axis off;
50 title('Foreground','fontsize',12);
51 mov = addframe(mov,getframe(1));
52 end
53 h = close(mov);
54
55 end
2. ObjDetection_DECOLOR.m
1 function [LowRank,Mask,tau,info] = ObjDetection_DECOLOR(ImData,opt)
2 % This function use DECOLOR to detect moving objects in sequence ImData
3 % http://arxiv.org/PS_cache/arxiv/pdf/1109/1109.0882v1.pdf
4 % eexwzhou@ust.hk
5 % Syntex: [LowRank,Mask,tau] = ObjDetection_DECOLOR(ImData).
6 % Input:
7 % ImData -- 3D array representing a image sequence.
8 % Each image is stored in ImData(:,:,i).
9 % opt -- options. Usually, default setting is good. No need to specify.
10 % opt.K: desired rank of the estimated low-rank component.
11 % Default: \sqrt(min(size(D))) is good generally.
12 % opt.lambda: a constant controls the strength of smoothness regularize
13 % lambda ~ [1 5] is recommended. Default: 5
14 % opt.sigma: STD of noise in the image. If not specified, computed online
15 % opt.flagAlign: whether alighment is needed or not.
16 % opt.tol: convergence precision. Default: 1e-4
17 % Output:
18 % LowRank -- Low-rank background component
19 % Mask -- Segmented object mask
20 % tau - transformation parameters to compensate for camera motion
21 % info -- other information
22
23 disp('^_^^_^^_^^_^^_^^_^ DECOLOR ^_^^_^^_^^_^^_^');
24 tic;
25
26 %% default parameter setting
27 if ~exist('opt','var'); opt = []; end
28 if ~isfield(opt,'tol'); opt.tol = 1e-4; end
29 if ~isfield(opt,'K'); opt.K = floor(sqrt(size(ImData,3))); end
30 if ~isfield(opt,'lambda'); opt.lambda = 5; end % gamma = opt.lambda * beta;
31 if ~isfield(opt,'sigma'); opt.sigma = []; end % sigma can be estimated online
32 if ~isfield(opt,'flagAlign'); opt.flagAlign = false; end % alignment or not
33
34 %% variable initialize
35 ImData = mat2gray(ImData); % 0~1
36 ImMean = mean(ImData(:));
37 ImData = ImData - ImMean; % subtract mean is recommended
38 numImg = size(ImData,3);
39 sizeImg = [size(ImData,1),size(ImData,2)];
40 if opt.flagAlign == true && sizeImg(2) > 1
41 disp('----------- Pre-alignment ----------');
42 [ImTrans,tau] = preAlign(ImData);
43 Dtau = reshape(ImTrans,prod(sizeImg),numImg);
44 else
45 tau = [];
46 Dtau = reshape(ImData,prod(sizeImg),numImg);
47 end
48 maxOuterIts = 20;
49 alpha = []; % Default setting by soft-impute
50 beta = 0.5*(std(Dtau(:,1)))^2; % Start from a big value
51 minbeta = 0.5*(3*std(Dtau(:,1))/20)^2; % lower bound: suppose SNR <= 20
52 sigma = opt.sigma; % if empty, will be estimated online
53 B = Dtau; % the low-rank matrix
54 Omega = true(size(Dtau)); % background support
55 OmegaOut = false(size(Dtau)); % OmegaOut is to record the extrapolated regions
56 ObjArea = sum(~Omega(:));
57 minObjArea = numel(Dtau(:,1))/1e4; % minimum number of outliers
58
59 % graph cuts initialization
60 % GCO toolbox is called
61 if opt.lambda > 0
62 hMRF = GCO_Create(prod(sizeImg),2);
63 GCO_SetSmoothCost( hMRF, [0 1;1 0] );
64 AdjMatrix = getAdj(sizeImg);
65 amplify = 10 * opt.lambda;
66 GCO_SetNeighbors( hMRF, amplify * AdjMatrix );
67 end
68
69 %% outer loop
70 energy_old = inf; % total energy
71 for outerIts = 1:maxOuterIts
72 disp(['---------------- Outer Loop: ' num2str(outerIts) ' ----------------']);
73
74 %% update tau
75 if opt.flagAlign == true && sizeImg(2) > 1
76 disp('*** Estimate Transformation ***');
77 for i = 1:numImg
78 % update once
79 [Iwarp,tau(:,i),dummy,Lout] = regImg(reshape(B(:,i),sizeImg),ImData(:,:,i),tau(:,i),double(reshape(Omega(:,i),sizeImg)),1);
80 Dtau(:,i) = reshape(Iwarp,prod(sizeImg),1);
81 OmegaOut(:,i) = reshape(Lout,prod(sizeImg),1); % extrapolated regions
82 end
83 end
84
85 %% update B
86 disp('*** Estimate Low-rank Matrix *** ');
87 [B,Bnorm,alpha] = softImpute(Dtau,B,~OmegaOut&Omega,alpha,opt.K);
88 E = Dtau - B;
89
90 %% estimate sigma
91 if isempty(opt.sigma)
92 sigma_old = sigma;
93 residue = sort(E(~OmegaOut(:)&Omega(:)));
94 truncate = 0.005;
95 idx1 = round(truncate*length(residue))+1;
96 idx2 = round((1-truncate)*length(residue));
97 sigma = std(residue(idx1:idx2));
98 if abs(sigma_old-sigma)/abs(sigma_old) < 0.01
99 sigma = sigma_old;
100 end
101 end
102 % update beta
103 if ObjArea < minObjArea
104 beta = beta/2;
105 else
106 beta = min(max([beta/2,0.5*(3*sigma)^2 minbeta]),beta);
107 end
108 gamma = opt.lambda * beta;
109
110 %% estimate S = ~Omega;
111 disp('*** Estimate Outlier Support *** ');
112 disp(['$$$ beta = ' num2str(beta) '; gamma = ' num2str(gamma) '; sigma = ' num2str(sigma)]);
113 if opt.lambda > 0
114 % call GCO to run graph cuts
115 energy_cut = 0;
116 for i = 1:size(Dtau,2)
117 GCO_SetDataCost( hMRF, (amplify/gamma)*[ 0.5*(E(:,i)).^2, ~OmegaOut(:,i)*beta + OmegaOut(:,i)*0.5*max(E(:)).^2]' );
118 GCO_Expansion(hMRF);
119 Omega(:,i) = ( GCO_GetLabeling(hMRF) == 1 )';
120 energy_cut = energy_cut + double( GCO_ComputeEnergy(hMRF) );
121 end
122 ObjArea = sum(Omega(:)==0);
123 energy_cut = (gamma/amplify) * energy_cut;
124 else
125 % direct hard thresholding if no smoothness
126 Omega = 0.5*E.^2 < beta;
127 ObjArea = sum(Omega(:)==0);
128 energy_cut = 0.5*norm(Dtau-B-E,'fro')^2 + beta*ObjArea;
129 end
130
131 %% display energy
132 energy = energy_cut + alpha * Bnorm;
133 disp(['>>> the object area is ' num2str(ObjArea)]);
134 disp(['>>> the objectvive energy is ' num2str(energy)]);
135
136 %% check termination condition
137 if ObjArea > minObjArea && abs(energy_old-energy)/energy < opt.tol; break; end
138 energy_old = energy;
139
140 end
141
142 LowRank = uint8(mat2gray(reshape(B,size(ImData))+ImMean)*256);
143 Mask = reshape(~Omega,size(ImData));
144
145 info.opt = opt;
146 info.time = toc;
147 info.outerIts = outerIts;
148 info.energy = energy;
149 info.rank = rank(B);
150 info.alpha = alpha;
151 info.beta = beta;
152 info.sigma = sigma;
153
154 if opt.lambda > 0
155 GCO_Delete(hMRF);
156 end
157
158 end
159
160
161
162 %% function to get the adjcent matirx of the graph
163 function W = getAdj(sizeData)
164 numSites = prod(sizeData);
165 id1 = [1:numSites, 1:numSites, 1:numSites];
166 id2 = [ 1+1:numSites+1,...
167 1+sizeData(1):numSites+sizeData(1),...
168 1+sizeData(1)*sizeData(2):numSites+sizeData(1)*sizeData(2)];
169 value = ones(1,3*numSites);
170 W = sparse(id1,id2,value);
171 W = W(1:numSites,1:numSites);
172 end
173
174
175 %% function for soft-impute
176 function [Z,Znorm,alpha] = softImpute(X,Z,Omega,alpha0,maxRank)
177 %
178 % This program implements the soft-impute algorithm followed by
179 % postprocessing in the Matrix completion paper Mazumder'10 IJML
180 % min || Z - X ||_Omega + \alpha || Z ||_Nulear
181 % \alpha is decrease from alpha0 to the minima value that makes rank(Z) <= maxRank
182
183 % X is the incomplete matrix
184 % maxRank is the desired rank in the constraint
185 % Omega is the mask with value 1 for data and 0 for missing part
186 if isempty(Z)
187 Z = X;
188 end
189 if isempty(Omega)
190 Omega = true(size(X));
191 end
192 if isempty(alpha0)
193 [U,D] = svd(X,'econ');
194 alpha0 = D(2,2);
195 end
196 if isempty(maxRank)
197 maxRank = -1;
198 end
199 % parameters
200 eta = 0.707;
201 epsilon = 1e-4;
202 maxInnerIts = 20;
203 %% trivial
204 % no rank constraint
205 if maxRank >= min(size(X))
206 Z = X;
207 [U,D] = svd(Z,'econ');
208 Znorm = sum(diag(D));
209 alpha = 0;
210 return;
211 end
212 % no observation
213 if sum(Omega(:)) == 0
214 % no data
215 Z = zeros(size(X));
216 Znorm = 0;
217 alpha = alpha0;
218 return;
219 end
220 %% soft-impute
221 % 1. initialize
222 outIts = 0;
223 alpha = alpha0;
224 % 2. Do for alpha = alpha0 > alpha_1 > alpha_2 > ... > alpha_maxRank
225 disp('begin soft-impute iterations');
226 while 1
227 outIts = outIts + 1;
228 energy = inf;
229 for innerIts = 1:maxInnerIts
230 % (a)i
231 C = X.*Omega + Z.*(1-Omega);
232 [U,D,V] = svd(C,'econ');
233 VT = V';
234 % soft impute
235 d = diag(D);
236 idx = find(d > alpha);
237 Z = U(:,idx) * diag( d(idx) - alpha ) * VT(idx,:);
238 % (a)ii
239 Znorm = sum(d(idx)-alpha);
240 energy_old = energy;
241 energy = alpha*Znorm + norm(Z(Omega(:))-X(Omega(:)),'fro')/2;
242 if abs(energy - energy_old) / energy_old < epsilon
243 break
244 end
245 end
246 % check termination condition of alpha
247 k = length(idx); % rank of Z
248 disp(['alpha = ' num2str(alpha) '; rank = ' num2str(k) '; number of iteration: ' num2str(innerIts)]);
249 if k <= maxRank && alpha > 1e-3
250 alpha = alpha*eta;
251 else
252 break;
253 end
254 end
255 end