From 6dc23cc82883c8fc62823b8ce86d54f0fdeb9fc6 Mon Sep 17 00:00:00 2001 From: Jean-Yves TINEVEZ <tinevez@pasteur.fr> Date: Mon, 6 Jul 2020 14:12:39 +0200 Subject: [PATCH] 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 | 122 +++--------------- src/{ => @deproj}/add_plot_euler_alpha.m | 5 +- src/{ => @deproj}/add_plot_euler_beta.m | 5 +- src/{ => @deproj}/add_plot_euler_gamma.m | 5 +- src/{ => @deproj}/add_plot_id.m | 9 +- src/{ => @deproj}/add_plot_scalebar.m | 7 +- src/{ => @deproj}/add_plot_variable.m | 6 +- src/@deproj/deproj.m | 74 +++++++++++ src/@deproj/from_heightmap.m | 153 +++++++++++++++++++++++ src/{ => @deproj}/get_z.m | 2 +- src/{ => @deproj}/mask_to_objects.m | 0 src/{ => @deproj}/plot_fit_ellipse.m | 11 +- src/{ => @deproj}/plot_fit_plane.m | 14 ++- src/@epicell/epicell.m | 5 +- src/@epicell/fit_ellipse_3d.m | 2 +- 15 files changed, 291 insertions(+), 129 deletions(-) rename src/{ => @deproj}/add_plot_euler_alpha.m (70%) rename src/{ => @deproj}/add_plot_euler_beta.m (70%) rename src/{ => @deproj}/add_plot_euler_gamma.m (58%) rename src/{ => @deproj}/add_plot_id.m (54%) rename src/{ => @deproj}/add_plot_scalebar.m (91%) rename src/{ => @deproj}/add_plot_variable.m (65%) create mode 100644 src/@deproj/deproj.m create mode 100644 src/@deproj/from_heightmap.m rename src/{ => @deproj}/get_z.m (91%) rename src/{ => @deproj}/mask_to_objects.m (100%) rename src/{ => @deproj}/plot_fit_ellipse.m (61%) rename src/{ => @deproj}/plot_fit_plane.m (72%) diff --git a/RunExample.m b/RunExample.m index 7cf67f2..511f016 100755 --- a/RunExample.m +++ b/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 ); diff --git a/src/add_plot_euler_alpha.m b/src/@deproj/add_plot_euler_alpha.m similarity index 70% rename from src/add_plot_euler_alpha.m rename to src/@deproj/add_plot_euler_alpha.m index 83ccfa7..f20cce6 100644 --- a/src/add_plot_euler_alpha.m +++ b/src/@deproj/add_plot_euler_alpha.m @@ -1,6 +1,7 @@ -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 diff --git a/src/add_plot_euler_beta.m b/src/@deproj/add_plot_euler_beta.m similarity index 70% rename from src/add_plot_euler_beta.m rename to src/@deproj/add_plot_euler_beta.m index a359f4b..0689723 100644 --- a/src/add_plot_euler_beta.m +++ b/src/@deproj/add_plot_euler_beta.m @@ -1,6 +1,7 @@ -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 diff --git a/src/add_plot_euler_gamma.m b/src/@deproj/add_plot_euler_gamma.m similarity index 58% rename from src/add_plot_euler_gamma.m rename to src/@deproj/add_plot_euler_gamma.m index fcab772..2965905 100644 --- a/src/add_plot_euler_gamma.m +++ b/src/@deproj/add_plot_euler_gamma.m @@ -1,8 +1,9 @@ -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 diff --git a/src/add_plot_id.m b/src/@deproj/add_plot_id.m similarity index 54% rename from src/add_plot_id.m rename to src/@deproj/add_plot_id.m index 7345bea..e4acb12 100644 --- a/src/add_plot_id.m +++ b/src/@deproj/add_plot_id.m @@ -1,12 +1,15 @@ -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 diff --git a/src/add_plot_scalebar.m b/src/@deproj/add_plot_scalebar.m similarity index 91% rename from src/add_plot_scalebar.m rename to src/@deproj/add_plot_scalebar.m index 0e11918..56a9d91 100644 --- a/src/add_plot_scalebar.m +++ b/src/@deproj/add_plot_scalebar.m @@ -1,10 +1,13 @@ -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 ) ); diff --git a/src/add_plot_variable.m b/src/@deproj/add_plot_variable.m similarity index 65% rename from src/add_plot_variable.m rename to src/@deproj/add_plot_variable.m index 6573e1f..1a89b60 100644 --- a/src/add_plot_variable.m +++ b/src/@deproj/add_plot_variable.m @@ -1,7 +1,9 @@ -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 diff --git a/src/@deproj/deproj.m b/src/@deproj/deproj.m new file mode 100644 index 0000000..caec79a --- /dev/null +++ b/src/@deproj/deproj.m @@ -0,0 +1,74 @@ +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 + diff --git a/src/@deproj/from_heightmap.m b/src/@deproj/from_heightmap.m new file mode 100644 index 0000000..54e8408 --- /dev/null +++ b/src/@deproj/from_heightmap.m @@ -0,0 +1,153 @@ +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 diff --git a/src/get_z.m b/src/@deproj/get_z.m similarity index 91% rename from src/get_z.m rename to src/@deproj/get_z.m index 0904e2f..9f23e3d 100644 --- a/src/get_z.m +++ b/src/@deproj/get_z.m @@ -1,5 +1,5 @@ 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. diff --git a/src/mask_to_objects.m b/src/@deproj/mask_to_objects.m similarity index 100% rename from src/mask_to_objects.m rename to src/@deproj/mask_to_objects.m diff --git a/src/plot_fit_ellipse.m b/src/@deproj/plot_fit_ellipse.m similarity index 61% rename from src/plot_fit_ellipse.m rename to src/@deproj/plot_fit_ellipse.m index e459166..7fe5d27 100644 --- a/src/plot_fit_ellipse.m +++ b/src/@deproj/plot_fit_ellipse.m @@ -1,5 +1,10 @@ -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 diff --git a/src/plot_fit_plane.m b/src/@deproj/plot_fit_plane.m similarity index 72% rename from src/plot_fit_plane.m rename to src/@deproj/plot_fit_plane.m index de98557..c3d488a 100644 --- a/src/plot_fit_plane.m +++ b/src/@deproj/plot_fit_plane.m @@ -1,29 +1,33 @@ -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' ) diff --git a/src/@epicell/epicell.m b/src/@epicell/epicell.m index 2040133..7a5fdd1 100644 --- a/src/@epicell/epicell.m +++ b/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 ) diff --git a/src/@epicell/fit_ellipse_3d.m b/src/@epicell/fit_ellipse_3d.m index fe153c4..5f54f98 100644 --- a/src/@epicell/fit_ellipse_3d.m +++ b/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. -- GitLab