Detecting edge/region in contour plot

Question

I am interested in improving this plot:

which I produced with the following command:

dat // ListContourPlot[#, ContourShading -> False, 
   ContourStyle -> ColorData[10] /@ Range[10], Contours -> Range[21]/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(?)

Answer 1

How about this:

{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[10] /@ Range[10], 
   Contours -> Range[21]/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.

Answer 2

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[10] /@ Range[10], 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[6]] & // 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[[1]]]]] // 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[10] /@ Range[10], Contours -> Range[20]/8, 
      ContourLabels -> None, InterpolationOrder -> 1, FrameLabel -> {"x", "y"}, 
      DataRange -> {{-7, 2}, {-15, 15}}, PlotLabel -> "Contour plot of ϕ(x,y)", 
      Epilog -> First@Graphics[{AbsoluteThickness[3], 
         BSplineCurve[pts[[1 ;; ;; 40]] ~Join~ {Last[pts]}]}]] &