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 算法对应的ｃｏｄｅ如下：
%% function for soft-impute
function [Z,Znorm,alpha] = softImpute(X,Z,Omega,alpha0,maxRank)
%
% This program implements the soft-impute algorithm followed by
% postprocessing in the Matrix completion paper Mazumder'10 IJML
% min || Z - X ||_Omega + \alpha || Z ||_Nulear
% \alpha is decrease from alpha0 to the minima value that makes rank(Z) <= maxRank % X is the incomplete matrix
% maxRank is the desired rank in the constraint
% Omega is the mask with value for data and for missing part
if isempty(Z)
Z = X;
end
if isempty(Omega)
Omega = true(size(X));
end
if isempty(alpha0)
[U,D] = svd(X,'econ');
alpha0 = D(,);
end
if isempty(maxRank)
maxRank = -;
end
% parameters
eta = 0.707;
epsilon = 1e-;
maxInnerIts = ;
%% trivial
% no rank constraint
if maxRank >= min(size(X))
Z = X;
[U,D] = svd(Z,'econ');
Znorm = sum(diag(D));
alpha = ;
return;
end
% no observation
if sum(Omega(:)) ==
% no data
Z = zeros(size(X));
Znorm = ;
alpha = alpha0;
return;
end
%% soft-impute
% . initialize
outIts = ;
alpha = alpha0;
% . Do for alpha = alpha0 > alpha_1 > alpha_2 > ... > alpha_maxRank
disp('begin soft-impute iterations');
while
outIts = outIts + ;
energy = inf;
for innerIts = :maxInnerIts
% (a)i
C = X.*Omega + Z.*(-Omega);
[U,D,V] = svd(C,'econ');
VT = V';
% soft impute
d = diag(D);
idx = find(d > alpha);
Z = U(:,idx) * diag( d(idx) - alpha ) * VT(idx,:);
% (a)ii
Znorm = sum(d(idx)-alpha);
energy_old = energy;
energy = alpha*Znorm + norm(Z(Omega(:))-X(Omega(:)),'fro')/;
if abs(energy - energy_old) / energy_old < epsilon
break
end
end
% check termination condition of alpha
k = length(idx); % rank of Z
disp(['alpha = ' num2str(alpha) '; rank = ' num2str(k) '; number of iteration: ' num2str(innerIts)]);
if k <= maxRank && alpha > 1e-
alpha = alpha*eta;
else
break;
end
end
end
3.3.2 Estimation of the Outlier Support S
在假设给定低秩矩阵B的情况下，如何求前景支持S ? 我们可以得到如下的公式：
从文章中得知，
当给定B的时候，C 是一个常数，那么可以忽略，Sij前面的那一大堆项，也是定量，于是就成了 Sij 和 ||A vec(S)||1两家独大，这两项即是马尔科夫随机场的两项．
于是，通过上述两种方式，就可以迭代的求解出前景S和背景B.
对应的code 为：
for i = :size(Dtau,)
GCO_SetDataCost( hMRF, (amplify/gamma)*[ 0.5*(E(:,i)).^, ~OmegaOut(:,i)*beta + OmegaOut(:,i)*0.5*max(E(:)).^]' );
GCO_Expansion(hMRF);
Omega(:,i) = ( GCO_GetLabeling(hMRF) == )';
energy_cut = energy_cut + double( GCO_ComputeEnergy(hMRF) );
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
clear all
close all addpath('internal');
addpath(genpath('gco-v3.0')); %% data
% static background
dataList = {'people1','people2','cars6','cars7'}; for dataID = : dataName = dataList{dataID};
load(['data\' dataName],'ImData'); %% run DECOLOR
opt.flagAlign = ;
opt.tol = 1e-;
[LowRank,Mask,tau,info] = ObjDetection_DECOLOR(ImData,opt); % warp masks to match the original images
for i = :size(ImData,)
% use [Iwarp,Omega] = warpImg(I,tau,mode,extrapval)
Mask(:,:,i) = warpImg(double(Mask(:,:,i)),tau(:,i),,)>0.5;
cropRatio = 0.01;
Mask([:round(cropRatio*size(Mask,)),round((-cropRatio)*size(Mask,)):end],:,i) = false;
Mask(:,[:round(cropRatio*size(Mask,)),round((-cropRatio)*size(Mask,)):end],i) = false;
end
save(['result\' dataName '_DECOLOR.mat'],'dataName','Mask','LowRank','tau','info'); %% displaying
load(['result\' dataName '_DECOLOR.mat'],'dataName','Mask','LowRank','tau');
moviename = ['result\' dataName,'_DECOLOR.avi']; fps = 12;
mov = avifile(moviename,'fps',fps,'compression','none');
for i = :size(ImData,)
figure(); clf;
subplot(,,);
imshow(ImData(:,:,i)), axis off, colormap gray; axis off;
title('Original image','fontsize',);
subplot(,,);
imshow(LowRank(:,:,i)), axis off,colormap gray; axis off;
title('Low Rank','fontsize',);
subplot(,,);
imshow(ImData(:,:,i)), axis off,colormap gray; axis off;
hold on; contour(Mask(:,:,i),[ ],'y','linewidth',);
title('Segmentation','fontsize',);
subplot(,,);
imshow(ImData(:,:,i).*uint8(Mask(:,:,i))), axis off, colormap gray; axis off;
title('Foreground','fontsize',);
mov = addframe(mov,getframe());
end
h = close(mov); end
2. ObjDetection_DECOLOR.m
function [LowRank,Mask,tau,info] = ObjDetection_DECOLOR(ImData,opt)
% This function use DECOLOR to detect moving objects in sequence ImData
% http://arxiv.org/PS_cache/arxiv/pdf/1109/1109.0882v1.pdf
% eexwzhou@ust.hk
% Syntex: [LowRank,Mask,tau] = ObjDetection_DECOLOR(ImData).
% Input:
% ImData -- 3D array representing a image sequence.
% Each image is stored in ImData(:,:,i).
% opt -- options. Usually, default setting is good. No need to specify.
% opt.K: desired rank of the estimated low-rank component.
% Default: \sqrt(min(size(D))) is good generally.
% opt.lambda: a constant controls the strength of smoothness regularize
% lambda ~ [ ] is recommended. Default:
% opt.sigma: STD of noise in the image. If not specified, computed online
% opt.flagAlign: whether alighment is needed or not.
% opt.tol: convergence precision. Default: 1e-
% Output:
% LowRank -- Low-rank background component
% Mask -- Segmented object mask
% tau - transformation parameters to compensate for camera motion
% info -- other information disp('^_^^_^^_^^_^^_^^_^ DECOLOR ^_^^_^^_^^_^^_^');
tic; %% default parameter setting
if ~exist('opt','var'); opt = []; end
if ~isfield(opt,'tol'); opt.tol = 1e-; end
if ~isfield(opt,'K'); opt.K = floor(sqrt(size(ImData,))); end
if ~isfield(opt,'lambda'); opt.lambda = ; end % gamma = opt.lambda * beta;
if ~isfield(opt,'sigma'); opt.sigma = []; end % sigma can be estimated online
if ~isfield(opt,'flagAlign'); opt.flagAlign = false; end % alignment or not %% variable initialize
ImData = mat2gray(ImData); % ~
ImMean = mean(ImData(:));
ImData = ImData - ImMean; % subtract mean is recommended
numImg = size(ImData,);
sizeImg = [size(ImData,),size(ImData,)];
if opt.flagAlign == true && sizeImg() >
disp('----------- Pre-alignment ----------');
[ImTrans,tau] = preAlign(ImData);
Dtau = reshape(ImTrans,prod(sizeImg),numImg);
else
tau = [];
Dtau = reshape(ImData,prod(sizeImg),numImg);
end
maxOuterIts = ;
alpha = []; % Default setting by soft-impute
beta = 0.5*(std(Dtau(:,)))^; % Start from a big value
minbeta = 0.5*(*std(Dtau(:,))/)^; % lower bound: suppose SNR <=
sigma = opt.sigma; % if empty, will be estimated online
B = Dtau; % the low-rank matrix
Omega = true(size(Dtau)); % background support
OmegaOut = false(size(Dtau)); % OmegaOut is to record the extrapolated regions
ObjArea = sum(~Omega(:));
minObjArea = numel(Dtau(:,))/1e4; % minimum number of outliers % graph cuts initialization
% GCO toolbox is called
if opt.lambda >
hMRF = GCO_Create(prod(sizeImg),);
GCO_SetSmoothCost( hMRF, [ ; ] );
AdjMatrix = getAdj(sizeImg);
amplify = * opt.lambda;
GCO_SetNeighbors( hMRF, amplify * AdjMatrix );
end %% outer loop
energy_old = inf; % total energy
for outerIts = :maxOuterIts
disp(['---------------- Outer Loop: ' num2str(outerIts) ' ----------------']); %% update tau
if opt.flagAlign == true && sizeImg() >
disp('*** Estimate Transformation ***');
for i = :numImg
% update once
[Iwarp,tau(:,i),dummy,Lout] = regImg(reshape(B(:,i),sizeImg),ImData(:,:,i),tau(:,i),double(reshape(Omega(:,i),sizeImg)),);
Dtau(:,i) = reshape(Iwarp,prod(sizeImg),);
OmegaOut(:,i) = reshape(Lout,prod(sizeImg),); % extrapolated regions
end
end %% update B
disp('*** Estimate Low-rank Matrix *** ');
[B,Bnorm,alpha] = softImpute(Dtau,B,~OmegaOut&Omega,alpha,opt.K);
E = Dtau - B; %% estimate sigma
if isempty(opt.sigma)
sigma_old = sigma;
residue = sort(E(~OmegaOut(:)&Omega(:)));
truncate = 0.005;
idx1 = round(truncate*length(residue))+;
idx2 = round((-truncate)*length(residue));
sigma = std(residue(idx1:idx2));
if abs(sigma_old-sigma)/abs(sigma_old) < 0.01
sigma = sigma_old;
end
end
% update beta
if ObjArea < minObjArea
beta = beta/;
else
beta = min(max([beta/,0.5*(*sigma)^ minbeta]),beta);
end
gamma = opt.lambda * beta; %% estimate S = ~Omega;
disp('*** Estimate Outlier Support *** ');
disp(['$$$ beta = ' num2str(beta) '; gamma = ' num2str(gamma) '; sigma = ' num2str(sigma)]);
if opt.lambda >
% call GCO to run graph cuts
energy_cut = ;
for i = :size(Dtau,)
GCO_SetDataCost( hMRF, (amplify/gamma)*[ 0.5*(E(:,i)).^, ~OmegaOut(:,i)*beta + OmegaOut(:,i)*0.5*max(E(:)).^]' );
GCO_Expansion(hMRF);
Omega(:,i) = ( GCO_GetLabeling(hMRF) == )';
energy_cut = energy_cut + double( GCO_ComputeEnergy(hMRF) );
end
ObjArea = sum(Omega(:)==);
energy_cut = (gamma/amplify) * energy_cut;
else
% direct hard thresholding if no smoothness
Omega = 0.5*E.^ < beta;
ObjArea = sum(Omega(:)==);
energy_cut = 0.5*norm(Dtau-B-E,'fro')^ + beta*ObjArea;
end %% display energy
energy = energy_cut + alpha * Bnorm;
disp(['>>> the object area is ' num2str(ObjArea)]);
disp(['>>> the objectvive energy is ' num2str(energy)]); %% check termination condition
if ObjArea > minObjArea && abs(energy_old-energy)/energy < opt.tol; break; end
energy_old = energy; end LowRank = uint8(mat2gray(reshape(B,size(ImData))+ImMean)*);
Mask = reshape(~Omega,size(ImData)); info.opt = opt;
info.time = toc;
info.outerIts = outerIts;
info.energy = energy;
info.rank = rank(B);
info.alpha = alpha;
info.beta = beta;
info.sigma = sigma; if opt.lambda >
GCO_Delete(hMRF);
end end %% function to get the adjcent matirx of the graph
function W = getAdj(sizeData)
numSites = prod(sizeData);
id1 = [:numSites, :numSites, :numSites];
id2 = [ +:numSites+,...
+sizeData():numSites+sizeData(),...
+sizeData()*sizeData():numSites+sizeData()*sizeData()];
value = ones(,*numSites);
W = sparse(id1,id2,value);
W = W(:numSites,:numSites);
end %% function for soft-impute
function [Z,Znorm,alpha] = softImpute(X,Z,Omega,alpha0,maxRank)
%
% This program implements the soft-impute algorithm followed by
% postprocessing in the Matrix completion paper Mazumder'10 IJML
% min || Z - X ||_Omega + \alpha || Z ||_Nulear
% \alpha is decrease from alpha0 to the minima value that makes rank(Z) <= maxRank % X is the incomplete matrix
% maxRank is the desired rank in the constraint
% Omega is the mask with value for data and for missing part
if isempty(Z)
Z = X;
end
if isempty(Omega)
Omega = true(size(X));
end
if isempty(alpha0)
[U,D] = svd(X,'econ');
alpha0 = D(,);
end
if isempty(maxRank)
maxRank = -;
end
% parameters
eta = 0.707;
epsilon = 1e-;
maxInnerIts = ;
%% trivial
% no rank constraint
if maxRank >= min(size(X))
Z = X;
[U,D] = svd(Z,'econ');
Znorm = sum(diag(D));
alpha = ;
return;
end
% no observation
if sum(Omega(:)) ==
% no data
Z = zeros(size(X));
Znorm = ;
alpha = alpha0;
return;
end
%% soft-impute
% . initialize
outIts = ;
alpha = alpha0;
% . Do for alpha = alpha0 > alpha_1 > alpha_2 > ... > alpha_maxRank
disp('begin soft-impute iterations');
while
outIts = outIts + ;
energy = inf;
for innerIts = :maxInnerIts
% (a)i
C = X.*Omega + Z.*(-Omega);
[U,D,V] = svd(C,'econ');
VT = V';
% soft impute
d = diag(D);
idx = find(d > alpha);
Z = U(:,idx) * diag( d(idx) - alpha ) * VT(idx,:);
% (a)ii
Znorm = sum(d(idx)-alpha);
energy_old = energy;
energy = alpha*Znorm + norm(Z(Omega(:))-X(Omega(:)),'fro')/;
if abs(energy - energy_old) / energy_old < epsilon
break
end
end
% check termination condition of alpha
k = length(idx); % rank of Z
disp(['alpha = ' num2str(alpha) '; rank = ' num2str(k) '; number of iteration: ' num2str(innerIts)]);
if k <= maxRank && alpha > 1e-
alpha = alpha*eta;
else
break;
end
end
end