Clear Filters
Clear Filters

How to find the point that resides on a contour that is closest to a known point?

7 views (last 30 days)
I got this from another post of mine.
Z = peaks(50)/10;
level = 0.04;
surf(Z)
hold on
% extract all isocline for a given level
[C,h] = contour(Z,level*[1 1]);
[m,n] = size(C);
ind = find(C(1,:)==level); % index of beginning of each isocline data in C
ind = [ind n+1]; % add end (+1)
for k = 1:numel(ind)-1
xc = C(1,ind(k)+1:ind(k+1)-1);
yc = C(2,ind(k)+1:ind(k+1)-1);
zc = level*ones(size(xc));
plot3(xc,yc,zc,'linewidth',5);
end
hold off
Now, suppose I know a point (16,35,0.05). I want to find the closest point on a contour whose height is zc = 0.056.
How do I do that? What sort of minimization problem needes to be solved? And how do I replicate the usual minimization that we use in differential calculus like dy/dx =0?

Accepted Answer

John D'Errico
John D'Errico on 28 Apr 2024
Edited: John D'Errico on 28 Apr 2024
You can't. That is, you cannot use calculus to do this, at least not as it is, not directly. That contour is a set of piecewise linear segments. Could you solve the problem in some way with calc? Well, yes. You could do it in the case of some simple functions.
You CAN use my distance2curve utility, found on the File Exchange for free download. (In fact, if you allow it to do so, it does use calculus to solve the problem, in a sense. So if the curve is a spline, then it needs to formulate and solve a problem, where it implicitly minimizes the distance to a smooth curve.)
It seems like you want to learn how to formulate the problem in an analytical form, so I'll do it that way. If all you want to do is solve the problem for a specific contour, then download and use my function. I'll give an example for a simple function though, so you can learn how to solve the problem.
syms x y
z(x,y) = x^2 + x*y + y^2;
And that would be just a simple ellipse, in terms of its contours. We can see one such contour, for example at z(x,y) == 1
fimplicit(z - 1)
axis equal
xlabel x
ylabel y
Now can we formulate the problem in terms of calculus, as you seem to want? The curve is z(x,y)-1==0. Ugh. I can think of at least 4 different ways to solve even this simple problem.
The most direct seems to be to use Lagrange multipliers.
syms lambda x0 y0
We want to mimimize the expression
L = (x - x0)^2 + (y - y0)^2 + lambda*(z(x,y)-1)
L = 
We can do that by differentiating, setting the derivatives to zero. The result will be a nonlinear system of three equations, in three unknowns, here x,y, and lambda.
I'll pick some specific point, say (x0,y0) = (2,3).
xylamb = solve(gradient(subs(L,[x0,y0],[2,3]),[x,y,lambda]) == 0,[x,y,lambda],returnconditions = true)
xylamb = struct with fields:
x: [4x1 sym] y: [4x1 sym] lambda: [4x1 sym] parameters: [1x0 sym] conditions: [4x1 sym]
vpa([xylamb.x,xylamb.y])
ans = 
In there, we can ignore the complex solutions. The one that actually mimizes the distance happens to be the 4th one. (The first is probably the point of maximum distance.)
fimplicit(z - 1)
hold on
plot(xylamb.x(1),xylamb.y(1),'rx')
plot([2,xylamb.x(4)],[3,xylamb.y(4)],'-go')
axis equal
hold off
Ok, can you solve this in a simpler way, that does not force you to use Lagrange multipliers? We can certainly use fmincon.
zlevel = 1;
f = @(xy) xy(1)^2 + xy(1)*xy(2) + xy(2)^2 - zlevel;
nonlcon = @(xy) deal([],f(xy));
xy0 = [2 3];
obj = @(xy) sum((xy - xy0).^2);
xystart = [0,0];
[xysol,fval,exitflag] = fmincon(obj,xystart,[],[],[],[],[],[],nonlcon)
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
xysol = 1x2
0.3290 0.7940
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
fval = 7.6585
exitflag = 1
And that is the same as we got from using Lagrange multipliers. But you might argue that it implicitly used Lagrange multipliers anyway.
So, can we do it in a simpler way yet? Probably, since the equation of that ellipse simplifies greatly in a parametric form. Or, you could use my distance2curve utility.
However, your question was on a more general function, such as peaks, where the contour would be less easy to define. For a completely general surface, we might start with contour or contourc.
[X,Y] = meshgrid(linspace(-5,5,101));
Z = peaks(X,Y);
C = contour(X,Y,Z,[1,1])
C = 2x206
1.0000 -1.1984 -1.1000 -1.0000 -0.9000 -0.8000 -0.7000 -0.6605 -0.6000 -0.5000 -0.4328 -0.4000 -0.3000 -0.2522 -0.2000 -0.1037 -0.1000 0 0.0140 0.1000 0.1006 0.1588 0.1806 0.1733 0.1436 0.1000 0.0977 0.0498 0 -0.0050 59.0000 -1.2000 -1.2607 -1.2766 -1.2691 -1.2476 -1.2165 -1.2000 -1.1778 -1.1337 -1.1000 -1.0840 -1.0285 -1.0000 -0.9658 -0.9000 -0.8968 -0.8112 -0.8000 -0.7006 -0.7000 -0.6000 -0.5000 -0.4000 -0.3000 -0.2054 -0.2000 -0.1000 -0.0082 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The result here are two distinct contours, both of which have z(x,y)=1. So we would then find the nearest point on each of those distinct curves, then choose the solution with the smaller distance. As an example, I'll pick the point (x,y)=(-2,1).
xy0 = [-2,1];
[xymin1,d1] = distance2curve(C(:,2:60)',xy0)
xymin1 = 1x2
-0.6000 0.0048
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
d1 = 1.7177
[xymin2,d2] = distance2curve(C(:,62:end)',xy0)
xymin2 = 1x2
-1.2655 1.4000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
d2 = 0.8364
So this point fell at a distance of 1.7177 units from the lower curve, but only 0.8364 units from the upper curve.
hold on
plot([xy0(1),xymin1(1)],[xy0(2),xymin1(2)],'r-x',[xy0(1),xymin2(1)],[xy0(2),xymin2(2)],'g-s')
axis equal

More Answers (1)

Torsten
Torsten on 28 Apr 2024
Moved: Torsten on 28 Apr 2024
format long
Z = peaks(50)/10;
surf(Z)
obj = @(x,y)(x-16.35)^2+(y-0.05)^2;
nonlcon = @(z)deal([],interp2(0:49,0:49,Z,z(1),z(2))-0.056);
lb = [0 0];
ub = [49 49];
xy0 = [16 0];
sol = fmincon(@(z)obj(z(1),z(2)),xy0,[],[],[],[],lb,ub,nonlcon)
Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.
sol = 1x2
36.000000001470539 15.766401131249982
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
obj(sol(1),sol(2))
ans =
6.331277645761478e+02
interp2(0:49,0:49,Z,sol(1),sol(2))-0.056
ans =
0

Tags

Community Treasure Hunt

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

Start Hunting!