Finding central point in a ContourPlot

Question

I am using ContourPlot for defining a 2-parameter space feasibility of a function. This function simply returns 1 if everything is True, or -1 if everything is False. My question is, without knowing the algebraic expression of the function, how can you find the point inside the feasibility region that is most distant from the border?

Just as an example, if my contour was a square, the point I would be looking for would be in the center.

ContourPlot[Sin[y] + Sin[x], {x, -4.7, 1.6}, {y, -4.7, 1.6}, 
 Contours -> 1,
 FrameLabel -> {x, y}]

Any ideas are appreciated! Thank you in advance

Answer 1

Pre-Note: For this I assume, that you change your function to return 1 inside and 0 outside the region. You can simply achieve this using Boole.

Pre-Note 2: This answer will probably not work in a real live example, because the used center of mass is not equivalent with the point most distant from the border. The purpose of my answer is to show, that it is not necessary to rely on ContourPlot.

If you have a function, and as turned out in chat you are dealing with an InterpolatingFunction, there is no need to make a ContourPlot first. One simple solution is to sample your region with Table and then calculate the center of mass. Note that the returned position is the position inside your raster, so you probably need to rescale it to your original region

data = Table[ Boole[Sin[x] + Sin[y] < 0], {y, -4.7, 1.5, 0.05}, {x, -4.7, 1.5, 0.05}];
Mean[Position[data, 1]] // N

(* {63.5, 63.5} *)

Basically, the solution through ContourPlot would use similar data, because for this plot your function is sampled as well. Although, ContourPlot uses a more fancy locally adapting sampling technique.

A solution with more possibilities is to use the data, create an image and use all the morphological measurements which already come with Mathematica

img = Image[data];
ComponentMeasurements[MorphologicalComponents[img], "Centroid"]

(* {63., 62.} *)

Again pixel coordinates. The difference in coordinates comes from the fact how morphological components treat pixels as objects and that the y axis is reversed in images. Both is easily fixable. The nice thing with this approach is, that you can use all the fancy image component measures, so e.g.

"Centroid"              center of mass coordinate
"Medoid"                coordinate of the closest element to the centroid
"MeanCentroidDistance"  mean distance of all elements from the centroid
"MaxCentroidDistance"   maximum distance of all elements from the centroid
"MinCentroidDistance"   minimum distance of all elements from the centroid

Finally, if you have already an interpolation function, nothing stops you from integrating. NIntegrate will probably complain with an interpolation function with such hard jumps from 0 to 1 and honestly, I wouldn't use this approach in real live, but it is surely possible. Assuming ip is your interpolation function, although I have here the exact expression

ip[x_, y_] := Boole[Sin[x] + Sin[y] < 0]
center = 1/NIntegrate[ip[x, y], {x, -4.7, 1.5}, {y, -4.7, 1.5}]*
  NIntegrate[ip[x, y]*{x, y}, {x, -4.7, 1.5}, {y, -4.7, 1.5}];

ContourPlot[Sin[y] + Sin[x], {x, -4.7, 1.5}, {y, -4.7, 1.5}, 
 Contours -> 1, Epilog -> {Red, PointSize[0.01], Point[center]}]

Here you get, when your interpolation function has the right region, real coordinates as you can see since I have put the calculated center directly in the plot. The above integration is equivalent to the formula from the center of mass Wikipedia site:

Answer 2

How about this...

Create the contour plot in black and white

cp = ContourPlot[Sin[y] + Sin[x], {x, -4.7, 1.6}, {y, -4.7, 1.6}, 
  Contours -> 1, ContourShading -> {White, Black}, Frame -> False, 
  PlotRangePadding -> 0, ImageSize -> 200]

Then do a DistanceTransform and locate the pixel with the maximum value.

indices = First@PixelValuePositions[#, Max@ImageData@#] &[DistanceTransform[cp]]
(* {100, 99} *)

Convert to coordinates using Rescale

centre = MapThread[Rescale[#1, {1, 200}, #2] &, {indices, PlotRange[cp]}]
(* {-1.56583, -1.59749} *)

Show[cp, Graphics[{Red, Point[centre]}]]

Answer 3

f must return only two distinct values: one for inside, the other for outside.
The first point checked, {xmin,ymin}, must be outside.

{outside,inside} = Last@Reap@Do[Sow[{x,y},f[x,y]],{x,xmin,xmax,dx},{y,ymin,ymax,dy}];
fun = Nearest[outside];
inside[[Ordering[(((fun[#,1]&/@inside)[[All, 1]]-inside)^2).{1.,1.},-1]]]

is the inside point with the largest minimum distance to a boundary.

For v7 and beyond, that have DistanceTransform, Simon's code will be faster than mine, which is for v6. The following will be much slower but will work in v5, which doesn't have Nearest.

<<Developer`
{outside,inside} = Last@Reap@Do[Sow[{x,y},f[x,y]],{x,xmin,xmax,dx},{y,ymin,ymax,dy}];
outside = -Transpose@ToPackedArray@outside;
inside[[Ordering[Min[{1.,1.}.(outside+ToPackedArray@#)^2]&/@inside,-1]]]

For the packing to help, all the coordinate pairs in inside and outside must be of the same type, Real or Integer. If they are Integer then use {1,1} instead of {1.,1.} in the last line.

Of course, the best way to speed things is to reduce the area of the region enclosing the figure. In particular, it need not be rectangular. The only requirement is that it must entirely enclose the figure.