function corr_calibrate(imageData, firstFrame, lastFrame, makeMovie)

% imageData: stack of camera frames
% firstFrame, lastFrame: frames image data that are used for tracking
% makeMovie = 1 makes a movie of the tracking result
%
% Uses Matlab function cpcorr for digital image correlation as e.g. described 
% in http://www.mathworks.com/matlabcentral/fileexchange/12413 by C. Eberl. 
% We typically set CORRSIZE = 5  in cpcorr. We typically use movies with a
% resolution of 100x100 pixels.
%
% created: Johannes Seelig, 2009


%--------------------------------------------------------------------------
% define the fields for correlation tracking 
% the center of the image has to coincide with the center of the ball
[size_y, size_x, size_z] = size(imageData);
centerBall_x = round(size_x/2);
centerBall_y = round(size_y/2);

%number of pixels that generate a field for correlation
fieldPixels_x = 10;
fieldPixels_y = 10;

noPoints_x = floor(size_x/fieldPixels_x);
noPoints_y = floor(size_y/fieldPixels_y);

grid_x = noPoints_x-mod(noPoints_x,2) - 2;
grid_y = noPoints_y-mod(noPoints_y,2) - 2;

gridVector_x = [centerBall_x - grid_x/2*fieldPixels_x - fieldPixels_x/2: ...
        fieldPixels_x:centerBall_x + grid_x/2*fieldPixels_x + fieldPixels_x/2];
gridVector_y = [centerBall_y - grid_y/2*fieldPixels_y - fieldPixels_y/2: ...
        fieldPixels_y:centerBall_y + grid_y/2*fieldPixels_y + fieldPixels_y/2];

[X,Y] = meshgrid(gridVector_x, gridVector_y);

X2 = reshape(X,1,[]);
Y2 = reshape(Y,1,[]);
inputPointsIn = [X2;Y2]';
basePointsIn = [X2;Y2]';
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
%lowpass filter
for p = 1:lastFrame - firstFrame
    
	lowCutoff = 0.4;
	h = fir1(32,lowCutoff);
	hh = ftrans2(h);
	imageData(:,:,p) = imfilter(imageData(:,:,p),hh,'replicate');
    
end
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
%correlation for each field for all frames
for p = 1:lastFrame-firstFrame
    
    base = imageData(:,:,p);
    input = imageData(:,:,p+1);
    inputPointsIn2 = cpcorr(inputPointsIn,basePointsIn,input,base);
    deltaPoints = inputPointsIn2-inputPointsIn;
    deltaPointsMat1(:,:,p) = deltaPoints;
    
    %remove the points at the border of the image that can't be tracked
    deltaField = reshape(deltaPoints,length(gridVector_y), length(gridVector_x), 2);
    noBorderPoints = 2;
    deltaField = deltaField(noBorderPoints + 1: length(gridVector_x) - noBorderPoints,:, :);
    deltaField = deltaField(:, noBorderPoints + 1: length(gridVector_y) - noBorderPoints, :);
    deltaPoints = reshape(deltaField,[],2);
    deltaPointsMat(:,:,p) = deltaPoints(:,:);

end


cumsum_x=cumsum(deltaPointsMat(:,1,:),3);
cumsum_x=squeeze(cumsum_x);

cumsum_y=cumsum(deltaPointsMat(:,2,:),3);
cumsum_y=squeeze(cumsum_y);
 
%mean displacement in x and y direction over all fields
meanCsX = mean(cumsum_x, 1);
meanCsY = mean(cumsum_y, 1);
%-------------------------------------------------------------------------- 


%--------------------------------------------------------------------------
%make a movie
if makeMovie == 1
    
    aviFilename = 'ballTracking.avi';
    ballMovieObj = avifile(aviFilename,'FPS', 5, 'Compression', 'None');

    f =  inline('mean2(A)');

    figure, 


    for r = 1:lastFrame - firstFrame

        U = squeeze(deltaPointsMat1(:,1,r));
        U = reshape(U,length(gridVector_y), length(gridVector_x));
        V = squeeze(deltaPointsMat1(:,2,r));
        V = reshape(V,length(gridVector_y), length(gridVector_x));
        imagesc(squeeze(imageData(:,:,r)));
        colormap gray
        hold on
        quiver(X,Y,U,V,'w');
        hold off
        h = gcf;
        ballMovieObj = addframe(ballMovieObj,h);

    end
    

    ballMovieObj = close(ballMovieObj);
    
end
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
%plot tracking results
figure, subplot(2,1,1)
plot(meanCsX)
ylabel('x displacement [pixels]')
axis tight

subplot(2,1,2)
plot(meanCsY)
xlabel('frame number')
ylabel('y displacement [pixels]')
axis tight
%--------------------------------------------------------------------------
