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. 算法实现

首先给出数据结构,PointPolygon

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. 漫话二维凸包

Fuzzy C Means 算法及其 Python 实现

[latexpage]

1. $\quicklatex{size=25}K-Means $ 算法向 $\quicklatex{size=25}FCM$ 算法的扩展

在 $K-Means$ 算法中,如果要将数据集合 $X=\{ X_1,X_2,X_3,\dots,X_n \}$ 划分为 $ k\,(1 \le k \le n)$ 个类,使得任意数据对象 $X_i$ 必须属于并且仅属于一个类,同时每一个类至少包含一个数据对象,那么可以用一个 $k\times n$ 的矩阵 $U$ 来表示,矩阵中的任意一个元素 $u_{ij}$ 可以表示为:
$$ {u_{ij}} = \left{
\begin{cases}
1& {X_i \in G_j} \\
0& {X_i \notin G_j}
\end{cases}
\right.
$$
其中 ${G_j}\left( {j = 1,2, \ldots ,k} \right)$ 表示第 $j$ 个类。并且 $U$ 需要满足如下条件 $(1 \le i \le k,\,1 \le j \le n)$:
$$
\left{
\begin{cases}
{u_{ij}} \in \{0,1\} \\
\sum\limits_{i = 1}^k {{u_{ij}}} = 1 \\
\sum\limits_{j = 1}^n {{u_{ij}}} > 0
\end{cases}
\right.
$$
如果上述矩阵 $U$ 中的元素 $u_{ij}$ 的取值范围不仅仅是 0 或者 1,那么就可以推广到模糊集合上的划分,$U$ 就变成了模糊判定矩阵。此时 ${u_{ij}}$ 需满足:
\begin{equation} \label{eq:restrict}
\left{
\begin{cases}
{u_{ij}} \in [ 0,1 ]} \\
\sum\limits_{i = 1}^k {{u_{ij}}} = 1 \\
\sum\limits_{j = 1}^n {{u_{ij}}} > 0
\end{cases}
\right.
\end{equation}

2. 目标函数与聚类中心

$K-Means$ 算法在度量数据对象的非相似性(或者说距离)时一般使用欧几里得距离,要求每个类的聚类中心与数据对象的距离平方之和最小,目标函数可以表示为:
$$J = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^n {s_{ij}^2} } $$
$${s_{ij}} = Eculid({C_i},{X_j})$$
其中 $C_i$ 表示任意聚类中心,而聚类中心一般取类内所有对象在各属性上的平均值,因此可以表示为:
$$
{C_i} = \frac{{\sum\limits_{j,{X_j} \in {G_i}} {{X_j}} }}{{\sum\limits_{j = 1}^n {{u_{ij}}} }}
$$
${G_i}{\kern 1pt} \left( {1 \le i \le k} \right)$ 表示任意一个类。

将算法推广到模糊集后,$Dunn$ 对样本与类中心之间的距离采用隶属度的平方来加权,$Bezdek$ 则进一步引入了隶属度的加权指数 $m$ 从而得到了新的目标函数:
\begin{equation} \label{eq:target}
J = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^n {{{\left( {{u_{ij}}} \right)}^m}s_{ij}^2} }
\end{equation}

要使得 (\ref{eq:target}) 式达到最小值则要求聚类中心 $C_i$ 和隶属度 $u_{ij}$ 满足如下条件:
\begin{equation} \label{eq:C_i}
{C_i} = \frac{{\sum\limits_{j = 1}^n {u_{ij}^m{X_j}} }}{{\sum\limits_{j = 1}^n {u_{ij}^m} }}
\end{equation}

\begin{equation} \label{eq:u_ij}
{u_{ij}} = \frac{1}{{\sum\limits_{l = 1}^k {{{\left( {\frac{{{u_{ij}}}}{{{u_{lj}}}}} \right)}^{2/\left( {m - 1} \right)}}} }}
\end{equation}

3. $\quicklatex{size=25}FCM$ 算法计算过程

由于聚类中心和隶属度两个必要条件的存在,$FCM$ 算法的可以表示为一个简单的迭代过程。将 $n$ 个拥有 $m$ 个属性的对象构建成一个 $n\times m$ 的矩阵记为 $D$,$FCM$ 的算法可以描述为:

