%This program loads in a sequence of images, determines particle locations
%and then makes links between the particles in each succeeding images.
%The linked list of particle locations are output in an appropriately named
%data file, which can then be used to determine velocities at different
%positions in the flow field.  We begin by inputing the name of the images
%in the image sequence.  The base name is assumed to be followed by a space,
%a two digit number, and by .jpg, the default output of iMovie.  The reduced
%data is saved in the file basenamelinks.mat, where basename is the input text field.

basename=input('Enter the base filename: ','s');

%There are a number of variables which may need to be changed for different
%series of images.  First, we have:
npics=92; %We operate on npics images.

%Note that this must be less than or equal to the number of images stored,
%and also that it needs to be less than 99 - if you have stored more than
%99 frames, the naming convention for imovie changes.

%Second, we have a number of image processing parameters.  These are:

%A. We use a convolution matrix to clean up the image.  This trims "p"
%pixels off of the image on each side.  Thus:
p=5;

%B. We specify a threshold used for removing background intensity, as
%well as a minimum intensity for registering a particle.  Thus:
thresh=15;

%C. We multiply the convolved image by some gain factor to improve display,
%and which also feeds back into the threshold value chosen.  Thus:
gain=35;

%D. We set a maximum number of particles we will look for in each image.
%Note that we get the brightest ones first.
maxpart=400;

%E. We set the maximum distance between particles in successive images for
%making a link.  Because the velocity in the x-direction is greater (usually)
%than that in the y-direction, we contract the displacement in the x-direction
%by some factor.  Thus:
xfactor=3;
maxdist=2.5; %max distance in pixels

%F. We set the minimum number of frames that we follow a particle before we
%keep the link in the final output:
minframes=min(npics/2-.1,7.5); %e.g., the lesser of 8 or half the total.

%G. We set the maximum number of frames that we follow a particle:
maxframes=15; %This means that we will break up our linked lists to get more points.


%OK, now we're ready to do the image analysis.

npt=zeros(1,npics); %The number of points found in each image.
ptall=zeros(maxpart,2,npics); %Where we will store all particle locations

for i=1:npics; %We work on each image at a time.
  if i<10
   fname=[basename,' 0',num2str(i),'.jpg'];
  else
   fname=[basename,' ',num2str(i),'.jpg'];
  end

  %Now we input the image we wish to operate on
  a=imread(fname,'jpg');
  at=mean(double(a),3); %We convert to double precision, and then gray scale
  figure(1)
  image(uint8(at));colormap('gray')
  title(['The original grayscaled image ',fname])
  drawnow

  %Next we want to filter the image to get a better image quality.  To do this
  %we construct a matrix which emphasizes the individual points and
  %subtracts off the local background:

  %We convolve with a 2p+1 square matrix.
  pt=-ones(2*p+1,2*p+1)/(2*p+1)/(2*p+1); %We normalize the matrix
  pt(p+1,p+1)=1;
 
  %We do the convolution:
  pm=max(conv2(at,pt,'valid')-thresh,0)*gain;

  %Note that the size of the new image is smaller than the old.  This is
  %because we are only using the "valid" part of the convolution, which
  %removes p pixels from all four sides.  You have to account for this in
  %determining the true x-y location after convolution - basically, just add
  %p pixels to each value!

  figure(2)
  image(uint8(pm));colormap('gray')
  title(['The filtered image ',fname])
  drawnow

  %Finally, we identify the particle positions.
  %We find the location of the maximum point:
  [temp,xptarray]=max(pm);
  [val,ypt]=max(temp);
  xpt=xptarray(ypt);

  %We can improve the accuracy of our estimated position by fitting a
  %parabola to the data around the maximum:
  [nw,mw]=size(pm);
  yl=[-1,-1,-1,0,0,0,1,1,1]'+min(max(2,ypt),mw-1);
  xl=[-1,0,1,-1,0,1,-1,0,1]'+min(max(2,xpt),nw-1);
  amax=[ones(9,1),yl,yl.^2,xl,xl.^2];
  bmax=diag(pm(xl,yl));
  af=amax\bmax;
  ypf=-af(2)/2/af(3);
  xpf=-af(4)/2/af(5);
  delt=1; %We shift by no more than one pixel from the maximum
  ypf=min(max(ypt-delt,ypf),ypt+delt);
  xpf=min(max(xpt-delt,xpf),xpt+delt);

  npt(i)=0;
  ptcoord=[0,0]; %We keep the coordinate values in this array.
  while val>thresh & npt(i)<maxpart
    npt(i)=npt(i)+1;
    ptcoord=[ptcoord;[xpf,ypf]];
    %We blank out a 2k+1 square array of values around the maximum:
    k=4;

    pm(max(1,xpt-k):min(nw,xpt+k),max(1,ypt-k):min(mw,ypt+k))=0;

    %We find the location of the new maximum point in the window:
    [temp,xptarray]=max(pm);
    [val,ypt]=max(temp);
    xpt=xptarray(ypt);
      
    %Again we do a parabolic interpolation
    yl=[-1,-1,-1,0,0,0,1,1,1]'+min(max(2,ypt),mw-1);
    xl=[-1,0,1,-1,0,1,-1,0,1]'+min(max(2,xpt),nw-1);
    amax=[ones(9,1),yl,yl.^2,xl,xl.^2];
    bmax=diag(pm(xl,yl));
    af=amax\bmax;
    ypf=-af(2)/2/af(3);
    xpf=-af(4)/2/af(5);
    delt=1; %We shift by no more than one pixel from the maximum
    ypf=min(max(ypt-delt,ypf),ypt+delt);
    xpf=min(max(xpt-delt,xpf),xpt+delt);

  end
  figure(2)
  ptcoord=ptcoord(2:npt(i)+1,:); %We trim off the junk first row.
  hold on
  plot(ptcoord(:,2),ptcoord(:,1),'.')
  hold off
  drawnow

  %OK, we need to keep this list of particles.  We use the 3-D array ptall:
  ptall(1:npt(i),:,i)=ptcoord;

  %And now it's on to the next image!

