CODEWITHSUNDEEP
matlabimage-processing2dleast-squares
pse18.10
pse18.10

Reputation: 15

Matlab: Least square fit to 2d data set

I have the following issue:

Model: centrally symmetric circle with a profile which is a combination of gaussian and lorentzian distribution. To get the plot of the model just insert the following code to Matlab:

N = 501;                            %Matrix size
R = zeros(N,N);                     %Initializing matrix R
x0 = ceil(N/2); y0 = x0;            %Barycenter coordinates

for i=1:N                           %Calculation of matrix R
    for j=1:N
        R(i,j) = sqrt((x0-j)^2 + (y0-i)^2);          
    end       
end

%Model z1 and the parameters:
peak = 0.275*N;                       %Peak location
m = 0.3;
sigma = 0.1*N;
gamma = 15;
A1 = 1000;
A2 = 50;
z1 = (1-m)*A1/(sigma*sqrt(pi))*exp(-(abs(R - peak)).^2/sigma^2) + m*A2/pi * (gamma./((abs(R - peak1)).^2 + gamma^2));

figure('name','Show Model')
surf(z1,'EdgeColor','none','LineStyle','none','FaceLighting','phong');

So this is the "idealistic" model. To simulate real data, I will add random noise to z1:

z2 = z1 + random('Normal',0,1,N,N);
figure('name','Show random noise data')
surf(z2,'EdgeColor','none','LineStyle','none','FaceLighting','phong');

Finally a plot of the intersecting plane through the barycenter:

figure('name','Show intersecting plane with model and random noise data')
xaxis = -floor(N/2):1:floor(N/2);
intersectionline1 = z1(ceil(N/2),:);
intersectionline2 = z2(ceil(N/2),:);
plot(xaxis,intersectionline1,xaxis,intersectionline2,'.r');

Z2 could be for example a real dataset of my measurements. The circle is anywhere in my image. I am actually able to find the barycenter and get a squared sector. Besides that I get a good approximation of the radius r (variable "peak").

My question now: Is it possible to make a least square fit to my dataset to get the parameter m, sigma, gamma, peak, A1 and A2 from my model??? I don´t know if Matlab is able to do so...

Thanks in advance!

Upvotes: 0

Views: 942

Answers (1)

user2271770
user2271770

Reputation:

You may use fminsearch to find the set of coefficients that give the minimal least squares:

% peak  = K(1)^2
% m     = K(2)^2
% sigma = K(3)^2
% gamma = K(4)^2
% A1    = K(5)^2
% A2    = K(6)^2

model = @(z0, K) (1-K(2)^2)*K(5)^2/(K(3)^2*sqrt(pi))*exp(-(abs(z0 - K(1)^2)).^2/K(3)^4) + K(2)^2*K(6)^2/pi * (K(4)^2./((abs(z0 - K(1)^2)).^2 + K(4)^4));
errsq = @(K) sum(sum((z2-model(R,K)).^2));
K0 = ones(6,1);
K  = fminsearch(errsq,K0);
figure('name','Reconstructed');
surf(model(R,K),'EdgeColor','none','LineStyle','none','FaceLighting','phong');

To make sure that the coefficients stay positive, the model uses squared parametrization. Of course, the initial guess K0 is very important for the final result.

Upvotes: 2

Related Questions