使用ArcGIS Desktop在多边形内填充多边形？

#1 楼

有很多方法可以解决此问题。数据的栅格格式建议使用基于栅格的方法。在回顾这些方法时，将问题表示为二进制整数线性程序看起来很有希望，因为它非常符合许多GIS站点选择分析的精神，并且可以轻松地应用于它们。

在此公式中，我们列举了填充多边形的所有可能的位置和方向，我将它们称为“平铺”。与每个图块相关联的是其“优”的度量。目的是找到其总优度尽可能大的非重叠图块的集合。在这里，我们可以将每个图块的优势作为其覆盖的区域。 （在数据更加丰富和复杂的决策环境中，我们可能会通过将每个图块中包含的单元格的属性，可能与可见性相关的属性，与其他事物的接近程度等等结合起来来计算优度。）
这个问题的约束条件是解决方案中不能有两个图块重叠。

通过枚举像元，可以更有效地进行抽象化处理在要填充的多边形（“区域”）1、2，...，M中。可以使用零和一的指示符矢量对任何图块的位置进行编码，让这些图块对应于图块所覆盖的单元格，而在其他位置为零。在这种编码中，可以通过对它们的指示符向量求和（通常是逐个分量）求和来找到有关瓦片集合所需的所有信息：在至少一个瓦片覆盖一个单元的确切位置上，总和将为非零值，并且总和将更大。多于两个或多个图块重叠的位置。 （总和有效地计算出重叠量。）

还有一点抽象：可以枚举可能的瓷砖放置位置集，例如1、2，...，N。任何一组瓷砖放置位置的选择本身都对应一个指示符矢量，其中指示符矢量指定要放置的瓷砖。

这里有一个很小的插图，可以解决这些想法。它与用于进行计算的Mathematica代码相伴随，因此编程困难（或缺乏编程困难）显而易见。 br />

region =  {{0, 0, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}};

如果从上到下从左到右编号其单元格，则该区域的指标向量有16个条目：

Flatten[region]

{0，0，1，1，1，1，1，1，1，1，1，1，1，1，1，1}


让我们使用以下磁贴以及所有90度的倍数旋转：


产生旋转（和反射）的代码：

tileSet = {{{1, 1}, {1, 0}}};

（这种不透明的计算在https://math.stackexchange.com上的回复中进行了解释。 / a / 159159，该图显示了它仅会生成图块的所有可能的旋转和反射，然后删除所有重复的结果。）

假设我们将图块放置在此处，如下所示：


单元格3、6和7在此放置。由指示符向量


{0，0，1，0，0，1，1，0，0，0，0，0，0，0，0， 0}


如果将该图块向右移动一列，则该指标向量将改为


{0，0，0， 1，0，0，1，1，0，0，0，0，0，0，0，0}


试图同时在这两个位置放置图块的组合由这些指标的总和确定，


{0，0，1，1，0，1，2，1，0，0，0，0，0，0， 0，0}


第七位置的2表示这些重叠在一个单元格中（向下第二行，左起第三列）。因为我们不希望重叠，所以我们将要求任何有效解中的向量之和不得超过1。

事实证明，对于此问题，方向和位置的29个组合是可能的瓷砖。 （这是通过一些涉及穷举搜索的简单编码而发现的。）我们可以通过将其指标绘制为列向量来描述所有29种可能性。 （通常使用列代替行。）这是结果数组的图片，该数组将有16行（矩形中每个可能的单元格一个）和29列：




（前两个指示符矢量显示在左侧的前两列。）敏锐的读者可能已经注意到了进行并行处理的几种机会：这些计算可能需要花费一些时间。

所有上述内容都可以使用矩阵表示法紧凑地重述：


F是此选项数组，具有M行和N列。 > X是长度为N的一组拼贴放置的指示符。
b是一个N的矢量向量。
R是该区域的指示符；它是一个M向量。 （如果希望解决方案偏爱或避开该地区的某些区域，可以对R进行加权。）这是要最大化的。因为我们可以将其写为（R.F）.X，所以它是X的线性函数：这很重要。 （在下面的代码中，变量c包含RF）

