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

Fix the ellipse fit in 3D and its plot.

parent d77f866e
......@@ -5,9 +5,6 @@ function [ f, Q ] = fit_ellipse_2d( p, method )
method = 'direct';
end
c = mean( p );
p = p - repmat( c, size( p, 1 ), 1 );
switch( lower( method ) )
case 'direct'
......@@ -17,9 +14,7 @@ function [ f, Q ] = fit_ellipse_2d( p, method )
end
f = quadratic_to_cartesian2( Q );
f( 1 ) = f( 1 ) + c( 1 );
f( 2 ) = f( 2 ) + c( 2 );
%% Subfunctions
function f = quadratic_to_cartesian2( Q )
......
function [ f, v ] = fit_ellipse_3d( p, E, method )
function [ f3d, v ] = fit_ellipse_3d( p, E, method )
%FIT_ELLIPSE_3D Fit a 2D ellipse on 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
% ellipse fit on the projected points. This turns to be much more robust
% than a 3D fit, and also closely match our configuration.
% This function returns f3d = [ x0 y0 z0 a b theta ]
% and v, the ransformation matrix that rotates the 3D points close to the
% XY plane. It will be used to plot the ellipse in 3D.
% Greatly inspired from https://stackoverflow.com/questions/29051168
% /data-fitting-an-ellipse-in-3d-space
......@@ -27,5 +30,9 @@ function [ f, v ] = fit_ellipse_3d( p, E, method )
f = fit_ellipse_2d( pt( :, 1:2 ), method );
f( 1 ) = f( 1 ) + c( 1 );
f( 2 ) = f( 2 ) + c( 2 );
f3d = [ f(1:2) c(3) f(3:5) ];
end
function h = plot_ellipse_3d( f, v, npoints)
function h = plot_ellipse_3d( f3d, v, npoints)
%PLOT_ELLIPSE_3D Plot a 2D ellipse in 3D.
if nargin < 3
npoints = 20;
end
x0 = f(1);
y0 = f(2);
a = f(3);
b = f(4);
theta = f(5);
x0 = f3d(1);
y0 = f3d(2);
z0 = f3d(3);
a = f3d(4);
b = f3d(5);
theta = f3d(6);
R = [ cos( theta ) sin( theta ) ;
-sin( theta ) cos( theta ) ] ;
......@@ -17,16 +18,16 @@ function h = plot_ellipse_3d( f, v, npoints)
t = linspace( 0 , 2 * pi, npoints )';
XY0 = [ a * sin(t), b * cos(t) ];
XY1 = XY0 * R;
xr = XY1( :, 1 ) + x0;
yr = XY1( :, 2 ) + y0;
xr = XY1( :, 1 );
yr = XY1( :, 2 );
zr = zeros( numel( xr ), 1 );
pr = [ xr yr zr ];
% Transform back
pb = pr * v';
xb = pb( :, 1 );
yb = pb( :, 2 );
zb = pb( :, 3 );
xb = pb( :, 1 ) + x0;
yb = pb( :, 2 ) + y0;
zb = pb( :, 3 ) + z0;
xb = [ xb ; xb(1,:) ];
yb = [ yb ; yb(1,:) ];
......
......@@ -12,12 +12,13 @@ function [ hf, hc, he ] = plot_fit_ellipse( epicells )
for i = 1 : n_obj
o = epicells( i );
[ f, v ] = fit_ellipse_3d( double( o.boundary ) );
[ f3d, v ] = fit_ellipse_3d( double( o.boundary ) );
hc( i ) = o.plot_contour_3d;
he( i ) = plot_ellipse_3d( f, v );
he( i ) = plot_ellipse_3d( f3d, v );
end
set( he, 'Color', 'k' )
add_plot_scalebar( gca, epicells, 10, 'µm' );
end
......
Markdown is supported
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