%Need Matlab version 9 or higher

%Convert the measurement of the reference point printed on the optic to
%meters and plug this number in as the variable "lengthmeters." 

%Plug the number of pixels across the reference point measured using
%the code "resolution.m" in as the variable "lengthpixels."

lengthpixels=2922; %length of reference point in pixels 
lengthmeters=0.0254; %length of reference point in meters  
diameterofoptic=0.0762; %diameter of optic in meters 

ratioperim= (lengthpixels/lengthmeters);  


%read image
I = imread('000200 (teleconverter).jpg');
%display image
imshow(I)

%measures the background illumination
background = imopen(I,strel('disk',15));

%Subtracts background illumination from original image to get rid of glare
I2 = I - background;
imshow(I2)


%measures image dimensions
[x,y,z]=size(I2);

% converts image to grayscale
I3 = rgb2gray(I2);
imshow(I3)


% crop circle
imageSize = size(I3);
ci = [((x+0)/2),((y+0)/2),((ratioperim*diameterofoptic)/2)];%edit x+0 and y+0 to move origin coordinates of crop circle a certain number of pixels
[xx,yy] = ndgrid((1:imageSize(1))-ci(1),(1:imageSize(2))-ci(2));
mask = uint8((xx.^2 + yy.^2)<ci(3)^2);
croppedImage = uint8(zeros(size(I3)));
croppedImage(:,:,1) = I3(:,:,1).*mask;
imshow(croppedImage);

%convert grayscale image to Binary
BW = im2bw(croppedImage, .14); 
figure
imshow(BW);


%Get rid of any remaining glare by eliminating particles over a desired size
lower = (ratioperim*(.00001)) %lower limit of 10 microns
Z = round(lower); %rounds the lower limit number to the nearest integer 
LB = Z;
upper = (ratioperim*(.02)) %upper limit of 20,000 microns
Y = round(upper); %rounds the upper limit number to the nearest integer
UB = Y;
Iout = xor(bwareaopen(BW,LB),  bwareaopen(BW,UB));
figure, imshow(Iout);
Inew = Iout;



hold on

%outline the particles
[B,L] = bwboundaries(Iout, 'noholes');
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)
end

%number the particles
s = regionprops(L, 'centroid');
centroids = cat(1, s.Centroid);

for K=1: length(centroids)
    x=centroids(K,1);
    y=centroids(K,2);
    a=num2str(K);
    text(x,y,sprintf('%d',K),'color','g','Fontsize',12)

end

hold off


%find the areas and perimeters of the objects detected in pixels
D1 = regionprops(Iout, 'area'); 
D2 = regionprops (Iout, 'perimeter');
%results are in structure array form
%will say perimeter=0 when particle is a single pixel

%converting from structure array to matrix form
areatemp = cell2mat(struct2cell(D1)) ;
perimtemp = cell2mat(struct2cell(D2)) ;

%converting areas of particles to microns
area_microns = (areatemp/ratioperim)*1e6; % pixels to microns^2

%converting perimeter of particles to microns
perim_microns = (perimtemp/ratioperim)*1e6; %pixels to microns

%calculating radii assuming particles are perfect circles
radius_microns = ((area_microns)/pi).^(1/2);

%calculating diameter assuming particles are perfect circles
diameter_microns = radius_microns.*2;