end

%Now we can plot up all the particle positions on one graph:

%First we set up a list of colors and symbols we want to use:
cols='bgrcm';
cols=[cols,cols,cols,cols,cols,cols]; %That gives us 30
cols=[cols,cols,cols,cols,cols,cols]; %That gives us 30*6
syms='ox+*sdv^<>ph';
syms=[syms,syms,syms]; %Which gives us 39 symbols
syms=[syms,syms,syms]; %Which gives us 39*3 symbols

labels=['image ',num2str(1)]; %We create an array of labels
for i=1:npics
  ptcoord=ptall(1:npt(i),:,i);
  figure(3)
  if i>1.5
    hold on
    labels=str2mat(labels,['image ',num2str(i)]);
  end
  plot(ptcoord(:,2),ptcoord(:,1),[cols(i),syms(i)])
  hold off
end
title(['Particle Position Identifications for ',basename])
xlabel('y')
ylabel('-x')
% legend(labels,'Location','eastoutside'); %This has been removed for
% clarity due to changes in matlab

%OK, now we want to determine the linkages between the identified particles
%in each image.  In the short time between images, particles don't move more
%than a couple of pixels in the x-direction, and less than that in the
%y-direction.  We shall set up a couple of arrays containing the link 
%information. The arrays will be a nparticles x nlinks each, one containing
%the x-positions and one containing the y-positions.

xlinks=zeros(npics*maxpart,npics); %We predimension the arrays to the maximum possible values
ylinks=zeros(npics*maxpart,npics);
nlinks=zeros(npics*maxpart,1);

i=1;
lastpts=ptall(1:npt(i),:,i);
xlinks(1:npt(i),1)=lastpts(:,1);
ylinks(1:npt(i),1)=lastpts(:,2);
nlinks(1:npt(i))=ones(npt(i),1);
lastident=[1:npt(i)]';
nextident=npt(i)+1;

