Finding a point from a set of contours such that it is nearest to a given point

Question

I have several contour lines and one point. How can I find a point in one of those contour lines which is nearest to the given point?

(*Create the implicit curves*)
Data={{10,20,1},{10,40,2},{10,60,3},{10,80,4},{20,25,2},{20,45,3},{20,65,4},{30,30,3},{30,50,4},{40,35,4},{40,55,5},{50,20,4},{50,40,5},{60,25,5}};
U=NonlinearModelFit[Data,a x^b (y^(1-b))+c,{a,b,c},{x,y}];
L={ContourPlot[U[x,y]=={1},{x,0,100},{y,0,100},ContourStyle->Red],ContourPlot[U[x,y]=={2},{x,0,100},{y,0,100},ContourStyle->Magenta],ContourPlot[U[x,y]=={3},{x,0,100},{y,0,100},ContourStyle->Brown],ContourPlot[U[x,y]=={4},{x,0,100},{y,0,100},ContourStyle->Blue],ContourPlot[U[x,y]=={5},{x,0,100},{y,0,100},ContourStyle->Green]};
(*Point nearest to which we need to find the points on the curves*)
pt={30,50};
(*Graphic*)
Show[L,Graphics[{PointSize[Large],Blue,Point[pt]}],FrameLabel->{"X","Y"}]

enter image description here

Would it be possible to approach this by looking for intersections of a circle (centered at pt, with variable radius) with each of the curves? As radius decreases, the circle would intersect with fewer curves, eventually only intersecting one of them. As it shrinks further, the pair of intersections with that curve will become a single point, the closest point to the center of the circle. (Apologies if this is naive.) — DavidC, Oct 18 '13 at 13:37
@c_susetyo I deleted the end part of your question containing my old code as I just adapted my answer to your updated question. — PlatoManiac, Oct 18 '13 at 20:08

PlatoManiac · Answer 1 · 2013-10-18T20:06:27.693

For a general purpose & more accurate solution I can suggest the following.

Update(the missing data part)

