# cbca **Repository Path**: zyflzz/cbca ## Basic Information - **Project Name**: cbca - **Description**: 基于 Zalik 和 Kolingerova 两位大佬于2001年提出的 Cell-Based Containment Algorithm (CBCA) 算法修改实现 - **Primary Language**: Java - **License**: Not specified - **Default Branch**: master - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 0 - **Forks**: 0 - **Created**: 2024-04-23 - **Last Updated**: 2024-04-23 ## Categories & Tags **Categories**: Uncategorized **Tags**: None ## README # 算法实现与测试 > 相关绘图资源（[geogebra.org](https://www.geogebra.org/classic)）源文件在目录 `geogebra_plot` 下，其中 > > - `points-example.ggb` 是图形的绝对位置 > - `points-example-innet.ggb` 是图形的相对位置 > - `points-example-rasterized.ggb` 是栅格化后的图形 ## 概述 点与多边形的相对位置问题（Relative Position of Point and Polygon，RPPP）有多种解决方案，大体上可以分为两类：一类是依靠**计算某些参数**得到判定结果，一般是借助一些几何原理进行计算；另一类是将多边形分解为一些简单的几何元素并进行有效组织，或者通过建立高效的**空间划分结构**来加速判定操作，即使用了空间换时间的思路。 本文采取的是网格法（Cell-based Containment Algorithm，CBCA），来自于 Zalik 和 Kolingerova [1] 在 2001 年的论文《A cell-based point-in-polygon algorithm suitable for large sets of points》。不过本文并没有完全采用论文中的一些方案，而是在 cell-based 的基础上进行一些改进。本算法适用于**非自交的任意单连通多边形**，所谓非自交就是多边形中的任意两条边都不会有交点（除了起止点），所谓单连通就是多边形内部没有“洞”。 CBCA 的思路如下 1. 构建网格 1. 根据多边形确定网格信息，包括：网格的大小、网格中每一个单元的大小等； 2. 确定多边形的边经过哪些单元； 3. 通过边界代数法确定哪些点在多边形外部，哪些点在多边形内部。 2. 点的判定 1. 如果目标点不在边经过的单元，那么可以直接判定，复杂度为 $O(1)$​； 2. 如果目标点在边经过的单元，由于单元左下角在多边形内外是已经确定了的，可以通过计算目标点和其所在单元左下角的点的连线，与在这个单元内的所有多边形的交点的个数的奇偶性，来判断目标点在多边形内外 按照这个思路，程序可以分为两部分，一部分是初始化，一部分是点的判定。初始化又分为三个阶段，第一个阶段是计算网格信息，第二阶段是边界确定，第三阶段是边界代数法扫描。判定部分也分为两个子分支，即根据是否在边经过的单元，进行一系列计算和判断。 但是在实现的时候，由于要知道边经过的单元有哪些，因此如果构建的网格仅仅是简单的数值型数组，操作起来就很麻烦，而且很容易造成 $O(n^2)$ 的复杂度，得不偿失。为此可以充分利用空间换时间，定义多个类型的数据结构来简化操作。例如，定义网格单元（Cell）数据结构、定义边界单元（EdgeCell）数据结构、定义边界列表（EdgeCellList）数据结构。 下面将详细阐述算法原理和代码实现逻辑。 ## 算法思路 > 具体代码实现可以直接看代码，里面有很详细的注释 ### STEP-1 网格信息计算 网格的大小称为网格分辨率，即 X 方向上多少个网格、Y 方向上多少个网格，计算公式如下 $$ NumOfCell_x = d\left[r\cdot\sqrt{n}\right]\\ NumOfCell_y = d\left[\frac{\sqrt{n}}{r}\right] $$ 其中，$d$ 表示网格的密度，是一个常数，一般为 2；$n$ 表示多边形的顶点个数；$r$ 表示 cell 的宽长比，表达式如下 $$ r=\frac{x_{max}-x_{min}}{y_{max}-y_{min}} $$ 基于网格的大小，可以得到每一个单元的单位长度 $$ SizeOfCell_x=\frac{x_{max}-x_{min}}{NumOfCell_x}\\ SizeOfCell_y=\frac{y_{max}-y_{min}}{NumOfCell_y} $$ 在有了网格基本信息的基础上，可以对多边形进行坐标变换，变换公式如下 $$ x_{new}=\frac{x-x_{min}}{x_{max}-x_{min}}\cdot (NumOfCell_x-1)\\ y_{new}=\frac{y-y_{min}}{y_{max}-y_{min}}\cdot (NumOfCell_y-1) $$ **变换坐标是为了方便操作和保证精度**。保证精度是因为，多边形的各个顶点的原始坐标可能只在小数点后面有差距，换言之就是**在原始的坐标系下的点可能是比较密集的**，如果还是用原始坐标，就很容易造成精度损失。方便操作是因为，把原始坐标变为网格中的相对坐标后，对其**只要向下取整就可以得到对应单元左下角的点**，对于程序的实现是很有用的（在程序实现时，使用一个左下角的点的坐标来表示这个单元的位置）。此外，网格一般都是定义为数组，而数组下标都是从 0 开始的，因此需要 $NumOfCell-1$ 作为乘数，而不是 $NumOfCell$​​ 作为乘数。 在网格信息计算中，有一个可以稍微优化的地方就是坐标最值的搜索，若采用二分法求最值，那么可以将局部复杂度降低到 $O(\log{n})$。 ### STEP-2 边界扫描 边界扫描的时候，需要根据需求来确定是否有序。在本算法中，是要求**边界扫描一定要有序**，一方面是程序所维护的数据结构所需，另一方面是方便边界扫描法的系列操作。基于这一要求，**输入的多边形的点也是有序的**，例如 `[0,0],[1,0],[1,1],[0,1]` 是正方形，但是如果传入 `[0,0],[1,1],[0,1],[1,0]` 就是一个沙漏形状。至于所给的多边形数组最后是否包含起点，无关紧要，因为默认最后一个点和起点相连。 边界扫描这一步的操作是：**将对应的单元设置为“边”，并保存该单元所经过的边**。即两个操作，设置操作和保存操作。 先说设置操作，对于一条有向边 $(x_1,y_1)\to(x_2,y_2)$，通过计算线段与网格的全部交点，来判断线段经过哪些单元。 首先对坐标点单元化，即向下取整，得到 $(cx_{1},cy_{1})\to(cx_{2},cy_{2})$，先考虑两种特别情况： 1. 线段在同一列网格中（$cx_1=cx_2$） 2. 线段在同一行网格中（$cy_1=cy_2$） **这两种情况下只需要依次设置起点到终点之间的所有单元即可**。需要注意的是，线段在同一列网格中，只需要 X 坐标向下取整的值相同，而无需严格的铅垂线；对于线段在同一行网格中的情况同理，无需严格的水平线。 这两种情况之外，就是一般的斜线，这一步首先要确定直线表达式，有两个参数，$k$ 和 $b$。 $$ k=\frac{y_2-y_1}{x_2-x_1},b=y_1-k\cdot x_1 $$ 然后分别计算 $y=kx+b$​ 在 X 方向上和 Y 方向上与网格线的交点，可以初步写出如下伪代码 ```text 边界点列表 points 对于 x_1 至 x_2 之间的所有点横坐标 x 计算 y = 向下取整(kx + b) points.add(x, y) 对于 y_1 至 y_2 之间的所有点纵坐标 y 计算 x = 向下取整((y - b) / k) points.add(x, y) ``` 这里需要说明一下，在**具体写代码的时候，是需要根据起止点很坐标、纵坐标的大小分成四种情况的**（或者说线段的方向有四种情况：向右上方、向左上方、向左下方、向右下方）。 而且这个伪代码的思路是有问题的，例如，若两个点为 $(1.5,1.5),(3.5,0.5)$​，可以计算出线段与网格线有如下交点（如下图） ![边界扫描的特殊情况](./README.assets/边界扫描的特殊情况.png) 在确定边在哪个单元的时候，交点也是需要向下取整的，因此点 $C,D$ 都会被判定在 $(2,1)$ 的单元中，$E$ 会被判定在 $(3,0)$ 的单元中。但是我们发现，$(2,0)$ 也应当是边，但是按照这个逻辑并没有判定为边单元，这是有问题的。问题的根源在于向下取整的逻辑。解决方法很简单，就是**在 $k<0$ 时，Y 方向扫描时，最后需要 `points.add(x,y-1)`**，而不是 `points.add(x,y)`。至于 X 方向上和 $k>0$ 的情况，向下取整的逻辑是不会有问题的。 由此可以得到如下逻辑 ```text 边界点列表 points 对于 x_1 至 x_2 之间的所有点横坐标 x 计算 y = 向下取整(kx + b) points.add(x, y) 对于 y_1 至 y_2 之间的所有点纵坐标 y 计算 x = 向下取整((y - b) / k) 如果 k > 0 points.add(x, y) 否则 points.add(x, y - 1) ``` 接下来说保存操作，保存操作一方面是把有向边 $(x_1,y_1)\to(x_2,y_2)$ 保存到对应的边单元中，例如上面的这个例子，$(1.5,1.5)\to(3.5,0.5)$ 这条边就需要保存到 $(1,1),(2,1),(2,0),(3,0)$ 对应的单元中；另一方面，是把对应的边单元保存到边列表中，这里要求边列表是有序的。因此，但是在保存之前，需要对得到的 `points` 数组进行重排序，即**按照直线方向的顺序进行排序**。 这里的 $(1,1),(2,1),(2,0),(3,0)$ 是我们观察判断得到的数据，但是按照上述逻辑得到的顺序是 $(1,1),(2,1),(3,0),(2,0)$，这是因为 X 方向和 Y 方向上的扫描是分开的，因此只能保证 X 方向有序、Y 方向有序，但是放在一起就可能无序了。进行排序的逻辑也很简单，只需要计算对应点到起点的距离，**依据距离的从小到大来确定哪些点在前、哪些点在后**。 需要注意的是，这里计算的距离不是 $AC,AD$ 之类的长度，因为都向下取整了，但是这并不妨碍 $(1,1)$ 的下一个单元是 $(2,1)$，再下一个单元才是 $(2,0)$​，所以这样的操作是可行的。 对于排序的具体逻辑，因为 X 方向有序、Y 方向有序，所以这就好像归并排序的最后一步，或者说仅使用一轮归并排序即可，时间复杂度是 $O(n)$。 ### STEP-3 边界填充 使用边界填充法确定哪些单元是内、哪些单元是外，具体来说就是每个单元的左下角的点，是在多边形的内部还是外部。操作方法很简单，基于有序的边列表进行遍历，如果是水平的，那么`continue`；**如果是上行的，那么就对当前单元左侧所有单元的值进行`+1`操作；如果是下行的，就是`-1`操作**。值得注意的是，这时候所谓的边仅仅是边单元，或者说，**这时候的多边形已经被“栅格化了”**。 ## 测试 使用如下点作为多边形 ```java double[][] polygon = new double[][]{{118.897299, 32.02599}, {118.893303, 32.026212}, {118.895235, 32.028178}, {118.894384, 32.032069}, {118.892003, 32.034575}, {118.885522, 32.034631}, {118.883062, 32.035544}, {118.885012, 32.036735}, {118.885706, 32.039075}, {118.882035, 32.041}, {118.882474, 32.042578}, {118.878943, 32.041401}, {118.877634, 32.045555}, {118.879646, 32.046496}, {118.877766, 32.048919}, {118.87954, 32.049736}, {118.881587, 32.047881}, {118.887313, 32.05076}, {118.888262, 32.054747}, {118.888244, 32.058637}, {118.885969, 32.06099}, {118.88806, 32.060962}, {118.884915, 32.067108}, {118.882658, 32.067648}, {118.881499, 32.06935}, {118.875904, 32.072395}, {118.877002, 32.075343}, {118.880225, 32.072617}, {118.883879, 32.070762}, {118.881604, 32.074167}, {118.8851, 32.075523}, {118.884959, 32.078319}, {118.888051, 32.076644}, {118.888657, 32.074707}, {118.892961, 32.07461}, {118.894173, 32.077862}, {118.893066, 32.078554}, {118.897177, 32.081045}, {118.89766, 32.082457}, {118.900918, 32.082429}, {118.90834, 32.086512}, {118.907883, 32.088103}, {118.904054, 32.090663}, {118.900637, 32.095326}, {118.897124, 32.094565}, {118.891538, 32.098066}, {118.889834, 32.098315}, {118.887006, 32.101802}, {118.883589, 32.103642}, {118.878197, 32.103255}, {118.877547, 32.100847}, {118.874323, 32.09974}, {118.871741, 32.101013}, {118.868685, 32.099671}, {118.867437, 32.100709}, {118.863854, 32.100169}, {118.859963, 32.098343}, {118.858619, 32.096627}, {118.854316, 32.09772}, {118.849037, 32.096973}, {118.847851, 32.097554}, {118.84281, 32.09707}, {118.841765, 32.099035}, {118.844698, 32.098869}, {118.843029, 32.101885}, {118.834466, 32.093998}, {118.833491, 32.095063}, {118.840307, 32.101041}, {118.839569, 32.103296}, {118.842318, 32.105275}, {118.841677, 32.106036}, {118.834431, 32.104279}, {118.831392, 32.105275}, {118.826438, 32.104901}, {118.82628, 32.103601}, {118.819149, 32.10367}, {118.818033, 32.095672}, {118.811674, 32.095575}, {118.809127, 32.09617}, {118.806606, 32.098066}, {118.802461, 32.097042}, {118.799773, 32.095133}, {118.799115, 32.093237}, {118.796831, 32.091839}, {118.797788, 32.089113}, {118.793476, 32.088408}, {118.784096, 32.092268}, {118.784104, 32.041761}, {118.827167, 32.038729}, {118.826676, 32.034174}, {118.828292, 32.0338}, {118.830066, 32.039034}, {118.831655, 32.039809}, {118.838067, 32.040058}, {118.836934, 32.035752}, {118.839745, 32.034921}, {118.838076, 32.029757}, {118.843126, 32.028801}, {118.844171, 32.027804}, {118.844285, 32.022445}, {118.845981, 32.022043}, {118.846472, 32.020008}, {118.850275, 32.018845}, {118.853367, 32.020132}, {118.85327, 32.023636}, {118.857416, 32.021573}, {118.862536, 32.023567}, {118.864838, 32.021143}, {118.875105, 32.024204}, {118.876572, 32.021116}, {118.876273, 32.017446}, {118.880638, 32.014842}, {118.884959, 32.018914}, {118.895051, 32.022847}, {118.897299, 32.02599}}; ``` 图像如下所示 ![测试案例](./README.assets/测试案例.png) 部分测试输出如下 ```text 初始化用时：31ms CBCA_BaseInfo{ n = 115 density = 2.0, numOfCellX = 30, numOfCellY = 16, ratio = 1.3624141939161363, sizeOfCellX = 0.004141466666666342, sizeOfCellY = 0.005699625000000097, xMax = 118.90834, xMin = 118.784096, yMax = 32.106036, yMin = 32.014842 } 栅格画像 X/Y | 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0 | - - - - @ @ @ @ @ @ @ @ @ - - - 1 | - - - - @ * * * * * * * @ - - - 2 | - - - - @ * * * * * * * @ - - - 3 | - - - - @ * * * * * * * @ @ - - 4 | - - - - @ * * * * * * * * @ - - 5 | - - - - @ * * * * * * * * @ - - 6 | - - - - @ * * * * * * * * @ - - 7 | - - - - @ * * * * * * * * @ - - 8 | - - - - @ * * * * * * * * * @ - 9 | - - - @ * * * * * * * * * * @ - 10 | - - - @ @ * * * * * * * * * @ - 11 | - - - - @ * * * * * * * * @ @ - 12 | - - @ @ @ * * * * * * * * @ @ - 13 | - - @ * * * * * * * * * * @ @ @ 14 | @ @ @ * * * * * * * * * * @ - - 15 | @ * * * * * * * * * * * * @ - - 16 | @ @ * * * * * * * * * * * @ - - 17 | - @ * * * * * * * * * * * @ - - 18 | - @ * * * * * * * * * * * * @ - 19 | - @ * * * * * * * * * * * @ @ - 20 | - @ * * * * * * * * * * * * @ - 21 | @ @ * * * @ * * * @ * * * @ @ - 22 | @ * * * @ @ * * @ @ * * * * @ - 23 | @ * * @ @ @ * @ @ @ @ * * * @ - 24 | @ @ * @ - @ @ @ - @ @ * * @ @ - 25 | - @ @ @ - - - - - @ @ * * @ - - 26 | - @ - - - - - - - - @ @ * @ - - 27 | - - - - - - - - - - - @ @ @ - - 28 | - - - - - - - - - - - @ @ - - - 29 | - - - - - - - - - - - @ - - - - [118.82231, 32.03611] : OUTER 用时：0ms [118.85482, 32.05551] : INNER 用时：0ms [118.82, 32.06] : INNER 用时：0ms [118.86824, 32.0135] : OUTER 用时：0ms [118.90767, 32.04096] : OUTER 用时：0ms [118.81612, 32.02166] : OUTER 用时：0ms [118.87, 32.03] : INNER 用时：0ms 1000000 个点用时: 161ms, 384939个点在多边形内部, 613641个点在多边形外部 ``` 总之如果网格单元数量级是 $1000\times1000$，初始化耗时也仅仅在 300-400 ms，判断速度远小于初始化速度。 ## 参考文献
[1] ZALIK B., KOLINGEROVA I.. A cell-based point-in-polygon algorithm suitable for large sets of points[J]. Computers & geosciences,2001,27(10):1135-1145. [doi: 10.1016/S0098-3004(01)00037-1]