diff --git a/RunExample.m b/RunExample.m index e317d9b80138323251b46d5c5b71b09b853f5454..d39626755a2c41447a40002fe9cadc4c383396c8 100755 --- a/RunExample.m +++ b/RunExample.m @@ -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; diff --git a/src/area3d.m b/src/area3d.m deleted file mode 100644 index b47b52d9e7cf91b781d68a7eeba94ebf5d52dd40..0000000000000000000000000000000000000000 --- a/src/area3d.m +++ /dev/null @@ -1,40 +0,0 @@ -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 - diff --git a/src/centered_points.m b/src/centered_points.m deleted file mode 100644 index 29f7cbb4f38b34478203ebdac208766e1c5302c8..0000000000000000000000000000000000000000 --- a/src/centered_points.m +++ /dev/null @@ -1,12 +0,0 @@ -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 - diff --git a/src/epicell.m b/src/epicell.m new file mode 100644 index 0000000000000000000000000000000000000000..c353d9d0f4fa4ea5f38b750d581e057b0fc067ba --- /dev/null +++ b/src/epicell.m @@ -0,0 +1,110 @@ +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 + diff --git a/src/fit_ellipse.m b/src/fit_ellipse.m index 2922dccd7f713476af12580c0a5b6e85c38c6f5d..4b60f7a5bf09a4856f04c98cb920403e558162fd 100644 --- a/src/fit_ellipse.m +++ b/src/fit_ellipse.m @@ -1,4 +1,4 @@ -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; diff --git a/src/perimeter3d.m b/src/perimeter3d.m deleted file mode 100644 index d026425f54ba2841f3759b3484957ea29542cfd1..0000000000000000000000000000000000000000 --- a/src/perimeter3d.m +++ /dev/null @@ -1,23 +0,0 @@ -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 -