% This code use the 'trackmate' dataset to perform the cell
% segmentation of the complete time-lapse acquisition saved in the folder
% 'folder_out'. 
%
% The variable 'folder_out' must contain the location of the folder containing the reconstructed phase image.
%
% The array 'trackmate' must first be initialized with the output of the cell tracking Tracmate Fiji plugin, 
% i.e. it must contain the complete data saved into the' Spots in tracks statistics' output file
% 
% After segmentation several metrics are computed for each segmented cell. 
% In the 'trackmate' array the following are modified as follows :
%8  POSITION_Z       >> becomes the max. intensity [grey level intensity]
%9  POSITION_T       unchanged
%10 FRAME            unchanged
%11 RADIUS           >> becomes cell area [µm2]
%12 VISIBILITY       >> becomes cell dry mass [pg]
%13 MANUAL_COLOR     >> becomes cell volume
%14 MEAN_INTENSITY   >> becomes cell diameter [µm]
%15 MEDIAN_INTENSITY >> becomes cell thickness [µm]
%16 MIN_INTENSITY    >> becomes cell major axis length [µm]
%17 MAX_INTENSITY    >> becomes cell aspect_ratio [a.u.]
%18 TOTAL_INTENSITY  >> becomes cell circularity [a.u.]
%
% The sort_nat function can be found in the Mathworks file exchange.

bVisu=false;  % boolean for visual inspection
display('Cell segmentation ...');
cell=struct('track',0,'N',0,'pixels',0,'mean',0,'std',0);
Nframe=max(trackmate(:,10));
dx=1.67;
factor_CDM=pi*0.00000065/0.00000018*(dx^2)*1E-6^2/0.000000000001/128/(2*pi);
factor_V=dx^2*1E-6^2/128*pi/2/pi/0.025*0.00000065/1E-15*1000;

liste_im=dir(strcat(folder_out,'*.',Pic_format));
[~,order] = sort_nat({liste_im.name});
liste_im = liste_im(order);

