I want to plane fit my sperical data points.

I want to correct this data with the radius, so that it creates a plane. The radius is known due measurements.

Thank you in advance.

Answer 1

[Bora Eryilmaz]

Assuming a spherical surface with unknown origin (and perhaps radius), you can run an optimization algorithm to estimate the model parameters (origin, etc.) and subtract the location of the points on the surface of the sphere from your data. That would give you, roughly, a plane-like view of your data at a distance equal to the radius of the sphere.

A better planar view might be to actually project the data into an r-theta-z plane, but this would require a more complicated algorithm.

% "Unknown" model parameters.
x0 = 45;
y0 = 30;
z0 = -198;
r = 200.0;

% Construct the data surface
[X,Y] = meshgrid(0:1:100, 0:1:60);
Z = sqrt((r+randn(size(X))).^2 - (X-x0).^2 - (Y-y0).^2) + z0;
surf(X,Y,Z)

% Estimate model parameters x0, y0, z0, and r.
% You can remove "r" from the estimation if it is already known.
P0 = [30, 30, -30, 10];
options = optimset('MaxFunEvals', 2000);
P = fminsearch(@(p) fcn(p,X,Y,Z), P0, options)

P = 1×4

44.1685 29.7778 -200.0555 201.9682

% Project into the plane at distance r.
x0 = P(1);
y0 = P(2);
z0 = P(3);
r = P(4);
Zplane = Z - sqrt(r^2 - (X-x0).^2 - (Y-y0).^2) - z0;
surf(X,Y,Zplane)

function cost = fcn(p, X, Y, Z)
x0 = p(1);
y0 = p(2);
z0 = p(3);
r = p(4);
z = sqrt(r^2 - (X-x0).^2 - (Y-y0).^2) + z0;
cost = norm(z-Z, 2); % Least squares cost.
end

[Chiel van Wanrooij]

Thank you for your answer.

In the code below you can see how I implemented it. My xco is a 768x1 double, my yco is a 576x1 double and my xco is a 768x576 double. I use the inpaint_nans to interpolate for the NaNs in my measured data. My X0, Y0 and Z0 are the coordinates of my origin. The area of the rectangle is a 86x64, so I chose the middle of this plane. the radius is 1500. Than I use the formula for Zplane. However, the arrays have incompatible sizes for this operation.

Is it possible to solve this problem?

% Raw data

N_rows = fread(fid,1,'int');  % Number of rows
disp(N_rows')
M_rows = fread(fid,1,'int');    % number of columns
disp(M_rows')
Data = fread(fid,N_rows*M_rows,'float');   % Data array
%disp(Data')
Data_max = fread(fid,1,'float');    % Maximum of data
disp(Data_max')
Data_min = fread(fid,1,'float');    % minimum of data
disp(Data_min')
fclose(fid);

%%% Make raw data file with Sensofar settings and measured data

wrong_dat=find(Data==1000001);
Data(wrong_dat)=NaN;

xco=[0 :mppx:(M_rows-1)*mppx]'; % x-coordinates with x=0 is the first x coordinate
yco=[y0 :-mppy:(N_rows-1)*(-mppy)+y0]'; % y-coordinates with y=0 is the middle
zco=[Data(1:M_rows)]';  %   measured z-coordinates in a matrix
for i=2:N_rows
    z=[Data(((i-1)*M_rows)+1:(i)*M_rows)]';
    zco(i,:)=z;
end

X0 = 43
Y0 = 32
Z0 = -1500
r = 1500
fzco = inpaint_nans(zco,4);

Zplane = fzco - sqrt(r^2 - (xco-X0).^2 - (yco-Y0).^2) - Z0
mesh(xco,yco,Zplane)


% figure(1)
% A = mesh(xco,yco,zco) 
% figure(2)
% B = mesh(xco,yco,fzco)
% C = surfarea(B)

[Bora Eryilmaz]

Yes, xco and yco need to be the same size for this to work. You may need to use ngrid() or meshgrid() commands to convert your data into the meshgrid format: https://www.mathworks.com/help/matlab/ref/meshgrid.html#mw_6ae7effe-9402-4974-b1a6-0391eba290d8.

[Chiel van Wanrooij]

Thank you for your help. I appreciate it.

Answer 2

[Matt J]

You can use sphericalFit() from this FEX download

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

to fit a sphere to your points (non-iteratively). This does not currently allow you to constrain the radius to a known value, however, if that is essential, it would at least give you a pretty good initial guess of (x0,y0,z0) for the iterative optimization proposed by @Bora Eryilmaz.