Finding points of zero slope from 3D function

3 次查看（过去 30 天）
James Akula
James Akula 2017-9-19
编辑： Karan Gill 2017-9-21

I will be fitting 3D differences of Gaussians to data, and from these fits I want to be able to identify troughs and peaks (points of zero slope) in my resultant functions. I can compute XYZ data points and look for local maxima, but the results are not satisfactory, at least the way I am doing it.

Below is a sample script showing my problem:

syms a1 a2 b1 b2 b3 b4 c x x0 y y0

DoG3D(a1, a2, b1, b2, b3, b4, c, x, x0, y, y0) = ...
    (a1 * exp(-1/2 * (((x - x0)/b1)^2 + ((y - y0)/b2)^2))) - ...   % Gaus 1
    (a2 * exp(-1/2 * (((x - x0)/b3)^2 + ((y - y0)/b4)^2))) + ...   % Gaus 2
    c;

% The constants below will differ, depending on a fit, but here are some
% convenient ones
%                      a1   a2   b1  b2   b3   b4   c   x  x0   y  y0
DoG3DFit(x, y) = DoG3D(-50, -25, 50, 150, 100, 300, 10, x, 220, y, 180);

pretty(DoG3D)

figure
subplot(2, 2, 1)
fsurf(DoG3DFit, [-1000, 1000])
daspect([50, 50, 1])

% Trying to solve it numerically is unsatisfactory (at least this way)
[X, Y] = meshgrid(-1000:1000, -1000:1000);
matDoG3DFit = matlabFunction(DoG3DFit);
Z = matDoG3DFit(X, Y);
Zpoints = or(imregionalmax(Z), imregionalmax(-Z));
subplot(2, 2, 2)
imagesc(Zpoints)
axis image xy

% In 2D, this is easy
DoG(a1, a2, b1, b2, c, x, x0) = ...
    a1 * exp(-1/2 * (((x - x0)/b1)^2)) - ...                       % Gaus 1
    a2 * exp(-1/2 * (((x - x0)/b2)^2)) + ...                       % Gaus 2
    c;
pretty(DoG)

zeroPoints = solve(diff(DoG(-50, -25, 50, 150, 10, x, 220), x) == 0, x);
subplot(2, 2, 3)
fplot(DoG(-50, -25, 50, 150, 10, x, 220), [-1000, 1000])
hold on
scatter(zeroPoints, DoG(-50, -25, 50, 150, 10, zeroPoints, 220), 'r')

% How do I do the same in 3D?

I am hoping there is a way, using the Symbolic Math or other toolbox, to generate formulae that will describe the locations where the slope is zero. These should be, to my understanding, at point (x0, y0), and then at points around the rim of the "caldera." The location of those points will depend upon what constants I supply DoG3D, above.

If this were a 2D difference of Gaussians, I would simply take the derivative and solve for 0, but I do not know what the equivalent approach is for a 3D function.

Can someone please help me to achieve this? Or is there no such solution?

  1 个评论
Tim Berk
Tim Berk 2017-9-19
  • I don't think this is what you want, but if you just want to plot the maximum you can use contourslice, i.e.
...
zeroPoints =...
% start added code
max_val = max(max(eval(DoG3DFit(X(1:100:end,1:100:end),Y(1:100:end,1:100:end)))));
[xmesh,ymesh,zmesh] = meshgrid(-1000:1000, -1000:1000,[0 max_val]);
V = matDoG3DFit(xmesh, ymesh);
hold on
cs = contourslice(xmesh,ymesh,zmesh,V,[],[],max_val,[max_val max_val]);
cs(1).EdgeColor = 'red';
cs(1).LineWidth = 3;
% end added code
subplot(2,2,3) ...
  • But from your question it looks as if you want to be able to generate an analytical solution.Doesn't the original (2D) method still work? I think at the 'caldera' both dz/dx = 0 and dz/dy = 0.So if you find either of these derivatives, you can find the maximum/minimum.Offcourse these derivatives are a bit more complicated..
Hope this helps.

请先登录，再进行评论。

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

采纳的回答

