diff --git a/src/epicell.m b/src/@epicell/epicell.m
similarity index 85%
rename from src/epicell.m
rename to src/@epicell/epicell.m
index 3f2d48f6a23ff49146bddbb21c9777136e7d94e7..204013355ee61537a660e576c08922592e034dc3 100644
--- a/src/epicell.m
+++ b/src/@epicell/epicell.m
@@ -75,7 +75,7 @@ classdef epicell
end
end
- %% Static methods: compute final properties value.
+ %% Private static methods: compute final properties value.
methods ( Access = private, Hidden = true, Static = true )
function p = centered_points( p )
@@ -142,8 +142,34 @@ classdef epicell
p = p - repmat( c, size( p, 1 ), 1 );
% Fit a plane to these points.
[ ~, ~, v ] = svd( p );
- E = rot2eulerZXZ( v );
+ E = epicell.rot2eulerZXZ( v );
end
+
+
+
+ end
+
+ %% Public static methods: utilities.
+
+ methods ( Access = public, Hidden = false, Static = true )
+ % Convert euler angles to rotation matrix.
+ R = euleurZXZ2rot( E )
+
+ % Convert rotation matrix to euler angles.
+ [ E, E_deg ] = rot2eulerZXZ( R )
+
+ % Fit a 2D ellipse to a set of 2D points.
+ [ f, Q ] = fit_ellipse_2d( p, method )
+
+ % Fit a 2D ellipse to a set of 3D points.
+ [ f3d, v ] = fit_ellipse_3d( p, E, method )
+
+ % Plot a 2D ellipse in 3D.
+ h = plot_ellipse_3d( f3d, v, npoints )
+
+ % Plot an ellipse in XY plane.
+ h = plot_ellipse_2d( f, npoints )
+
end
end
diff --git a/src/euleurZXZ2rot.m b/src/@epicell/euleurZXZ2rot.m
similarity index 100%
rename from src/euleurZXZ2rot.m
rename to src/@epicell/euleurZXZ2rot.m
diff --git a/src/fit_ellipse_2d.m b/src/@epicell/fit_ellipse_2d.m
similarity index 100%
rename from src/fit_ellipse_2d.m
rename to src/@epicell/fit_ellipse_2d.m
diff --git a/src/@epicell/fit_ellipse_3d.m b/src/@epicell/fit_ellipse_3d.m
new file mode 100644
index 0000000000000000000000000000000000000000..fe153c4e77dc465255c1a1af02e13bf3f5492228
--- /dev/null
+++ b/src/@epicell/fit_ellipse_3d.m
@@ -0,0 +1,44 @@
+function [ f3d, v ] = fit_ellipse_3d( p, E, method )
+%FIT_ELLIPSE_3D Fit a 2D ellipse to a set of 3D points.
+% The fit requires (or compute) 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 ]. These are the
+% cartesian parameters of the ellipse in the rotated plane. The ellipse
+% semi-major axis and semi-minor axis are always so that a >= b and theta
+% measure the angle of the semi-major axis with respect to the X axis (in
+% the rotated plane).
+%
+% v is the transformation 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
+
+ if nargin < 3
+ method = 'direct';
+ end
+
+ c = mean( p );
+ p = p - repmat( c, size( p, 1 ), 1 );
+
+ % Fit a plane to these points.
+ if nargin < 2
+ [ ~, ~, v ] = svd( p );
+ else
+ v = euleurZXZ2rot( E );
+ end
+
+ % Rotate the points into the principal axes frame.
+ pt = p * v;
+
+ f = epicell.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
+
diff --git a/src/plot_ellipse_2d.m b/src/@epicell/plot_ellipse_2d.m
similarity index 100%
rename from src/plot_ellipse_2d.m
rename to src/@epicell/plot_ellipse_2d.m
diff --git a/src/plot_ellipse_3d.m b/src/@epicell/plot_ellipse_3d.m
similarity index 100%
rename from src/plot_ellipse_3d.m
rename to src/@epicell/plot_ellipse_3d.m
diff --git a/src/rot2eulerZXZ.m b/src/@epicell/rot2eulerZXZ.m
similarity index 94%
rename from src/rot2eulerZXZ.m
rename to src/@epicell/rot2eulerZXZ.m
index 7ef58205c0fd9af352de9c1f5dac058a4865cb61..1e4c3e712ddf32ecf37c649fea554036a0b51ffe 100644
--- a/src/rot2eulerZXZ.m
+++ b/src/@epicell/rot2eulerZXZ.m
@@ -1,4 +1,4 @@
-function[E,E_deg]=rot2eulerZXZ(R)
+function[ E, E_deg ] = rot2eulerZXZ( R )
%ROT2EULERZXZ Convert rotation matrix to ZX'Z' Euler angles.
r13 = R( 1, 3 );
diff --git a/src/fit_ellipse.m b/src/fit_ellipse.m
deleted file mode 100644
index 4b60f7a5bf09a4856f04c98cb920403e558162fd..0000000000000000000000000000000000000000
--- a/src/fit_ellipse.m
+++ /dev/null
@@ -1,137 +0,0 @@
-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
-% ellipse fit on the projected points. This turns to be much more robust
-% than a 3D fit, and also closely match our configuration.
-
-% Greatly inspired from
-% https://stackoverflow.com/questions/29051168/data-fitting-an-ellipse-in-3d-space
-
- % Fit a plane to these points.
- if nargin < 2
- [ ~, ~, v ] = svd( p );
- else
- v = euleurZXZ2rot( E );
- end
-
- % Rotate the points into the principal axes frame.
- p = p * v;
-
- % Direct ellipse fit.
-% A = direct_ellipse_fit( p( :, 1:2 ) );
- A = taubin_ellipse_fit( p( :, 1:2 ) );
- f = quadratic_to_cartesian( A );
-
- %% Subfunctions
-
- function f = quadratic_to_cartesian( A )
- % Equations taken from Wolfram website.
-
- a = A(1);
- b = A(2);
- c = A(3);
- d = A(4);
- f = A(5);
- g = A(6);
-
- x0 = ( c * d - b * f ) / ( b^2 - a * c );
- y0 = ( a * f - b * d ) / ( b^2 - a * c );
-
- l1 = sqrt( 2*(a*f^2+c*d^2+g*b^2-2*b*d*f-a*c*g) / ((b^2-a*c)*(sqrt((a-c)^2+4*b^2)-(a+c))));
- l2 = sqrt( 2*(a*f^2+c*d^2+g*b^2-2*b*d*f-a*c*g) / ((b^2-a*c)*(-sqrt((a-c)^2+4*b^2)-(a+c))));
-
- if b == 0 && a < c
- phi = 0;
- elseif b == 0 && a > c
- phi = 0.5*pi;
- elseif b ~= 0 && a < c
- phi = 0.5* acot((a-c)/(2*b));
- else
- phi = 0.5*pi + 0.5* acot((a-c)/(2*b));
- end
-
- f = [ x0 y0 l1 l2 phi ];
-
- end
-
- function A = direct_ellipse_fit(XY) %#ok<DEFNU>
- % Direct ellipse fit, proposed in article
- % A. W. Fitzgibbon, M. Pilu, R. B. Fisher
- % "Direct Least Squares Fitting of Ellipses"
- % IEEE Trans. PAMI, Vol. 21, pages 476-480 (1999)
- %
- % Adapted from https://fr.mathworks.com/matlabcentral/fileexchange/22684-ellipse-fit-direct-method
-
- centroid = mean(XY); % the centroid of the data set
- D1 = [(XY(:,1)-centroid(1)).^2, (XY(:,1)-centroid(1)).*(XY(:,2)-centroid(2)),...
- (XY(:,2)-centroid(2)).^2];
- D2 = [XY(:,1)-centroid(1), XY(:,2)-centroid(2), ones(size(XY,1),1)];
- S1 = D1'*D1;
- S2 = D1'*D2;
- S3 = D2'*D2;
- T = -inv(S3)*S2';
- M = S1 + S2*T;
- M = [M(3,:)./2; -M(2,:); M(1,:)./2];
- [evec,eval] = eig(M); %#ok<ASGLU>
- cond = 4*evec(1,:).*evec(3,:)-evec(2,:).^2;
- A1 = evec(:,find(cond>0)); %#ok<FNDSB>
- A = [A1; T*A1];
- A4 = A(4)-2*A(1)*centroid(1)-A(2)*centroid(2);
- A5 = A(5)-2*A(3)*centroid(2)-A(2)*centroid(1);
- A6 = A(6)+A(1)*centroid(1)^2+A(3)*centroid(2)^2+...
- A(2)*centroid(1)*centroid(2)-A(4)*centroid(1)-A(5)*centroid(2);
- A(4) = A4; A(5) = A5; A(6) = A6;
- A = A/norm(A);
- end
-
- function A = taubin_ellipse_fit(XY)
- % Ellipse fit by Taubin's Method published in
- % G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
- % Space Curves Defined By Implicit Equations, With
- % Applications To Edge And Range Image Segmentation",
- % IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
- %
- % Input: XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)
- %
- % Output: A = [a b c d e f]' is the vector of algebraic
- % parameters of the fitting ellipse:
- % ax^2 + bxy + cy^2 +dx + ey + f = 0
- % the vector A is normed, so that ||A||=1
- %
- % Among fast non-iterative ellipse fitting methods,
- % this is perhaps the most accurate and robust
- %
- % Note: this method fits a quadratic curve (conic) to a set of points;
- % if points are better approximated by a hyperbola, this fit will
- % return a hyperbola. To fit ellipses only, use "Direct Ellipse Fit".
-
- centroid = mean(XY); % the centroid of the data set
- Z = [(XY(:,1)-centroid(1)).^2, (XY(:,1)-centroid(1)).*(XY(:,2)-centroid(2)),...
- (XY(:,2)-centroid(2)).^2, XY(:,1)-centroid(1), XY(:,2)-centroid(2), ones(size(XY,1),1)];
- M = Z'*Z/size(XY,1);
- P = [M(1,1)-M(1,6)^2, M(1,2)-M(1,6)*M(2,6), M(1,3)-M(1,6)*M(3,6), M(1,4), M(1,5);
- M(1,2)-M(1,6)*M(2,6), M(2,2)-M(2,6)^2, M(2,3)-M(2,6)*M(3,6), M(2,4), M(2,5);
- M(1,3)-M(1,6)*M(3,6), M(2,3)-M(2,6)*M(3,6), M(3,3)-M(3,6)^2, M(3,4), M(3,5);
- M(1,4), M(2,4), M(3,4), M(4,4), M(4,5);
- M(1,5), M(2,5), M(3,5), M(4,5), M(5,5)];
- Q = [4*M(1,6), 2*M(2,6), 0, 0, 0;
- 2*M(2,6), M(1,6)+M(3,6), 2*M(2,6), 0, 0;
- 0, 2*M(2,6), 4*M(3,6), 0, 0;
- 0, 0, 0, 1, 0;
- 0, 0, 0, 0, 1];
- [V,D] = eig(P,Q);
- [Dsort,ID] = sort(diag(D)); %#ok<ASGLU>
- A = V(:,ID(1));
- A = [A; -A(1:3)'*M(1:3,6)];
- A4 = A(4)-2*A(1)*centroid(1)-A(2)*centroid(2);
- A5 = A(5)-2*A(3)*centroid(2)-A(2)*centroid(1);
- A6 = A(6)+A(1)*centroid(1)^2+A(3)*centroid(2)^2+...
- A(2)*centroid(1)*centroid(2)-A(4)*centroid(1)-A(5)*centroid(2);
- A(4) = A4; A(5) = A5; A(6) = A6;
- A = A/norm(A);
- end % Taubin
-
-
-end
-
diff --git a/src/fit_ellipse_3d.m b/src/fit_ellipse_3d.m
deleted file mode 100644
index 78ece117e5a6a6732fb9613809addf3c05e357fa..0000000000000000000000000000000000000000
--- a/src/fit_ellipse_3d.m
+++ /dev/null
@@ -1,38 +0,0 @@
-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
-
- if nargin < 3
- method = 'direct';
- end
-
- c = mean( p );
- p = p - repmat( c, size( p, 1 ), 1 );
-
- % Fit a plane to these points.
- if nargin < 2
- [ ~, ~, v ] = svd( p );
- else
- v = euleurZXZ2rot( E );
- end
-
- % Rotate the points into the principal axes frame.
- pt = p * v;
-
- 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
-