Compute and store tissue local curvatures.

src/@deproj/deproj.m

@@ -31,6 +31,9 @@ classdef deproj
         % Plot the 2D ellipses on the tissue surface.
         [ hf, hc, he ] = plot_fit_ellipse( obj, scale_bar_length )
         
+        % Figure with the local curvaure for a collection of epicells.
+        [ hf, ax1, ax2, ax3 ] = plot_curvatures( obj, scale_bar_length )
+        
         %% Helpers.
         % They are public in case of.
         
@@ -54,7 +57,19 @@ classdef deproj
         
         %Plots the ellipses, colored by the specified values.
         hts = add_ellipse_variable( obj, values, cmap, ax )
-
+        
+        % Add the tissue plot colored with the mean curvature.
+        hts = add_plot_curvature_mean( obj, ax )
+        
+        % Add the tissue plot colored with the Gaussian curvature.
+        hts = add_plot_curvature_gauss( obj, ax )
+        
+        % Add the tissue plot colored with the first principal curvature.
+        hts = add_plot_curvature_k1( obj, ax )
+        
+        % Add the tissue plot colored with the second principal curvature.
+        hts = add_plot_curvature_k2( obj, ax )
+                
     end
     
      %% Public static methods: builders.

src/@deproj/from_heightmap.m

@@ -82,14 +82,14 @@ function obj = from_heightmap( ...
     end
     
     % Smooth the height-map over a scale smaller than a cell.
-    H = imgaussfilt( H, object_scale );
+    Hs1 = imgaussfilt( H, object_scale );
     
     % For junction.
-    z_junction = deproj.get_z( junction_graph.Nodes.Centroid, H, pixel_size, voxel_depth );
+    z_junction = deproj.get_z( junction_graph.Nodes.Centroid, Hs1, pixel_size, voxel_depth );
     
     % For objects.
     for i = 1 : n_objects
-        z_obj = deproj.get_z( objects( i ).boundary, H, pixel_size, voxel_depth );
+        z_obj = deproj.get_z( objects( i ).boundary, Hs1, pixel_size, voxel_depth );
         objects( i ).boundary = [ objects( i ).boundary z_obj ];
         objects( i ).center(3) = mean( z_obj );
     end
@@ -135,6 +135,33 @@ function obj = from_heightmap( ...
         end
     end
     
+    %% Compute local curvature from the smoothed height-map.
+    % Ref: http://users.vcnet.com/simonp/curvature.pdf
+    
+    fprintf('Compuing tissue local curvature.\n' )
+    
+    % We need to smooth the height-map over the scale of several cells.
+    Hs2 = imgaussfilt( Hs1, 3 * object_scale );
+    
+    if invert_z
+        Hs2 = -Hs2;
+    end
+    
+    [ Hx , Hy  ] = gradient( Hs2 );
+    [ Hxx, Hxy ] = gradient( Hx );
+    [  ~ , Hyy ] = gradient( Hy );
+    % Gaussian curvature.
+    Nk = ( 1. + Hx .^ 2 + Hy .^ 2 );
+    curvGauss = ( Hxx .* Hyy - Hxy .^ 2 ) ./ ( Nk .^ 2 );
+    % Mean curvature.
+    Dk1 = ( 1. + Hx .^ 2 ) .* Hyy;
+    Dk2 = ( 1. + Hy .^ 2 ) .* Hxx;
+    Dk3 = - 2. * Hx .* Hy .* Hxy;
+    curvMean = ( Dk1 + Dk2 + Dk3 ) ./ ( 2 * Nk .^ 1.5 );
+    % Principle curvatures.
+    curvK1 = curvMean + sqrt( curvMean .* curvMean - curvGauss );
+    curvK2 = curvMean - sqrt( curvMean .* curvMean - curvGauss );
+    
     %% Create epicell instances.
     
     fprintf('Computing morphological descriptors of all objects.\n' )
@@ -143,7 +170,14 @@ function obj = from_heightmap( ...
     epicells = repmat( epicell, n_objects, 1 );
     for i = 1 : n_objects
         o = objects( i );
-        epicells( i ) = epicell( o.boundary, o.junctions, i  );
+        
+        % Get center position in pixel again.
+        xc = round( o.center( 1 ) / pixel_size );
+        yc = round( o.center( 2 ) / pixel_size );
+        
+        curvatures = [ curvMean( yc, xc ), curvGauss( yc, xc ), curvK1( yc, xc), curvK2( yc, xc ) ];
+        
+        epicells( i ) = epicell( o.boundary, o.junctions, i, curvatures );
     end
     
     obj = deproj( epicells, junction_graph, units );

src/@epicell/epicell.m

@@ -9,21 +9,26 @@ classdef epicell
         area
         perimeter
         euler_angles
+        curvatures
         ellipse_fit
         eccentricity
         proj_direction
         uncorrected_area
         uncorrected_perimeter
-        id
+        id        
     end
     
     methods
-        function obj = epicell( boundary, junction_ids, id )
-            %EPICELL Construct an epicell from a N x 3 list of points.
+        
+        function obj = epicell( boundary, junction_ids, id , curvatures )
+        %EPICELL Construct an epicell from a N x 3 list of points.
             
             if nargin == 0
                 return
             end
+            if nargin >= 4
+                obj.curvatures = curvatures;
+            end
             
             % Reduce polygon in 2D then put back Z. Takes time.
             p2d = boundary( :, 1 : 2 );