Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am interested in improving this plot:

contour 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:

badly-sorted points

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(?)

share|improve this question
    
"On the other hand these points are not sorted correctly..." - you've already tried FindCurvePath[], I presume? –  J. M. 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. Oct 1 '12 at 13:58
add comment

2 Answers 2

up vote 4 down vote accepted

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

enter image description here

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.

share|improve this answer
add comment

Here's my solution. Step-by-step explanation follows.

final solution

  • 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)"] &
    

    contour

  • 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
    

    edge

  • 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]}]}]] & 
    
share|improve this answer
add comment

Your Answer

 
discard

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.