Open modenaxe opened 3 years ago
@renaultJB I did not check this on GIBOC-core before, but the SphereFit is actually a linear least-square sphere fitting. See page 15 of this report: https://eprintspublications.npl.co.uk/5050/1/DITC140.pdf If data are accurate, the fitting will be accurate, otherwise it is better be used to initialise a proper non-linear model. This is ok for now but must be improved in the future.
Yes indeed, I have used it while developing because it was faster and forgot to change it afterwards, but there is an implementation in the non linear least square fit library provided by the functions associated with the document you mention, the code is already in the repo non linear least square sphere fit We could wrap it inside the SphereFit.m to avoid touching other parts of the code. Also, as you suggested we could use the linear least square sphere for initialisation, Here is an untested proposition (I'am not sure about the tolerance parameters and defaulting back to linear least squares in case convergence is not reached) :
function [Center,Radius,ErrorDist] = sphereFit(X)
% This fits a sphere in two steps
%
% Step 1 : Fit a linear least square sphere
%---------------------------------------------------------------------------
% this fits a sphere to a collection of data using a closed form for the
% solution (opposed to using an array the size of the data set).
% Minimizes Sum((x-xc)^2+(y-yc)^2+(z-zc)^2-r^2)^2
% x,y,z are the data, xc,yc,zc are the sphere's center, and r is the radius
% Assumes that points are not in a singular configuration, real numbers, ...
% if you have coplanar data, use a circle fit with svd for determining the
% plane, recommended Circle Fit (Pratt method), by Nikolai Chernov
% http://www.mathworks.com/matlabcentral/fileexchange/22643
% Input:
% X: n x 3 matrix of cartesian data
% Outputs:
% Center: Center of sphere
% Radius: Radius of sphere
% Author:
% Alan Jennings, University of Dayton
% Modified to add distance to sphere -> ErrorDist
% ---------------------------------------------------------------------------
% Step 2 : Fit a non linear least square sphere
%---------------------------------------------------------------------------
%
% LSSPHERE.M Least-squares sphere using Gauss-Newton.
%
% Version 1.0
% Last amended I M Smith 27 May 2002.
% Created I M Smith 08 Mar 2002
%
% ---------------------------------------------------------------------
% See A B Forbes: Least-squares best-fit geometric elements,
% NPL Report DITC 140/89.
%
% ---------------------------------------------------------------------
% Author A B Forbes
% National Physical Laboratory
% England
%
% Created November 1988
% Version 2.0 93/07/09
% Crown Copyright
% ---------------------------------------------------------------------------
%% STEP 1
A=[mean(X(:,1).*(X(:,1)-mean(X(:,1)))), ...
2*mean(X(:,1).*(X(:,2)-mean(X(:,2)))), ...
2*mean(X(:,1).*(X(:,3)-mean(X(:,3)))); ...
0, ...
mean(X(:,2).*(X(:,2)-mean(X(:,2)))), ...
2*mean(X(:,2).*(X(:,3)-mean(X(:,3)))); ...
0, ...
0, ...
mean(X(:,3).*(X(:,3)-mean(X(:,3))))];
A=A+A.';
B=[mean((X(:,1).^2 + X(:,2).^2 + X(:,3).^2).*(X(:,1) - mean(X(:,1))));...
mean((X(:,1).^2 + X(:,2).^2 + X(:,3).^2).*(X(:,2) - mean(X(:,2))));...
mean((X(:,1).^2 + X(:,2).^2 + X(:,3).^2).*(X(:,3) - mean(X(:,3))))];
Center0=(A\B).';
Radius0 = sqrt(mean(sum([X(:,1)-Center0(1), X(:,2)-Center0(2), X(:,3)-Center0(3)].^2,2)));
% (P-Centre)^2 - Radius^2
ErrorDist0 = sum([X(:,1)-Center0(1), X(:,2)-Center0(2), X(:,3)-Center0(3)].^2,2) - Radius0^2;
%% STEP 2
% Fit a full non linear least square sphere
[Center, Radius, ErrorDist, ~, conv, ~, ~, ~, ~, ~] = lssphere(X, Center0, Radius0, 1e-5, 1e-5)
%% If STEP 2 does not converge use STEP 1 results
if conv == 0 :
warning('*** Gauss-Newton algorithm has not converged for non linear least square sphere fit***');
warning('*** Using linear least square method for sphere fit instead***');
Center = Center0;
Raduis = Radius0;
ErrorDist = ErrorDist0;
@renaultJB good idea! I did notice that LSGE included a sphere fit but did not think about merging the two. I will give it a try.
@renaultJB still have to do this - I will deal with it this week