Commit d0fe0d60 authored by Jean-Yves TINEVEZ's avatar Jean-Yves TINEVEZ
Browse files

Use a data class to store morphological parameters of a cell.

parent b37c53f0
......@@ -135,32 +135,17 @@ 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
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
%% Compute object de-projected area.
%% Create epicell instances.
for i = 1 : n_objects
epicells = repmat( epicell, n_objects, 1 );
for i = 1 : n_objects
o = objects( i );
[ area, uncorr_area ] = area3d( o );
[ perim, uncorr_perim ] = perimeter3d( o );
[ f, E ] = fit_ellipse( o );
objects( i ).id = i;
objects( i ).area = area;
objects( i ).perimeter = perim;
objects( i ).euler_angles = E;
objects( i ).ellipse_fit = f;
uncorr = struct();
uncorr.area = uncorr_area;
uncorr.perimeter = uncorr_perim;
objects( i ).uncorr = uncorr;
epicells( i ) = epicell( o.boundary, o.junctions, i );
end
......@@ -200,7 +185,7 @@ axis equal
for i = 1 : n_objects
o = objects( i );
o = epicells( i );
P = o.boundary;
% err = o.perimeter / o.uncorr.perimeter - 1;
......
function [ area, uncorr_area ] = area3d( o )
%AREA3D Computes the area of the object.
%% Deprojected 3D version.
% Put all vertex coordinates with respect to center.
p = centered_points( o );
n_vertices = size( p, 1 );
% Build small triangles.
index = [ 2 : n_vertices 1 ];
p1 = p;
p2 = p( index, : );
% Cross product.
cp = cross( p1, p2 );
% Norm of each vector.
vn = euclidean_norm( cp );
% Positive area.
area_triangle = abs( vn );
% Total positive area.
area = sum( area_triangle ) / 2;
%% 2D area.
uncorr_area = polyarea( o.boundary(:,1), o.boundary(:,2) );
%% Subfunction.
function n = euclidean_norm( v )
n = sqrt( sum( v .* v, ndims( v ) ) );
end
end
function p = centered_points( o )
%CENTERED_POINTS Returns the 3D coordinates of the object bounds, with
%respect to its center.
p = o.boundary;
n_vertices = size( p ,1 );
center = mean( p );
center = repmat( center, [ n_vertices, 1 ] );
p = p - center;
end
classdef epicell
%EPICELL Data class that store data resulting from tissue cell segmentation.
% Immutable class: all properties are read-only.
properties (SetAccess = private)
boundary
center
junction_ids
area
perimeter
euler_angles
ellipse_fit
uncorrected_area
uncorrected_perimeter
id
end
methods
function obj = epicell( boundary, junction_ids, id )
%EPICELL Construct an epicell from a N x 3 list of points.
if nargin == 0
return
end
% Base properties.
obj.boundary = boundary;
obj.junction_ids = junction_ids;
obj.id = id;
obj.center = mean( boundary );
% Morphological descriptors.
p = epicell.centered_points( boundary );
[ obj.area, obj.uncorrected_area ] = epicell.area3d( p );
[ obj.perimeter, obj.uncorrected_perimeter ] = epicell.perimeter3d( p );
obj.euler_angles = epicell.fit_plane( p );
obj.ellipse_fit = fit_ellipse( p, obj.euler_angles );
end
end
%% Static methods: compute final properties value.
methods ( Access = private, Hidden = true, Static = true )
function p = centered_points( p )
%CENTERED_POINTS Returns the 3D coordinates of the object bounds, with
%respect to its center.
n_vertices = size( p ,1 );
center = mean( p );
center = repmat( center, [ n_vertices, 1 ] );
p = double( p - center );
end
function [ area, uncorr_area ] = area3d( p )
%AREA3D Computes the area of the object.
n_vertices = size( p, 1 );
% Build small triangles.
index = [ 2 : n_vertices 1 ];
p1 = p;
p2 = p( index, : );
% Cross product.
cp = cross( p1, p2 );
% Norm of each vector.
vn = euclidean_norm( cp );
% Positive area.
area_triangle = abs( vn );
% Total positive area.
area = sum( area_triangle ) / 2;
% 2D area.
uncorr_area = polyarea( p(:,1), p(:,2) );
function n = euclidean_norm( v )
n = sqrt( sum( v .* v, ndims( v ) ) );
end
end
function [ perim, uncorr_perim ] = perimeter3d( p )
%PERIMETER3D Perimeter of a closed N-dimensional polygon.
perim = compute_perim( p );
uncorr_perim = compute_perim( p( : , 1:2 ) );
function l_perim = compute_perim( p )
% p can be a N x d matrix, with d being the dimensionality.
p2 = [ p ; p( 1, : ) ];
p_diff = diff( p2 );
p_diff_2 = p_diff .* p_diff;
p_diff_2_sum = sum( p_diff_2, 2 );
sls = sqrt( p_diff_2_sum );
l_perim = sum( sls );
end
end
function E = fit_plane( p )
% Fit a plane to these points.
[ ~, ~, v ] = svd( p );
E = rot2eulerZXZ( v );
end
end
end
function [ f, E ] = fit_ellipse( o )
function f = fit_ellipse( p, E )
%FIT_ELLIPSE Fit a 2D ellipse to the 3D points.
% The fit requires the Euler angles of the plane fitted through the
% opints, so that we can project them on this plane. We then make a 2D
......@@ -8,11 +8,12 @@ function [ f, E ] = fit_ellipse( o )
% Greatly inspired from
% https://stackoverflow.com/questions/29051168/data-fitting-an-ellipse-in-3d-space
p = double( centered_points( o ) );
% Fit a plane to these points.
[ ~, ~, v ] = svd( p );
E = rot2eulerZXZ( v );
if nargin < 2
[ ~, ~, v ] = svd( p );
else
v = euleurZXZ2rot( E );
end
% Rotate the points into the principal axes frame.
p = p * v;
......
function [ perim, uncorr_perim ] = perimeter3d( o )
%PERIMETER3D Perimeter of a closed N-dimensional polygon.
perim = compute_perim( o.boundary );
uncorr_perim = compute_perim( o.boundary( : , 1:2 ) );
%% Subfunction
function l_perim = compute_perim( p )
% p can be a N x d matrix, with d being the dimensionality.
p2 = [ p ; p( 1, : ) ];
p_diff = diff( p2 );
p_diff_2 = p_diff .* p_diff;
p_diff_2_sum = sum( p_diff_2, 2 );
sls = sqrt( p_diff_2_sum );
l_perim = sum( sls );
end
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment