%
% Copyright (c) 2011, Kai Guo, Prakash Ishwar, and Janusz Konrad. All
% rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1) Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%
% 2) Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
% A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
% HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%

% this function is to calculate the feature covariance matrix given the
% intensity and optical flow of a video segment

function CC=optical_flow_feature(flow,seq_gray)
thresh=5; 
% threshold that roughly tells moving pixels, can be 10, 15 for
% stationary camera. for moving camera, using global mv to compensate
% the frame offset

thresh_cov = 100; 
% threshold that decides if current cov matrix is reliable or not,
% theoretically should be at least 13, but larger num is prefered
% (tunable)

% spatial low-pass filtering
fl=fspecial('gaussian',[3 3],1);

seq=double(seq_gray);
[H,V,T]=size(seq);

for k=1:T
    seq(:,:,k)=imfilter(seq(:,:,k),fl,'same');
end

% compute optical flow gradients ux,uy,ut,vx,vy,vt and intensity temporal
% graident dt
ux=zeros(H,V,T);
uy=zeros(H,V,T);
ut=zeros(H,V,T);
vx=zeros(H,V,T);
vy=zeros(H,V,T);
vt=zeros(H,V,T);
dt=zeros(H,V,T);

u=flow(:,:,:,1);
v=flow(:,:,:,2);
for k=1:T
    [ux(:,:,k),uy(:,:,k)]=gradient(u(:,:,k));
    [vx(:,:,k),vy(:,:,k)]=gradient(v(:,:,k));    
end

for i=1:H
    [ut(i,:,:),temp]=gradient(squeeze(u(i,:,:)));
    [vt(i,:,:),temp]=gradient(squeeze(v(i,:,:)));
end

for i=1:H
    [dt(i,:,:),temp]=gradient(squeeze(seq(i,:,:)));
end

[xcoor,ycoor,tcoor]=meshgrid(1:V,1:H,1:T);

%============compute kinametic features======================

% divergence
f_div=ux+vy;
% vorticity
f_vor=vx-uy;
    
% Gradient tensor feature: Q,S_Q
f_Q=1/2*(f_div.^2-(ux.^2+2*uy.*vx+vy.^2));
f_S_Q=1/2*(f_div.^2-(ux.^2+1/2*(uy+vx).^2+vy.^2));


%========compute covariance matrix based on kinametic features==========

feature=[xcoor(:),ycoor(:),tcoor(:),dt(:),u(:),v(:),ut(:),vt(:),...
         f_div(:),f_vor(:),f_Q(:),f_S_Q(:)];
feature(find(abs(feature(:,4))<=thresh),:)=[];
if size(feature,1)<=thresh_cov
   CC=[];
else
   CC=cov(feature);
end

%========compute covariance matrix based on gradients of optical flow===
%
% feature=[xcoor(:),ycoor(:),tcoor(:),dt(:),u(:),v(:),...
%          ux(:),uy(:),ut(:),vx(:),vy(:),vt(:)];
% feature(find(abs(feature(:,4))<=threshold),:)=[];
% 
% if size(feature,1)<=thresh_cov
%    CC=[];
% else
%    CC=cov(feature);
% end
%

%%%%%%%%
