# Detecting edge/region in contour plot

I am interested in improving this plot: which I produced with the following command:

dat // ListContourPlot[#, ContourShading -> False,
ContourStyle -> ColorData /@ Range, Contours -> Range/8,
ContourLabels -> None, FrameLabel -> {x, y},
DataRange -> {{-7, 2}, {-15, 15}}, PlotLabel -> "Contour plot of ϕ(x,y)"] &


The actual dat is available here; the edge corresponds to zero values (though, it could be set to something else).

Question

I would like Mathematica to identify the edge of the contours and do a region plot to avoid the ragged line, or alternatively, to draw a thick line over it to make it publication-friendly.

Attempt

I guess I know how to find the edge:

dat2 = dat // Image // EdgeDetect // ImageData // Position[#, 1] & // Sort;


On the other hand, these points are not sorted correctly: I supposed I could use 2D neighbours as a criterion for sorting, but I feel there is a smarter way to achieve my overall goal(?)

• "On the other hand these points are not sorted correctly..." - you've already tried FindCurvePath[], I presume? – J. M. will be back soon Oct 1 '12 at 13:41
• @J.M. no I didn't know about FindCurvePath... arg! – chris Oct 1 '12 at 13:48
• @J.M. FindCurvePath is not producing miracles, probably because the curve is too jagged? – chris Oct 1 '12 at 13:54
• Yes, it doesn't always work (but when it does, it's great); that's why I was asking if you've seen it... – J. M. will be back soon Oct 1 '12 at 13:58

{d1, d2} = Dimensions[dat];
xvals = Range[-15, 15, (15 + 15)/(d1 - 1)];
yvals = Range[-7, 2, (2 + 7)/(d2 - 1)];
dat2 = Flatten[Join[Table[{yvals[[j]], xvals[[i]]}, {i, d1}, {j, d2}],
Table[Partition[dat[[i]], 1], {i, d1}], 3], 1];
dat2 = Extract[dat2, Position[dat2[[All, 3]], x_ /; x != 0]];
lcp = ListContourPlot[dat2, ContourShading -> False,
ContourStyle -> ColorData /@ Range,
Contours -> Range/8, ContourLabels -> None,
FrameLabel -> {x, y}, DataRange -> Automatic,
PlotLabel -> "Contour plot of \[Phi](x,y)",
BoundaryStyle -> Directive[Black, Thick]] The "tricks" are to pass ListContourPlot a list of tuples where the zero valued tuples are removed, and to set DataRange -> Automatic so that it interprets the list as a list of {x, y, f} tuples instead of as many datasets.

Here's my solution. Step-by-step explanation follows. • First, import the data

dat = ToExpression@Import["http://pastebin.com/raw.php?i=XWyb7jFJ"];

• Instead of using EdgeDetect, I'll just use the outermost contour (on a fine mesh) as the "edge":

contour =  dat // ListContourPlot[#, ContourShading -> False,
ContourStyle -> ColorData /@ Range, Contours -> {1}/20,
ContourLabels -> None, Frame -> False, InterpolationOrder -> 1,
FrameLabel :> {x, y}, DataRange -> {{-7, 2}, {-15, 15}},
PlotLabel -> "Contour plot of ϕ(x,y)"] & • Next, I use the image processing functions to close the gap, get rid of the egg-shaped contour and thin the edge.

edge = contour // Image // Binarize // ColorNegate //
SelectComponents[#, "Elongation", 0.5 < # &] & // DeleteSmallComponents //
Closing[#, BoxMatrix] & // Thinning • Get the positions of the 1s (this forms the curve in the binary image) and then use FindCurvePath to get the proper ordering for the curve

pos = Position[ImageData@edge, 1 | 1.];
curve = FindCurvePath@pos;

• Next, convert the positions and ordering from filled curve to ListPlot coordinates.

pts = N@With[{rescale = Rescale[#, Through[{Min, Max}@#], #2] &},
{rescale[#2, {-7, 2.2}], rescale[-#1, {-14.5, 12.5}]} & @@
Transpose[pos[[curve[]]]] // Transpose];


Here, I've eyeballed the rescaling from the extents of the original image. It is possible to get these programmatically, but from my trials, some amount of fine tuning is eventually necessary to get a nice fit.

• Finally, downsample the points and form a BSplineCurve and overlay on the original image

dat // ListContourPlot[#, ContourShading -> False,
ContourStyle -> ColorData /@ Range, Contours -> Range/8,
ContourLabels -> None, InterpolationOrder -> 1, FrameLabel -> {"x", "y"},
DataRange -> {{-7, 2}, {-15, 15}}, PlotLabel -> "Contour plot of ϕ(x,y)",
Epilog -> First@Graphics[{AbsoluteThickness,
BSplineCurve[pts[[1 ;; ;; 40]] ~Join~ {Last[pts]}]}]] &