Reorganize processing in a mother class 'deproj'.

This class mainly manages a collection of epicells, has
static methods for its creation from images, and has some
methods for plotting things.

RunExample.m

@@ -31,129 +31,45 @@ heightmap_filename   = 'HeightMap-2.tif';
 % Physical size of a pixel in X and Y. This will be used to report sizes in
 % µm.
 pixel_size = 0.183; % µm
+units = 'µm';
 
 % Z spacing.
 % How many µm bewtween each Z-slices. We need this to convert the
 % height-map values, that stores the plane of interest, into µm.
 voxel_depth = 1.; % µm
 
+
 % Try to remove objects that have a Z position equal to 0. Normally this
 % value reports objects that are out of the region-of-interest.
 prune_zeros = true;
+inpaint_zeros = true;
+
+% Invert z for plotting.
+invert_z = true;
 
-%% Read a mask file and convert it to objects.
+%% Read files.
 
 % ImageJ mask.
 fprintf('Opening mask image: %s\n', mask_filename )
 I = imread( fullfile( root_folder, mask_filename ) );
 
-% The conversion can take up to 30s for a 1000 x 1000 mask image.
-fprintf('Converting mask to objects.\n' )
-tic 
-% We will downsample later: for now we will need all the pixels to have a
-% robust estimate of the ellipse and plane fit.
-downsample = false;
-[ objects, junction_graph ] = mask_to_objects( I, downsample );
-n_objects = numel( objects );
-n_junctions = numel( junction_graph.Nodes );
-fprintf('Converted a %d x %d mask in %.1f seconds. Found %d objects and %d junctions.\n', ...
-    size(I,1), size(I,2), toc, n_objects, n_junctions )
-
-%% Compute some scale of an object radius.
-
-% Median area (pixel units)
-med_area = median(arrayfun( @(s) polyarea( s.boundary(:,1), s.boundary(:,2) ),  ...
-    objects ) ); % pixels^2
-
-object_scale = 2 * sqrt( med_area )/ pi; % pixel units
-
-fprintf('Typical object scale: %.1f pixels or %.2f µm.\n', ...
-    object_scale, object_scale * pixel_size )
-
-
-%% Scale to physical coordinates.
-
-% Scale junction XY coordinates to µm.
-junction_graph.Nodes.Centroid( :, 1:2 ) = junction_graph.Nodes.Centroid( :, 1:2 ) * pixel_size;
-
-% Scale object XY boundary to µm.
-for i = 1 : n_objects
-    objects( i ).boundary   = objects( i ).boundary * pixel_size;
-    objects( i ).center     = objects( i ).center * pixel_size;
-end
-
-%% Load the height-map and add the Z coordinates to objects.
-
 fprintf('Opening height-map image: %s\n', heightmap_filename )
 H = imread( fullfile( root_folder, heightmap_filename ) );
 
-if prune_zeros
-    mask = H == 0;
-    H = regionfill( H, mask );
-    H( H == 0 ) = NaN;
-end
-
-% Smooth the height-map over a scale smaller than a cell.
-H = imgaussfilt( H, object_scale );
-
-% For junction.
-z_junction = get_z( junction_graph.Nodes.Centroid, H, pixel_size, voxel_depth );
-
-% For objects.
-for i = 1 : n_objects
-    z_obj = get_z( objects( i ).boundary, H, pixel_size, voxel_depth );
-    objects( i ).boundary = [ objects( i ).boundary z_obj ];
-    objects( i ).center(3) = mean( z_obj );
-end
-
-%% Possibly remove the junctions and cells found at z = 0.
-
-if prune_zeros
-    
-    bad_junctions = find( z_junction == 0 | isnan( z_junction ) );
-    junction_graph = junction_graph.rmnode( bad_junctions );
-    z_junction( bad_junctions ) = [];
-
-    fprintf( 'Removed %d junctions at Z=0.\n', numel( bad_junctions ) )
-    
-    objects_to_prune = false( n_objects, 1 );
-    for i = 1 : numel( objects )
-        o = objects( i );
-        if any( isnan( o.boundary(:) ) ) || ~isempty( intersect( bad_junctions, o.junctions ) )
-            objects_to_prune( i ) = true;
-        end        
-    end
-    
-    objects( objects_to_prune ) = [];
-    fprintf( 'Removed %d objects touching these junctions.\n', sum( objects_to_prune ) )
-    
-   n_objects = numel( objects );
-   n_junctions = numel( junction_graph.Nodes );
-    
-end
-
-%% Invert Z position for plotting.
-
-max_z = max( z_junction );
-z_junction = max_z -  z_junction;
-junction_graph.Nodes.Centroid = [ junction_graph.Nodes.Centroid z_junction ];
-
-for i = 1 : n_objects
-    objects( i ).boundary( :, 3 ) = max_z - objects( i ).boundary( :, 3 );
-    objects( i ).center( 3 ) = max_z - objects( i ).center( 3 );
-end
-
-%% Create epicell instances.
-
-epicells = repmat( epicell, n_objects, 1 );
-for i = 1 : n_objects    
-    o = objects( i );
-    epicells( i ) = epicell( o.boundary, o.junctions, i  );
-end
+%% Create deproj instance.
+
+dpr = deproj.from_heightmap( ...
+    I, ...
+    H, ...
+    pixel_size, ...
+    voxel_depth, ...
+    units, ...
+    invert_z, ...    
+    inpaint_zeros, ...
+    prune_zeros );
 
 %% Plot euler angles.
 
 close all
