function out=pkfnd(im,th,sz)
% finds local maxima in an image to pixel level accuracy.
% this provides a rough guess of particle
% centers to be used by cntrd.m. Inspired by the lmx subroutine of Grier
% and Crocker's feature.pro
% INPUTS:
% im: image to process, particle should be bright spots on dark background with little noise
% ofen an bandpass filtered brightfield image (fbps.m, fflt.m or bpass.m) or a nice
% fluorescent image
% th: the minimum brightness of a pixel that might be local maxima.
% (NOTE: Make it big and the code runs faster
% but you might miss some particles. Make it small and you'll get
% everything and it'll be slow.)
% sz: if your data's noisy, (e.g. a single particle has multiple local
% maxima), then set this optional keyword to a value slightly larger than the diameter of your blob. if
% multiple peaks are found withing a radius of sz/2 then the code will keep
% only the brightest. Also gets rid of all peaks within sz of boundary
%OUTPUT: a N x 2 array containing, [row,column] coordinates of local maxima
% out(:,1) are the x-coordinates of the maxima
% out(:,2) are the y-coordinates of the maxima
%CREATED: Eric R. Dufresne, Yale University, Feb 4 2005
%MODIFIED: ERD, 5/2005, got rid of ind2rc.m to reduce overhead on tip by
% Dan Blair; added sz keyword
% ERD, 6/2005: modified to work with one and zero peaks, removed automatic
% normalization of image
% ERD, 6/2005: due to popular demand, altered output to give x and y
% instead of row and column
% ERD, 8/24/2005: pkfnd now exits politely if there's nothing above
% threshold instead of crashing rudely
% ERD, 6/14/2006: now exits politely if no maxima found
% ERD, 10/5/2006: fixed bug that threw away particles with maxima
% consisting of more than two adjacent points
%find all the pixels above threshold
%im=im./max(max(im));
ind=find(im > th);
[nr,nc]=size(im);
tst=zeros(nr,nc);
n=length(ind);
if n==0
out=[];
display('nothing above threshold');
return;
end
mx=[];
%convert index from find to row and column
rc=[mod(ind,nr),floor(ind/nr)+1];
for i=1:n
r=rc(i,1);c=rc(i,2);
%check each pixel above threshold to see if it's brighter than it's neighbors
% THERE'S GOT TO BE A FASTER WAY OF DOING THIS. I'M CHECKING SOME MULTIPLE TIMES,
% BUT THIS DOESN'T SEEM THAT SLOW COMPARED TO THE OTHER ROUTINES, ANYWAY.
if r>1 & r<nr & c>1 & c<nc
if im(r,c)>=im(r-1,c-1) & im(r,c)>=im(r,c-1) & im(r,c)>=im(r+1,c-1) & ...
im(r,c)>=im(r-1,c) & im(r,c)>=im(r+1,c) & ...
im(r,c)>=im(r-1,c+1) & im(r,c)>=im(r,c+1) & im(r,c)>=im(r+1,c+1)
mx=[mx,[r,c]'];
%tst(ind(i))=im(ind(i));
end
end
end
%out=tst;
mx=mx';
[npks,crap]=size(mx);
%if size is specified, then get ride of pks within size of boundary
if nargin==3 & npks>0
%throw out all pks within sz of boundary;
ind=find(mx(:,1)>sz & mx(:,1)<(nr-sz) & mx(:,2)>sz & mx(:,2)<(nc-sz));
mx=mx(ind,:);
end
%prevent from finding peaks within size of each other
[npks,crap]=size(mx);
if npks > 1
%CREATE AN IMAGE WITH ONLY PEAKS
nmx=npks;
tmp=0.*im;
for i=1:nmx
tmp(mx(i,1),mx(i,2))=im(mx(i,1),mx(i,2));
end
%LOOK IN NEIGHBORHOOD AROUND EACH PEAK, PICK THE BRIGHTEST
for i=1:nmx
roi=tmp( (mx(i,1)-floor(sz/2)):(mx(i,1)+(floor(sz/2)+1)),(mx(i,2)-floor(sz/2)):(mx(i,2)+(floor(sz/2)+1))) ;
[mv,indi]=max(roi);
[mv,indj]=max(mv);
tmp( (mx(i,1)-floor(sz/2)):(mx(i,1)+(floor(sz/2)+1)),(mx(i,2)-floor(sz/2)):(mx(i,2)+(floor(sz/2)+1)))=0;
tmp(mx(i,1)-floor(sz/2)+indi(indj)-1,mx(i,2)-floor(sz/2)+indj-1)=mv;
end
ind=find(tmp>0);
mx=[mod(ind,nr),floor(ind/nr)+1];
end
if size(mx)==[0,0]
out=[];
else
out(:,2)=mx(:,1);
out(:,1)=mx(:,2);
end
end
function [out,out1,out2]=cntrd(im,mx,sz,interactive)
%cntrd: calculates the centroid of bright spots to sub-pixel accuracy.
% Inspired by Grier & Crocker's feature for IDL, but greatly simplified and optimized
% for matlab
% INPUTS:
% im: image to process, particle should be bright spots on dark background with little noise
% ofen an bandpass filtered brightfield image or a nice fluorescent image
%
% mx: locations of local maxima to pixel-level accuracy from pkfnd.m
%
% sz: diamter of the window over which to average to calculate the centroid.
% should be big enough
% to capture the whole particle but not so big that it captures others.
% if initial guess of center (from pkfnd) is far from the centroid, the
% window will need to be larger than the particle size. RECCOMMENDED
% size is the long lengthscale used in bpass plus 2.
%
%
% interactive: OPTIONAL INPUT set this variable to one and it will show you the image used to calculate
% each centroid, the pixel-level peak and the centroid
%
% NOTE:
% - if pkfnd, and cntrd return more then one location per particle then
% you should try to filter your input more carefully. If you still get
% more than one peak for particle, use the optional sz parameter in pkfnd
% - If you want sub-pixel accuracy, you need to have a lot of pixels in your window (sz>>1).
% To check for pixel bias, plot a histogram of the fractional parts of the resulting locations
% - It is HIGHLY recommended to run in interactive mode to adjust the parameters before you
% analyze a bunch of images.
%
%OUTPUT: a N x 3 array containing, x, y and brightness for each feature
% out(:,1) is the x-coordinates
% out(:,2) is the y-coordinates
% out(:,3) is the brightnesses
% out(:,4) is the sqare of the radius of gyration
%
%CREATED: Eric R. Dufresne, Yale University, Feb 4 2005
% 5/2005 inputs diamter instead of radius
% Modifications:
% D.B. (6/05) Added code from imdist/dist to make this stand alone.
% ERD (6/05) Increased frame of reject locations around edge to 1.5*sz
% ERD 6/2005 By popular demand, 1. altered input to be formatted in x,y
% space instead of row, column space 2. added forth column of output,
% rg^2
% ERD 8/05 Outputs had been shifted by [0.5,0.5] pixels. No more!
% ERD 8/24/05 Woops! That last one was a red herring. The real problem
% is the "ringing" from the output of bpass. I fixed bpass (see note),
% and no longer need this kludge. Also, made it quite nice if mx=[];
if nargin==3
interactive=0;
end
if sz/2 == floor(sz/2)
warning('sz must be odd, like bpass');
end
if isempty(mx)
warning('there were no positions inputted into cntrd. check your pkfnd theshold')
out=[];
return;
end
r=(sz+1)/2;
%create mask - window around trial location over which to calculate the centroid
m = 2*r;
x = 0:(m-1) ;
cent = (m-1)/2;
x2 = (x-cent).^2;
dst=zeros(m,m);
for i=1:m
dst(i,:)=sqrt((i-1-cent)^2+x2);
end
ind=find(dst < r);
msk=zeros([2*r,2*r]);
msk(ind)=1.0;
%msk=circshift(msk,[-r,-r]);
dst2=msk.*(dst.^2);
ndst2=sum(sum(dst2));
[nr,nc]=size(im);
%remove all potential locations within distance sz from edges of image
ind=find(mx(:,2) > 1.5*sz & mx(:,2) < nr-1.5*sz);
mx=mx(ind,:);
ind=find(mx(:,1) > 1.5*sz & mx(:,1) < nc-1.5*sz);
mx=mx(ind,:);
[nmx,crap] = size(mx);
%inside of the window, assign an x and y coordinate for each pixel
xl=zeros(2*r,2*r);
for i=1:2*r
xl(i,:)=(1:2*r);
end
yl=xl';
pts=[];
%loop through all of the candidate positions
for i=1:nmx
%create a small working array around each candidate location, and apply the window function
tmp=msk.*im((mx(i,2)-r+1:mx(i,2)+r),(mx(i,1)-r+1:mx(i,1)+r));
%calculate the total brightness
norm=sum(sum(tmp));
%calculate the weigthed average x location
xavg=sum(sum(tmp.*xl))./norm;
%calculate the weighted average y location
yavg=sum(sum(tmp.*yl))./norm;
%calculate the radius of gyration^2
rg=(sum(sum(tmp.*dst2))/ndst2);
%concatenate it up
%pts=[pts,[mx(i,1)+xavg-r,mx(i,2)+yavg-r,norm,rg]']; if u wanna find
%the total brightness and raius of gyration^2.
x_final = mx(i,1)+xavg-r ;
y_final = mx(i,2)+yavg-r ;
pts=[pts,[x_final,y_final]'];
%OPTIONAL plot things up if you're in interactive mode
if interactive==1
imagesc(tmp)
axis image
hold on;
plot(xavg,yavg,'x')
plot(xavg,yavg,'o')
plot(r,r,'.')
hold off
pause
end
end
out=pts';
end
