% Bead tracking algorithm using center of mass
% Written by Heejin for BE309 lab. Oct 16, 2006

clear all; close all;

NumberFrames = 5349;
ImageName = 'beadone';   % 10x10 image
BeadNumber = input('Bead # (1 or 2):');
Sum = input('Input the number of frames to sum:');
SumFrames = floor(NumberFrames/Sum);  % Number of 10 sumed image frames 

TempImage = imread([ImageName, num2str(1)],'jpg');
Xsize = size(TempImage,1);   % Number of rows
Ysize = size(TempImage,2);   % Number of columns
SumA = zeros(Xsize,Ysize,SumFrames);  
BeadPos = zeros(2,SumFrames);
TempA = zeros(Xsize,Ysize);

Threshold = input('Contour threshold: ');
ThresholdMtx = Threshold .* ones(Xsize, Ysize);

j = 1;   % index for SumFrame number
for i = 2 : NumberFrames
    TempA = TempA + double(imread([ImageName, num2str(i)],'jpg'));
    if (mod(i,Sum) == 1)
        SumA(:,:,j) = TempA;
        TempA = zeros(Xsize,Ysize);
        j = j + 1;
    end
end
'SumA matrix is filled'

for k = 1 : SumFrames
    BinaryMtx = ( SumA(:,:,k) > ThresholdMtx );
    LabeledMtx = double(BinaryMtx);
    LabeledSumA = double(BinaryMtx) .* double(SumA(:,:,k)); 
   
    %Find the centroid
    xsum = 0; ysum = 0;
    for p = 1: Xsize  % row: y-coordinate
        for q = 1: Ysize  % column: x-cooridinate
            xsum = xsum + q*LabeledSumA(p,q);
            ysum = ysum + p*LabeledSumA(p,q);    
        end
    end
    xc = xsum/sum(sum(LabeledSumA));
    yc = ysum/sum(sum(LabeledSumA));

    BeadPos(:,k) = [xc; yc];
    % BeadPos(i,:) = [stats.Centroid];
    % BeadPos is a matrix of 2 columns (X and Y) and NumberFrames rows.
end

y = [BeadPos(1,1:SumFrames)*0.141 ; BeadPos(2,1:SumFrames)*0.141];
fid = fopen(['BeadPos',num2str(BeadNumber),'.txt'],'w');
fprintf(fid,'%6.4f \t %6.4f \n',y);
fclose(fid)
['Position of the Bead #',num2str(BeadNumber),' in \mum is saved in BeadPos',num2str(BeadNumber),'.txt']

figure(1)
plot(BeadPos(1,1:SumFrames)*0.141,BeadPos(2,1:SumFrames)*0.141,'o-')
title('Bead Position using Centroid Method')
xlabel('\mum'); ylabel('\mum');