-plot_fit_plane( epicells );
-
+plot_fit_plane( dpr );
 

src/add_plot_euler_alpha.m → src/@deproj/add_plot_euler_alpha.m

-function hts = add_plot_euler_alpha( ax, epicells )
+function hts = add_plot_euler_alpha( obj, ax )
 %ADD_PLOT_EULER_ALPHA Add the tissue plot colored with the 1st euler angle.
     
+    epicells = obj.epicells;
     angles = vertcat( epicells.euler_angles );
     alphas = rad2deg( angles( :, 1 ) );
     
@@ -8,6 +9,6 @@ function hts = add_plot_euler_alpha( ax, epicells )
     neg_alphas = alphas < 0;
     alphas( neg_alphas ) = 180 + alphas( neg_alphas );
 
-    hts = add_plot_variable( ax, { epicells.boundary }, alphas );
+    hts = add_plot_variable( obj, alphas, ax );
 end
 

src/add_plot_euler_beta.m → src/@deproj/add_plot_euler_beta.m

-function hts = add_plot_euler_beta( ax, epicells )
+function hts = add_plot_euler_beta( obj, ax )
 %ADD_PLOT_EULER_BETA Add the tissue plot colored with the 2nd euler angle.
     
+    epicells = obj.epicells;
     angles = vertcat( epicells.euler_angles );
     betas = rad2deg( angles( :, 2 ) );
     
@@ -8,6 +9,6 @@ function hts = add_plot_euler_beta( ax, epicells )
     large_betas = betas > 90;
     betas( large_betas ) = 180 - betas( large_betas );
 
-    hts = add_plot_variable( ax, { epicells.boundary }, betas );
+    hts = add_plot_variable( obj, betas, ax );
 end
 

src/add_plot_euler_gamma.m → src/@deproj/add_plot_euler_gamma.m

-function hts = add_plot_euler_gamma( ax, epicells )
+function hts = add_plot_euler_gamma( obj, ax )
 %ADD_PLOT_EULER_GAMMA Add the tissue plot colored with the 3rd euler angle.
     
+    epicells = obj.epicells;
     angles = vertcat( epicells.euler_angles );
     gammas = rad2deg( angles( :, 3 ) );
-    hts = add_plot_variable( ax, { epicells.boundary }, gammas );
+    hts = add_plot_variable( obj, gammas, ax );
 end
 

src/add_plot_id.m → src/@deproj/add_plot_id.m

-function hts = add_plot_id( epicells, ax )
+function hts = add_plot_id( obj, ax )
 %ADD_PLOT_ID Add the epicell ids to the specified plot axes.
 
+    epicells = obj.epicells;
     n_objects = numel( epicells );
     hts = NaN( n_objects, 1 );
     
     for i = 1 :  n_objects
