%%%%%%%%%

% computeMorphoMaps_colorlimits('file directory', -25, 25);

%%%%%%%
%%
function computeMorphoMaps_colorlimits(cell_mat_file,min_val,max_val)

    if( nargin == 1 )
        min_val= -1;
        max_val = -1;
    end
    [PathName,filename,~] = fileparts(cell_mat_file);
    tmp = strsplit(filename,'_');
    acquisition_name = [tmp{1} '_' tmp{2} '_']; % to find the prefix of the acquisition ex: c012_pos4

    cell_file = dir([PathName '\' acquisition_name  '*res_cell.mat']);
    golgi_file = dir([PathName '\' acquisition_name  '*res_golgi.mat']);
    nucl_file = dir([PathName '\' acquisition_name  '*res_nucleus.mat']);
    if( isempty(cell_file) || isempty(golgi_file) || isempty(nucl_file) )
        msgbox('You should first segment independently the cell/golgi/nuclei of your data set','Incorrect file','error');
    else
        data_cell = load(fullfile(cell_file.folder,cell_file.name));
        data_golgi = load(fullfile(golgi_file.folder,golgi_file.name));
        data_nucleus = load(fullfile(nucl_file.folder,nucl_file.name));
    end
    
    
    if ( isfield(data_cell,'cell_centroid') ) % old version, put the result in "cell_centroid_imgRef"
        cell_centroid_imgRef = data_cell.cell_centroid;
    else
        cell_centroid_imgRef = data_cell.cell_centroid_imgRef;
    end
    BW_cell = data_cell.BW_cell;
    [~, Hg, Nf] = size(BW_cell);
    
    if ( isfield(data_golgi,'golgi_centroid') ) % old version, put the result in "golgi_centroid_imgRef"
        golgi_centroid_imgRef = data_golgi.golgi_centroid;
    else
        golgi_centroid_imgRef = data_golgi.golgi_centroid_imgRef;
    end
    
    if ( isfield(data_nucleus,'nucl_centroid') ) % old version, put the result in "golgi_centroid_imgRef"
        nucl_centroid_imgRef = data_nucleus.nucl_centroid;
    else
        nucl_centroid_imgRef = data_nucleus.nucl_centroid_imgRef;
    end

    %read position file
    pos_file = dir([PathName '\' acquisition_name '*Positions.mat']);
    if ( ~isempty(pos_file) )
        load(fullfile(PathName,pos_file.name));
    else
        warning = 'No positions file found, impossible to compute the translations';
    end
    % keep only some positions
    p = p(1:2,:); % taking on point every pt_interval, starting with firstpoint
    X = p(1,:)-p(1,1);
    Y = p(2,:)-p(2,1);
       
    scale_param_pos = 1/0.2167;
    trX = X*scale_param_pos; % carefull there seems to be more positions than the number of images
    trY = Y*scale_param_pos;
  
    
    coord_cell_refOrtho = [cell_centroid_imgRef(2,:)+trX(1:Nf); Hg-cell_centroid_imgRef(1,:)+trY(1:Nf)];
    smoothX_cell = sgolayfilt(coord_cell_refOrtho(1,:), 2, 35);
    smoothY_cell = sgolayfilt(coord_cell_refOrtho(2,:), 2, 35);
    coord_cell_smooth = [smoothX_cell; smoothY_cell];
    
    coord_golgi_refOrtho = [golgi_centroid_imgRef(2,:)+trX(1:Nf); Hg-golgi_centroid_imgRef(1,:)+trY(1:Nf)];
    coord_nucl_refOrtho = [nucl_centroid_imgRef(2,:)+trX(1:Nf); Hg-nucl_centroid_imgRef(1,:)+trY(1:Nf)];
        
    numPoints = 100;
    
    MD = zeros(numPoints,Nf-1);
    MDn = zeros(numPoints,Nf-1);
    MDns = zeros(numPoints,Nf-1);
    ind_pt_contours_golgi = zeros(Nf,1);
    ind_pt_contours_traj = zeros(Nf,1);
    ind_pt_contours_centroid = zeros(Nf,1);
    
    
    boundary=bwboundaries(BW_cell(:,:,1));
    boundaryX_refOrtho = boundary{1}(:,2) + trX(1);
    boundaryY_refOrtho = (Hg-boundary{1}(:,1)) + trY(1); 
    boundaryLength=size(boundaryX_refOrtho,1)-1;
    
    indexCol=[0:(numPoints/boundaryLength):numPoints];
    previousEdgeCoors(:,1)=interp1(indexCol,boundaryX_refOrtho,0.5:numPoints);
    previousEdgeCoors(:,2)=interp1(indexCol,boundaryY_refOrtho,0.5:numPoints);
    
    for k=1:Nf-1
        frame = k;
        boundary=bwboundaries(BW_cell(:,:,frame+1)); % in image referential -> resultat en (Y,X) %% for looking: figure; imshow(BW_cell(:,:,frame+1)); hold on; plot(boundary{1}(:,2),boundary{1}(:,1),'r');

        boundaryX_refOrtho = boundary{1}(:,2) + trX(frame+1);
        boundaryY_refOrtho = (Hg-boundary{1}(:,1)) + trY(frame+1);

        boundaryLength=size(boundaryX_refOrtho,1)-1;

        indexCol=0:(numPoints/boundaryLength):numPoints;
        edgeCoors(:,1)=interp1(indexCol,boundaryX_refOrtho,0.5:numPoints); %X-direction
        edgeCoors(:,2)=interp1(indexCol,boundaryY_refOrtho,0.5:numPoints); %Y-direction

        distSum=zeros(numPoints-1,1);

        distSum(1) = sum(sum((edgeCoors-previousEdgeCoors).^2));
        for i_pt=1:(numPoints-1)
            temp=matrixCircularlyPermuteRows(edgeCoors,i_pt);
            distSum(i_pt)=sum(sum((temp-previousEdgeCoors).^2));
        end

        [~,phaseVal]=min(distSum);
        edgeCoors=matrixCircularlyPermuteRows(edgeCoors,phaseVal); 
        % record morphodynamics
        [a,b] = findnorm(previousEdgeCoors);

        ind_Nan_a = find(isnan(a)); 

        norm = sqrt((edgeCoors(:,1)-previousEdgeCoors(:,1)).^2+(edgeCoors(:,2)-previousEdgeCoors(:,2)).^2);
        dproj = abs(a.*edgeCoors(:,1)+b-edgeCoors(:,2))./sqrt(a.^2+1);
        dproj(ind_Nan_a) = abs(edgeCoors(ind_Nan_a,1)-previousEdgeCoors(ind_Nan_a,1));

        dnorm = sqrt(norm.^2-dproj.^2);
        dnorm(dproj>norm) = 0;

        signC = inpolygon(edgeCoors(:,1),edgeCoors(:,2),[previousEdgeCoors(:,1); previousEdgeCoors(1,1)],[previousEdgeCoors(:,2); previousEdgeCoors(1,2)]);
        signC = double(signC); 
        signC = 1-signC;
        %signC(signC == 1) = -1;
        signC(signC == 0) = -1;
        MD(:,k) = signC.*norm;
        MDn(:,k) = signC.*dnorm;
        ind_pt_contours_golgi(k) = findClosestPoint(coord_golgi_refOrtho(:,k+1),coord_nucl_refOrtho(:,k+1),edgeCoors);
        ind_pt_contours_traj(k) = findClosestPoint(coord_cell_smooth(:,k+1),coord_cell_smooth(:,k),edgeCoors);
        ind_pt_contours_centroid(k) = findClosestPoint(coord_cell_refOrtho(:,k+1),coord_cell_refOrtho(:,k+1)-[1; 0],edgeCoors);
        previousEdgeCoors = edgeCoors; 
    end
        
       
        
    for i_pt=1:numPoints
        MDns(i_pt,:) = smooth(MDn(i_pt,:),10,'rloess');
    end

    
    %MDn = circshift(MDn,shift);
    %MD = circshift(MD,shift);

    displayImageAndPoints(MDns,ind_pt_contours_golgi,numPoints,min_val,max_val);
    title('Nucleus/golgi axis');
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_noshift_nuclgolg.fig']));
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_noshift_nuclgolg.png']));
    
    displayImageAndPoints(MDns,ind_pt_contours_traj,numPoints,min_val,max_val);
    title('Cell trajectory');
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_noshift_celltraj.fig']));
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_noshift_celltraj.png']));
    
    displayImageAndPoints(MDns,ind_pt_contours_centroid,numPoints,min_val,max_val);
    title('X axis');
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_noshift_Xaxis.fig']));
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_noshift_Xaxis.png']));
    
    ShiftMorphoMap(MDn,numPoints,length(ind_pt_contours_golgi),ind_pt_contours_golgi,min_val,max_val);
    title('Golgi centered');
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_shift_golgcenter.fig']));
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_shift_golgcenter.png']));
    
    ShiftMorphoMap(MDn,numPoints,length(ind_pt_contours_traj),ind_pt_contours_traj,min_val,max_val);
    title('Cell trajectory centered');
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_shift_cellTrajcenter.fig']));
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_shift_cellTrajcenter.png']));
    ShiftMorphoMap(MDn,numPoints,length(ind_pt_contours_centroid),ind_pt_contours_centroid,min_val,max_val);
    title('X axis centered');
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_shift_Xaxis.fig']));
    saveas(gcf,fullfile(PathName,[ acquisition_name 'morpho_map_shift_Xaxis.png']));
    
    function ShiftMorphoMap(MDn,numPoints,nbTimePt, ind_pt_cont_plot,val_min,val_max)
        % shift
        MDcs = [];
        for i=1:nbTimePt-1
            shift = ind_pt_cont_plot(i)-numPoints/2;
            v = circshift(MDn(:,i),-shift);
            MDcs = [MDcs,v];
        end

        for i=1:numPoints
                    MDns(i,:) = smooth(MDcs(i,:),10,'rloess');
        end

        figure('units','normalized','outerposition',[0 0 1 1]);
        f1 = fspecial('Gaussian',5,1);
        MDns= imfilter(MDns,f1,'replicate');
        maxM = max(max(MDns));
        minM = min(min(MDns));
        imshow(MDns,[0.7*minM 0.7*maxM],'initialmagnification','fit'); hold on; axis on;
        ax = gca();
        ax.YTickLabel = [];
        ax.YTick = [];
        plot(1:nbTimePt,repmat(50,1,nbTimePt),'k','linewidth',3);
        colormap('jet'); 
         if( val_min ~= val_max )
            caxis manual
            caxis([val_min val_max]);
        end
        colorbar
        
    end
    
%     saveas(gcf,[PathName,acquisition_name,'Morpho.jpg']);
%     saveas(gcf,[PathName,acquisition_name,'Morpho.fig']);
%     save([PathName,acquisition_name,'Morpho.mat'],'ind_pt_cont_plot','MDn');
    save([acquisition_name,'Morpho.mat'],'MDn');
    
    function displayImageAndPoints(MDns,ind_pt_contours,numPoints,val_min,val_max)
        maxM = max(max(MDns));
        minM = min(min(MDns));
        shift = numPoints/2-ind_pt_contours(1);
        MDns_plot = circshift(MDns,shift);
        
        figure('units','normalized','outerposition',[0 0 1 1]);
        f1 = fspecial('Gaussian',5,1);
        MDns_plot= imfilter(MDns_plot,f1,'replicate');
        
        imshow(MDns_plot,[0.7*minM 0.7*maxM],'initialmagnification','fit'); hold on; axis on;
        ax = gca();
        ax.YTickLabel = [];
        ax.YTick = [];
        ind_pt_contours = ind_pt_contours(1:Nf-1)+shift;
        ind_pt_contours(ind_pt_contours>numPoints) = ind_pt_contours(ind_pt_contours>numPoints)-numPoints;
        ind_pt_contours(ind_pt_contours<1) = numPoints+ind_pt_contours(ind_pt_contours<1);        
        plot(1:Nf-1,ind_pt_contours,'k','linewidth',3);
        colormap('jet'); 
        
        if( val_min ~= val_max )
            caxis manual
            caxis([val_min val_max]);
        end
        colorbar;
    end
      
    function ind_pt_cont = findClosestPoint(pt1, pt2, edgeCoors)
        a_d = (pt1(2)-pt2(2))/( pt1(1) - pt2(1));
        b_d = pt1(2)-a_d*pt1(1) ;
        sign_inc = sign(pt1(1) - pt2(1));
        xp = pt1(1);

        xp = xp+sign_inc*[0.5:0.5:300];
        yp = a_d*xp+b_d;

        pos_inCell = inpolygon(xp,yp,[edgeCoors(:,1); edgeCoors(1,1)],[edgeCoors(:,2); edgeCoors(1,2)]);
        I_out = find(pos_inCell == 0);
        coordX = xp(min(I_out)); 
        coordY = yp(min(I_out));      

        [~,ind_pt_cont] = min(sqrt( (edgeCoors(:,1)-coordX).^2+(edgeCoors(:,2)-coordY).^2 )); 
    end



     function newedges = matrixCircularlyPermuteRows(edges,phase)
            tot = length(edges);
            newedges = zeros(size(edges));
            newedges(1:tot-phase,:) = edges(phase+1:end,:);
            newedges(tot-phase+1:end,:) = edges(1:phase,:);
     end

 
    function [a,b] = findnorm(pt)
        npt = length(pt);
        a = zeros(npt,1);
        b = zeros(npt,1);
        for i_pt=2:npt-1
            x1 = pt(i_pt-1,1);
            y1 = pt(i_pt-1,2);
            x2 = pt(i_pt+1,1);
            y2 = pt(i_pt+1,2);
            if y1==y2
                a(i_pt) = nan;
                b(i_pt) = nan;

            else
                a(i_pt) = -(x2-x1)/(y2-y1);
                b(i_pt) = pt(i_pt,2)-a(i_pt)*pt(i_pt,1);
            end
        end
        
        x1 = pt(end,1);
        y1 = pt(end,2);
        x2 = pt(2,1);
        y2 = pt(2,2);
        if y1==y2
            a(1) = nan;
            b(1) = nan;
        else
            a(1) = -(x2-x1)/(y2-y1);
            b(1) = pt(1,2)-a(1)*pt(1,1);
        end
        
        x1 = pt(end-1,1);
        y1 = -pt(end-1,2);
        x2 = pt(1,1);
        y2 = -pt(1,2);
        if y1==y2
            a(end) = nan;
            b(end) = nan;
        else
            a(end) = -(x2-x1)/(y2-y1);
            b(end) = pt(end,2)-a(end)*pt(end,1);

        end
    end


end