Adding ellipse angle and refactoring v0p12.

Added the angular transformation of the ellipses to create a 3D orientation of the ellipses
Refactoring the code for new variable names easier to handle.
Last non nuclear commit!

Signed-off-by: sebherbert <sherbert@pasteur.fr>

MainScripts/displayCombinedMap.m

@@ -3,12 +3,14 @@
 function displayCombinedMap(tableOutputDeproj,tableOutputBell,contour3D,outputFolder)
 % If not specified, the default display is now in real space (not projected)
 % For shape descriptor see publication Zdilla et al. 2016
-% Usual call is 
+% Usual call is
 % displayCombinedMap(tableOutputDeproj,tableOutputBell,dataCells.cellContour3D)
 
 % close all
 
-scriptVersion = 'displayCombinedMap_v0p3';
+scriptVersion = 'displayCombinedMap_v0p5';
+
+nSec = 5; % give 5 sec to the software to save the image before going to the next
 
 if nargin==3
     outputFolder = uigetdir(pwd,'Select the output folder');
@@ -24,16 +26,18 @@ dataBool.doPerimeter = true;
 dataBool.doAnisotropy = true;
 dataBool.doCircularity = true;
 dataBool.doRoundness = true;
-% 
+dataBool.doOrientation = true;
+
+
 % dataBool.doApical = false;
 % dataBool.doNeighbours = false;
 % dataBool.doAreaRatioRP = false;
 % dataBool.doAreaRatioPB = false;
 % dataBool.doPerimeter = false;
 % dataBool.doAnisotropy = false;
-% dataBool.doCircularity = true;
+% dataBool.doCircularity = false;
 % dataBool.doRoundness = false;
-
+% dataBool.doOrientation = true;
 
 save('displayParameters','scriptVersion','dataBool');
 
@@ -54,7 +58,7 @@ if dataBool.doApical
     dataVal = tableOutputDeproj.AreaReal;
     titleFig = 'Cell apical surface area (Real)';
     figName = 'mapApicalArea';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
 
 % Nbr of neighbours
@@ -62,7 +66,7 @@ if dataBool.doNeighbours
     dataVal = tableOutputDeproj.NbrNeighbours;
     titleFig = 'Nbr of neighbours';
     figName = 'mapNbrNeighbours';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
 
 % Area ratio (real vs proj)
@@ -70,15 +74,15 @@ if dataBool.doAreaRatioRP
     dataVal = tableOutputDeproj.AreaReal./tableOutputDeproj.AreaProj;
     titleFig = 'Cell apical surface area ratio (Real/Proj)';
     figName = 'mapAreaRatio_RP';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
-    
+
 % Area ratio (Proj vs Bell)
 if dataBool.doAreaRatioPB
     dataVal = tableOutputDeproj.AreaProj./tableOutputBell.AreaBell;
     titleFig = 'Cell apical surface area ratio (Proj/Bell)';
     figName = 'mapAreaRatio_PB';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
 
 % Perimeter
@@ -86,7 +90,7 @@ if dataBool.doPerimeter
     dataVal = tableOutputDeproj.PerimeterReal;
     titleFig = 'Cell perimeter (Real)';
     figName = 'mapPerimeter';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
 
 % Anisotropy (Real) (as calculated by Bellaiche lab: 1-1/elongation)
@@ -94,65 +98,100 @@ if dataBool.doAnisotropy
     dataVal = tableOutputDeproj.AnisotropyReal;
     titleFig = 'Cell anisotropy (Real)';
     figName = 'mapAnisotropy';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
 
 % Circularity is a shape descriptor that can mathematically indicate the degree of similarity to a perfect circle
-if dataBool.doCircularity % Error in the method   
+if dataBool.doCircularity % Error in the method
     dataVal = (4*pi*(tableOutputDeproj.AreaEllipseReal) ./ (tableOutputDeproj.PerimeterEllipseReal.^2));
     titleFig = 'Cell circularity (Real)';
     figName = 'mapCircularity';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
 end
 
 % Roundness is similar to circularity but is insensitive to irregular borders along the perimeter
 if dataBool.doRoundness
     dataVal = 4*tableOutputDeproj.AreaEllipseReal./(pi*(tableOutputDeproj.semiMajAxReal*2).^2);
-    titleFig = 'Cell roundness (Real)';
+    titleFig = 'Cell roundness (Real)'; 
     figName = 'mapRoundness';
