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):
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.
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 .
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: 
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:
Take a look at the following code:
% plotting some data:
data = [0 0 0 0 0 0 1 1
0 0 0 0 0 1 1 1
0 0 0 0 1 1 1 1
0 0 0 0 0 1 1 1
0 0 0 0 1 1 1 1
0 0 0 0 1 1 1 1
0 0 0 0 1 1 1 1];
p = pcolor(data);
axis ij
% compute the contour
x = size(data,2)-cumsum(data,2)+1;
x = x(:,end);
y = (1:size(data,1));
% compute the edges shift
Y = get(gca,'YTick');
y_shift = (Y(2)-Y(1))/2;
% plot it:
hold on
plot(x,y+y_shift,'g','LineWidth',3,'Marker','o',...
'MarkerFaceColor','r','MarkerEdgeColor','none')
It produces this:
Is this what you look for?
The most important lines above is:
x = size(data,2)-cumsum(data,2)+1;
x = x(:,end);
which finds the place of shifting between 0 to 1 for every row ( assuming there is only one in a row ).
Then, within the plot I shift y by half of the distance between two adjacent y-axis tick, so they will be placed at the center of the edge.
EDIT:
After some trials with this kind of data, I have got this result:
imagesc(data);
axis ij
b = bwboundaries(data.','noholes');
x = b{1}(:,1);
y = b{1}(:,2);
X = reshape(bsxfun(@plus,x,[0 -0.5 0.5]),[],1);
Y = reshape(bsxfun(@plus,y,[0 0.5 -0.5]),[],1);
k = boundary(X,Y,1);
hold on
plot(X(k),Y(k),'g','LineWidth',3,'Marker','o',...
'MarkerFaceColor','r','MarkerEdgeColor','none')
It's not perfect, but may get you closer to what you want in a more simple approach:
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);
OK, I think I've solved it... well close enough to be happy.
First I take the original data (which I call mask_rho and use this to make masks mask_u , mask_v , which is similar to mask_rho but is shifted slightly in the horizontal and vertical directions, respectively.
%make mask_u and mask_v
for i = 2:size(mask_rho,2)
for j = 1:size(mask_rho,1)
mask_u(j, i-1) = mask_rho(j, i) * mask_rho(j, i-1);
end
end
for i = 1:size(mask_rho,2)
for j = 2:size(mask_rho,1)
mask_v(j-1, i) = mask_rho(j, i) * mask_rho(j-1, i);
end
end
I then make modified masks mask_u1 and mask_v1 which are the same as mask_rho but averaged with the neighbouring points in the horizontal and vertical directions, respectively.
%make mask which is shifted E/W (u) and N/S (v)
mask_u1 = (mask_rho(1:end-1,:)+mask_rho(2:end,:))/2;
mask_v1 = (mask_rho(:,1:end-1)+mask_rho(:,2:end))/2;
Then I use the difference between the masks to locate places where the masks change from 0 to 1 and 1 to 0 in the horizontal direction (in the u mask) and in the vertical direction (in the v mask).
% mask_u-mask_u1 gives the NEXT row with a change from 0-1.
diff_mask_u=logical(mask_u-mask_u1);
lon_u_bnds=lon_u.*double(diff_mask_u);
lon_u_bnds(lon_u_bnds==0)=NaN;
lat_u_bnds=lat_u.*double(diff_mask_u);
lat_u_bnds(lat_u_bnds==0)=NaN;
lon_u_bnds(isnan(lon_u_bnds))=[];
lat_u_bnds(isnan(lat_u_bnds))=[];
%now same for changes in mask_v
diff_mask_v=logical(mask_v-mask_v1);
lon_v_bnds=lon_v.*double(diff_mask_v);
lon_v_bnds(lon_v_bnds==0)=NaN;
lat_v_bnds=lat_v.*double(diff_mask_v);
lat_v_bnds(lat_v_bnds==0)=NaN;
lon_v_bnds(isnan(lon_v_bnds))=[];
lat_v_bnds(isnan(lat_v_bnds))=[];
bnd_coords_cat = [lon_u_bnds,lon_v_bnds;lat_u_bnds,lat_v_bnds]'; %make into 2 cols, many rows
And the result grabs all the coordinates at the edges of the boundary:
Now my answer goes a bit awry. If I plot the above vector as points plot(bnd_coords_cat(:,1),bnd_coords_cat(:,2),'kx' I get the above image, which is fine. However, if I join the line, as in: plot(bnd_coords_cat(:,1),bnd_coords_cat(:,2),'-' then the line jumps around, as the points aren't sorted. When I do the sort (using sort and pdist2 ) to sort by closest points, Matlab sometimes chooses odd points... nevertheless I figured I'd include this code as an appendix, and optional extra. Someone may know a better way to sort the output vector bnds_coords_cat :
% now attempt to sort
[~,I]=sort([lon_u_bnds,lon_v_bnds]);
bnd_coords_inc1 = bnd_coords_cat(I,1);
bnd_coords_inc2 = bnd_coords_cat(I,2);
bnd_coords = [bnd_coords_inc1,bnd_coords_inc2];
bnd_coords_dist = pdist2(bnd_coords,bnd_coords);
bnd_coords_sort = nan(1,size(bnd_coords,1));
bnd_coords_sort(1)=1;
for ii=2:size(bnd_coords,1)
bnd_coords_dist(:,bnd_coords_sort(ii-1)) = Inf; %don't go backwards?
[~,closest_idx] = min(bnd_coords_dist(bnd_coords_sort(ii-1),:));
bnd_coords_sort(ii)=closest_idx;
end
bnd_coords_final(:,1)=bnd_coords(bnd_coords_sort,1);
bnd_coords_final(:,2)=bnd_coords(bnd_coords_sort,2);
Note that the pdist2 method was suggested by a colleague and also from this SO answer, Sort coordinates points in matlab. This is the final result:
To be honest, plotting without the line is fine. So as far as I'm concerned this is close enough to be answered!
链接地址: http://www.djcxy.com/p/67302.html