-        o = objects( i );
-        hts(i ) = text( ax,  o.center(1), o.center(2), o.center(3) + 0.5, num2str( o.id ), ...
+        o = epicells( i );
+        hts(i ) = text( ax,  ...
+            o.center(1), o.center(2), o.center(3) + 0.5, ...
+            num2str( o.id ), ...
             'HorizontalAlignment', 'center', ...
             'VerticalAlignment', 'middle' );
     end

src/add_plot_scalebar.m → src/@deproj/add_plot_scalebar.m

-function [ hsb, ht ] = add_plot_scalebar( ax, epicells, length, units )
+function [ hsb, ht ] = add_plot_scalebar( obj, length, ax )
 %ADD_PLOT_SCALEBAR Add a scale-bar to the plot.
 %   We used the collection of epicells to determine a sensible position.
 
     POS_RATIO = 0.10; % of width to position the bar
     THICKNESS_RATIO = 0.02;
-
+    
+    epicells = obj.epicells;
+    units = obj.units;
+    
     p = vertcat( epicells.boundary );
     minp = double( min( p ) );
     maxp = double( max( p ) );

src/add_plot_variable.m → src/@deproj/add_plot_variable.m

-function hts = add_plot_variable( ax, boundaries, values )
+function hts = add_plot_variable( obj, values, ax )
 %ADD_PLOT_VARIABLE Plots the boundaries as patches, colored by the specified values.
     
-    n_objects = numel( boundaries );
+    epicells = obj.epicells;
+    boundaries = { epicells.boundary };
+    n_objects = numel( boundaries );    
     hts = NaN( n_objects, 1 );
     
     for i = 1 :  n_objects

src/@deproj/deproj.m

+classdef deproj
+    %DEPROJ Manage a collection of epicells.
+    
+    properties
+        epicells
+        units
+    end
+    
+    methods
+        function obj = deproj( epicells, units )
+            obj.epicells = epicells;
+            obj.units = units;
+        end
+    end
+    
+    %% Plotting routines.
+    methods
+        
+        %% Generate figures.
+        
+        % Figure with the local plan orientation for a collection of epicells.
+        [ hf, ax1, ax2, ax3 ] = plot_fit_plane( obj, scale_bar_length )
+        
+        % Plot the 2D ellipses on the tissue surface.
+        [ hf, hc, he ] = plot_fit_ellipse( obj, scale_bar_length )
+        
+        %% Helpers.
+        % They are public in case of.
+        
+        % Add the tissue plot colored with the 1st euler angle.
+        hts = add_plot_euler_alpha( obj, ax )
+        
+        % Add the tissue plot colored with the 2nd euler angle.
+        hts = add_plot_euler_beta( obj, ax )
+        
+        % Add the tissue plot colored with the 3rd euler angle.
+        hts = add_plot_euler_gamma( obj, ax )
+        
+        % Add the epicell ids to the specified plot axes.
+        hts = add_plot_id( obj,  ax )
+        
+        % Add a scale-bar to the plot.
+        [ hsb, ht ] = add_plot_scalebar( obj, length, ax )
+        
+        % Plots the boundaries as patches, colored by the specified values.
+        hts = add_plot_variable( obj, values, ax )
+    end
+    
+     %% Public static methods: builders.
+    methods ( Access = public, Hidden = false, Static = true )
+        
+        % Returns the Z position of points taken from a height-map.
+        obj = from_heightmap( I, ...
+            H, ...
+            pixel_size, ...
+            voxel_depth, ...
+            units, ...
+            invert_z, ...
+            inpaint_zeros, ...
+            prune_zeros );
+    end
+    
+    %% Private static methods: utilities.
+    methods ( Access = private, Hidden = true, Static = true )
+        
+        % Returns the Z position of points taken from a height-map.
+        z_coords = get_z( P, H, pixel_size, voxel_depth )
+        
+        % Returns the cells from a BW image with ridges.
+        [ objects, junction_graph ] = mask_to_objects( I, downsample )
+        
+    end
+end
+

src/@deproj/from_heightmap.m

+function obj = from_heightmap( ...
+    I, ...
+    H, ...
+    pixel_size, ...
+    voxel_depth, ...
+    units, ...
+    invert_z, ...    
+    inpaint_zeros, ...
+    prune_zeros )
+%FROM_HEIGTHMAP Returns a deproj oebjct built from segmentatio and height-map.
+
+    %% Default values for some parameters.
+
+    if nargin < 8
+        prune_zeros = true;
+    end
+    if nargin < 7
+        inpaint_zeros = true;
+    end
+    if nargin < 6
+        invert_z = false;
+    end
+    if nargin < 5
+        units = 'pixels';
+    end
+    if nargin < 4
+        voxel_depth = 1.;
+    end
+    if nargin < 3
+        pixel_size = 1.;
+    end
+    
+    %% Load objects from segmentation mask.
+    
+    fprintf('Converting mask to objects.\n' )
+    tic
+    % We will downsample later: for now we will need all the pixels to have
+    % a robust estimate of the ellipse and plane fit.
+    downsample = false;
+    [ objects, junction_graph ] = deproj.mask_to_objects( I, downsample );
+    n_objects = numel( objects );
+    n_junctions = numel( junction_graph.Nodes );
+    fprintf('Converted a %d x %d mask in %.1f seconds. Found %d objects and %d junctions.\n', ...
+        size(I,1), size(I,2), toc, n_objects, n_junctions )
+    
+    
+    %% Compute some scale of an object radius.
+    
+    % Median area (pixel units)
+    med_area = median(arrayfun( @(s) polyarea( s.boundary(:,1), s.boundary(:,2) ),  ...
+        objects ) ); % pixels^2
+    
+    object_scale = 2 * sqrt( med_area )/ pi; % pixel units
+    
+    fprintf('Typical object scale: %.1f pixels or %.2f µm.\n', ...
+        object_scale, object_scale * pixel_size )
+    
+    
+    %% Scale to physical coordinates.
+    
+    % Scale junction XY coordinates to µm.
+    junction_graph.Nodes.Centroid( :, 1:2 ) = junction_graph.Nodes.Centroid( :, 1:2 ) * pixel_size;
+    
+    % Scale object XY boundary to µm.
+    for i = 1 : n_objects
+        objects( i ).boundary   = objects( i ).boundary * pixel_size;
+        objects( i ).center     = objects( i ).center * pixel_size;
+    end
+    
+    %% Load the height-map and add the Z coordinates to objects.
+    
+    fprintf('Collecting Z coordinates.\n' )
+    tic
+    
+    if inpaint_zeros
+        mask = H == 0;
+        H = regionfill( H, mask );
+    end
+    
+    if prune_zeros
+        H( H == 0 ) = NaN;
+    end
+    
+    % Smooth the height-map over a scale smaller than a cell.
+    H = imgaussfilt( H, object_scale );
+    
+    % For junction.
+    z_junction = deproj.get_z( junction_graph.Nodes.Centroid, H, pixel_size, voxel_depth );
+    
+    % For objects.
+    for i = 1 : n_objects
+        z_obj = deproj.get_z( objects( i ).boundary, H, pixel_size, voxel_depth );
+        objects( i ).boundary = [ objects( i ).boundary z_obj ];
+        objects( i ).center(3) = mean( z_obj );
+    end
+    
+    fprintf('Done in %.1f seconds.\n', toc )
+    
+    
+    %% Possibly remove the junctions and cells found at z = 0.
+    
+    if prune_zeros
+        
+        bad_junctions = find( z_junction == 0 | isnan( z_junction ) );
+        junction_graph = junction_graph.rmnode( bad_junctions );
+        z_junction( bad_junctions ) = [];
+        
+        fprintf( 'Removed %d junctions at Z=0.\n', numel( bad_junctions ) )
+        
+        objects_to_prune = false( n_objects, 1 );
+        for i = 1 : numel( objects )
+            o = objects( i );
+            if any( isnan( o.boundary(:) ) ) || ~isempty( intersect( bad_junctions, o.junctions ) )
+                objects_to_prune( i ) = true;
+            end
+        end
+        
+        objects( objects_to_prune ) = [];
+        fprintf( 'Removed %d objects touching these junctions.\n', sum( objects_to_prune ) )
+        
+        n_objects = numel( objects );
+        
+    end
+    
+    %% Invert Z position for plotting.
+    
+    if invert_z
+        max_z = max( z_junction );
+        z_junction = max_z -  z_junction;
+        junction_graph.Nodes.Centroid = [ junction_graph.Nodes.Centroid z_junction ];
+        
+        for i = 1 : n_objects
+            objects( i ).boundary( :, 3 ) = max_z - objects( i ).boundary( :, 3 );
+            objects( i ).center( 3 ) = max_z - objects( i ).center( 3 );
+        end
+    end
+    
+    %% Create epicell instances.
+    
+    fprintf('Computing morphological descriptors of all objects.\n' )
+    tic
+    
+    epicells = repmat( epicell, n_objects, 1 );
+    for i = 1 : n_objects
+        o = objects( i );
+        epicells( i ) = epicell( o.boundary, o.junctions, i  );
+    end
+    
+    obj = deproj( epicells, units );
+    
+    fprintf('Done in %.1f seconds.\n', toc )
+    
+end

src/get_z.m → src/@deproj/get_z.m

 function z_coords = get_z( P, H, pixel_size, voxel_depth )
-%GET_Z Returns the Z position of points taken from a height-map
+%GET_Z Returns the Z position of points taken from a height-map.
 %   - P is a Nx2 list of points, in physical coordinates.
 %   - H is the height map, encoding the z plane of interest for all X & Y.
 %   - pixel_size: convert pixel coordinates to physical coordinates.

src/plot_fit_ellipse.m → src/@deproj/plot_fit_ellipse.m

-function [ hf, hc, he ] = plot_fit_ellipse( epicells )
-%PLOT_FIT_ELLIPSE Summary of this function goes here
+function [ hf, hc, he ] = plot_fit_ellipse( obj, scale_bar_length )
+%PLOT_FIT_ELLIPSE Plot the 2D ellipses on the tissue surface.
+
+    epicells = obj.epicells;
+    if nargin < 2
+        scale_bar_length = 10;
+    end
 
     hf = figure( 'Position', [ 1204 20 600 500 ] );
     hold on
@@ -19,7 +24,7 @@ function [ hf, hc, he ] = plot_fit_ellipse( epicells )
     end
 
     set( he, 'Color', 'k' )
-    add_plot_scalebar( gca,  epicells, 10, 'µm' );
+    add_plot_scalebar(  obj, scale_bar_length, gca );
 
 end
 

src/plot_fit_plane.m → src/@deproj/plot_fit_plane.m

-function [ hf, ax1, ax2, ax3 ] = plot_fit_plane( epicells )
+function [ hf, ax1, ax2, ax3 ] = plot_fit_plane( obj, scale_bar_length )
 %PLOT_FIT_PLANE Figure with the local plan orientation for a collection of epicells.
 
+    if nargin < 2
+        scale_bar_length = 10;
+    end
+
     hf = figure( 'Position', [ 1204 20 600 1000 ] );
     
     ax1 = subplot( 3, 1, 1 );
     hold on
     axis equal
-    add_plot_euler_alpha( ax1, epicells );
+    add_plot_euler_alpha( obj, ax1 );
     colormap( ax1, 'hsv' )
     colorbar(ax1)
 
     ax2 = subplot( 3, 1, 2 );
     hold on
     axis equal
-    add_plot_euler_beta( ax2, epicells );
+    add_plot_euler_beta( obj, ax2 );
     colorbar(ax2)
 
     ax3 = subplot( 3, 1, 3 );
     hold on
     axis equal
-    add_plot_euler_gamma( ax3, epicells );
+    add_plot_euler_gamma( obj, ax3 );
     colormap( ax3, 'hsv' )
     colorbar(ax3)   
     
-    add_plot_scalebar( ax3,  epicells, 10, 'µm' );
+    add_plot_scalebar( obj, scale_bar_length, ax3 );
     linkaxes( [ ax3 ax2 ax1 ] )
     
     axis( ax1, 'off' )

src/@epicell/epicell.m

@@ -43,7 +43,7 @@ classdef epicell
             % Morphological descriptors on non-downsampled boundary.
             p = epicell.centered_points( boundary );
             obj.euler_angles = epicell.fit_plane( p );
-            obj.ellipse_fit = fit_ellipse( p, obj.euler_angles );            
+            obj.ellipse_fit = epicell.fit_ellipse_3d( p, obj.euler_angles );            
         end
         
         function h = plot_patch_2d( obj, val )
@@ -149,8 +149,7 @@ classdef epicell
         
     end
     
-    %% Public static methods: utilities.
-    
+    %% Public static methods: utilities.    
     methods ( Access = public, Hidden = false, Static = true )
         % Convert euler angles to rotation matrix.
         R = euleurZXZ2rot( E )

src/@epicell/fit_ellipse_3d.m

@@ -28,7 +28,7 @@ function [ f3d, v ] = fit_ellipse_3d( p, E, method )
     if nargin < 2
         [ ~, ~, v ] = svd( p );
     else
-        v = euleurZXZ2rot( E );
+        v = epicell.euleurZXZ2rot( E );
     end
     
     % Rotate the points into the principal axes frame.