Nicolas Schmit
Nicolas Schmit 2017-9-20
To find the solution, you search the x and y for which the partial derivatives are null. This is tricky for an automated solver because the solution is complicated. A single solution in the center of the caldera, and an infinite number of solutions along the rim.
Below is an example of how you can calculate the symbolic solutions (R2017a).
syms a1 a2 b1 b2 b3 b4 c x x0 y y0
DoG3D(a1, a2, b1, b2, b3, b4, c, x, x0, y, y0) = ...
(a1 * exp(-1/2 * (((x - x0)/b1)^2 + ((y - y0)/b2)^2))) - ... % Gaus 1
(a2 * exp(-1/2 * (((x - x0)/b3)^2 + ((y - y0)/b4)^2))) + ... % Gaus 2
c;
% The constants below will differ, depending on a fit, but here are some
% convenient ones
% a1 a2 b1 b2 b3 b4 c x x0 y y0
DoG3DFit(x, y) = DoG3D(-50, -25, 50, 150, 100, 300, 10, x, 220, y, 180);
assume(x, 'real')
assume(y, 'real')
% Calculate the partial derivatives
dfdx = simplify(diff(DoG3DFit(x, y), x));
dfdy = simplify(diff(DoG3DFit(x, y), y));
% Solve with respect to y (3 solutions)
soly = simplify(solve(dfdy == 0, y));
% Substitute
dfdy = simplify(subs(dfdy, y, soly));
dfdx = simplify(subs(dfdx, y, soly));
% Only solution 1 has dependency in x
solx = simplify(solve(dfdx(1) == 0, x));
% Point solutions
pointSolutions = [solx repmat(soly(1), numel(solx), 1)]
% Solutions 2 and 3 are functions of x
functionSolutions = soly(2:3)
%%Plot the solutions
f2 = matlabFunction(soly(2), 'Vars', x);
f3 = matlabFunction(soly(3), 'Vars', x);
X = linspace(-500, 500, 3000);
Y2 = f2(X);
i2 = arrayfun(@(x) imag(x) == 0, Y2);
Y3 = f3(X);
i3 = arrayfun(@(x) imag(x) == 0, Y3);
figure;
fsurf(DoG3DFit, [-1000, 1000])
hold on
for k=1:size(pointSolutions, 1)
plot3(pointSolutions(k, 1), pointSolutions(k, 2), DoG3DFit(pointSolutions(k, 1), pointSolutions(k, 2)), 'or', 'MarkerSize', 10, 'LineWidth', 5);
end
plot3(X(i2), Y2(i2), DoG3DFit(X(i2), Y2(i2)), 'm', 'LineWidth', 5);
plot3(X(i3), Y3(i3), DoG3DFit(X(i3), Y3(i3)), 'm', 'LineWidth', 5);
hold off
xlabel('x');
ylabel('y')
zlabel('z');
  3 个评论
显示 1更早的评论
Karan Gill
Karan Gill 2017-9-21
编辑：Karan Gill 2017-9-21
If you just want the single minima in the caldera centre, you can solve the equations simultaneously with vpasolve with initial guesses. Providing the initial guesses is necessary here. But as Nicholas says, for the curve along the rim, this won't work.
clear all, reset(symengine)
syms a1 a2 b1 b2 b3 b4 c x x0 y y0
assume([x y],'real')
DoG3D(a1, a2, b1, b2, b3, b4, c, x, x0, y, y0) = ...
(a1 * exp(-1/2 * (((x - x0)/b1)^2 + ((y - y0)/b2)^2))) - ... % Gaus 1
(a2 * exp(-1/2 * (((x - x0)/b3)^2 + ((y - y0)/b4)^2))) + ... % Gaus 2
c;
% The constants below will differ, depending on a fit, but here are some
% convenient ones
% a1 a2 b1 b2 b3 b4 c x x0 y y0
DoG3DFit(x, y) = DoG3D(-50, -25, 50, 150, 100, 300, 10, x, 220, y, 180);
f = DoG3DFit;
dfdx = diff(f,x);
dfdy = diff(f,y);
eqns = [dfdx==0, dfdy==0];
vars = [x y];
init_guesses = [250 150];
[xSol ySol] = vpasolve(eqns,vars,init_guesses)
xSol =
220.0
ySol =
180.0

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Probability Distributions 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by