(* Create the implicit curves analytically *)
conts = ((a x^b  (y^(1 - b)) + c /. U["BestFitParameters"]) == #) & /@Range[5];
(* Given x val function to find y=f(x) values per curve *)
fun[exp_] := {#,y /. Quiet@FindRoot[Evaluate@(exp /. x -> #), {y, 1}]} &;
(* Find discrete (x,f(x)) for 0<x<100 *)
data = (fun[#] /@ Range[.1, 101, .25]) & /@ conts;

Continuing the old code

(* Interpolate the data *)
    ints = Interpolation /@ data;
    (* Given a point or a list of points following function finds the minimizer
    (x,f(x)) on each curve *)
    xpts[pt_] := With[{int = #},
    NMinimize[{Evaluate@Total[Sqrt[(x - #1)^2 + (int[x] - #2)^2] &@@@ pt], .1 <x<100},
    {x}]] & /@ ints;
    {#2[[2]], #3[#2[[2]]]} & @@ Flatten@First@Sort@Transpose@{xpts[pt], ints}

{36.9385, 52.1471}

For fun lets allow some dynamics! Note that I did not update the animations here. The code is self-contained now so one can reproduce similar things easily.

enter image description here

Note:

The above function is written in such a way that one can have more than one point say $X=(x_1,...,x_n)$. Then it will find the unique minimizer $P=(a,b)$ on the given contours such that $dist=\sum_{i=1}^{n}d(P,x_i)$ is minimized where $d(*,*)$ is the Euclidean distance metric. Such a case follows in the figure. Such an example follows below. Here the number or points $\#(X)=n$ is increasing at every step of the animation. They are the black dots on the plot.

One should also note given a large number of points one can try FindMinimum in place of NMinimize but for up to 500 points NMinimize seems to do the job well enough.

enter image description here

Code:

xptsSum[pt_] := 
 With[{int = #}, 
    x /. Last@
      NMinimize[{Evaluate@
         Total[Sqrt[(x - #1)^2 + (int[x] - #2)^2] & @@@ pt], .1 < x < 
         100}, {x}]] & /@ ints; CptSum = 
 Table[{50 + Sin[u^.6] u , 40 - Cos[u^.4] u}, {u, 0, 50, 50./89}];
clus = Take[CptSum, #] & /@ Range[Length@CptSum];
imsSum = ParallelMap[With[{pt = #},
     xVals = xptsSum[pt];
     Plot[Evaluate[(#[x]) & /@ ints], {x, .1, 100}, 
      PlotRange -> {{0, 100}, {-10, 110}}, Frame -> True, 
      AspectRatio -> 1, ImageSize -> 500, 
      Epilog -> {{Directive[Opacity[.7], Gray], 
         Line[pt]}, {Directive[Opacity[.7], Black], PointSize@Small, 
         Point /@ pt}, {Red, PointSize@Large, 
         MapThread[Point@{#2, #1[#2]} &, {ints, xVals}]}}]] &, clus];
ListAnimate[imsSum, DefaultDuration -> 10]

I put together this answer and my data. It works very well. Great job, PlatoManiac. — C_S, Oct 18 '13 at 14:05

cormullion · Answer 2 · 2013-10-18T09:36:35.273

8

One way you could try is to examine the points of the plot. If cp is the plot:

points = Cases[cp, {x_?NumericQ, y_?NumericQ} :> {x, y}, Infinity];
nearest = Last@Nearest[points, First[pt], 2]

{36.7666, 51.0523}

This is not the nearest point to pt - that's pt itself - but the next nearest:

Show[
 cp,
 Graphics[{Red, Disk[nearest]}]
 ]

plot

As yves points out in the comments, this only works as well as the original plot allows:

blow up

so technically I haven't found the nearest possible point, just the nearest sampled point.

edited Oct 18 '13 at 09:36

answered Oct 18 '13 at 07:50

cormullion

24,243
4
64
133

1

Nice and concise. The only thing to be aware of is that this entirely depends on the discrete sampling of the contour lines by Plot (and associated options). – Yves Klett Oct 18 '13 at 08:55
@yves good point, I'll add a note. – cormullion Oct 18 '13 at 09:36
@cormullion very nice visual treat! – PlatoManiac Oct 18 '13 at 09:44
@PlatoManiac Well, not as cool as another answer here...! :) – cormullion Oct 18 '13 at 09:47

score 3 · Answer 3 · edited Apr 13 '17 at 12:55

My original idea was the same as cormullion's, and then followed by FindMinimum. But here another, still related way using MeshFunctions. It's not as efficient in general as other ways, perhaps, but as a method, it is sometimes convenient.

The idea is to create a mesh function that has the same sign as the angle between the vector from the pt to the curve and the normal to the curve. Then use Mesh to include the points where the angle is zero (thereby overcoming the limitation on cormullion's ContourPlot/Nearest method). To avoid including the mesh points in Nearest, I set MeshStyle -> None. The mesh is still computed, I have to select those points that are on the contours. These appear in Lines.

clevels = {3, 6, 9, 12, 15, 18, 24, 35};
test[a_, b_, c_] = a x^b (y^(1 - b)) + c;
meshfn = Numerator @ Together[
           Cross[(pt - {x, y})].D[Rationalize@test[.42, .55, .1], {{x, y}}]]
plot = ContourPlot[Evaluate@test[.42, .55, .1],
   {x, 0, 100}, {y, 0, 100},
   Contours -> clevels, 
   ContourStyle -> ColorData[1] /@ Range@Length@clevels, 
   ContourShading -> None,
   Mesh -> {{0}}, MeshStyle -> None,
   MeshFunctions -> {Function[{x, y, z}, Evaluate[meshfn]]}];
cpts = First @ Cases[plot,
                     GraphicsComplex[pts_, __] :>
                       pts[[Flatten @ Cases[plot, Line[p_] :> p, Infinity]]],
                     Infinity];
nf = Nearest @ cpts;

Show[
 plot,
 Graphics[{Blue, PointSize@Large, Point@pt, Red, Point@nf[pt]}]
 ]

Mathematica graphics

You can make a nice plot by including the mesh line.

nf2 = Nearest[Function[{x, y}, Evaluate[meshfn]] @@@ cpts -> cpts];
ContourPlot[Evaluate@test[.42, .55, .1],
 {x, 0, 100}, {y, 0, 100},
 Contours -> clevels, 
 ContourStyle -> ColorData[1] /@ Range @ Length @ clevels, 
 ContourShading -> None,
 Mesh -> {{0}},
 MeshFunctions -> {Function[{x, y, z}, Evaluate[meshfn]]},
 Epilog -> {Red, PointSize[Large], Point@pt, PointSize[Medium], 
   Thin, {Point[#], Line[{pt, #} & /@ #]} &@
    nf2[0., 2 + Length @ clevels]}]

Mathematica graphics

C_S · Answer 4 · 2013-10-18T14:33:50.610

I combined my data and answer from PlatoManiac, as follows;

(*Input Data*)
Data={{10,20,1},{10,40,2},{10,60,3},{10,80,4},{20,25,2},{20,45,3},{20,65,4},{30,30,3},{30,50,4},{40,35,4},{40,55,5},{50,20,4},{50,40,5},{60,25,5}};
U = NonlinearModelFit[Data, a x^b (y^(1 - b)) + c, {a, b, c}, {x, y}];
(*Create the implicit curves analytically*)
conts=(U[x,y]==#==#)&/@{1,2,3,4,5};
(*Given x val following function finds y=f(x) values per curve*)
fun[exp_]:={#,y/.Quiet@FindRoot[Evaluate@(exp/.x->#),{y,1}]}&;
(*Find discrete (x,f(x)) for 0<x<100*)
data=(fun[#]/@Range[.1,101,10])&/@conts;
(*point closest to which we need to find the points on the curves*)
pt={{34,52}};
(*Interpolate the data*)
ints=Interpolation/@data;
(*Given a point or a list of points following function finds the minimizer (x,f(x)) on each curve*)
xpts[pt_]:=With[{int=#},NMinimize[{Evaluate@Total[Sqrt[(x-#1)^2+(int[x]-#2)^2]&@@@pt],.1<x<100},{x}]]&/@ints;
{#2[[2]],#3[#2[[2]]]}&@@Flatten@First@Sort@Transpose@{xpts[pt],ints}
{35.9226, 53.3513}

The result can be evaluated in code;

ContourPlot[Evaluate@conts,{x,0,100},{y,0,100},ContourStyle->{Red,Magenta,Brown,Blue,Green},FrameLabel->{"X","Y"},Epilog->{{Blue,PointSize@Large,Point[pt]},{Red,PointSize@Large,Point[{35.9226,53.3513}]}}]

enter image description here

Thank you all for your answers.

Finding a point from a set of contours such that it is nearest to a given point

4 Answers4

Linked