Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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}]

enter image description here

Any ideas are appreciated! Thank you in advance

share|improve this question
up vote 6 down vote accepted

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]}]

Mathematica graphics

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:


share|improve this answer
+1 for a thorough answer, though I think it needs a caveat that the centroid is not necessarily the most distant point from the border, and may even lie outside the feasibility region. – Simon Woods Aug 27 '13 at 8:49
@SimonWoods Yes, you are right. My answer is only the long version of the chat discussion I had with the OP. What I wanted to stress out is, that it is not required to make this through ContourPlot when you have a function. The center of mass will not even be inside the region for certain forms. Maybe I should point this out upfront. – halirutan Aug 27 '13 at 9:18

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]

enter image description here

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]}]]

enter image description here

share|improve this answer

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.

{outside,inside} = Last@Reap@Do[Sow[{x,y},f[x,y]],{x,xmin,xmax,dx},{y,ymin,ymax,dy}];
outside = -Transpose@ToPackedArray@outside;

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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.