%
% video_total_motion_v2( strFilename, rcCrop, bShow, bWrite )
%
%       Estimation of total motion based on absolute differences between
%       frames, without respect to the direction of motion that causes 
%       differences.
%
%       Test program for video motion estimation and quantification.
%
%   USAGE: arrMotion = video_total_motion_v2( 'Cosmo_torso_no-instruction_640x480.avi', 0, true, true )
%           OR
%          arrMotion = video_total_motion_v2( 'Cosmo_torso_no-instruction_640x480.avi', [100 1 470 256], true, true )
%
%   ARGUMENTS:
%
%       strFilename:    Video file name. 
%                       .avi,.mpg and.wmv recognized on Windows systems.
%                       .avi and .mov on recognized on Mac.
%
%       rcCrop:     	0 for full frame (no cropping)
%                       [ left top right bottom ] Cropping rectangele
%                       boundary with respect to top-left corner at ( 1, 1 ).
%                           Use get_video_frame( 'Cosmo_torso_no-instruction_640x480.avi', 1, true, false );
%                           and data cursor tool to find crop corner coordinates.
%
%       bShow:          Show progress, plot result, display, total motion image.
%
%       bWrite:         Write the total motion array to a .mat file, the
%                       plot to a .fig file, and the total motion to a .jpg
%                       file. The crop bounds are encoded in the file name.
%
%   RETURN VALUES:
%
%       arrMotion:     An array with time courses of estimated motion,
%                      one time course for each of tsMax time-scale.
%
%   HARDCODED:
%
%       nDiffThreshold:     [1,255] Difference threshold for inclusion, 
%                           independently for each 8-bit channel.
%                           It might be better use a single luminance channel
%                           particularly for speed (not yet implemented).
%
%       tsMax:          Number of temporal differences found, 
%                       from nSkip*2^1 to nSkip*2^(tsMax-1)
%
%       nSkip:          Frame index increment.
%                       1 for every frame
%
%   CALLS:          (none)
%
%
% University of Oregon Brain Development Laboratory
% Mark Dow, http://lcni.uoregon.edu/~mark/
% Created     February 2, 2009 
% Modified    March    5, 2009	(optional crop rectangle, total motion
%                                image, show and write options) 
%

function arrMotion = video_total_motion_v2( strFilename, rcCrop, bShow, bWrite )

%%%%%%%%%%%%%%%%%%%%%%%%
% Hardcoded information:

nDiffThreshold  = 20    % [1,255] Difference threshold for inclusion, 
                        % independently for each 8-bit channel.
                        % It might be better use a single luminance channel
                        % particularly for speed (not yet implemented).
tsMax           =  7    % Number of temporal differences found, 
                        % from nSkip*2^1 to nSkip*2^(tsMax-1)
nSkip           =  1    % Frame index increment.
                        % 1 for every frame

%
%%%%%%%%%%%%%%%%%%%%%%%%

strFilestem = [ strFilename( 1 : length(strFilename) - 4 ) ];

videoIn = mmreader( strFilename )

nframes = get( videoIn, 'NumberOfFrames' );

% Get a full frame.
imFrameFull = read( videoIn, 1 );

% Default to no cropping
if rcCrop == 0
    rcCrop = [ 1 1 size(imFrameFull,2) size(imFrameFull,1) ];
else
    % To Do: sanity check on crop boundaries
    
    strFilestem = [ strFilestem '_' num2str( rcCrop(1) ) '-' num2str( rcCrop(2) )  '-' num2str( rcCrop(3) ) '-' num2str( rcCrop(4) ) ];
end

imFrame = imFrameFull( rcCrop(2):rcCrop(4), rcCrop(1):rcCrop(3), : );

% Allocate memory for reference and difference frames.
imFrame0    = zeros( [rcCrop(4)-rcCrop(2)+1 rcCrop(3)-rcCrop(1)+1 3 3 ], class(imFrame) ); 
imDiff      = zeros( [rcCrop(4)-rcCrop(2)+1 rcCrop(3)-rcCrop(1)+1 3 ] );
imSumDiff   = zeros( [rcCrop(4)-rcCrop(2)+1 rcCrop(3)-rcCrop(1)+1 3 ] );

nthDiff = 0
nftDiff( 1 : nframes - 1, 1:2, 1:tsMax-1 ) = 0; % currently only (:, 2, :) are used

