Search code examples
matlabplotmatlab-figurecontouredge-detection

Find contour/edge in pcolor in Matlab


I'm trying to make a contour that follows the edges of the 'pixels' in a pcolor plot in Matlab. This is probably best explained in pictures. Here is a plot of my data. There is a distinct boundary between the yellow data (data==1) and the blue data (data==0):

Pcolor plot

Note that this is a pcolor plot so each 'square' is essentially a pixel. I want to return a contour that follows the faces of the yellow data pixels, not just the edge of the yellow data.

i want this output

So the output contour (green line) passes through the mid-points of the face (red dots) of the pixels.

Note that I don't want the contour to follow the centre points of the data (black dots), which would do something like this green line. This could be achieved easily with contour.

I don't want this outcome

Also, if it's any help, I have a few grids which may be useful. I have the points in the middle of the pixels (obviously, as that's what I've plotted here), I also have the points on the corners, AND I have the points on the west/east faces and the north/south faces. IF you're familiar with Arakawa grids, this is an Arakawa-C grid, so I have the rho-, u-, v- and psi- points.

I've tried interpolation, interweaving grids, and a few other things but I'm not having any luck. Any help would be HUGELY appreciated and would stop me going crazy.

Cheers, Dave

EDIT:

Sorry, I simplified the images to make what I was trying to explain more obvious, but here is a larger (zoomed out) image of the region I'm trying to separate: zoomed out image

As you can see, it's a complex outline which heads in a "southwest" direction before wrapping around and moving back "northeast". And here is the red line that I'd like to draw, through the black points:

intended output line


Solution

  • You can solve this with a couple of modifications to a solution I posted to a related question. I used a section of the sample image mask in the question for data. First, you will need to fill the holes in the mask, which you can do using imfill from the the Image Processing Toolbox:

    x = 1:15;  % X coordinates for pixels
    y = 1:17;  % Y coordinates for pixels
    mask = imfill(data, 'holes');
    

    Next, apply the method from my other answer to compute an ordered set of outline coordinates (positioned on the pixel corners):

    % Create raw triangulation data:
    [cx, cy] = meshgrid(x, y);
    xTri = bsxfun(@plus, [0; 1; 1; 0], cx(mask).');
    yTri = bsxfun(@plus, [0; 0; 1; 1], cy(mask).');
    V = [xTri(:) yTri(:)];
    F = reshape(bsxfun(@plus, [1; 2; 3; 1; 3; 4], 0:4:(4*nnz(mask)-4)), 3, []).';
    
    % Trim triangulation data:
    [V, ~, Vindex] = unique(V, 'rows');
    V = V-0.5;
    F = Vindex(F);
    
    % Create triangulation and find free edge coordinates:
    TR = triangulation(F, V);
    freeEdges = freeBoundary(TR).';
    xOutline = V(freeEdges(1, [1:end 1]), 1);  % Ordered edge x coordinates
    yOutline = V(freeEdges(1, [1:end 1]), 2);  % Ordered edge y coordinates
    

    Finally, you can get the desired coordinates at the centers of the pixel edges like so:

    ex = xOutline(1:(end-1))+diff(xOutline)./2;
    ey = yOutline(1:(end-1))+diff(yOutline)./2;
    

    And here's a plot showing the results:

    imagesc(x, y, data);
    axis equal
    set(gca, 'XLim', [0.5 0.5+size(mask, 2)], 'YLim', [0.5 0.5+size(mask, 1)]);
    hold on;
    plot(ex([1:end 1]), ey([1:end 1]), 'r', 'LineWidth', 2);
    plot(ex, ey, 'k.', 'LineWidth', 2);
    

    enter image description here