% 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 = .5;

figure(1)
clf
set(gcf,'doublebuffer','on');



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

%======================================
while(0)
set(gcf,'position',[30 500 60 80])
axis([0 1 0 1]);
nT = 60;
rsize=.25;
for i=1:nT
    cla
    rectangle('position',[.2 .25+.2*cos(4*pi*i/nT) rsize rsize], 'facecolor','k');
    axis off
    im = getframe(gcf);
    im = rgb2gray(im.cdata);
    vol(:,:,i) = double(im) ;
end
end
%======================================

vol = smooth3(vol,'gaussian',[1,1,3]) ;
slicenum = 1;

[a,b,c]=size(vol);
resRange = fix(log2(min(a,b)))-1;


for res=0:resRange
    
    vol = smooth3(vol,'box',[3,3,1]) ;
    [gx,gy,gt]=gradient(vol) ;
    
    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 ;
    
    for i=1:nT;
       %avgSpeed = sum(sum(speed(:,:,i)))/(nX*nY/4^(slicenum-1));
        s(i,slicenum)=sum(sum(speed(:,:,i)))/(nX*nY/4^(slicenum-1)) ;
        imagesc(vol(:,:,i))
        axis off
        drawnow
    end
    slicenum = slicenum + 1 ;
    vol = reducevolume(vol,[2,2,1]);
    clear gx gy gt vx vy speed;
    
end

figure(2)
clf
h = surfl(s)

%p = patch(surf2patch(h));
%delete(h);
%light
%set(p,'edgecolor','none')
%set(p,'facecolor',[.6 .6 .3])

set(gca,'xtick',1:resRange);
tickstr = min(a,b) * 2.^(-(1:resRange)+1);
set(gca, 'xticklabel',tickstr')
zlabel('Total Pixel Motion')
ylabel('Time (frame #)');
xlabel('resolution (pixels)')
colormap('pink')
rotate3d on