-    dispMap(dataVal,titleFig,figName,contour3D);
+    dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec);
+end
+
+% Orientation of the cell is representing the orientation of the fitted
+% ellipse on the cell contour on the mesh. Vector length is weighted by the
+if dataBool.doOrientation
+    % Arrow length is function of the anisotropy
+    dataVec = tableOutputDeproj.OrientationEllipseReal.*...
+        tableOutputDeproj.AnisotropyReal;
+    % Set the origin of the arrow
+    dataOri = tableOutputDeproj.centerEllipseReal;
+    titleFig = 'Cell orientation (Real)';
+    figName = 'mapOrientation';
+    dispQuiverMap(dataOri,dataVec,titleFig,figName,nSec);
 end
 
+
 end
 
 
-function dispMap(dataVal,titleFig,figName,contour3D)
-    
-    %% Preparing the colormap
-    dataRange = max(dataVal)-min(dataVal);
-    greyLvlNbr = 200;
-    color2plot = round(((dataVal-min(dataVal))/dataRange)*(greyLvlNbr-1)+1);
-    fprintf('displaying: %s\n',titleFig);
-    dataColor = parula(greyLvlNbr);
-    
-    %% Plot the figure
-    figure
-    hold on
-    for bioCell = 1:numel(dataVal)
-        if isnan(dataVal(bioCell))
-            fprintf('Warning: A value was not set properly: skipping cell %d\n',bioCell);
-            continue
-        else
-            fill3(contour3D{bioCell}.VerticesOrdered(:,1),...
-                contour3D{bioCell}.VerticesOrdered(:,2),...
-                contour3D{bioCell}.VerticesOrdered(:,3),...
-                dataColor(color2plot(bioCell),:), ...
-                'LineStyle', 'none');
-        end
+function dispPolygonMap(dataVal,titleFig,figName,contour3D,nSec)
+
+%% Preparing the colormap
+dataRange = max(dataVal)-min(dataVal);
+greyLvlNbr = 200;
+color2plot = round(((dataVal-min(dataVal))/dataRange)*(greyLvlNbr-1)+1);
+fprintf('displaying: %s\n',titleFig);
+dataColor = parula(greyLvlNbr);
+
+%% Plot the figure
+figure
+hold on
+for bioCell = 1:numel(dataVal)
+    if isnan(dataVal(bioCell))
+        fprintf('Warning: A value was not set properly: skipping cell %d\n',bioCell);
+        continue
+    else
+        fill3(contour3D{bioCell}.cellCt(:,1),...
+            contour3D{bioCell}.cellCt(:,2),...
+            contour3D{bioCell}.cellCt(:,3),...
+            dataColor(color2plot(bioCell),:), ...
+            'LineStyle', 'none');
     end
-    
-    %% Set axes and colorbar
-    h = colorbar;
-    caxis([min(dataVal) max(dataVal)])
-    axis equal
-    title(titleFig);
-    xlabel('X position (µm)'); ylabel('Y position (µm)');
-    zlabel('Z position (µm)');
-
-    savefig(figName)
-    
 end
 
+%% Set axes and colorbar
+h = colorbar;
+caxis([min(dataVal) max(dataVal)])
+axis equal
+title(titleFig);
+xlabel('X position ({\mu}m)'); ylabel('Y position ({\mu}m)');
+zlabel('Z position ({\mu}m)');
+
+savefig(figName)
+export_fig(figName,'-png','-m5');
+pause(nSec);
+end
+
+
+function dispQuiverMap(dataOri,dataVec,titleFig,figName,nSec)
+
+fprintf('displaying: %s\n',titleFig);
+
+% plot the figure
+figure
+quiver3(dataOri(:,1),dataOri(:,2),dataOri(:,3),dataVec(:,1),dataVec(:,2),dataVec(:,3),'ShowArrowHead','off')
+axis equal
+
+%% Set axes and colorbar
+axis equal
+title(titleFig);
+xlabel('X position ({\mu}m)'); ylabel('Y position ({\mu}m)');
+zlabel('Z position ({\mu}m)');
+
+savefig(figName)
+export_fig(figName,'-png','-m5');
+pause(nSec);
+end
 
 
 

MainScripts/findContour.m

@@ -20,7 +20,7 @@ function cellContour2D = findContour(contour2Dpos,bmethod,PARAMS)
 % - cellContour2D.Line: paired vertices matrix to infer the pixel
 % connections -> deleted
 % - cellContour2D.Vertices: XY sorted vertices -> deleted
