在 Matlab 中查找 pcolor 中的轮廓/边缘

Question

我正在尝试在 Matlab 的 pcolor 图中创建一个遵循“像素”边缘的轮廓。这可能在图片中得到了最好的解释。这是我的数据图。黄色数据（data==1）和蓝色数据（data==0）之间有明显的界限：

请注意，这是一个pcolor 图，因此每个“正方形”本质上都是一个像素。我想返回一个跟随黄色数据像素的面的轮廓，不只是黄色数据的边缘。

所以输出轮廓（绿线）通过像素的面部（红点）中点。

请注意，我不希望轮廓跟随数据的中心点（黑点），这会像这条绿线一样。这可以通过contour 轻松实现。

另外，如果有任何帮助，我有一些网格可能有用。我有像素中间的点（显然，这就是我在这里绘制的），我也有角落上的点，并且我有西/东面和北/南面的点。如果你熟悉Arakawa grids，这是一个 Arakawa-C 网格，所以我有 rho-、u-、v- 和 psi- 点。

我尝试过插值、交织网格和其他一些东西，但我没有任何运气。任何帮助将不胜感激，并会阻止我发疯。

干杯，戴夫

编辑：

抱歉，我简化了图像以使我试图解释的内容更明显，但这是我试图分离的区域的更大（缩小）图像：

如您所见，它是一个复杂的轮廓，先向“西南”方向前进，然后环绕并向后“东北”移动。这是我想通过黑点画的红线：

Answer 1

您可以通过将a solution I posted 修改为related question 来解决此问题。我在data 的问题中使用了示例图像掩码的一部分。首先，您需要填充遮罩中的孔，您可以使用Image Processing Toolbox 中的imfill 来完成：

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

接下来，应用我的其他答案中的方法来计算一组有序的轮廓坐标（位于像素角上）：

% 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

最后，您可以像这样在像素边缘的中心获得所需的坐标：

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

这是一个显示结果的图表：

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);

Answer 2

看看下面的代码：

% 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')

它产生这个：

这是你要找的吗？

上面最重要的几行是：

x = size(data,2)-cumsum(data,2)+1;
x = x(:,end);

找到每一行在 0 到 1 之间移动的位置（假设一行中只有一个）。

然后，在plot 内，我将y 移动两个相邻 y 轴刻度之间距离的一半，因此它们将放置在边缘的中心。

---

编辑：

在对这类数据进行了一些试验后，我得到了这样的结果：

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')

它并不完美，但可以通过更简单的方法让您更接近您想要的：

Answer 3

好的，我想我已经解决了……已经很接近了，可以开心了。

首先我获取原始数据（我称之为mask_rho 并使用它来制作掩码mask_u、mask_v，这类似于mask_rho，但分别在水平和垂直方向上略有偏移。

%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

然后我制作修改后的掩码mask_u1 和mask_v1，它们与mask_rho 相同，但分别与水平和垂直方向上的相邻点进行平均。

%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;

然后我使用掩码之间的差异来定位掩码在水平方向（在 u 掩码中）和在垂直方向（在 v 掩码中）从 0 变为 1 和 1 变为 0 的位置。

% 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

并且结果抓取了边界边缘的所有坐标：

现在我的回答有点不对劲。如果我将上面的向量绘制为点plot(bnd_coords_cat(:,1),bnd_coords_cat(:,2),'kx'，我会得到上面的图像，这很好。但是，如果我加入这条线，例如：plot(bnd_coords_cat(:,1),bnd_coords_cat(:,2),'-'，那么这条线就会跳来跳去，因为这些点没有被排序。当我进行排序（使用sort 和pdist2）按最近点排序时，Matlab 有时会选择奇数点......不过我想我会将这段代码作为附录包含在内，并且是可选的额外内容。可能有人知道更好的对输出向量进行排序的方法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);

请注意，pdist2 方法是由一位同事建议的，也来自这个 SO 答案Sort coordinates points in matlab。这是最终结果：

说实话，没有线的情节很好。因此，就我而言，这已经足够接近可以回答了！