Fitting a 3D gaussian to a 3D matrix

I have a 3D matrix that I need to fit with a 3D gaussian function:

I need to get A, and all three σ's as the output after fitting.

I have tried to do it using Least Square fitting as:

[xx,yy,zz]=meshgrid(x,y,z);
Mat(:,:,:,1)=xx;Mat(:,:,:,2)=yy;Mat(:,:,:,3)=zz;
F=@(param,Mat) param(1)*exp(-(Mat(:,:,:,1).^2/(param(2))^2+Mat(:,:,:,2).^2/(param(3))^2+Mat(:,:,:,3).^2/(param(4))^2));
param0=[80000,10,10,10]; % initial parameters
param=lsqcurvefit(F,param0,Mat,D); % D is the 3D matrix that needs to be fit.

But, I get various errors. I don't think my approach is correct. Is there any other way?

Answer 1

[Are Mjaavatten]

meshgrid arranges the xx and yy output matrices to fit the convention where the surface plots have x on the horizontal axis and y on the vertical, while you want to have x vary in the vertical direction in your matrices.

[yy,xx,zz]=meshgrid(y,x,z);

should fix the problem.

[Yatish Yatish]

Thanks Are! I guess that's not a problem. The code gives such message:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point 
is less than the default value of the optimality tolerance.

<stopping criteria details>

[Are Mjaavatten]

In order to say anything meaningful about this I need the data for your problem. Can you submit a mat file containing x, y, z, and D?

[Yatish Yatish]

Hi @Are MjaavattenI have attached the data file.

D=load('data_25x17x17.mat');
x=1:1:25;y=1:1:17;z=1:1:17;
% x=-12:1:12;y=-8:1:8;z=-8:1:8; % can be tried
[x1,y1,z1]=meshgrid(y,x,z);
Mat(:,:,:,1)=x1;Mat(:,:,:,2)=y1;Mat(:,:,:,3)=z1;
F=@(param,Mat) param(1)*exp(-(Mat(:,:,:,1).^2/(param(2))^2+Mat(:,:,:,2).^2/(param(3))^2+Mat(:,:,:,3).^2/(param(4))^2));
param0=[max(max(max(D))),8,8,8]; % at the centre of this volume, lies the maxima/amplitude of the fit fuction
param=lsqcurvefit(F,param0,Mat,D);

Thanks :-)

[Are Mjaavatten]

I made a few minor changes to your code, so now it runs and converges, However, it seems that your Gaussian model does not give a very good fit to the data. The plot shows the D surface, together with the model results, for z = 10. These are typical for other values of z as well. The model results are the relatively flat surface, while the wavy surface is your D data.

D=load('data_25x17x17.mat');
D = D.D;
x=1:1:25;y=1:1:17;z=1:1:17;
[y1,x1,z1]=meshgrid(y,x,z);
Mat(:,:,:,1)=x1;Mat(:,:,:,2)=y1;Mat(:,:,:,3)=z1;
F=@(param,Mat) param(1)*exp(-(Mat(:,:,:,1).^2/(param(2))^2+Mat(:,:,:,2).^2/(param(3))^2+Mat(:,:,:,3).^2/(param(4))^2));
param0=[3e6,8,8,8];
param=lsqcurvefit(F,param0,Mat,D);

DD = F(param,Mat);   % Values calculated from converged parameters

figure;
surf(y,x,D(:,:,10));hold on;surf(y,x,DD(:,:,10));xlabel y;ylabel x;title('z = 10 plane')

You may also try other slices, e.g.:

k=25;figure;
surf(y,z,permute(D(k,:,:),[2,3,1]));hold on;surf(y,z,permute(DD(k,:,:),[2,3,1]));
xlabel z;ylabel y;title('x = 25 plane') 

It is obvious that the Gaussian model is unable to follow the wavy nature of the data. Maybe you are able to find a better model from you knowledge of the process that produced D.

Good luck.

[Yatish Yatish]

The problem was with the matrix that I had attached. It works perfectly for the right 3D data. Thanks a lot for your help! :-)

Answer 2

[Matt J]

You can also use gaussfitn,

https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn

[xx,yy,zz]=meshgrid(x,y,z);

lbCell={0,[],[], zeros(3)};
ubCell={0,[],[], +diag(inf(1,3))}
param0={0,80000,[0,0,0], (1/10)*eye(3) };

param = gaussfitn([xx(:),yy(:),zz(:)], D(:),param0, lbCell, ubCell)