for m=1:Nframe+1
    ifile=m;
    display(strcat(num2str(ifile),'/',num2str(length(liste_im))));
    filename=strcat(folder_out,liste_im(ifile).name);
    I=double(imread(filename));
    if (size(I,3)==3)
        I_RGB=I;
        I=I(:,:,1);
        I_RGB_contrast=1-(I_RGB-100)*2/255;
        I_RGB_contrast=imsharpen(I_RGB_contrast,'Radius',1.5,'Amount',1);
        mode_RGB=true;
    else
        mode_RGB=false;
    end
    
    phase=128*ones(size(I)+2*crop_size);phase=double(phase);
    phase(crop_size:crop_size+size(I,1)-1,crop_size:crop_size+size(I,2)-1)=I;
    mask=zeros(size(I));
    track_pos=find(trackmate(:,10)==m-1);
    if not(isempty(track_pos))
        xcg=max(uint32(trackmate(track_pos,6)/dx),1);
        ycg=max(uint32(trackmate(track_pos,7)/dx),1);
        list_xy=ycg+xcg*size(mask,1);
        % Seeding
        mask(ycg+(xcg-1)*size(mask,1))=trackmate(track_pos,4);
        track_pos_k=trackmate(track_pos,4);
        size_mask=size(mask,1);
        end_mask=size(mask,1)*size(mask,2);
        flag=true;n_loop=0;
        % Filling
        bw=uint16(imerode(~bwdist(I>128),strel('disk',1))); 
        mask=uint16(mask);
        flag=true;
        while((flag)&(n_loop<100))
            n_loop=n_loop+1;
            new_mask=imdilate(mask,strel('disk',2))-mask;
            new_mask=(uint16(new_mask).*bw).*uint16(mask==0);
            mask=max(mask,new_mask);
            if(sum(sum(new_mask))==0)
                flag=false;
            end
        end
        mask2=0*ones(size(mask)+2*crop_size);
        mask2(crop_size:crop_size+size(mask,1)-1,crop_size:crop_size+size(mask,2)-1)=mask;
        % Metrics calculation
        for k=track_pos'
            kk=trackmate(k,4);
            cell(k).track=kk;
            cell(k).N=k;
            cx=max(1,uint16(trackmate(k,6)/dx));
            cy=max(1,uint16(trackmate(k,7)/dx));
            phase_crop=imcrop(phase,[cx cy crop_size*2 crop_size*2]);
            support_crop=imcrop(mask2,[cx cy crop_size*2 crop_size*2]);
            support_crop=(support_crop==kk);
            if (sum(sum(support_crop))==0)
                trackmate(k,8)=128;
                trackmate(k,11:18)=0;
                cell(k).pixels=0;
                cell(k).std=0;
            else
                pointdata=regionprops(support_crop,phase_crop,'Area','Eccentricity','MaxIntensity','MinorAxisLength','MajorAxisLength','PixelValues'); % g=pointdata.WeightedCentroid;gx=uint8(g(2));gy=uint8(g(1));
                trackmate(k,8)=pointdata(1).MaxIntensity;                                                             % max intensity [grey level]
                trackmate(k,16)=pointdata(1).MinorAxisLength/pointdata(1).MajorAxisLength;                            % aspect ratio [a.u.]
                trackmate(k,17)=pointdata(1).MajorAxisLength*dx;
                trackmate(k,18)=pointdata(1).Eccentricity;
                trackmate(k,11)=round(sum(sum(support_crop))*dx*dx);                                                 % area [µm2]
                trackmate(k,12)=1/10*round(round(sum(sum(double(support_crop).*(phase_crop-128))))*factor_CDM*10);   % dry mass [pg]
                trackmate(k,13)=1/100*round(round(sum(sum(double(support_crop).*(phase_crop-128))))*factor_V*100);   % volume [fl]
                trackmate(k,14)=1/100*round(100*sqrt(trackmate(k,11)/pi)*2);                                         % mean diameter [µm]
                trackmate(k,15)=1/100*round(100*trackmate(k,13)/trackmate(k,11));                                    % mean thickness [µm]
                %                     cell(k).pixels=pointdata(1).PixelValues;
                cell(k).mean=mean(pointdata(1).PixelValues);
                cell(k).std=std(pointdata(1).PixelValues);
            end
        end
    end
    
    if ((bVisu)&(not(isempty(track_pos))))
        if(mode_RGB)
            c_I_RGB=I_RGB_contrast;
        else
            c_I_RGB=zeros(size(I,1),size(I,2),3);
            c_I_RGB(:,:,1)=1-(I-100)*2/255;
            c_I_RGB(:,:,2)=1-(I-100)*2/255;
            c_I_RGB(:,:,3)=1-(I-100)*2/255;
        end
        for k=track_pos_k'
            new_mask=zeros(size(mask));new_mask(mask==k)=255;
            R=c_I_RGB(:,:,1);R(new_mask>0)=R(new_mask>0)*cc(mod(k,128)+1,1);c_I_RGB(:,:,1)=R;
            G=c_I_RGB(:,:,2);G(new_mask>0)=G(new_mask>0)*cc(mod(k,128)+1,2);c_I_RGB(:,:,2)=G;
            B=c_I_RGB(:,:,3);B(new_mask>0)=B(new_mask>0)*cc(mod(k,128)+1,3);c_I_RGB(:,:,3)=B;
            new_mask=zeros(size(mask));new_mask(mask==k)=255;new_mask=edge(new_mask);
            R=c_I_RGB(:,:,1);R(new_mask>0)=cc(mod(k,128)+1,1);c_I_RGB(:,:,1)=R;
            G=c_I_RGB(:,:,2);G(new_mask>0)=cc(mod(k,128)+1,2);c_I_RGB(:,:,2)=G;
            B=c_I_RGB(:,:,3);B(new_mask>0)=cc(mod(k,128)+1,3);c_I_RGB(:,:,3)=B;
        end
        
        figure(1);clf;subplot(1,2,1);imagesc(c_I_RGB,[0 1]);
        set(1,'units','normalized','position',[0.05 0.2 0.9 0.62]);set(1,'Color','w');
        for k=1:max(trackmate(:,4))
            track_pos=find(trackmate(:,4)==k);
            track=trackmate(track_pos,:);
            track_pos=find(track(:,10)==m-1);
            if not(isempty(track_pos))
                xcg=track(1:track_pos,6)/dx;
                ycg=track(1:track_pos,7)/dx;
                imin=max(numel(xcg)-49,1);
                hold on;plot(xcg(imin:end),ycg(imin:end),'-','color',cc(mod(k,128)+1,:),'LineWidth',1);
                hold off;
            end
        end
        subplot(1,2,2);imagesc(c_I_RGB,[0 1]);
        for k=1:max(trackmate(:,4))
            track_pos=find(trackmate(:,4)==k);
            track=trackmate(track_pos,:);
            track_pos=find(track(:,10)==m-1);
            if not(isempty(track_pos))
                xcg=track(1:track_pos,6)/dx;
                ycg=track(1:track_pos,7)/dx;
                hold on;
                text(xcg(track_pos)+10,ycg(track_pos),num2str(k),'color','r','FontSize',5);
                text(xcg(track_pos)+10,ycg(track_pos)+8,strcat(num2str(track(track_pos,12)),'pg'),'color','r','FontSize',5);
                hold off;
            end
        end
        drawnow;
    end
end
save(filename,'trackmate','-v7.3');