filename = 'cell32_c31_s7580.tif';
framenumber = 1;

for framenumber=1:32
    I = imread(filename,framenumber);
    I = imadjust(I);
    imshow(I,[])
    %% first mask: cell border
    th = graythresh(I);
    bw = im2bw(I,th-0.52*th);
    %bw2 = uint16(bw);
    %imshowpair(I,bw2, 'blend')
    % cleaning first mask
    bw = bwareaopen(bw,5000);
    se = strel('disk',20);
    bw = imclose(bw,se);
    bw = imfill(bw,'holes');
    
    %% for some specific cells (in frame with others)
%     bw=imclearborder(bw,4);
    
    %% added for some specific cells
    BWclearborder=imclearborder(bw,4);
    
    %bw2 = uint16(bw);
    %imshowpair(I,bw2, 'blend')
    % analysis
    stats = regionprops(bw,'Centroid','MajorAxisLength','MinorAxisLength','Orientation');
    if ~isempty(stats)
        ss = length(stats);
        if ss>1
            d=[];
           for j=1:ss
               d(j)=(stats(j).Centroid(1)-512).^2+(stats(j).Centroid(2)-512).^2;
           end
        [~,ig] = min(d)
        else
            ig=1;
        end
        cell(:,framenumber) = [stats(ig).Centroid(1),stats(ig).Centroid(2)];
        % get the contour
        bwp = bwperim(bw);
        [r,c] = find(bwp==1,1);
        cellcontour = bwtraceboundary(bw,[r c],'W',8,Inf,'counterclockwise');
    else
        cell(:,framenumber) = [1,1];
    end
    %% second mask: golgi
    I = imread(filename,framenumber);
%     I = imadjust(I);
    th = graythresh(I);
%     bw = im2bw(I,th);
    
    bw = im2bw(I,th-0.98*th);
    
    %bw2 = uint16(bw);
    %imshowpair(I,bw2, 'blend')
    % cleaning second mask
    bw = bwareaopen(bw,3000);
    se = strel('disk',20);
    bw = imclose(bw,se);
    bw = imfill(bw,'holes');
    if min(min(bw))==1
        I = imread(filename,framenumber);
        I = imadjust(I);
        th = graythresh(I);
        bw = im2bw(I,0.99);
        %bw2 = uint16(bw);
        %imshowpair(I,bw2, 'blend')
        % cleaning second mask
        bw = bwareaopen(bw,3000);
        se = strel('disk',20);
        bw = imclose(bw,se);
        bw = imfill(bw,'holes');
    end
    %bw2 = uint16(bw);
    %imshowpair(I,bw2, 'blend')
    % anlaysis
    stats = regionprops(bw,'Centroid');
    if ~isempty(stats)
        ss = length(stats);
        if ss>1
            d=[];
           for j=1:ss
               d(j)=(stats(j).Centroid(1)-512).^2+(stats(j).Centroid(2)-512).^2;
           end
        [~,ig] = min(d)
        else
            ig=1;
        end
        golgi(:,framenumber) = [stats(ig).Centroid(1),stats(ig).Centroid(2)];
        % get the contour
        bwp = bwperim(bw);
        [r,c] = find(bwp==1,1);
        golgicontour = bwtraceboundary(bw,[r c],'W',8,Inf,'counterclockwise');
    else
       golgi(:,framenumber) = [1,1];
    end
    %% final figure
    I = imread(filename,framenumber);
    figure(1)
    imshow(I,[]);
    hold on
    plot(cell(1,framenumber),cell(2,framenumber),'oc','markersize',6)
    plot(golgi(1,framenumber),golgi(2,framenumber),'og','markersize',6)
    plot(cellcontour(:,2),cellcontour(:,1),'c','LineWidth',0.5);
    plot(golgicontour(:,2),golgicontour(:,1),'g','LineWidth',0.5);
    %pause(0.5)
    title(['frame ',num2str(framenumber)])
    drawnow
end
%% for saving
save('cell32.mat','cell','golgi');