Fit gaussian surface on 3D data

Using the following script and attached image tif file, I am trying to fit a gaussian surface on the image data. Is there something that I skip? The result soes not seem to be gaussian. I also tried using the following but couldnt obtain a meaningful result.

https://www.mathworks.com/matlabcentral/fileexchange/41938-fit-2d-gaussian-with-optimization-toolbox

I appreciate any help/comment in advence.

pixelsize = 0.65512E-6; %pixel conversion [m]

folder = pwd;
baseFileName = 'imageprocessing.tif';
fullFileName = fullfile(folder, baseFileName); 
% info = imfinfo(baseFileName);
info = imfinfo(fullFileName);

grayImage = imread(fullFileName);
Z = grayImage./max(grayImage);
xx =  [1:size(grayImage,1)]*pixelsize;
yy =  [1:size(grayImage,2)]*pixelsize;
[X,Y]=meshgrid(xx,yy);

surf(X,Y,Z');shading interp
xlabel('x [\mu m]')
ylabel('y [\mu m]')
zlabel('Signal')
hold on 
xpdf = X(:); ypdf = Y(:); zpdf = Z(:);

gaussianfit = fit([xpdf, ypdf], zpdf, 'a*exp(-( ((x-b)/c)^2 + ((y-d)/e)^2 ) ) + f', 'StartPoint', [400, 3, 50, 0, 40, 450])
plot(gaussianfit, [xpdf, ypdf], zpdf)

[Matt J]

The image itself doesn't look very Gaussian. It seems to consist of a bright and a dark lobe. Which lobe are you trying to fit?

load grayImage

imshow(grayImage,[])

[Ali Aykut]

I am trying to fit the dark lob, The white region is an impurity. However, the fit seem to be something more plane than a gaussian surface. The data forms a surface attached below.

[Matt J]

a*exp(-( ((x-b)/c)^2 + ((y-d)/e)^2 ) )

It seems inadvisable to model the Gaussian term as separable in x and y, with only 5 parameters, instead of the more general 6-parameter model which has cross-terms,

Even if the object truly follows a separable Gaussian, a slight rotation of the imaging device would render the separable model invalid.

Answer 1

[Matt J]

Using gaussfitn from,

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

load Image %cropped to include only the dark lobe

Z=double(grayImage);
Z=Z(1:2:end,1:2:end);

[X,Y]=ndgrid(1:size(Z,1),1:size(Z,2));


b=diag([inf,inf]);
p=gaussfitn([X(:),Y(:)],Z(:),[],...
            {[],[],[],-b}, {[],[],[],+b});  %force covariance matrix to be diagonal

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

[D,A,mu,sig]=deal(p{:});


delta=([X(:),Y(:)]-mu')';

Zfit= D + A*exp( -0.5 * sum( (sig\delta).*delta,1) ); %gaussian fit
 Zfit=reshape(Zfit,size(X));


 %Display
 f=@(q) reshape(q(1:10:end,1:10:end),[],1);
 scatter3(f(X),f(Y),f(Z)); hold on
 surf(X,Y,Zfit,'FaceColor','r','FaceAlpha',0.5,'EdgeColor','none'); hold off
 xlabel X, ylabel Y; view(45,35)
 legend Data Fit

[Ali Aykut]

I appreciate your help Matt. Do you have any idea on how to find w (the waist of the light source at 1/e^2(13.5%) of the peak intensity of the Gaussian) using the fit function you wrote? Can it be found by using one of the parameters of the function output?

[Matt J]

Your diagram in IlluminationWaist.jpg only makes sense in the context of a 1D Gaussian. What does the "waist" mean in 2D?

[Ali Aykut]

It is a similar concept in both 1D and 2D. It is the half width at the approximately 14% of bleaching depth from the top.

[Matt J]

The isocontours of a 2D Gaussian surface are ellipses. What do you consider to be the "width" of an ellipse? Is it the semi-major axis? The semi-minor axis? Something in between?

Answer 2

[Image Analyst]

That's not 3-D data. That is 2-D data. How about you find the centroid with code I gave you before, then get the average radial profile, then fit a Guassian to that?

[Ali Aykut]

You are right, I meant a surface in 3D as attached(datasurface.png).

In fact, I found this dark regions radius change with time thanks to your help. It has a profile as in the attached file(radiuswtime.png). I found it by finding the area of the dark region with your script. But it would be more appropriate to locate the center and use accumarray as you suggest. This process might evolve into the way you expressed.

However, if I do so I am afraid I would loss some of the data or the quality of it. Wouldn't it less prone to error to directly fit a gaussian surface to the profile provided in 'datasurface.png'?

I appreciate your effort and help.

[Image Analyst]

I don't think I recommended accumarray. What is the profile for any one point in time? Does it look like a Gaussian? If the profile is gaussian, then I think the histogram would also be Gaussian and you could also fit the histogram to a Gaussian. The surface shown in 'datasurface.png' looks slightly tilted so you might want to flatten/level it before, though I think my demo had a tilt allowed. If you had two different Gaussians, one for x and one for y, then you'd have to account for that with additional terms in the model.

I don't know what this dark thing represents but I think youi may be overthinking it by trying for more accuracy that the situation really needs. Let's say you were somehow able to get the two Gaussian parameter sets, one along each principal direction, then what would you do differently than if you just had the standard deviation of the thing assuming a single radially symmetric Gaussian like I suggested? Probably nothing different, so why make it more complicated than it needs to be?