输入:隶属度加权指数 $m$,数据矩阵 $D$,聚类的个数 $k$,最大迭代次数 $iterMax$,阈值 $threshold$。
输出:隶属度矩阵 $U$
算法过程:
Step1: 用区间 [0,1] 内的满足 (\ref{eq:restrict}) 式的随机数初始化隶属度矩阵,记为 $D^{(0)}$。初始化迭代计数器 $iterCounter$。
Step2: 根据 (\ref{eq:C_i}) 式和 $U^{(iterCounter)}$ 计算 $k$ 个聚类中心,记做 $C^{(iterCounter)}$。
Step3: 根据 (\ref{eq:target}) 式计算目标函数的值。如果本次迭代的目标函数值与上一次计算的目标函数值之差小于 $threshold$ 或者迭代计数器 $iterCounter$ 大于等于 $iterMax$ 那么迭代结束,输出结果,算法停止。
Step4: 依据 $C^{(iterCounter)}$ 和 (\ref{eq:u_ij}) 式计算 $U^{(iterCounter+1)}$,$iterCounter++$ ,返回Step2。

4. Python 实现代码

import numpy as np

def loadDataFromTxt(fileName):
    dataSet=np.mat(np.loadtxt(fileName,delimiter='\t'))
    return dataSet

def normalization(dataSet):
    colNum=np.shape(dataSet)[1]
    for index in range(colNum):
        colMax=np.max(dataSet[:,index])
        dataSet[:,index]=dataSet[:,index]/colMax
    return dataSet

def initWithFuzzyMat(n,k):
    fuzzyMat=np.mat(np.zeros((k,n)))
    for colIndex in range(n):
        memDegreeSum=0
        randoms=np.random.rand(k-1,1)
        for rowIndex in range(k-1):
            fuzzyMat[rowIndex,colIndex]=randoms[rowIndex,0]*(1-memDegreeSum)
            memDegreeSum+=fuzzyMat[rowIndex,colIndex]
        fuzzyMat[-1,colIndex]=1-memDegreeSum
    return fuzzyMat

def eculidDistance(vectA,vectB):
    return np.sqrt(np.sum(np.power(vectA-vectB,2)))

def calCentWithFuzzyMat(dataSet,fuzzyMat,p):
    n,m=dataSet.shape
    k=fuzzyMat.shape[0]
    centroids=np.mat(np.zeros((k,m)))
    for rowIndex in range(k):
        degExpArray=np.power(fuzzyMat[rowIndex,:],p)
        denominator=np.sum(degExpArray)
        numerator=np.array(np.zeros((1,m)))
        for colIndex in range(n):
            numerator+=dataSet[colIndex]*degExpArray[0,colIndex]
        centroids[rowIndex,:]=numerator/denominator
    return centroids

def calFuzzyMatWithCent(dataSet,centroids,p):
    n,m=dataSet.shape
    c=centroids.shape[0]
    fuzzyMat=np.mat(np.zeros((c,n)))
    for rowIndex in range(c):
        for colIndex in range(n):
            d_ij=eculidDistance(centroids[rowIndex,:],dataSet[colIndex,:])
            fuzzyMat[rowIndex,colIndex]=1/np.sum([np.power(d_ij/eculidDistance(centroid,dataSet[colIndex,:]),2/(p-1)) for centroid in centroids])
    return fuzzyMat

def calTargetFunc(dataSet,fuzzyMat,centroids,k,p):
    n,m=dataSet.shape
    c=fuzzyMat.shape[0]
    targetFunc=0
    for rowIndex in range(c):
        for colIndex in range(n):
            targetFunc+=eculidDistance(centroids[rowIndex,:],dataSet[colIndex,:])**2*np.power(fuzzyMat[rowIndex,colIndex],p)
    return targetFunc

def fuzzyCMean(dataSet,k,p,initMethod=initWithFuzzyMat):
    n,m=dataSet.shape
    fuzzyMat=initWithFuzzyMat(n,k)
    centroids=calCentWithFuzzyMat(dataSet,fuzzyMat,p)
    lastTargetFunc=calTargetFunc(dataSet,fuzzyMat,centroids,k,p)

    fuzzyMat=calFuzzyMatWithCent(dataSet,centroids,p)
    centroids=calCentWithFuzzyMat(dataSet,fuzzyMat,p)
    targetFunc=calTargetFunc(dataSet,fuzzyMat,centroids,k,p)
    while lastTargetFunc*0.99>targetFunc:
        lastTargetFunc=targetFunc
        fuzzyMat=calFuzzyMatWithCent(dataSet,centroids,p)
        centroids=calCentWithFuzzyMat(dataSet,fuzzyMat,p)
        targetFunc=calTargetFunc(dataSet,fuzzyMat,centroids,k,p)
    return fuzzyMat,centroids


if __name__=='__main__':
    dataSet=loadDataFromTxt('data.txt')
    dataSet=normalization(dataSet)
    fuzzyMat,centroids=fuzzyCMean(dataSet,3,2)
    print 'fuzzyMat=\n',fuzzyMat
    print np.sum(fuzzyMat,axis=0)
    print 'centroids=\n',centroids

