function P=centers_list_img(image_file_name)
%P - Four column list of (1)Numeric label of individual cell 
%                        (2&3)x,y coordinates of centers of cells
%                        (4)Estimate of size of individual cell - crude
%                        approximation, by counting pixels in binary image.
%
%There are three primary parameters in this program that need modification,
%depending on the images being used. (1)Image file name (2)Crop-region
%(3)Fourier filtering mask parameters. These are indicated when they are
%defined in the program, and are currently defaulted to values suitable to
%the image file 'polydisperse.jpg'
%
%DISPLAY RESULTS ONLY section may be commented out for efficiency.
%For bugs, comments, etc. please contact Kapil Krishan :
%kkrishan@uci.edu
%www.physics.uci.edu/~foams
%version: I


%Image file read as matrix 'a'
%(1)IMAGE FILE NAME
%image_file_name='./imgs/image0.bmp';
a=double(imread(image_file_name));
%Cropping the image..
%Parameters describing region to crop from raw image along y and x
%(2)CROP REGION
Cy=[50 320];Cx=[1 640];
img=a(Cy(1):Cy(2),Cx(1):Cx(2));

%...and fourier filtering image file
f=fftshift(fft2(img));
%Parameters for fourier filtering of image. Elipsoidal masks are used for
%the high and low frequency cuttoffs. 'rad' contains the major and minor
%radii for these filters
%(3)FOURIER FILTERING MASK PARAMETERS
rad=[35 75;15 35];
sz=size(f);
f_msk=mask_ellipse(rad(1,:),sz)-mask_ellipse(rad(2,:),sz);%msk: mask for fourier filter

af=real(ifft2(ifftshift(f.*f_msk)));%'af' is now fourier filtered

%Center of mass of bubbles
%Making binary image (Black and White) from filtered image by simple
%thresholding. The choice of threshold is a free parameter, but should be
%applied consistantly across a number of images.
BW=(af<mean(af(:))-std(af(:)));
%Labeling each connected component with a different integer.
[L,num] = bwlabel(BW,8);

%Locating center of mass and size of each connected component
P=zeros(num,4);
P(:,1)=[1:num]';
for ic=P(:,1)'
    [I,J]=find(L==ic);
    P(ic,[2 3 4])=[mean(I) mean(J) size(I,1)];
end

%DISPLAY RESULTS ONLY
%The following two lines display the primary results of the program. They may
%be commmented out.
figure(1);imagesc(img);axis image;colormap gray;
figure(2);imagesc(L);axis image;colormap([1 1 1;jet]);hold on;plot(P(:,3),P(:,2),'k.','Markersize',10);drawnow;

function msk=mask_ellipse(rad,sz)
%function that makes a binary mask in the shape of an ellipse: Matrix of
%size 'sz', ellipse centered at origin with radii rad(1) and rad(2).

ctr=floor(sz/2+1);
msk=zeros(sz(1),sz(2));
prm=sqrt(2)*pi*norm(rad);
for theta=0:pi/prm:2*pi
    x=rad(1)*cos(theta);
    y=rad(2)*sin(theta);
    msk(ctr(1)+ceil(x),ctr(2)+ceil(y))=1;
end

msk=imfill(msk,'holes');