function [x,y,u,v,SnR]=simplepiv(im1,im2,winsize,Dt,overlap)
%SIMPLEPIV - basic Particle Image Velocimetry.
%
% Image read
if ischar(im1)
    A=readmyimage(im1); B=readmyimage(im2);
else
    A=im1; B=im2;
end
A=double(A); B=double(B);

% Various declarations and initialisations
M=winsize; N=M; cj1=1; ci1=1;
x=zeros(ceil((size(A,1)-winsize)/((1-overlap)*winsize)), ...
    ceil((size(A,2)-winsize)/((1-overlap)*winsize)));
y=x; u=x; v=x; SnR=x; Pkh=x; u2=x;

% Create bias correction matrix
BiCor=xcorrf2(ones(winsize),ones(winsize))/(winsize*winsize);
SA=std(A(:));

tic
% MAIN LOOP
for jj=1:(1-overlap)*winsize:size(A,1)-winsize+1
    for ii=1:(1-overlap)*winsize:size(A,2)-winsize+1
        C=A(jj:jj+winsize-1,ii:ii+winsize-1);
        D=B(jj:jj+winsize-1,ii:ii+winsize-1);
        % Calculate the standard deviation of each of the subwindows
        stad1=std(C(:));stad2=std(D(:));
        % the following test removes some potential breakpoint
        if stad1<0.2*SA, stad1=NaN; end
        if stad2<0.2*SA, stad2=NaN; end
        % Subtract the mean of each window to avoid correlation of
        % the mean background intensities.
        C=C-mean(C(:));
        D=D-mean(D(:));
        
        % Calculate the correlation using the xcorrf2 (found at the
        % Mathworks web site). Divide by N (winsize) and the stad's
        % to normalize the peakheights.
        if isnan(stad1)~=1 && isnan(stad2)~=1
            R=xcorrf2(C,D)./(winsize*winsize*stad1*stad2);
            % Correct for displacement bias
            R=R./BiCor;
            % Locate the highest point
            [y1,x1]=find(R==max(max(R(0.5*winsize+2:1.5*winsize-3,...
                0.5*winsize+2:1.5*winsize-3))));
            if size(x1,1)>1 || size(y1,1)>1
                x1=round(sum(x1.*(1:length(x1))')./sum(x1));
                y1=round(sum(y1.*(1:length(y1))')./sum(y1));
            end
            % Interpolate to find the peak position at subpixel resolution,
            % X0,Y0 now denotes the displacements.
            x01=x1 + ( (log(R(y1,x1-1))-log(R(y1,x1+1)))/( (2*log(R(y1,x1-1)))-(4*log(R(y1,x1)))+(2*log(R(y1,x1+1)))) );
            y01=y1 + ( (log(R(y1-1,x1))-log(R(y1+1,x1)))/( (2*log(R(y1-1,x1)))-(4*log(R(y1,x1)))+(2*log(R(y1+1,x1)))) );
            x0=x01-(M);
            y0=y01-(N);
            
            %estimate signal to noise ratio
            R2=R;
            R2(y1-3:y1+3,x1-3:x1+3)=NaN;
            [p2_y2,p2_x2]=find(R2==max(max(R2( 0.5*N+2:1.5*N-3,0.5*M+2:1.5*M-3))));
            if length(p2_x2)>1
                p2_x2=p2_x2(round(length(p2_x2)/2));
                p2_y2=p2_y2(round(length(p2_y2)/2));
            end
            % Store the data
            x(cj1,ci1)=(winsize/2)+ii-1;
            y(cj1,ci1)=(winsize/2)+jj-1;
            u(cj1,ci1)=-x0/Dt;
            v(cj1,ci1)=-y0/Dt;
            SnR(cj1,ci1)=R(y1,x1)/R2(p2_y2,p2_x2);
        else
            %if the standard deviations were rubbish we jump directly
            %here to reduce calculation time and avoid crashes.
            u(cj1,ci1)=NaN; v(cj1,ci1)=NaN; SnR(cj1,ci1)=NaN; Pkh(cj1,ci1)=NaN; u2(cj1,ci1)=nan;
            x(cj1,ci1)=(ii+(winsize/2)-1);
            y(cj1,ci1)=(jj+(winsize/2)-1);
        end
        % Update counters
        ci1=ci1+1;
        
    end
    % Display calculation time
    fprintf('\r No. of vectors: %d', (cj1-1)*(ci1)+ci1-1 -sum(isnan(u(:))))
    fprintf(' , Seconds taken: %f', toc); %etime(clock,t0));
    ci1=1;
    cj1=cj1+1;
end

fprintf('.\n');

% now we inline the function XCORRF2 to shave off some time.

function c = xcorrf2(a,b)
%  c = xcorrf2(a,b)
%   Two-dimensional cross-correlation using Fourier transforms.
%       XCORRF2(A,B) computes the crosscorrelation of matrices A and B.
%       XCORRF2(A) is the autocorrelation function.
%       This routine is functionally equivalent to xcorr2 but usually faster.
%       See also XCORR2.

%       Author(s): R. Johnson
%       $Revision: 1.0 $  $Date: 1995/11/27 $

[ma,na] = size(a);
% if nargin == 1
%     %       for autocorrelation
%     b = a;
% end
[mb,nb] = size(b);
%       make reverse conjugate of one array
b = conj(b(mb:-1:1,nb:-1:1));

%       use power of 2 transform lengths
mf = 2^nextpow2(ma+mb);
nf = 2^nextpow2(na+nb);
at = fft2(b,mf,nf);
bt = fft2(a,mf,nf);
%       multiply transforms then inverse transform
c = ifft2(at.*bt);
%       make real output for real input
if ~any(any(imag(a))) && ~any(any(imag(b)))
    c = real(c);
end
%  trim to standard size
c(ma+mb:mf,:) = [];
c(:,na+nb:nf) = [];
