Robust Estimation of Sphere

Hi,

I have a point cloud file which represents a sphere. I want to find out the center and the radius of that sphere using some kind of robust estimation (M-estimation, for example). Can someone please help me regarding this?

Thanks in advance.

Zereen

[Walter Roberson]

When you say "estimation", do you mean you want to do it without examining all of the data points?

[Zereen]

Sorry, I am not sure whether I got your question or not. Still let me try. What I want to do is - using some kind of weight, I want to get rid of those points that are not part of sphere (outliers). As ordinary least squares may act severely in the presence of noise, I want to use robust estimation. Does it make sense? Please feel free to ask if I haven't answered your question properly.

Answer 1

[the cyclist]

Zereen, I am not aware of any way to do this "out of the box", and I did not find anything in the File Exchange. Do you have the Statistics Toolbox? If so, I think you could probably do something close to what you want by using the nlinfit() function, along with the 'Robust' option. Read "doc nlinfit" for details.

[Zereen]

Thanks. I am going to give it a try and will let you know if I have any further question.

[Zereen]

I was trying "nlinfit". As I was trying sphere fitting, my model function looks like:

y = sqrt((xi-x0)^2 + (yi-y0) + (zi - z0))-r

where,

xi, yi, zi (points on a sphere)= predictor variables

x0, y0, z0, r = variables to be determined

y = response variable

But when I tried to design the model function, I got an error - "??? Error using ==> mpower

Inputs must be a scalar and a square matrix."

Can you please help me to find out the problem?

[Sean de Wolski]

you probably have something like

(magic(5))^rand(5);

and actually want

magic(5).^rand(5)

note the .^ means element by element power as opposed to m(atrix)power

[Zereen]

Thanks Sean. But I solved the problem. Yes, you were right. I missed the . operator in one place. Thanks anyway.

:)

[Zereen]

I am trying to solve the problem using "nlinfit". But I am facing some problems. One of them its always produces the same values of beta as was set initially. Need suggestions. I am attaching my code here.

close all

clear all

clc;

[X, Y, Z] = File_Read ('20100204-000000-000-Ball-6.stl.txt');

X = X'; Y = Y'; Z = Z';

x = [X, Y, Z];

xo = 114; yo = 348; zo = 535; r = 17;

beta0 = [xo, yo, zo, r];

y = hModfun(beta0, x);

options = statset('Display', 'iter', 'MaxIter', 50, 'Robust', 'on', 'WgtFun', 'huber', 'TolX', 1e-20);

[par_out, residual, J, sigma] = nlinfit (x, y, @hModfun, beta0, options);

function y = hModfun (par, x)

y = sqrt(((x(:,1)-par(1)).^2)+((x(:,2)-par(2)).^2)+((x(:,3)-par(3)).^2))-par(4);

format long

disp (par)

end

[the cyclist]

Are you able to upload the txt file somewhere?

[Zereen]

Here is the link:

http://www.esnips.com/web/rushdecoder-Matlab

Initial values of beta was set using ordinary least-squares. The exact values were:

center = (114.4409432169675, 350.3414752720500, 535.3976447143736)

radius = 19.401690103987931

[the cyclist]

If you breakpoint immediately after your call to nlinfit(), you'll notice that the residual is all zeros at the first iteration, so it has found a solution. The reason it has found a solution is that you have defined the second and third input arguments (y and @hModfun) to be identical. In other words, your fitting function and your empirical data are the same, so the fit is perfect. You have misspecified the inputs.

The basic problem is that your definition of "y" should be purely from your data, but you have included the fitting parameter in there as well. I haven't thought carefully enough about this yet to know how to fix it.

[Zereen]

I think I got your point. But I am kind of confused what to do to solve it.

[the cyclist]

I think your problem boils down to this. This is going to be a bit pedantic, but that is mostly so I get the logic straight in my own head.

You have a series of (x,y,z), which lie on the surface of a sphere, except there is noise so they are not perfect. Given that series, you want to estimate the center and radius of the sphere. Theoretically (i.e. the fit), the radius is a fixed value, r_fit = r0. Empirically, the radius is r_empirical = sqrt((x-x0).^2+(y-y0).^2+(z-z0).^2).

So, in MATLAB you would LIKE to be able to write that as

>> nlinfit (x, sqrt((x-x0).^2+(y-y0).^2+(z-z0).^2), r0, beta0, options); % Won't work.

But you can't (I think), because that mixes the fitting parameters into both the empirical and the fitted part of nlinfit(). However, I believe that it is OK to instead do that fit by subtracting the "r_empirical" term from both pieces, getting

>> nlinfit (x, 0,r0-sqrt((x-x0).^2+(y-y0).^2+(z-z0).^2), beta0, options); % Almost there.

Getting the syntax right, and using your function handle, I think the following does that:

>> [par_out, residual, J, sigma] = nlinfit (x, zeros(size(X)),@hModfun, beta0, options)

I am not 100% sure that that is theoretically correct, but it does give a solution that is very close to your least-squares answer, so I have some confidence in its correctness.

[Zereen]

Thanks a lot. Let me give it a try. I will let you know the result.

[the cyclist]

One other comment. One thing that made this tricky is that your (x,y,z) data do not contain clear explanatory and response variables, which would normally be the case when fitting data. It really is more of an optimization problem, but we mashed it into a "fitting framework" to make use of the ability to use robust fitting. (Maybe there was a way to do that via something in the Optimization Toolbox, but I don't know.)

An alternative way to have done the fit would be to say that you have a z_fit and a z_empirical, with (x,y) being the explanatory. However, that puts z on a different footing than (x,y), which I didn't think made sense. But it is also a little odd that "r" is the response variable, given (x,y,z). I am not sure how theoretically strong my suggestions were. Let the buyer beware.