-% - cellContour2D.VerticesOrdered: XY sorted vertices in "CW" order
+% - cellContour2D.cellCt: XY sorted vertices in "CW" order
 
 fprintf('Recreating the cellular contours\n');
 
@@ -69,7 +69,7 @@ for bioCell = 1:length(contour2Dpos) % for all cells
             
         case 'bwboundary'
             [boundaries,~] = bwboundaries(miniImg,'noholes');
-            cellContour2D{bioCell}.VerticesOrdered = boundaries{1}+...
+            cellContour2D{bioCell}.cellCt = boundaries{1}+...
                 min(contour2Dpos{bioCell})-2; % replace in full image % -1 for start table at 1 and -1 for appened lines
             
         case 'mask2poly'
@@ -79,7 +79,7 @@ for bioCell = 1:length(contour2Dpos) % for all cells
             polyCell(:,1) = polyCell(:,2);
             polyCell(:,2) = polyCell(:,3);
             polyCell(:,3) = [];
-            cellContour2D{bioCell}.VerticesOrdered = polyCell+...
+            cellContour2D{bioCell}.cellCt = polyCell+...
                 min(contour2Dpos{bioCell})-2; % replace in full image % -1 for start table at 1 and -1 for appened lines
         
         otherwise
@@ -100,10 +100,10 @@ for bioCell = 1:length(contour2Dpos) % for all cells
                 end
                 Ori = find(cellContour2D{bioCell}.Lines(:,1)==Dest);
             end
-            cellContour2D{bioCell}.VerticesOrdered = tempVert;
+            cellContour2D{bioCell}.cellCt = tempVert;
             clear tempVert
     end
-    cellContour2D{bioCell}.VerticesOrdered = cellContour2D{bioCell}.VerticesOrdered*PARAMS.imSettings.latPixSize; % scale back to µm
+    cellContour2D{bioCell}.cellCt = cellContour2D{bioCell}.cellCt*PARAMS.imSettings.latPixSize; % scale back to µm
 end
 cellContour2D = cellContour2D';
 

MainScripts/formatTableOuputSurfaceCombine.m

@@ -10,10 +10,16 @@ tableOutputBell = table;
 tableOutputDeproj.Numbers = dataCells.numbers;
 tableOutputBell.Numbers = dataSeg.CELLS.numbers(dataCells.numbers);
 
+%% Cell center (ellipse real)
+for bioCell = 1:numel(dataCells.numbers)
+    EllipseRealCenter(bioCell,:) = dataCells.cellContour3D{bioCell}.ellipseReal.center;
+end
+tableOutputDeproj.EllipseRealCenter = EllipseRealCenter;
+
 %% Cell areas
 for bioCell = 1:numel(dataCells.numbers)
-    areaEllReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseFlatten.areaEllipse;
-    areaEllProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseFlatten.areaEllipse;
+    areaEllReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseRotViaOri.areaEllipse;
+    areaEllProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseProj.areaEllipse;
 end
 tableOutputDeproj.AreaReal = cell2mat(dataCells.area.areaRealTot);
 tableOutputDeproj.AreaProj = cell2mat(dataCells.area.areaProjTot);
@@ -25,12 +31,12 @@ tableOutputBell.AreaBell = dataSeg.CELLS.areas(dataCells.numbers);
 tableOutputDeproj.nbrNeighbours = dataSeg.CELLS.n_neighbors(dataCells.numbers);
 tableOutputBell.nbrNeighbours = dataSeg.CELLS.n_neighbors(dataCells.numbers);
 
-%% Cell anisotropy (Bellaiche style)
+%% Cell anisotropy (Bellaiche style) % Change for precalculated values !
 for bioCell = 1:numel(dataCells.numbers)
-    elongationReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseFlatten.semiMajAx / ...
-        dataCells.cellContour3D{bioCell}.ellipseFlatten.semiMinAx;
-    elongationProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseFlatten.semiMajAx / ...
-        dataCells.cellContour2D{bioCell}.ellipseFlatten.semiMinAx;
+    elongationReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseRotViaOri.semiMajAx / ...
+        dataCells.cellContour3D{bioCell}.ellipseRotViaOri.semiMinAx;
+    elongationProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseProj.semiMajAx / ...
+        dataCells.cellContour2D{bioCell}.ellipseProj.semiMinAx;
     elongationBell(bioCell) = dataSeg.CELLS.anisotropies(dataCells.numbers(bioCell));
 end
 tableOutputDeproj.anisostropyReal = 1-1./elongationReal';
