% motion analysis
%from a rendered frame

clear all

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

%=======================================
% 
% load 'c:\my documents\bigdata\mov.mat'
% [nX,nY,junk] = size(mov(1).cdata);
% [junk,nT] = size(mov);
% rfactor = .25;
% for i=1:nT
%     im = imresize(rgb2gray(mov(i).cdata),rfactor,'bicubic');
%     imagesc(im);
%     colormap(gray)
%     drawnow
%     vol(:,:,i) = double(im) ;
% end

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

%vol = smooth3(vol,'box',[1,1,3]) ;
[a,b,c]=size(vol);

%compute the 1/f phase correction
 for i=1:a
     for j=1:b
         rd(i,j) = sqrt( (i-a/2).^2 + (j-b/2).^2 );
         if rd(i,j) == 0
             rd(i,j) = 1;
         end
     end
 end

for i=1:nT
    fvol(:,:,i) = fft2(vol(:,:,i));  
end

fvol = fftshift(fvol,1);
fvol = fftshift(fvol,2);

vol = angle(fvol);
for i=1:nT
    vol(:,:,i) = vol(:,:,i) ./ rd ;
end

%get (spectral strength) .* (frame phase shift)
for i=1:nT-1
    dvol(:,:,i) = (vol(:,:,i+1) - vol(:,:,i)) .* abs(fvol(:,:,i)) ;  
end


%=================================================================

figure(2)
clf
%dvol = smooth3(dvol,'gaussian',[3,3,1]);
r=2;
rvol = reducevolume(dvol,[r,r,1]);

level = max(max(max(rvol)))/200 ;  % 2*mean(mean(mean(dvol)));
p=patch(isosurface(rvol,level));
reducepatch(p,.3)
set(p,'facecolor',[1,1,.5],'edgecolor','none','FaceLighting','phong')

level = -level ;  % 2*mean(mean(mean(dvol)));
p=patch(isosurface(rvol,level));
reducepatch(p,.3)
set(p,'facecolor',[1,.5, 1],'edgecolor','none','FaceLighting','phong')

camlight
light
set(gca,'xlim',[0 b/r],'ylim',[0 a/r],'zlim',[0 c])
daspect([1,1,1])
view(30,30)
zlabel('time')
xlabel('x freq scale')
ylabel('y freq scale')
box on
rotate3d on


%================================================================
figure(3)
clf
%dvol = smooth3(dvol,'box',[1,1,5]);
%subplot(2,1,1)
cx=4;
for i=1:nT-1;
    s(i)=sum(sum(abs(dvol([1:fix(a/2)-cx,fix(a/2)+cx:end],[1:fix(b/2)-cx,fix(b/2)+cx:end],i))));
end
s = medfilt1(s,3);
plot(s,'-o')
ylabel('Total Spectral Motion')




