查找二维数据集的区域

Question

我有一个.txt文件,在2-D平面上有大约100,000个点.当我绘制点时,有一个明确定义的二维区域(想象一下已经变形了一点的二维光盘).

计算该地区面积的最简单方法是什么？在Matlab中轻松做任何事情？

我通过在区域的边界上找到一堆(如40个)点并在Matlab中计算多边形区域的区域来进行多边形近似,但我想知道是否存在另一种不那么繁琐的方法而不是在边界上找到40个点.

Answer 1

考虑这个例子:

%# random points
x = randn(300,1);
y = randn(300,1);

%# convex hull
dt = DelaunayTri(x,y);
k = convexHull(dt);

%# area of convex hull
ar = polyarea(dt.X(k,1),dt.X(k,2))

%# plot
plot(dt.X(:,1), dt.X(:,2), '.'), hold on
fill(dt.X(k,1),dt.X(k,2), 'r', 'facealpha', 0.2);
hold off
title( sprintf('area = %g',ar) )

有一个简短的截屏视频由道格赫尔解决这个问题.

编辑:

我发布了第二个答案,灵感来自@ Jean-FrançoisCorbett提出的解决方案.

首先,我创建随机数据,并使用交互式 画笔工具,删除一些点,使其看起来像所需的"肾"形状...

要使用基线进行比较,我们可以使用IMFREEHAND函数手动跟踪封闭区域(我使用笔记本电脑的触摸板进行此操作,因此不是最精确的绘图!).然后我们使用POLYAREA找到该多边形的区域.就像我之前的回答一样,我也计算了凸包:

现在,基于之前回答的SO问题(2D直方图),我们的想法是在数据上铺设网格.网格分辨率的选择非常重要,我的是numBins = [20 30];使用的数据.

接下来我们计算包含足够点数的正方形(我至少使用1点作为阈值,但你可以尝试更高的值).最后,我们将此计数乘以一个网格的面积,以获得近似的总面积.

%### DATA ###
%# some random data
X = randn(100000,1)*1;
Y = randn(100000,1)*2;

%# HACK: remove some point to make data look like a kidney
idx = (X<-1 & -4<Y & Y<4 ); X(idx) = []; Y(idx) = [];
%# or use the brush tool
%#brush on

%### imfreehand ###
figure
line('XData',X, 'YData',Y, 'LineStyle','none', ...
    'Color','b', 'Marker','.', 'MarkerSize',1);
daspect([1 1 1])
hROI = imfreehand('Closed',true);
pos = getPosition(hROI);        %# pos = wait(hROI);
delete(hROI)

%# total area
ar1 = polyarea(pos(:,1), pos(:,2));

%# plot
hold on, plot(pos(:,1), pos(:,2), 'Color','m', 'LineWidth',2)
title('Freehand')

%### 2D histogram ###
%# center of bins
numBins = [20 30];
xbins = linspace(min(X), max(X), numBins(1));
ybins = linspace(min(Y), max(Y), numBins(2));

%# map X/Y values to bin-indices
Xi = round( interp1(xbins, 1:numBins(1), X, 'linear', 'extrap') );
Yi = round( interp1(ybins, 1:numBins(2), Y, 'linear', 'extrap') );

%# limit indices to the range [1,numBins]
Xi = max( min(Xi,numBins(1)), 1);
Yi = max( min(Yi,numBins(2)), 1);

%# count number of elements in each bin
H = accumarray([Yi(:), Xi(:)], 1, [numBins(2) numBins(1)]);

%# total area
THRESH = 0;
sqNum = sum(H(:)>THRESH);
sqArea = (xbins(2)-xbins(1)) * (ybins(2)-ybins(1));
ar2 = sqNum*sqArea;

%# plot 2D histogram/thresholded_histogram
figure, imagesc(xbins, ybins, H)
axis on, axis image, colormap hot; colorbar; %#caxis([0 500])
title( sprintf('2D Histogram, bins=[%d %d]',numBins) )
figure, imagesc(xbins, ybins, H>THRESH)
axis on, axis image, colormap gray
title( sprintf('H > %d',THRESH) )

