当前位置：文江博客话题详情

Matlab algorithm computational-geometry

求二维数据集的面积

发布于 2024-11-27 22:51:27 字数 221 浏览 1 评论 0 原文

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

计算该区域面积的最简单方法是什么？有什么方法可以在Matlab中轻松实现吗？

我通过在区域边界上找到一堆（例如 40 个）点并在 Matlab 中计算多边形区域的面积来进行多边形近似，但我想知道是否有另一种比在边界上找到 40 个点更简单的方法。

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

睡美人的小仙女 2024-12-04 22:51:27

考虑这个例子：

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

screenshot

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

这是所有三种方法叠加的最终结果，并显示面积近似值：

Consider this example:

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

screenshot

There is a short screencast By Doug Hull which solves this exact problem.

EDIT:

I am posting a second answer inspired by the solution proposed by @Jean-FrançoisCorbett.

First I create random data, and using the interactive brush tool, I remove some points to make it look like the desired "kidney" shape...

To have a baseline to compare against, we can manually trace the enclosing region using the IMFREEHAND function (I'm doing this using my laptop's touchpad, so not the most accurate drawing!). Then we find the area of this polygon using POLYAREA. Just like my previous answer, I compute the convex hull as well:

freehand
convexhull

Now, and based on a previous SO question I had answered (2D histogram), the idea is to lay a grid over the data. The choice of the grid resolution is very important, mine was numBins = [20 30]; for the data used.

Next we count the number of squares containing enough points (I used at least 1 point as threshold, but you could try a higher value). Finally we multiply this count by the area of one grid square to obtain the approximated total area.

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

Here is the final result of all three methods overlayed, with the area approximations displayed:

final

回复收藏 0 原文

吃不饱 2024-12-04 22:51:27

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

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

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

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

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

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

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

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

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

My answer is the simplest and perhaps the least elegant and precise. But first, a comment on previous answers:

Since your shape is usually kidney-shaped (not convex), calculating the area of its convex hull won't do, and an alternative is to determine its concave hull (see e.g. http://www.concavehull.com/home.php?main_menu=1) and calculate the area of that. But determining a concave hull is far more difficult than a convex hull. Plus, straggler points will cause trouble in both he convex and concave hull.

Delaunay triangulation followed by pruning, as suggested in @Ed Staub's answer, may a bit be more straightforward.

My own suggestion is this: How precise does your surface area calculation have to be? My guess is, not very. With either concave hull or pruned Delaunay triangulation, you'll have to make an arbitrary choice anyway as to where the "boundary" of your shape is (the edge isn't knife-sharp, and I see there are some straggler points sprinkled around it).

Therefore a simpler algorithm may be just as good for your application.

Divide your image in an orthogonal grid. Loop through all grid "pixels" or squares; if a given square contains at least one point (or perhaps two points?), mark the square as full, else empty. Finally, add the area of all full squares. Bingo.

The only parameter is the resolution length (size of the squares). Its value should be set to something similar to the pruning length in the case of Delaunay triangulation, i.e. "points within my shape are closer to each other than this length, and points further apart than this length should be ignored".

Perhaps an additional parameter is the number of points threshold for a square to be considered full. Maybe 2 would be good to ignore straggler points, but that may define the main shape a bit too tightly for your taste... Try both 1 and 2, and perhaps take an average of both. Or, use 1 and prune away the squares that have no neighbours (game-of-life-style). Simlarly, empty squares whose 8 neighbours are full should be considered full, to avoid holes in the middle of the shape.

There is no end to how much this algorithm can be refined, but due to the arbitrariness intrinsic to the problem definition in your particular application, any refinement is probably the algorithm equivalent of "polishing a turd".

回复收藏 0 原文

你是年少的欢喜 2024-12-04 22:51:27

我几乎一无所知，所以不要对此投入太多信心......考虑进行德劳内三角测量。然后删除任何长于某个最大值的船体（外）边缘。重复直到没有任何可删除的内容。填充剩余的三角形。

这将孤立一些异常点。

回复收藏 0 原文

↙温凉少女 2024-12-04 22:51:27

我建议使用空间填充曲线，例如 z 曲线或更好的摩尔曲线。 sfc 填充了整个空间，并且可以很好地对每个点进行索引。例如，对于所有 f(x)=y，您可以按升序对曲线点进行排序，并根据该结果获取尽可能多的点，直到获得完整的往返。然后您可以使用这些点来计算面积。因为您有很多点，所以您可能想使用更少的点并使用聚类，这会使结果不太准确。