function digitize_cavendish_movie
%video analysis program for Cavendish torsional pendulum
%this program will process a movie of the pendulum and determine the angular deflection from the position of the
%laser dot.
%The movie gets analyzed and the angular deflection of the laser beam is
%determined. For this the user is required to calibrate the camera in
%pixels_per_cm, and provide the distance between the projected screen and
%mirror.
%
%The user inputs include
%1) name of movie: file_name
%2) which frames to analyze in terms of: start_frame, step and end_frame
%3) calibration of the camera in terms of: pixels_per_cm
%4) distance between the projected screen and mirror: L
%5) The region of interest which contains the laser dot, ROI ( this should
%be small to increase the accuracy )
%6) Threshold value for discriminating the intensity of the laser vs. the
%background.
%
%Output is
%1) a .MAT file containing the angular deflection and time: file_name.MAT
%2) saved figure: filename.FIG
%
%Input paragraphs are delimited using %************************
%
% PS. I suggest you run cell segments of the code using "control+enter"
%sak@wpi.edu 1/10/12
%inputs are here%****************************************
file_name = 'Cavendish2.mov';

vid=VideoReader(file_name);

k=1;
while hasFrame(vid)
    

start_frame = 10;
end_frame = inf;
% end_frame = 500;
step = 10;
threshold = 220;
pixel_per_cm = (608-91)/(2.54*11);
%calibration of the camera. Note that I used the paper width

ROI = [100 260 150 225]; %Region of interest in pixels - the format is
% ROI = [200 250 150 225]; %Region of interest - the format is
%[horizontal_left, horizontal_right, vertical_lowest, vertical_highest]

show_frames_flag = 0;
%set this flag = 1 to show the frame with calculated position of laser dot
%set this flag = 0 for faster processing

L = 57*2.54;
%distance from pendulum to screen
%****************************************
%end of input
%%
obj = VideoReader(file_name);
end_frame = min( [get( obj, 'NumberOfFrames' ), end_frame] );
figure( sum( file_name ) ); set( clf, 'name', file_name );
clf;
to_do = ceil( linspace( start_frame, end_frame, 5 ) );
fprintf( 'loading %i frames', numel(to_do ) );
for i=1:numel(to_do)
    if i==1
        tmp = mean( read( obj, to_do(i) ), 3 );
        im = repmat( tmp, [1 1 numel(to_do)] );
    end
    tmp = mean( read( obj, to_do(i) ),3 );
    im(:,:,i) = tmp;
    fprintf( '.' );
end;
fprintf( 'done \n' )
%%
clf;
sbplt1 = subplot( 1,2,1);
imagesc( mean(im,3) );
colormap( 'gray' );
hold on;
plot( ROI([1 1 2 2 1]), ROI( [3 4 4 3 3] ), 'r' );
axis image;
title( { sprintf( 'superposition of %i frames', numel( to_do ) ); ''; ''; ''} );
xlabel( 'pixels' );
ylabel( 'pixels' );
axes( 'OuterPosition', get( sbplt1, 'OuterPosition' ), 'Position', get( sbplt1, 'Position' ), ...
    'PlotBoxAspectRatio', get( sbplt1, 'PlotBoxAspectRatio' ), ...
    'xaxislocation', 'top', 'yaxislocation', 'right', 'color', 'none', 'DataAspectRatio', [1 1 1] );
xlim( gca, xlim( sbplt1 )/pixel_per_cm );
ylim( gca, ylim( sbplt1 )/pixel_per_cm );
xlabel( 'cm' );
ylabel( 'cm' );
drawnow;
%%
to_do = start_frame:step:end_frame;
xy = single(zeros( numel(to_do), 2 ) );
intensity = nan( numel( to_do ), 1 );
t = [to_do]/get( obj, 'framerate' );
[x, y] = meshgrid( [ROI(1):ROI(2)]/pixel_per_cm, [ROI(3):ROI(4)]/pixel_per_cm );
%%
sbplt2 = subplot(1,2,2); hold on;
clc
fprintf( 'analyzing selected frames \n' );
fprintf( '\n%05i/%05i: %4.2E', 0, end_frame, 0 );
for j=1:numel(to_do)
    im = mean( read( obj, to_do(j) ), 3);
    im = im( ROI(3):ROI(4), ROI(1):ROI(2) );
    im(im<threshold) = 0;
    intensity(j) = sum(sum(im));
    xy(j,:) = [sum( sum( x.*im ) ), sum( sum( y.*im ) )]/intensity(j);
    if show_frames_flag
        plot( sbplt2, t(j), xy(j,1)-xy(1,1), 'or', t(j), xy(j,2)-xy(1,2), 'sk' );
        title( sprintf( '%04i/%04i: %4.2E', j, numel(to_do), intensity(j) ) );
        drawnow;
        if j==1
            ylabel( 'displacement (cm)' );
            xlabel( 't (sec)' );
            legend( {'x', 'y'} );
        end
    else
        fprintf( '\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%05i/%05i: %4.2E', j, numel(to_do), intensity(j) );
    end
end
%%
center = mean( xy );
theta = atan2( xy(:,1)-center(1), L );
cla;
plot( t, theta );
xlabel( 't (sec)' );
ylabel( '\theta (rad)' );
title( 'oscillations' );
%%
if any( intensity < max(intensity)*.3 )
    warning( 'unusual fluctuation in the intensity - check ROI' );
    figure(1); set( clf, 'name', 'bad ROI' );;
    subplot(2,1,1);
    plot( t, xy(:,1), '.' );
    xlabel( 't' );
    ylabel( 'x' );
    subplot( 2,1,2);
    plot( t, intensity );
    xlabel( 't' );
    ylabel( 'intensity' );
end
%%
saveas( gcf, strcat( file_name, '.fig' ) );
%%
fprintf( '\nsaving as a MAT-file to %s \n', strcat( file_name, '.mat' ) );
save( strcat( file_name, '.mat' ), 't', 'xy', 'theta', 'ROI', 'L' );