for i = 1 : tsMax-1
%     size(imFrame0)
%     size(imFrame)
    imFrame0( :, :, :, i ) = imFrame;
end

nFrame0(1:tsMax-1) = 1;

for k = 1 : nSkip : nframes - 1
    
    nthDiff = nthDiff + 1
    
    imFrameFull	= read( videoIn, k );
    imFrame     = imFrameFull( rcCrop(2):rcCrop(4), rcCrop(1):rcCrop(3), : );
    imFrameNextFull = read( videoIn, k + nSkip );
    imFrameNext     = imFrameNextFull( rcCrop(2):rcCrop(4), rcCrop(1):rcCrop(3), : );

    if k == 1
        
        % Auto shift first frame difference by one pixel for scaling.
        x1Diff = abs( imFrame( :, 1:rcCrop(3)-rcCrop(1)-1, : ) - imFrame( :, 2:rcCrop(3)-rcCrop(1), : ) );
        x1Scale = length( find( x1Diff > nDiffThreshold ) )
%         size( x1Diff )
%         figure
%         imshow(16*x1Diff)
       
        y1Diff = abs( imFrame( 1:rcCrop(4)-rcCrop(2)-1,:, : ) - imFrame( 2:rcCrop(4)-rcCrop(2),:, : ) );
        y1Scale = length( find( y1Diff > nDiffThreshold ) )
%         size( y1Diff )
%         figure
%         imshow(16*y1Diff)
        
    end
    
    imDiff = double( abs( imFrameNext - imFrame ) );
    imSumDiff( find( imDiff > nDiffThreshold ) ) = imSumDiff( find( imDiff > nDiffThreshold ) ) + imDiff( find( imDiff > nDiffThreshold ) );
   
%     nftDiff( nthDiff, 1, 1 ) = length( find( diff > nDiffThreshold +  0 ) );
    nftDiff( nthDiff, 2, 1 ) = length( find( imDiff > nDiffThreshold ) );
%     nftDiff( nthDiff, 3, 1 ) = length( find( imDiff > nDiffThreshold + 20 ) );
%     nftDiff( nthDiff, 4, 1 ) = length( find( imDiff > nDiffThreshold + 30 ) );
    
    
    % Lower temporal frequency sampling.
    for ts = 2:tsMax

        if    ( k - nFrame0(ts-1) ) > nSkip*(2^(ts-1)) %...
         % &&  nftDiff( nthDiff, 2, 1 ) < 1000    % no large single-step difference

            imDiff = abs( imFrameNext - imFrame0( :, :,  :, ts-1 ) );
%             figure
%             imagesc(imDiff)
            nftDiff( nFrame0(ts-1):nthDiff, 2, ts ) = length( find( imDiff > nDiffThreshold + 10 ) )/( k - nFrame0(ts-1) );
            
            % Reset reference image to current image.
            imFrame0( :, :, :, ts-1 ) = imFrameNext;
            nFrame0(ts-1) = k;
       end
        
    end

%     figure
%     hist( double( imDiff(:) ), 50 );
    
end

arrMotion = squeeze( nftDiff( :, 2, : ) );

figure
%bar( squeeze( nftDiff( :, 2, : ) ) )
plot( arrMotion )

if bWrite
    saveas( gcf, [ strFilestem '_total-motion.fig' ] )
end

dsScale = sqrt( x1Scale*x1Scale + y1Scale*y1Scale );
fprintf( [ '\nApproximate one pixel deviation scale: ' num2str( dsScale ) ' \n\n' ] );

imSumDiff = 3*imSumDiff/(max(imSumDiff(:)));
imSumDiff( find( imSumDiff > 1 ) ) = 1;
for jj = 1 : rcCrop(4)-rcCrop(2) 
for ii = 1 : rcCrop(3)-rcCrop(1) 
for cc = 1 : 3 
    if imSumDiff( jj, ii, cc ) > 0
        imSumDiff( jj, ii, cc ) = 1 - ( .8*( 1 - imSumDiff( jj, ii, cc ) ) );
    end
end
end
end

if bShow
    
    figure
    imagesc( imSumDiff )
end

if bWrite

    imwrite( imSumDiff, [ strFilestem '_sumdiff.jpg' ] );
    strFilenameOut = [ strFilestem '_total-motion.mat' ];
    save( strFilenameOut, 'arrMotion' )
end