Surface data cloud fitting to even asphere model

I've got a set of (r,z) data which represent a surface section. I would like to fit such data to the Even asphere expression:

Where is the surface sagita, r is the radial coordinate, R is the radius of curvatre at the vertex, κ is the conic constant and are the coefficients describing the deviation of the surface from a pure conic section. Can anyone help to address the problem?

I have calculated the value of R by first fitting the data to a sphere with the code attached below, but I can't find a way to fit the data to the term once R has been calculated. Thank you, your help is very much appreciated!

-----------------------------------------------------------------

Function to fit circle to data set:

function Par = CircleFitting(XY)

centroid = mean(XY); % the centroid of the data set

X = XY(:,1) - centroid(1); % centering data

Y = XY(:,2) - centroid(2); % centering data

Z = X.*X + Y.*Y;

Zmean = mean(Z);

Z0 = (Z-Zmean)/(2*sqrt(Zmean));

ZXY = [Z0 X Y];

[U,S,V]=svd(ZXY,0);

A = V(:,3);

A(1) = A(1)/(2*sqrt(Zmean));

A = [A ; -Zmean*A(1)];

Par = [-(A(2:3))'/A(1)/2+centroid , sqrt(A(2)*A(2)+A(3)*A(3)-4*A(1)*A(4))/abs(A(1))/2];

end

Answer 1

[Torsten]

% your data
x = ...;
y = ...;
z = ...;
p0 = [...];  % initial guess for R and kappa
p = lsqnonlin(@(p)fun(p,x,y,z),p0)
    
function res = fun(p,x,y,z)
  rho = 1/p(1);
  kappa = p(2);
  r = sqrt(x.^2+y.^2);
  res = z - r.^2*rho./(1+sqrt(1-(1+kappa)*(r*rho).^2));
end

[Cl ES]

Thank you Torsten for your quick response.

I have tried the method you proposed both for the pure conic and a fourth order asphere, but yet I cannot reproduce the data. I attatch my data and paste here the script, in case you detect any mistake... thanks!

CurrentFolder = pwd;

FittingDir = '\FittingFiles';

myData = 'data.txt'

filename = fullfile(CurrentFolder,FittingDir,myData);

s = load(filename); % First column is radial coordinate, 2nd column is sagita value

T=TaubinSVD(s);

BaseR = T(1,3);

SphericalSag = BaseR - sqrt(BaseR.^2 - s(:,1).^2);

figure;hold on; grid on;

plot(s(:,2),s(:,1),'b-');hold on;

plot(SphericalSag,s(:,1),'r.');

xlabel('sag (mm)');ylabel('radial coordinate (mm)');

legend('Original data', 'FittedSphere')

% Fit to conic section

fun = @(p)rdata.^2./(BaseR*(1+sqrt(1-(1+p)*rdata.^2/BaseR^2)))-zdata;

x0 = 0;

constantk = lsqnonlin(fun,x0);

ConicSag = rdata.^2./(BaseR*(1+sqrt(1-(1+constantk)*rdata.^2/BaseR^2)));

figure;hold on; grid on

plot(s(:,2),s(:,1),'b-');

plot(SphericalSag,s(:,1),'r.');

plot(ConicSag,s(:,1),'m.');

xlabel('sag (mm)');ylabel('radial coordinate (mm)');

legend('Original data', 'Fitted sphere', 'Fitted conic');

%Fit to 4th order even asphere (no conic)

fun = @(q)rdata.^2./(BaseR*(1+sqrt(1-rdata.^2/BaseR^2))) + q*rdata.^2 - zdata;

x0 = 0;

a4 = lsqnonlin(fun,x0);

[Cl ES]

Sorry in the last part I made a typing mistake, the correct expression for fun would be:

fun = @(q)(...) + q*rdata.^4

[Torsten]

Sorry, I cannot recognize "my" code from what you programmed.

I tried to fit both R and kappa simultaneously.

So in order to test my approach, I need to know what the 2 columns of your data file contain. Is it correct to assume that the first column is r and the second column is z ?

[Cl ES]

Hi Trosten,

Thank you for your help. You're right, the first column ir r and the second is z.

[Torsten]

And r can be negative ? I thought it is sqrt(x^2+y^2) ?

[Bruno Luong]

it doesn't matter aspheric equation is even and depends on r^2 only

[Torsten]

data = readmatrix('https://de.mathworks.com/matlabcentral/answers/uploaded_files/1119340/data.txt');
r = data(:,1);
z = data(:,2);
p0 = [50 0.01];  % initial guess for R and kappa
lb = [0 0];
ub = [Inf Inf];
format long
%p = lsqnonlin(@(p)fun(p,r,z),p0)
p = fmincon(@(p)fun(p,r,z),p0,[],[],[],[],lb,ub,@(p)nonlcon(p,r,z))

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

p = 1×2

1.243395899999452 0.005562271920655

fun(p,r,z)

ans = 
     3.551188401724851e-05

    
function res = fun(p,r,z)
  rho = 1/p(1);
  kappa = p(2);
  %res = z - r.^2*rho./(1+sqrt(1-(1+kappa)*(r*rho).^2));
  res = sum((z - r.^2*rho./(1+sqrt(1-(1+kappa)*(r*rho).^2))).^2);
end

function [c,ceq] = nonlcon(p,r,z)
  rho = 1/p(1);
  kappa = p(2);
  c = -(1-(1+kappa)*(r*rho).^2);
  ceq = [];
end

[Cl ES]

Hi Trosten,

thank you for your help. I've got a question, though. Why are you squaring the denominator ?

[Torsten]

So many brackets ... I corrected the code.