@@ -39,8 +45,8 @@ tableOutputBell.anisotropyBell = elongationBell';
 
 %% Cell longer axis
 for bioCell = 1:numel(dataCells.numbers)
-    semiMajorAxisReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseFlatten.semiMajAx;
-    semiMajorAxisProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseFlatten.semiMajAx;
+    semiMajorAxisReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseRotViaOri.semiMajAx;
+    semiMajorAxisProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseProj.semiMajAx;
 end
 tableOutputDeproj.semiMajAxisReal = semiMajorAxisReal';
 tableOutputDeproj.semiMajAxisProj = semiMajorAxisProj';
@@ -48,20 +54,31 @@ tableOutputDeproj.semiMajAxisProj = semiMajorAxisProj';
 %% Cell angle (Projected only)
 for bioCell = 1:length(dataCells.numbers)
     %     angleReal(bioCell) =
-    %     dataCells.cellContour3D{bioCell}.ellipseFlatten.alpha; => Not
+    %     dataCells.cellContour3D{bioCell}.ellipseRotViaOri.alpha; => Not
     %     calculated yet
-    angleProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseFlatten.alpha;
+    angleProj(bioCell) = dataCells.cellContour2D{bioCell}.ellipseProj.alpha;
     angleBell(bioCell) = dataSeg.CELLS.orientations(dataCells.numbers(bioCell));
 end
 % tableOutputDeproj.angleReal = angleReal';
 tableOutputDeproj.angleProj = angleProj';
 tableOutputBell.angleBell = angleBell';
 
+%% Cell orientation real
+for bioCell = 1:numel(dataCells.numbers) % Only keep the end point of the vector. The beginning is the ellipse center
+    if isnan(dataCells.cellContour3D{bioCell}.ellipseReal.normOrientation(1,1))
+        EllipseRealOrient(bioCell,:) = dataCells.cellContour3D{bioCell}.ellipseReal.normOrientation(1,1);
+    else
+        EllipseRealOrient(bioCell,:) = dataCells.cellContour3D{bioCell}.ellipseReal.normOrientation(:,2)-...
+            dataCells.cellContour3D{bioCell}.ellipseReal.normOrientation(:,1); % Centered on the origin
+    end
+end
+tableOutputDeproj.EllipseRealOrient = EllipseRealOrient;
+
 %% Cell perimeters
 for bioCell = 1:length(dataCells.numbers)
-    perimeterReal(bioCell) = dataCells.cellContour3D{bioCell}.perimeter3D;
+    perimeterReal(bioCell) = dataCells.cellContour3D{bioCell}.perimeterReal;
     perimeterProj(bioCell) = dataCells.cellContour2D{bioCell}.perimeterProj;
-    perimeterEllipseReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseFlatten.perimeterEllipse;
+    perimeterEllipseReal(bioCell) = dataCells.cellContour3D{bioCell}.ellipseRotViaOri.perimeterEllipse;
 end
 tableOutputDeproj.perimeterReal = perimeterReal';
 tableOutputDeproj.perimeterProj = perimeterProj';
@@ -69,8 +86,8 @@ tableOutputDeproj.perimeterEllipseReal = perimeterEllipseReal';
 tableOutputBell.perimeterBell = dataSeg.CELLS.perimeters(dataCells.numbers);
 
 %% Field naming
-tableOutputDeproj.Properties.VariableNames = {'cellID' 'AreaReal' 'AreaProj' 'AreaEllipseReal' 'AreaEllipseProj'...
-    'NbrNeighbours' 'AnisotropyReal' 'AnisotropyProj' 'semiMajAxReal' 'semiMajAxProj' 'AngleProj'...
+tableOutputDeproj.Properties.VariableNames = {'cellID' 'centerEllipseReal' 'AreaReal' 'AreaProj' 'AreaEllipseReal' 'AreaEllipseProj'...
+    'NbrNeighbours' 'AnisotropyReal' 'AnisotropyProj' 'semiMajAxReal' 'semiMajAxProj' 'AngleProj' 'OrientationEllipseReal'...
     'PerimeterReal' 'PerimeterProj' 'PerimeterEllipseReal'};
 tableOutputBell.Properties.VariableNames = {'cellID' 'AreaBell' 'NbrNeighbours' 'AnisotropyBell' 'AngleBell'...
     'PerimeterBell'};