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

Question

我正在尝试制作一个轮廓，其跟随Matlab中pcolor绘图中的像素边缘。这最好通过图片来解释。这是我数据的一个绘图。黄色数据（data==1）和蓝色数据（data==0）之间有一个明显的边界：请注意，这是一个pcolor绘图，因此每个“方块”本质上都是一个像素。我想要返回一个轮廓，该轮廓遵循黄色数据像素的面，而不仅仅是黄色数据的边缘。因此，输出轮廓（绿线）穿过像素面的中点（红点）。请注意，我不希望轮廓沿着数据中心点（黑点）走，这将形成类似于这条绿线的形状。使用contour可以轻松实现此目标。如果有帮助的话，我有一些网格可能会有用。我有像素中心点的坐标（显然，因为这是我在此处绘制的），我还有角上的点，并且我有西/东面和南/北面的点。如果您熟悉Arakawa网格，那么这是一个Arakawa-C网格，因此我有rho、u、v和psi点。

我尝试过插值、交错网格和其他一些方法，但没有成功。非常感谢任何帮助，这将使我避免变得疯狂。

编辑：

对不起，我简化了图片，以使我试图解释更加明显，但这里是我正在分离的区域的更大（缩小）图像：如您所见，它是一个复杂的轮廓，在“向西南”方向前进，然后绕回并向“向东北”移动。这里是我想要绘制的红线通过黑点：

- David_G

Answer 1

您可以通过对我之前发布的解决方案进行一些修改来解决这个问题，该解决方案是针对一个相关问题的。我在问题中使用了样本图像掩模的一部分作为数据。首先，您需要填补掩模中的空洞，您可以使用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);

- gnovice

---

我喜欢这个答案，它的方向与我原本想法很不一样。我的答案解决了我的问题，但你的答案似乎更清晰易懂 - 所以我接受了你的答案！ - David_G

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

虽然不是完美的，但可能以更简单的方式让您更接近所需：

- EBH

---

非常好的建议，但不幸的是输入数据有点更复杂 :(（请参见编辑）。但是您的优秀建议可能无法在数据点之间的线是水平的情况下起作用，就像我添加的新图片中一样？ - David_G

@David_G 很抱歉我提供了一个太具体的解决方案，现在我需要澄清一些事情 - 新图像是全部图像吗？只有一个被蓝色包围的黄色区域吗？换句话说，您如何定义应该绘制红线的位置？为什么中间的蓝洞里没有红线呢？ - EBH

新图像是完整矩阵的放大版。黄色区域有效地定义了陆地/水域存在的掩模。因此，黄色（==1）表示水，蓝色（==0）表示陆地。我想要边界线。有一个大的黄色区域，南部边缘有一些蓝色，还有一些孤立的蓝色区域。捕捉到你描述的那个孤立的蓝洞并不是很关键。如果你也在那里画一个红线，那就可以了。希望这能帮到你。 - David_G

@David_G 抱歉耽搁了，我没有时间回来处理这个问题。请看我的编辑，又尝试了一部分。 - EBH

你的图像工具箱方法没有给我足够接近的结果。话虽如此，非常感谢您让我认识它！我以前从未使用过bwboundaries等工具，但看起来非常强大。至少我想选中您的答案或为您提供额外的积分！ - David_G

@David_G 不用谢。我现在没有时间继续处理这个问题，但是如果你将两种方法结合起来，可能会得到最佳解决方案。使用你的方法识别所有点（在你的答案中标记为红色x），然后使用“boundary”将它们连接起来。唯一剩下的问题是丢弃蓝色的“岛屿”。如果我很快有更多时间，我会尝试回到这个问题上。 - EBH

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。这是最终结果：

说实话，没有线条的绘图也可以。就我而言，这已经足够接近答案了！

- David_G

---

这种方法对我很有效 - 所以我认为这个答案不错。然而，它有点笨拙，而@gnovice的答案更清晰，所以我会接受那个答案。 - David_G