约束是


X的所有元素必须为非负数；
X的所有元素必须小于1（这是b中的对应项）；
X的所有元素必须是整数。

约束（1）和（2）使它成为线性程序，而第三个要求使它成为整数线性程序。

存在许多用于求解以这种形式表示的整数线性程序的程序包。它们能够处理成千上万的M和N值。对于某些实际应用来说，这可能已经足够了。



作为我们的第一个插图，我使用Mathematica 8的LinearProgramming命令为上述示例计算了一个解决方案。 （这将使线性目标函数最小化。通过将目标函数取反，可以将最小值轻松地变为最大化。）它在0.011秒内返回了一个解决方案（作为图块及其位置的列表）：

apply[s_List, alpha] := Reverse /@ s;
apply[s_List, beta] := Transpose[s];
apply[s_List, g_List] := Fold[apply, s, g];
group = FoldList[Append, {}, Riffle[ConstantArray[alpha, 4], beta]];
tiles = Union[Flatten[Outer[apply[#1, #2] &, tileSet, group, 1], 1]];

灰色单元格根本不在该区域；该解决方案未涵盖白细胞。

您可以（手动）计算出许多其他与该瓷砖一样好的瓷砖-但找不到更好的瓷砖。这是这种方法的潜在局限性：即使有多个解决方案，它也可以为您提供最佳解决方案。 （有一些解决方法：如果我们对X的列进行重新排序，问题仍然没有改变，但是结果是软件经常选择其他解决方案。但是，这种行为是不可预测的。）

为了更加现实，让我们考虑问题中的区域。通过导入图像并重新采样，我用69 x 81网格表示它：



该区域包含2156个该网格的单元。 />为了使事情有趣，并说明线性编程设置的一般性，让我们尝试使用两种矩形来尽可能覆盖该区域：



一个是17 x 9（153个单元），另一个是15 x 11（165个单元）。我们可能更喜欢使用第二个，因为它较大，但是第一个更紧身，可以放在更狭窄的地方。让我们看看吧！

该程序现在涉及N = 5589个可能的图块放置。相当大！经过6.3秒的计算，Mathematica提出了以下10块解决方案：最多左边的四列），显然还有一些其他解决方案与该解决方案略有不同。

评论

---

此解决方案的早期版本（但效果不佳）出现在Mathematica网站mathematica.stackexchange.com/a/6888上。也可能值得注意的是，可以使用较小的配方变化来解决用尽可能少的图块完全覆盖区域的问题（当然允许一些重叠）：这将解决“坑洞修补”问题问题。

– hu
2012-06-18 18:09

---

为了节省空间，此答案未描述某些可能有用的改进。例如，在找到所有可能的图块位置（作为指示符向量）之后，您可以将它们全部加起来，以找出哪些图块实际上可以覆盖哪些单元格。在第二个示例中，此类单元格的集合分为两个独立的连接组件。这意味着可以在两个组件中独立解决该问题，从而大大减少了它的大小（从而减少了计算时间）。这种最初的简化对于解决现实问题往往很重要。

– hu
2012年6月19日下午0:58

---

很大的努力和答案。克里斯的回答也很有帮助。谢谢大家的帮助！有效，又让我再次朝着正确的方向前进。

– Thad
2012年6月19日12:57

---

哇！我对类似的问题感兴趣，这篇文章给了我新的观点。谢谢。如果R较大（例如140x140≈20000）怎么办，有什么方法可以降低计算成本？您知道与此问题有关的论文吗？我的搜索关键字没有以正确的方式引导我（直到现在）。

–nimcap
13年8月29日在9:13

---

@nimcap这是一类重要的问题，因此需要进行大量研究。要搜索的关键字将从“混合整数线性程序”开始，然后根据找到的内容从那里分支出来。

– hu
13年8月29日在15:29

---