% motion analysis
%from a rendered frame

clear all

load 'c:\my documents\bigdata\mov.mat'

[nX,nY,junk] = size(mov(1).cdata);
[junk,nT] = size(mov);
rfactor = .25;

figure(1)
clf
for i=1:nT
    im = imresize(rgb2gray(mov(i).cdata),rfactor,'bicubic');
    imagesc(im);
    colormap(gray)
    drawnow
    vol(:,:,i) = double(im) ;
end

%!!!!
%vol = smooth3(vol,'box',[3,3,3]);

[gx,gy,gt]=gradient(vol);

figure(2)
clf
%set(gcf,'doublebuffer','on');
for i=1:nT
    vx(:,:,i) = gt(:,:,i).*gx(:,:,i);
    vy(:,:,i) = gt(:,:,i).*gy(:,:,i);
    speed(:,:,i) = sqrt( (vx(:,:,i).^2 + vy(:,:,i).^2) ./ (gx(:,:,i).^2 + gy(:,:,i).^2) );   
end

speed(isnan(speed) | isinf(speed)) = 0;
speed = speed /(nX*nY);
%speed = smooth3(speed,'box',[1,1,1]);

%while(0)
level = 2*mean(mean(mean(speed)));
p=patch(isosurface(speed,level));
reducepatch(p,.2)
set(p,'facecolor',[1,1,.5],'edgecolor','none','FaceLighting','phong')
camlight
light
set(gca,'xlim',[0 nY*rfactor],'ylim',[0 nX*rfactor],'zlim',[0 nT])
daspect([1,1,1])
view(30,30)
zlabel('time')
xlabel('x')
ylabel('y')
box on
rotate3d on
%end

figure(3)
clf

subplot(2,1,1)
for i=1:nT;
    s(i)=sum(sum(speed(:,:,i)));
end
plot(s)
ylabel('Total Pixel Motion')

subplot(2,1,2)
for i=1:nT
    sx(i)=sum(sum(vx(:,:,i)))/(nX*nY);
    sy(i)=sum(sum(vy(:,:,i)))/(nX*nY);
    ang(i)=(atan(sy(i)./sx(i)));
end
%plot(1:nT,ang)
plot3(1:nT,sx,sy)
box on
%set(gca,'xlim',[0 nT],'ylim',[-10 10],'zlim',[-10 10])
xlabel('time')
ylabel('sx')
zlabel('sy')
%ylabel('Average Angle')
rotate3d on