寻找一个二维数据集的面积

Question

寻找一个二维数据集的面积

algorithmmatlabcomputational-geometry

6

我有一个包含大约 100,000 个在二维平面上的点的 .txt 文件。当我绘制这些点时，有一个明确定义的二维区域（想象一下略微变形的二维圆盘）。

有什么最简单的方法来计算这个区域的面积？有没有在 Matlab 中轻松完成的方法？

我通过找到边界上大约 40 个点并计算 Matlab 中多边形区域的面积来进行了多边形逼近，但我想知道是否有另一种比找到边界上的 40 个点更少费力的方法。

- db1234

图片类似于：http://h1svr.anu.edu.au/static/configdb/poincare/kh0.20-phi0_small.png因此，数据集清楚地定义了一个二维区域。我的区域也不是凸的，并且没有链接中那几个小白点。此外，我认为这些点在该区域上分布并不均匀。 - db1234

我还无法回答自己的问题，因此在此发布我的答案：该区域非常小，适合于x =[0,7]和y =[0,0.5]的网格内。在矩形[0,7]x [0,0.5]上定义一个细密的网格，例如7000x500个点。在该网格上定义一个函数，如果一个数据点位于该方框内，则将该函数值设为1，否则为零。这个函数的积分就是面积。我假设Matlab能够集成一个函数，只要提供一个细密的网格？如果不能，那么实现一个黎曼和应该很简单，并且可以给出足够好的估计。 - db1234

5个回答

2

我对此知之甚少，所以不要过分相信...考虑进行Delaunay三角剖分。然后删除任何长于某个最大值的外壳（外部）边缘。重复直到没有可删除的内容。填充剩余的三角形。

这将使一些离群点成为孤立点。

- Ed Staub

看起来太繁琐了...我不需要非常精确的估计，Jean-Francois Corbett的答案已经足够好，并且易于实现。 - db1234

2

我的回答可能是最简单、最少精确和优美的。但先说一下之前的回答：

由于你的形状通常是肾脏形状（而不是凸形），计算其凸包的面积是不可行的，另一种选择是确定它的凹包（例如见http://www.concavehull.com/home.php?main_menu=1）并计算其面积。但确定凹包比确定凸包更困难。此外，流浪者点在凸包和凹包中都会造成麻烦。

Delaunay三角剖分后修剪，如@Ed Staub所建议的那样，可能会更加直截了当。

我自己的建议是：你需要多精确的表面积计算？我猜不太需要。无论是凹包还是修剪的Delaunay三角剖分，都必须对“形状”的“边界”做出任意选择（边缘不是刀锋，我看到有些流浪者点散落在周围）。因此，一个简单的算法也许对于你的应用来说同样好。

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

唯一的参数是分辨率长度（方块的大小）。它的值应设置为类似于Delaunay三角剖分的修剪长度，即“我的形状中的点彼此之间距离比这个长度近，而相互分开的点应被忽略”。也许额外的参数是一个方格被认为是满的点数阈值。也许2个会好一些来忽略流浪者点，但可能太紧密地定义了主要形状... 试着同时使用1和2，并取平均值。或者，使用1并修剪掉那些没有邻居的方块（生命游戏风格）。类似地，那些8个邻居都是满的空方格应被视为满的，以避免在形状中间出现空洞。

这个算法的改进没有止境，但由于问题定义本身的任意性，任何改进可能等价于“擦亮一坨屎”的算法。

- Jean-François Corbett

我同意这是正确的方法。如果您的区域中有任何空洞，这可能会破坏总面积，但可以通过形态学闭合操作轻松解决。 - Jonathan

这基本上就是我说的，你实际上是在伪装中做一个黎曼和...请看问题下面的评论。由于我没有足够的“声望点数”，所以无法将其发布为答案。 - db1234

0

我建议使用空间填充曲线，例如z曲线或更好的moore曲线。 SFC填满整个空间，并且适合索引每个点。例如，对于所有f（x）= y，您可以按升序对曲线上的点进行排序，并从该结果中取出尽可能多的点，直到获得完整的往返行程。然后，您可以使用这些点来计算面积。由于有许多点，您可能希望使用较少的点并使用聚类使结果不太准确。

- Micromega

0

我认为你可以使用凸包算法来获取边界点，限制边长（应该按垂直轴对点进行排序）。这样它将遵循您区域的非凸性。我建议长度为0.02。无论如何，您可以尝试使用不同的长度绘制结果并在视觉上检查它。

- Sergey Bankevich

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Amro · Accepted Answer

考虑以下例子：

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

有一个由Doug Hull制作的短视频，解决了这个精确问题。

编辑：

我发表了第二篇答案，受到@Jean-FrançoisCorbett提出的解决方案的启发。

首先，我创建了随机数据，并使用交互式刷子工具，删除了一些点，使其看起来像所需的“肾”形状...

为了有一个基准比较，我们可以手动跟踪包围区域，使用IMFREEHAND函数（我使用的是笔记本电脑的触摸板，所以绘画不太准确！）。然后我们使用POLYAREA计算此多边形的面积。就像我的上一个答案一样，我也计算了凸包：

freehand convexhull

现在，基于我之前回答的另一个SO问题（2D直方图），想法是在数据上放置网格。网格分辨率的选择非常重要，我的数据使用numBins = [20 30];。

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

hist2d hist2d_threshold

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

以下是三种方法的最终结果，同时显示面积近似值：

final