发布于

Melkman算法

[latexpage]

1. 问题描述

给定一组二维坐标的集合，求从这个集合中找出这些二维坐标的外围边界。经过查询，确定了这是一个二维凸包的问题，属于计算几何的范畴。

2. Melkman 算法

经过查询，发现二维凸包的解法有很多种，其中 Melkman 算法是其中比较好的解决方案。Melkman 算法在点排序后可以以 O(n) 的时间复杂度计算凸包。同时 Melkman 还是一种在线的算法，虽然我这里用不到这个特性。

Melkman 算法的大致思想是这样的：首先找出所有点集合中某一个维度上的最大或者最小值，然后按照顺时针或者逆时针的方向逐个遍历集合中剩余的点。要找出集合中最外围的点，有一个不变的特性必须要保持住：

假设 $P_i,P_j,P_k$ 是凸包上逆时针方向上的连续 3 个点，那么它们必然保持一种左转的趋势，即：
$$
\overrightarrow{P_iP_j}\times\overrightarrow{P_jP_k}>0
$$
如果用 $(x_i,y_i),(x_j,y_j),(x_k,y_k)$ 来表示这三个点，亦即：
$$
\overrightarrow{P_iP_j}=(x_j-x_i,y_j-y_i)
$$

$$
\overrightarrow{P_jP_k}=(x_k-x_j,y_k-y_j)
$$

那么根据向量叉积公式必须满足：
$$
(x_j-x_i)(y_k-y_j)-(x_k-x_j)(y_j-y_i)>0
$$

3. 算法描述

Melkman 算法使用双端队列来实现这个原理，假设双端队列为 $D$。所有队列操作可以用入队首、出队首、入队尾、出队尾来描述，实际操作过程可以描述为：

  1. 假设点集合的坐标可以表示为 $ (x,y) $，那么首先找出 $x$ 值最小的点记为 $P_0$。

  2. 然后选定一个方向，比如逆时针方向。然后计算所有剩余的每一个点 $P_0$ 与 $P_x$ 形成的向量 $\overrightarrow{P_0P_x}$ 与 $y$ 轴负方向的夹角。根据向量的点积公式计算出夹角之后根据夹角从小到大就行排序得到有序集合 $S={P_0,P_1,P_2…P_{n-1}}$。

  3. 记某一时刻 $D$ 的状态为：$P_tP_{t-1}…P_0…P_{b-1}P_b$ ，对于 $S$ 中的每一个点进行遍历：

    3.1 如果是 $P_0$ 则首先将 $P_0$ 入队尾，如果是 $P_1$ 则入队尾，如果是 $P_2$ 则入队首并且入队尾。

    3.2 假设遍历到当前节点 $P_i$ ：
    ​ 3.2.1 如果 $P_{b-1}P_bP_i$ 能保持左转特性则继续，否则 $P_b$ 出队尾，如此往复直到 $P_{b-m-1}P_{b-m}P_i$ 能够满足左转特性，$P_i$ 入队尾。

    ​ 3.2.2 如果 $P_iP_tP_{t-1}$ 能保持左转特性则继续，否则 $P_t$ 出队首，如此往复直到 $P_iP_{t-n}P_{t-n-1}$ 能够满足左转特性，$P_i$ 入队首。

  4. 返回 $D$。

4. 算法实现

首先给出数据结构，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 List<Point> points;

  /* 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++) {
      List<Point> points = DataExampleConfig.getPoints(String.format(fileName, i));
      Polygon polygon = pointService.getConvexHullWithMelkman(points);

      Boundary boundaryOfPolygon = GeoService.getBoundaryOfPolygon(polygon);
      Set<String> gridsOfPolygon = GridService.getGrids(boundaryOfPolygon);

      gridsOfPolygon = GridService.getGridsInPolygon(gridsOfPolygon, polygon);

      List<Point> edgeOfGrids = GridService.getEdgeOfGrids(gridsOfPolygon);

      Polygon newPolygon = new Polygon(edgeOfGrids);

      GridService.outputBoundaryGrids(boundaryOfPolygon);
      GridService.outputPoints(edgeOfGrids);

      Feature polygonFeature = GeoService.getPolygonFeature(newPolygon);

      collection.add(polygonFeature);

      List<Point> 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<Point> pointList) {
    if (pointList.size() < 3) {
      return null;
    }

    sortByAngle(pointList);

    Deque<Point> 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<Point> pointList) {
    Point minLngPoint = pointList.stream()
        .min(Comparator.comparing(Point::getLng)).get();
    pointList.remove(minLngPoint);
    pointList.sort(new Comparator<Point>() {
      @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);
  }
}

上文的代码引用了一个依赖：

<!-- geojson -->
<dependency>
    <groupId>de.grundid.opendatalab</groupId>
    <artifactId>geojson-jackson</artifactId>
    <version>1.7</version>
</dependency>

输出的 json 结果可以在 geojson 做可视化展示。

5. 参考文献

  1. http://maxgoldste.in/melkman/
  2. http://w3.impa.br/~rdcastan/Cgeometry/
  3. 漫话二维凸包