使用Numpy实现K-Means算法

[latexpage]

1. $ \quicklatex{size=25}K-Means $算法简介

1967年,MacQueen 首次提出了 $ K$ 均值聚类算法( $ K-Means $算法 )。迄今为止,很多聚类任务都选择该经典算法。$ K-Means $算法的优点是快速高效,其算法复杂度为 $ O(tKmn)$,其中 $ t$ 表示算法迭代的次数,$ K$ 表示聚类的数目,$ m$ 表示每个对象拥有的属性个数,$ n$ 表示待聚类的对象的个数。

该算法的核心思想是找出 $ k$ 个聚类中心,使得每一个数据点 $X_i$ 和与其最近的聚类中心 $ C_l$ 的平方距离和被最小化( 该平方距离和被称为代价函数值 $ E $ )。

$$ D(X_i,C_l) = \sum\limits_{j = 1}^{ m } {Euclid^2(X_{ij},C_{lj})} $$
$$ E = \sum\limits_{l = 1}^k {\sum\limits_{i = 1}^{{n}} {D({X_i},{C_l})} } $$

$ K-Means$算法的 $ mean$ 为类簇的中心,是一个类簇中所有对象在所有属性上的平均。在聚类的初始阶段,我们一般随机指定待聚类对象中的 $ K$ 个对象作为 $ mean$。

但是该方法也存在着一些缺陷:(1)$ K-Means$ 算法往往只能收敛到一个局部最优值。(2)聚类结果对聚类数目 $ K $ 和聚类开始时选取的 $ K$ 个初始对象非常敏感。

2. $ \quicklatex{size=25}K-Means$算法描述

输入:$ n$ 个带聚类对象,聚类中心个数 $ K$
输出:$ K$ 个类簇(包括类簇中心和类簇中的对象)
算法过程:

  1. 将 $ n$ 个对象的数据进行标准化。
  2. 随机选取 $ K$ 个对象作为初始的聚类中心。
  3. 计算每个对象到 $ K$ 个聚类中心的距离,将对象加入距离最近的类簇中。
  4. 根据类簇内的对象更新类簇的中心。
  5. 计算代价函数值 $ E $。
  6. 迭代第3步到第5步直至 $ E$收敛。

3. Python代码

下文的代码实现使用了Python的Numpy库,关于Numpy的安装和使用请参照这里

import numpy as np

def loadDataFromTxt(fileName):
    dataMat=[]
    fr=open(fileName)
    for line in fr.readlines():
        curLine=line.strip().split('\t')
        fltLine=map(float,curLine)
        dataMat.append(fltLine)
    dataMat=np.mat(dataMat)
    return dataMat

def autoNormal(dataSet):
    dataShape=np.shape(dataSet)
    n=dataShape[0]
    m=dataShape[1]
    for j in range(m):
        colMax=np.max(dataSet[:,j])
        colMin=np.min(dataSet[:,j])
        colRange=colMax-colMin
        dataSet[:,j]=(dataSet[:,j]-colMin)/colRange
    return dataSet

def randomInitCentroids(dataSet,k):
    m=np.shape(dataSet)[1]
    centroids=np.mat(np.zeros((k,m)))
    for j in range(m):
        colMin=np.min(dataSet[:,j])
        colMax=np.max(dataSet[:,j])
        rangeJ=float(colMax-colMin)
        centroids[:,j]=colMin+np.random.rand(k,1)*rangeJ
    return centroids

def euclidDistance(vectA,vectB):
    return np.sqrt(np.sum(np.power(vectA-vectB,2)))

def KMean(dataSet,k,normalMethod=autoNormal,intialMethod=randomInitCentroids,distMethod=euclidDistance):
    n=dataSet.shape[0]
    dataSet=normalMethod(dataSet)
    centroids=intialMethod(dataSet,k)
    clusterAssignment=np.mat(np.zeros((n,2)))
    clusterChanged=True
    while clusterChanged:
        clusterChanged=False
        for i in range(n):
            minIndex=-1;minDist=np.inf
            for j in range(k):
                distance=distMethod(dataSet[i,:],centroids[j,:])
                if distance<minDist:
                    minIndex=j
                    minDist=distance
            if clusterAssignment[i,0] != minIndex :
                clusterChanged=True
                clusterAssignment[i,:]=minIndex,minDist        
        for cent in range(k):
            ptsInCluster=dataSet[np.where(clusterAssignment[:,0].A==cent)[0]]
            centroids[cent,:]=np.mean(ptsInCluster,axis=0)
    return centroids,clusterAssignment

dataSet=loadDataFromTxt('data.txt')
centroids,clusterAssignment=KMean(dataSet,2)
print clusterAssignment
print centroids