Commit 811f1ddb authored by Krzysztof Gonciarz's avatar Krzysztof Gonciarz
Browse files

Initial commit

parents
% fitness function for bilinear fit used by moments.m
function f=bilinearfit(x,t,m)
f = 0;
for i=1:length(t),
if t(i)<x(4),
f=f+((m(i)-(x(1)*t(i)+x(2)))^2);
else
d = x(1)*x(4)+x(2)-x(3)*x(4);
f=f+((m(i)-(x(3)*t(i)+d))^2);
end
end
if x(4) >= t(length(t)),
f = 1e8*f;
end
return
%
% CHECKTRACKER
% Check tracker performance by plotting the setected particles and
% trajectories in overlay with the original movie frames.
%
% SYNTAX:
% checkTracker(instub,first,last,infmt,ext,framewise,rescale,resultfile,
% cutlength,numtraj)
%
% INPUT ARGUMENTS:
% instub location and file name stub of input files
% first number of first image to process
% last number of last image to process
% infmt numbering format of input file names (see below)
% ext file extension of input and output
% framewise wait for keypress after every frame?
% rescale normalize images before display?
% resultfile name of the tracker result file containing the
% trajectories and particle detections
% cutlength only trajectories longer than cutlength frames are
% considered
% numtraj the ntraj longest trajectories will be taken
%
% RETURN VALUE:
% none.
%
% FORMATS:
% infmt specifies the file numbering format for input and
% output as follows: if the parameter is equal to an integer number >0
% is specifies the number of digits in a fixed numbering (e.g. infmt=3
% -> 001, 002, ..., 010, 011, ..., 100, 101, ..., 999) if the parameter
% is equal to the string 'free', the numbering is free format (i.e.
% infmt='free' -> 1, 2, 3, ..., 10, 11, ..., 100, 101, ...).
%
% AUTHOR:
% Ivo Sbalzarini, ETHZ, August 15, 2003
%
% VERSION:
% 0.1 rev.B
%
%-------------------------------------------------------------------
% No user-adjustable parameters below this line
%-------------------------------------------------------------------
function checkTracker(instub,first,last,infmt,ext,framewise,rescale,resultfile,cutlength,numtraj)
% argument validity checks
if or(nargin<10, nargin>10),
disp('ERROR: wrong number of input arguments. See "help checkTracker" !');
return;
end;
if length(instub)<2,
disp('ERROR: instub must be a non-empty string !');
return;
end;
if length(resultfile)<2,
disp('ERROR: resultfile must be a non-empty string !');
return;
end;
if or(first<0, length(first)>1),
disp('ERROR: first must be a non-negative integer !');
return;
end;
if cutlength < 1
disp('cutlength must be larger than 0. Reset to 1.');
end;
if numtraj < 1
disp('numtraj must be larger than 0. Reset to 1.');
end;
if or(last<0, length(last)>1),
disp('ERROR: last must be a non-negative integer !');
return;
end;
if last<first,
disp('ERROR: last must be larger than first !');
return;
end;
if or(and(length(round(infmt))>1, ~strcmp(infmt,'free')),round(infmt)<1),
if ~strcmp(infmt,'free'), infmt = round(infmt); end;
disp('ERROR: infmt must be a positive integer or "free"');
return;
end;
if and(framewise ~= 0, framewise ~= 1),
disp('ERROR: framewise must be either 0 or 1 !');
return;
end;
if and(rescale ~= 0, rescale ~= 1),
disp('ERROR: rescale must be either 0 or 1 !');
return;
end;
if sum(strcmp(ext,{'tif','jpg','bmp','png','hdf','pcx','xwd'})) ~= 1,
disp('ERROR: ext must be one of the following: tif,jpg,bmp,png,hdf,pcx,xwd');
return;
end;
if rescale,
% read files and determine global extrema
[imin,imax] = getGlobalExtrema(instub,first,last,infmt,ext);
end;
data = loadTraj(resultfile,cutlength);
Ntraj = length(data);
if Ntraj < 1,
disp('No trajectories read.')
return;
end;
% sort trajectories by length
trajlen = zeros(Ntraj,1);
for itraj=1:Ntraj,
trajlen(itraj) = size(data{itraj},1);
end
[lensort,sortidx] = sort(trajlen); % sort on ascending order
figure(10)
hist(lensort,Ntraj)
title('Histogram of trajectory lengths')
if numtraj > Ntraj,
numtraj = Ntraj;
end;
disp(sprintf('%d longest trajectories will be included in analysis',numtraj))
disp(' ');
% save trajectories in files
for it=1:numtraj,
% get it-longest trajectory
t = sortidx(Ntraj-it+1);
% save trajectory to file
outfile = sprintf('SPT_onlyone_%3.3d.xy',it);
pdet = data{t};
save(outfile,'pdet','-ASCII')
end
disp(sprintf('Extracted %d trajectories to xy files',numtraj))
for img=round(first):round(last),
file = getFileName(img,infmt,instub,ext);
a = double(imread(file));
if rescale,
b = (a-imin)/(imax-imin);
else
b = a;
end;
figure(1)
clf
imshow(b);
% find and plot all particle detections and trajectories in this frame
for it=1:numtraj,
% get it-longest trajectory
t = sortidx(Ntraj-it+1);
% current particle detections
pdet = data{t}(find(data{t}(:,1)==img),2:3);
figure(1)
hold on
plot(pdet(:,2)+1,pdet(:,1)+1,'ro')
% trajectories
if data{t}(length(data{t}(:,1)),1) < img,
% trajectory t is completed by now => add to plot in green
plot(data{t}(:,3)+1,data{t}(:,2)+1,'g-');
elseif and(data{t}(1,1)<=img, data{t}(length(data{t}(:,1)),1)>=img),
% trajectory is active => add to plot using red lines
bp = find(data{t}(:,1)==img);
plot(data{t}(1:bp,3)+1,data{t}(1:bp,2)+1,'r-');
end;
end;
figure(1)
hold off
if framewise,
disp(sprintf('Frame %d of %d. Press any key to proceed...',img,last))
pause;
end;
end;
%======================================================================
%
% DETECT_PARTICLES: detect particle-shaped features in frame images
%
% SYNTAX: peak = detect_particles(orig,w,cutoff,pth,v)
%
% INPUTS: orig original image to detect features in
% w global size parameter, >particle radius and
% <interparticle spacing
% cutoff probability cutoff for non-particle discrimination
% pth percentile threshold for maxima selection
% v visualization parameters [viz, nfig] as:
% viz=1 if intermediate visualization is needed
% nfig: figure number for first image
%
% After detection, a 1-cell list of a matrix is returned in
% "peak". peak{1} contains the particle information
% for the present frame stored in a matrix:
%
% peak{1}(:,1) x (col)-positions of particles
% peak{1}(:,2) y (row)-positions of particles
% peak{1}(:,3) zero order intensity moments
% peak{1}(:,4) second order intensity moments
% peak{1}(:,5) *** unused ***
% peak{1}(:,6) *** empty *** to be filled by linker
%
% USES: saturate.m
%
% Ivo Sbalzarini, 12.2.2003
% Institute of Computational Science, Swiss Federal
% Institute of Technology (ETH) Zurich.
% E-mail: sbalzarini@inf.ethz.ch
%
% based on an algorithm by Crocker & Grier:
% Crocker, J.C. & Grier, D.G., Methods of digital video microscopy
% for colloidal studies, J. colloid interface sci., 1996,
% 179: 298-310.
%======================================================================
function peak = detect_particles(orig,w,cutoff,pth,v)
viz = v(1);
nfig = v(2);
% correlation length of camera noise (usu. set to unity)
lambdan = 1;
% some often used quantities
idx = [-w:1:w]; % index vector
dm = 2*w+1; % diameter
im = repmat(idx',1,dm);
jm = repmat(idx,dm,1);
imjm2 = im.^2+jm.^2;
siz = size(orig); % image size
%======================================================================
% STEP 1: Image restoration
%======================================================================
% build kernel K for background extraction and noise removal
% (eq. [4])
B = sum(exp(-(idx.^2/(4*lambdan^2))));
B = B^2;
K0 = 1/B*sum(exp(-(idx.^2/(2*lambdan^2))))^2-(B/(dm^2));
K = (exp(-(imjm2/(4*lambdan^2)))/B-(1/(dm^2)))/K0;
% apply convolution filter
filtered = conv2(orig,K,'same');
if viz == 1,
figure(nfig)
nfig = nfig + 1;
imshow(orig)
title('original image')
figure(21);
mfig = nfig + 1;
imshow(filtered)
title('after convolution filter')
end;
%======================================================================
% STEP 2: Locating particles
%======================================================================
% determining upper pth-th percentile of intensity values
pth = 0.01*pth;
[cnts,bins] = imhist(filtered);
l = length(cnts);
k = 1;
while sum(cnts(l-k:l))/sum(cnts) < pth,
k = k + 1;
end;
thresh = bins(l-k+1);
% generate circular mask of radius w
mask = zeros(dm,dm);
mask(find(imjm2 <= w*w)) = 1;
% identify individual particles as local maxima in a
% w-neighborhood that are larger than thresh
dil = imdilate(filtered,mask);
[Rp,Cp] = find((dil-filtered)==0);
particles = zeros(siz);
V = find(filtered(sub2ind(siz,Rp,Cp))>thresh);
R = Rp(V);
C = Cp(V);
particles(sub2ind(siz,R,C)) = 1;
npart = length(R);
if viz == 1,
figure(nfig)
nfig = nfig + 1;
imshow(particles)
title('intensity maxima of particles');
end;
%======================================================================
% STEP 3: Refining location estimates
%======================================================================
% zero and second order intensity moments of all particles
m0 = zeros(npart,1);
m2 = zeros(npart,1);
% for each particle: compute zero and second order moments
% and position corrections epsx, epsy
for ipart=1:npart,
epsx = 1; epsy = 1;
while or(abs(epsx)>0.5,abs(epsy)>0.5),
% lower and upper index bounds for all particle neighborhoods
% in local coordinates. Recalculate after every change in R,C
li = 1-(R-w-saturate(R-w,1,siz(1)));
lj = 1-(C-w-saturate(C-w,1,siz(2)));
ui = dm-(R+w-saturate(R+w,1,siz(1)));
uj = dm-(C+w-saturate(C+w,1,siz(2)));
% masked image part containing the particle
Aij = filtered(R(ipart)+li(ipart)-w-1:R(ipart)+ui(ipart)-w-1,...
C(ipart)+lj(ipart)-w-1:C(ipart)+uj(ipart)-w-1).* ...
mask(li(ipart):ui(ipart),lj(ipart):uj(ipart));
% moments
m0(ipart) = sum(sum(Aij)); % eq. [6]
% eq. [7]
m2(ipart) = sum(sum(imjm2(li(ipart):ui(ipart),lj(ipart):uj(ipart))...
.*Aij))/m0(ipart);
% position correction
epsx = sum(sum(im(li(ipart):ui(ipart),lj(ipart):uj(ipart))...
.*Aij))/m0(ipart);
epsy = sum(idx(lj(ipart):uj(ipart)).*sum(Aij))/m0(ipart);
% if correction is > 0.5, move candidate location
if abs(epsx)>0.5,
R(ipart) = R(ipart)+sign(epsx);
end;
if abs(epsy)>0.5,
C(ipart) = C(ipart)+sign(epsy);
end;
end;
% correct positions (eq. [5])
R(ipart) = R(ipart)+epsx;
C(ipart) = C(ipart)+epsy;
end;
%======================================================================
% STEP 4: Non-particle discrimination
%======================================================================
sigx = 0.1;
sigy = 0.1;
prob = zeros(size(m0));
Nm = length(m0);
for i=1:Nm,
prob(i)=sum(exp(-((m0(i)-m0).^2./(2*sigx))-((m2(i)-m2).^2./...
(2*sigy)))/(2*pi*sigx*sigy*Nm));
end;
if viz == 1,
figure(20)
clf
figure(nfig)
nfig = nfig + 1;
subplot(2,2,1)
hold on
m0in = m0(find(prob >= cutoff));
m2in = m2(find(prob >= cutoff));
plot(m0in,m2in,'go')
m0in = m0(find(prob < cutoff));
m2in = m2(find(prob < cutoff));
plot(m0in,m2in,'ro')
hold off
xlabel('m0')
ylabel('m2')
subplot(2,2,2)
hist(m0,50)
xlabel('m0')
subplot(2,2,3)
hist(m2,50)
xlabel('m2')
end;
% indices of valid particles
tmp = find(prob>=cutoff);
% pack data into return value
npart = length(tmp);
peak = zeros(npart,6);
peak(:,2) = R(tmp); % row position
peak(:,1) = C(tmp); % col position
peak(:,3) = m0(tmp); % zero order moment
peak(:,4) = m2(tmp); % second order moment
% field 5: unused
% field 6: used by linker to store linked list indices
%======================================================================
% STEP 5: Visualization
%======================================================================
if viz == 1,
% plot crosses at particle positions
C = peak(:,1);
R = peak(:,2);
X = [[C'-2; C'+2], [C'; C']];
Y = [[R'; R'], [R'-2; R'+2]];
figure(nfig)
nfig = nfig + 1
imshow(orig)
hold on
hand = line(X,Y);
set(hand(:),'Color',[1 0 0]);
set(hand(:),'LineWidth',[1.1]);
hold off
end;
peak = {peak};
return
%
% fname = getFileName(index,fmt,stub,ext)
%
% Construct formated file name.
%
% inputs: index index of image
% fmt numbering format in file name (integer giving the
% number of fixed digit positions or 'free' for
% floating length numbering)
% stub file name stub
% ext the file extension (without dot)
%
% outputs: fname the formated file name as a string
%
% Ivo Sbalzarini, July 29, 2003
%
function fname = getFileName(index,fmt,stub,ext)
if length(ext) > 0,
ext = ['.', ext];
end;
if fmt=='free',
if index < 10,
fname = sprintf('%s%1.1d%s',stub,index,ext);
elseif index < 100,
fname = sprintf('%s%2.2d%s',stub,index,ext);
elseif index < 1000,
fname = sprintf('%s%3.3d%s',stub,index,ext);
elseif index < 10000,
fname = sprintf('%s%4.4d%s',stub,index,ext);
elseif index < 100000,
fname = sprintf('%s%5.5d%s',stub,index,ext);
else
fname = sprintf('%s%6.6d%s',stub,index,ext);
end;
else
fname = sprintf(strcat('%s%',sprintf('%d',fmt),'.',sprintf('%d',fmt),'d%s'),stub,index,ext);
end;
return;
% end; function getFileName
%
% [imin,imax] = getGlobalExtrema(instub,first,last,infmt,ext)
%
% determines the global intensity extrema in a series of images.
% For multichannel images (color) the extrema are taken over all
% channels.
%
% inputs: instub file stub for the input files
% first first file index to be read
% last last file index to be read
% infmt numbering format of input files (an integer giving
% the number of fixed digits or 'free' for floating length
% numbering)
% ext the file extension (without dot)
%
% outputs: [l,u] with l being the lowest intensity value found and u
% the highest
%
% Ivo Sbalzarini, February 10, 2004
%
function [imin,imax] = getGlobalExtrema(instub,first,last,infmt,ext)
imax = 0;
imin = Inf;
disp('Scaning images for global extrema. Please wait...');
for img=round(first):round(last),
file = getFileName(img,infmt,instub,ext);
a = double(imread(file));
u = max(max(max(a)));
l = min(min(min(a)));
if u>imax, imax=u; end;
if l<imin, imin=l; end;
end;
return;
%======================================================================
%
% LINK_TRAJECTORIES: links particle positions into trajectories
%
% SYNTAX: peaks = link_trajectories(peaks, L, viz, nfig)
%
% INPUTS: peaks peaks cell list of particle positions in all
% frames and zero and second order intensity
% moments as produced by detect_particles
% L maximum allowed particle displacement between
% two subsequent frames (cutoff)
% viz =1 if intermediate visualization is needed
% nfig figure number for first viz image
%
% The particle matching is done such that the sum of squared
% displacements between two frames in minimized and also the
% quadratic differences in the zero and second order intensity
% moments of the particles are minimizes (more precise: between
% the particle properties peaks{i}(:,3) and peaks{i}(:,4)).
%
% After linkage, the same cell list "peaks" is returned but with
% the fields peaks{i}(:,6) now filled. All other fields are left
% untouched. peaks{i}(j,6) contains the index of the one
% particle in frame i+1 corresponding to particle j in frame i.
% If no correspondence was found (i.e. the trajectory terminates)
% -1 is given.
%
% Ivo Sbalzarini, 12.2.2003
% Institute of Computational Science, Swiss Federal
% Institute of Technology (ETH) Zurich.
% E-mail: sbalzarini@inf.ethz.ch
%
% Based on the logistic transportation algorithm described in:
% http://www.damtp.cam.ac.uk/user/fdl/people/sd/digimage/...
% document/track2d/matching.htm#L_2_PARTICLE_MATCHING
%======================================================================
%======================================================================
% STEP 5: Linking locations into trajectories
%======================================================================
function peaks = link_trajectories(peaks, L, viz, nfig)
nframe = length(peaks);
for iframe = 2:nframe,
% initialize all linked lists to -1
peaks{iframe-1}(:,6) = -1;
disp(sprintf('Linking paths in frame %d of %d',iframe,nframe))
m = length(peaks{iframe-1}(:,1)); % previous frame: p_i, i=1,..,m
n = length(peaks{iframe}(:,1)); % this frame: q_j, j=1,..,n
% empty association matrix
A = zeros(m+1,n+1); % a_{ij}=1 <=> p_i = q_j
% dummy particles a_{(m+1)j},
% a_{i(n+1)}
% delta(i,j): quadratic distance between p_i and q_j
xrep = repmat(peaks{iframe}(:,1)',m,1);
yrep = repmat(peaks{iframe}(:,2)',m,1);
xdiff = xrep - repmat(peaks{iframe-1}(:,1),1,n);
ydiff = yrep - repmat(peaks{iframe-1}(:,2),1,n);
delta = (xdiff.^2)+(ydiff.^2);
% dm0(i,j): quadratic difference between m0 moments of p_i and q_j
xrep = repmat(peaks{iframe}(:,3)',m,1);
xdiff = xrep - repmat(peaks{iframe-1}(:,3),1,n);
dm0 = xdiff.^2;
% dm2(i,j): quadratic difference between m2 moments of p_i and q_j
xrep = repmat(peaks{iframe}(:,4)',m,1);
xdiff = xrep - repmat(peaks{iframe-1}(:,4),1,n);
dm2 = xdiff.^2;
% C(i,j): cost function for link p_i, q_j
C = L*L*ones(m+1,n+1); % set broken dummy links to L^2
C(1:m,1:n) = (delta+dm0+dm2);
% set cost of matchings that will never occur to Inf
C1 = C(1:m,1:n) - repmat(C(m+1,1:n),m,1);
C2 = C(1:m,1:n) - repmat(C(1:m,n+1),1,n);
set1 = find(C1 > 0);
set2 = find(C2 > 0);
s = union(set1,set2);
[i,j] = ind2sub([m n],s);
C(sub2ind(size(C),i,j)) = Inf;
% initialize link matrix A
for i=1:m,
% sort costs of real particles
[srtcst,srtidx] = sort(C(i,:));
% append index of dummy particle
iidx = 1;
dumidx = find(srtidx==(n+1));
% search for available particle of smallest cost or dummy
while and(sum(A(:,srtidx(iidx)))~=0, iidx<dumidx),