%### convex hull ###
dt = DelaunayTri(X,Y);
k = convexHull(dt);

%# total area
ar3 = polyarea(dt.X(k,1), dt.X(k,2));

%# plot
figure, plot(X, Y, 'b.', 'MarkerSize',1), daspect([1 1 1])
hold on, fill(dt.X(k,1),dt.X(k,2), 'r', 'facealpha',0.2); hold off
title('Convex Hull')

%### plot ###
figure, hold on

%# plot histogram
imagesc(xbins, ybins, H>=1)
axis on, axis image, colormap gray

%# plot grid lines
xoff = diff(xbins(1:2))/2; yoff = diff(ybins(1:2))/2;
xv1 = repmat(xbins+xoff,[2 1]); xv1(end+1,:) = NaN;
yv1 = repmat([ybins(1)-yoff;ybins(end)+yoff;NaN],[1 size(xv1,2)]);
yv2 = repmat(ybins+yoff,[2 1]); yv2(end+1,:) = NaN;
xv2 = repmat([xbins(1)-xoff;xbins(end)+xoff;NaN],[1 size(yv2,2)]);
xgrid = [xv1(:);NaN;xv2(:)]; ygrid = [yv1(:);NaN;yv2(:)];
line(xgrid, ygrid, 'Color',[0.8 0.8 0.8], 'HandleVisibility','off')

%# plot points
h(1) = line('XData',X, 'YData',Y, 'LineStyle','none', ...
    'Color','b', 'Marker','.', 'MarkerSize',1);

%# plot convex hull
h(2) = patch('XData',dt.X(k,1), 'YData',dt.X(k,2), ...
    'LineWidth',2, 'LineStyle','-', ...
    'EdgeColor','r', 'FaceColor','r', 'FaceAlpha',0.5);

%# plot freehand polygon
h(3) = plot(pos(:,1), pos(:,2), 'g-', 'LineWidth',2);

%# compare results
title(sprintf('area_{freehand} = %g, area_{grid} = %g, area_{convex} = %g', ...
    ar1,ar2,ar3))
legend(h, {'Points' 'Convex Jull','FreeHand'})
hold off

以下是叠加的所有三种方法的最终结果,显示了面积近似值:

Answer 2

我的答案是最简单的，也许也是最不优雅和精确的。但首先，对之前的答案进行评论：

由于您的形状通常是肾形（不是凸形），因此计算其凸包的面积是行不通的，另一种方法是确定其凹包（参见例如http://www.concavehull.com/home.php ?main_menu=1）并计算其面积。但确定凹壳比确定凸壳要困难得多。另外，散乱的点都会在凸壳和凹壳中造成麻烦。

正如 @Ed Staub 的答案中所建议的，Delaunay 三角剖分后进行修剪可能会更直接一些。

我自己的建议是：表面积计算必须有多精确？我的猜测是，不是很。无论是凹形船体还是修剪过的 Delaunay 三角剖分，您都必须任意选择形状的“边界”在哪里（边缘不是刀锋利的，而且我看到周围散布着一些散布的点）它）。

因此，更简单的算法可能同样适合您的应用程序。

将图像划分为正交网格。循环遍历所有网格“像素”或方块；如果给定的正方形包含至少一个点（或者可能是两个点？），则将该正方形标记为满，否则为空。最后，添加所有完整正方形的面积。答对了。

唯一的参数是分辨率长度（方块的大小）。它的值应该设置为类似于 Delaunay 三角剖分情况下的修剪长度，即“形状内的点彼此之间的距离比该长度更近，并且距离比该长度更远的点应该被忽略”。

也许一个附加参数是一个正方形被认为已满的点数阈值。也许 2 会很好地忽略掉队点，但这可能会根据您的口味定义主要形状有点太紧...尝试 1 和 2，也许取两者的平均值。或者，使用 1 并修剪掉没有邻居的方块（生活游戏风格）。同样，8 个邻居已满的空方块应被视为已满，以避免形状中间出现洞。

该算法的改进程度是没有止境的，但由于特定应用程序中问题定义固有的任意性，任何改进都可能相当于“抛光粪便”的算法。