Finding points of zero slope from 3D function

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?

[Tim Berk]

• 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.

Answer 1

[Nicolas Schmit]

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');

[James Akula]

This is excellent. Thank you.

[Nicolas Schmit]

You're welcome.

[Karan Gill]

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