Fit gaussian surface on 3D data

14 次查看（过去 30 天）
Ali Aykut
Ali Aykut 2023-5-26
编辑： Matt J 2023-5-30
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)
  3 个评论
显示 1更早的评论
Ali Aykut
Ali Aykut 2023-5-27
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
Matt J 2023-5-27
编辑：Matt J 2023-5-27
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-5-27
编辑：Matt J 2023-5-27
Ran in:
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
  4 个评论
显示 2更早的评论
Ali Aykut
Ali Aykut 2023-5-30
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
Matt J 2023-5-30
编辑：Matt J 2023-5-30
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?

请先登录，再进行评论。

更多回答（1 个）

Image Analyst
Image Analyst 2023-5-27
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?
  2 个评论
Ali Aykut
Ali Aykut 2023-5-27
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
Image Analyst 2023-5-27
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by