本文共 12315 字,大约阅读时间需要 41 分钟。
首先,什么是凸包?
假设平面上有p0~p12共13个点,过某些点作一个多边形,使这个多边形能把所有点都“包”起来。当这个多边形是凸多边形的时候,我们就叫它“凸包”。如下图:
然后,什么是凸包问题?
我们把这些点放在二维坐标系里面,那么每个点都能用 (x,y) 来表示。
现给出点的数目13,和各个点的坐标。求构成凸包的点?
时间复杂度:O(n³)。
思路:两点确定一条直线,如果剩余的其它点都在这条直线的同一侧,则这两个点是凸包上的点,否则就不是。
步骤:
如何判断一个点 p3 是在直线 p1p2 的左边还是右边呢?(坐标:p1(x1,y1),p2(x2,y2),p3(x3,y3))
当上式结果为正时,p3在直线 p1p2 的左侧;当结果为负时,p3在直线 p1p2 的右边。
时间复杂度:O(n㏒n)。
思路:应用分治法思想,把一个大问题分成几个结构相同的子问题,把子问题再分成几个更小的子问题……。然后我们就能用递归的方法,分别求这些子问题的解。最后把每个子问题的解“组装”成原来大问题的解。
步骤:
然而怎么求距离某直线最远的点呢?我们还是用到解一中的公式:
设有一个点 P3 和直线 P1P2 。(坐标:p1(x1,y1),p2(x2,y2),p3(x3,y3))
对上式的结果取绝对值,绝对值越大,则距离直线越远。
注意:在步骤一,如果横坐标最小的点不止一个,那么这几个点都是凸包上的点,此时上包和下包的划分就有点不同了,需要注意。
时间复杂度:O(nH)。(其中 n 是点的总个数,H 是凸包上的点的个数)
思路:
注意:
时间复杂度:O(n㏒n)
思路:Graham扫描的思想和Jarris步进法类似,也是先找到凸包上的一个点,然后从那个点开始按逆时针方向逐个找凸包上的点,但它不是利用夹角。
步骤:
最后,栈中的元素就是凸包上的点了。
以下为用Graham扫描法动态求解的过程:
给定一组二维坐标的集合,求从这个集合中找出这些二维坐标的外围边界。经过查询,确定了这是一个二维凸包的问题,属于计算几何的范畴。
Melkman
算法 经过查询,发现二维凸包的解法有很多种,其中 Melkman
算法是其中比较好的解决方案。Melkman
算法在点排序后可以以 O(n)
的时间复杂度计算凸包。同时 Melkman
还是一种在线的算法,虽然我这里用不到这个特性。
Melkman
算法的大致思想是这样的:首先找出所有点集合中某一个维度上的最大或者最小值,然后按照顺时针或者逆时针的方向逐个遍历集合中剩余的点。要找出集合中最外围的点,有一个不变的特性必须要保持住:
假设 是凸包上逆时针方向上的连续 3 个点,那么它们必然保持一种左转的趋势,即:
如果用 来表示这三个点,亦即:
那么根据向量叉积公式必须满足:
Melkman
算法使用双端队列来实现这个原理,假设双端队列为 。所有队列操作可以用入队首、出队首、入队尾、出队尾来描述,实际操作过程可以描述为:
然后选定一个方向,比如逆时针方向。然后计算所有剩余的每一个点 与 形成的向量 与 轴负方向的夹角。根据向量的点积公式计算出夹角之后根据夹角从小到大就行排序得到有序集合 。
记某一时刻 的状态为: ,对于 中的每一个点进行遍历:
3.1 如果是 则首先将 入队尾,如果是 则入队尾,如果是 则入队首并且入队尾。
3.2 假设遍历到当前节点 :
3.2.1 如果 能保持左转特性则继续,否则 出队尾,如此往复直到 能够满足左转特性, 入队尾。 3.2.2 如果 能保持左转特性则继续,否则 出队首,如此往复直到 能够满足左转特性, 入队首。
返回 。
首先给出数据结构,Point
和 Polygon
:
Point.java
:
package com.xxx.dvp.rest.domain.model;/** * User: wuzhiyu. * Date: 16/12/21. * Time: 15:18. */public class Point { private double lng; private double lat; /* Constructors */ /* Setters and Getters */}
Polygon.java
:
package com.xxx.dvp.rest.domain.model;import com.google.common.collect.Lists;import java.util.ArrayList;import java.util.Arrays;import java.util.Collections;import java.util.List;/** * The type Melkman polygon. */public class Polygon { private Listpoints; /* Constructors */ /* Setter and Getter */ /** * Contains boolean. * 源代码来自 .NET Framework。多边形边界上的点不算多边形内 * * @param point the melkman point * @return the boolean */ public boolean contains(Point point) { int crossings = 0; for (int index = 0; index < getPoints().size() - 1; index++) { double slope = getPoints().get(index + 1).getSlope( getPoints().get(index)); boolean cond1 = (getPoints().get(index).getLng() <= point.getLng()) && (point.getLng() < getPoints().get(index + 1).getLng()); boolean cond2 = (getPoints().get(index + 1).getLng() <= point.getLng()) && (point.getLng() < getPoints().get(index).getLng()); boolean above = ((point.getLng() - getPoints().get(index).getLng()) * slope + getPoints().get(index).getLat()) > point.getLat(); if ((cond1 || cond2) && above) { crossings++; } } return (crossings % 2 != 0); } /** * On border boolean. * * @param point the point * @return the boolean */ public boolean onBorder(Point point) { for (int index = 0; index < getPoints().size() - 1; index++) { double slope = getPoints().get(index + 1).getSlope(getPoints().get(index)); double latOnLine = (point.getLng() - this.points.get(index).getLng()) * slope + this.points.get(index).getLat(); if (latOnLine == point.getLat()) { return true; } } return false; } @Override public String toString() { return "Polygon{" + "points=" + points + '}'; }}
Melkman
算法:
package com.xxx.dvp.rest.domain.service;import com.xxx.dvp.rest.domain.model.Boundary;import com.xxx.dvp.rest.domain.model.Point;import com.xxx.dvp.rest.domain.model.Polygon;import com.xxx.dvp.rest.infra.persistence.conf.DataExampleConfig;import com.fasterxml.jackson.databind.ObjectMapper;import org.geojson.Feature;import org.geojson.FeatureCollection;import java.io.IOException;import java.util.ArrayDeque;import java.util.ArrayList;import java.util.Comparator;import java.util.Deque;import java.util.List;import java.util.Set;import java.util.stream.Collectors;/** * User: wuzhiyu. * Date: 16/12/26. * Time: 10:45. */public class PointService { public static void main(String[] argv) throws IOException { PointService pointService = new PointService(); FeatureCollection collection = new FeatureCollection(); String fileName = "cluster_%d.txt"; for (int i = 0; i < 8; i++) { Listpoints = DataExampleConfig.getPoints(String.format(fileName, i)); Polygon polygon = pointService.getConvexHullWithMelkman(points); Boundary boundaryOfPolygon = GeoService.getBoundaryOfPolygon(polygon); Set gridsOfPolygon = GridService.getGrids(boundaryOfPolygon); gridsOfPolygon = GridService.getGridsInPolygon(gridsOfPolygon, polygon); List edgeOfGrids = GridService.getEdgeOfGrids(gridsOfPolygon); Polygon newPolygon = new Polygon(edgeOfGrids); GridService.outputBoundaryGrids(boundaryOfPolygon); GridService.outputPoints(edgeOfGrids); Feature polygonFeature = GeoService.getPolygonFeature(newPolygon); collection.add(polygonFeature); List errorPoints = points.stream() .filter(point -> !polygon.contains(point)) .filter(point -> !polygon.onBorder(point)) .collect(Collectors.toList()); System.out.println("errorPoints = " + errorPoints); } ObjectMapper objectMapper = new ObjectMapper(); System.out.println(objectMapper.writeValueAsString(collection)); } public Polygon getConvexHullWithMelkman(List pointList) { if (pointList.size() < 3) { return null; } sortByAngle(pointList); Deque deque = new ArrayDeque<>(); pointList.forEach(point -> { switch (deque.size()) { case 0: case 1: deque.offerLast(point); return; case 2: deque.offerLast(point); deque.offerFirst(point); return; default: Point lastRightVertex = deque.pollLast(); Point lastLeftVertex = deque.pollFirst(); if (this.isLeftTurn(deque.peekLast(), lastRightVertex, point) && this.isLeftTurn(point, lastLeftVertex, deque.peekFirst())) { return; } while (!this.isLeftTurn(deque.peekLast(), lastRightVertex, point)) { lastRightVertex = deque.pollLast(); } deque.offerLast(lastRightVertex); deque.offerLast(point); while (!this.isLeftTurn(point, lastLeftVertex, deque.peekFirst())) { lastLeftVertex = deque.pollFirst(); } deque.offerFirst(lastLeftVertex); deque.offerFirst(point); return; } }); return new Polygon(new ArrayList<>(deque)); } private void sortByAngle(List pointList) { Point minLngPoint = pointList.stream() .min(Comparator.comparing(Point::getLng)).get(); pointList.remove(minLngPoint); pointList.sort(new Comparator () { @Override public int compare(Point o1, Point o2) { return Double.compare(angleWithSouth(minLngPoint, o1), angleWithSouth(minLngPoint, o2)); } }); pointList.add(0, minLngPoint); } public boolean isLeftTurn(Point point1, Point point2, Point point3) { return crossProduct(point1, point2, point3) > 0; } public double crossProduct(Point point1, Point point2, Point point3) { double x1 = point2.getLng() - point1.getLng(); double y1 = point2.getLat() - point1.getLat(); double x2 = point3.getLng() - point2.getLng(); double y2 = point3.getLat() - point2.getLat(); return x1 * y2 - x2 * y1; } public double dotProduct(Point point1, Point point2, Point point3) { double x1 = point2.getLng() - point1.getLng(); double y1 = point2.getLat() - point1.getLat(); double x2 = point3.getLng() - point2.getLng(); double y2 = point3.getLat() - point2.getLat(); return x1 * x2 + y1 * y2; } public double angleWithSouth(Point point1, Point point2) { Point point = new Point(point1.getLng(), point1.getLat() + 1); return Math.acos(this.dotProduct(point, point1, point2) / (this.norm(point, point1) * this.norm(point1, point2))); } public double norm(Point point1, Point point2) { double deltaLat = point2.getLat() - point1.getLat(); double deltaLng = point2.getLng() - point1.getLng(); return Math.sqrt(deltaLat * deltaLat + deltaLng * deltaLng); }}
上文的代码引用了一个依赖:
de.grundid.opendatalab geojson-jackson 1.7
输出的 json
结果可以在 做可视化展示。
以上讨论的只是二维的凸包,如果延生为三维、多维的凸包问题呢?如何求解?
不过首先,二维凸包可以用来解决围栏问题、城市规划问题、聚类分析等等。但是三维、多维的凸包可能的使用范畴有哪些?