%OK, now for linking with the succeeding images:
for i=2:npics
  currentpts=ptall(1:npt(i),:,i);
  currentident=zeros(npt(i),1);
  mind=zeros(npt(i),1); %We keep an array of minimum distances.
  for j=1:npt(i)
    ndist=(((currentpts(j,1)-lastpts(:,1))/xfactor).^2+(currentpts(j,2)-lastpts(:,2)).^2).^.5;
    %Note that we underweight x variation by a factor of xfactor.
    [mind(j) pt]=min(ndist);
    if mind(j) < maxdist %OK, we're close enough to make a link.
      %First we check and see if some other particle has already been linked into the target
      k=find(currentident==lastident(pt));
      if k %This will only be true if lastident(pt) is already used up.
        if mind(k)>mind(j)  %We want the new link!
          currentident(j)=lastident(pt);
          currentident(k)=nextident; %We fix up the old one
          nlinks(currentident(k))=1;
          xlinks(currentident(k),1)=currentpts(k,1);
          ylinks(currentident(k),1)=currentpts(k,2);
          nextident=nextident+1;          
        else  %We want the old one, so we start a new chain
          currentident(j)=nextident;
          nextident=nextident+1;
        end
      else  
        currentident(j)=lastident(pt); %We make the link.
      end
    else
      currentident(j)=nextident; %We start a new chain of links
      nextident=nextident+1;
    end
    nlinks(currentident(j))=nlinks(currentident(j))+1;
    xlinks(currentident(j),nlinks(currentident(j)))=currentpts(j,1);
    ylinks(currentident(j),nlinks(currentident(j)))=currentpts(j,2);
  end

  %Now we prepare to skip to the next image:
  lastident=currentident;
  lastpts=currentpts;
end

%We wish to break up really long linked series so that particles are
%tracked only for a short period of time (e.g., that they don't move all
%that much over each period).  Thus, we look for the longer linked lists:
xlinks=xlinks(:,1:npics);
ylinks=ylinks(:,1:npics);
k=find(nlinks>maxframes);
while k %We do it recursively
    xkeep=xlinks(k,:);
    ykeep=ylinks(k,:);
    nkeep=nlinks(k);
    nlinks(k)=maxframes;
    xlinks(k,:)=[xlinks(k,1:maxframes),zeros(length(k),npics-maxframes)];
    ylinks(k,:)=[ylinks(k,1:maxframes),zeros(length(k),npics-maxframes)];
    newlinks=nkeep-maxframes; %These are the lengths of the leftovers
    newxlinks=[xkeep(:,maxframes+1:npics),zeros(length(k),maxframes)];
    newylinks=[ykeep(:,maxframes+1:npics),zeros(length(k),maxframes)];
    nlinks=[nlinks;newlinks];
    xlinks=[xlinks;newxlinks];
    ylinks=[ylinks;newylinks];
    
    k=find(nlinks>maxframes);
end

%We wish to keep linked series such that particles are linked through, say, at 
%least 8 frames or so.
k=find(nlinks>minframes);
nlinks=nlinks(k);
xlinks=xlinks(k,:);
ylinks=ylinks(k,:);

%Finally, we want to trim out false identifications.  We can do this by
%examining the residual - see if the velocities are constant for each set
%of linked particle identifications.
resid=zeros(size(nlinks));
for i=1:length(nlinks);
    a=[ones(nlinks(i),1),[1:nlinks(i)]']; %The matrix for fitting the velocity
    ufit=a\xlinks(i,1:nlinks(i))';
    residu=norm(xlinks(i,1:nlinks(i))'-a*ufit)/(nlinks(i)-2)^.5;
    vfit=a\ylinks(i,1:nlinks(i))';
    residv=norm(ylinks(i,1:nlinks(i))'-a*vfit)/(nlinks(i)-2)^.5;
    resid(i)=norm([residu,residv]);
end

%We get rid of links for which the average deviation is greater than 1.2
%pixels:
k=find(resid<1.2);
nlinks=nlinks(k);
xlinks=xlinks(k,:);
ylinks=ylinks(k,:);

filename=[basename,'links.mat'];
save(filename, 'nlinks', 'xlinks', 'ylinks', 'basename')

%We can add connecting lines to the plot of particle identifications too:
figure(3)
hold on
for k=1:length(nlinks)
 plot(ylinks(k,1:nlinks(k)),xlinks(k,1:nlinks(k)))
end
hold off
zoom on %The zoom command allows you to inspect things close up by clicking on them.