常用计算几何算法

gistfile1.txt

整理一下工作学习中用到的计算几何算法，大多使用JavaScript实现。

判断点是否在多边形内部

1. Ray Casting algorithm
   该算法通过判断由待判断点为起点的一条射线与构成多边形的每个线段的交点的总数的奇偶性来判断。

/**
 * 判断点是否在面内
 * 使用RayCasting算法
 * See: https://rosettacode.org/wiki/Ray-casting_algorithm
 * @param {[number, number][]}} polygonRing
 * @param {number} px
 * @param {number} py
 * @returns
 */
export function isPointInPolygon(polygonRing, px, py) {
  var count = 0
  for (var i = 0; i < polygonRing.length; i++) {
    if (intersects(polygonRing[i], polygonRing[(i + 1) % polygonRing.length])) ++count
  }
  return !!(count % 2)

  function intersects(a, b) {
    if (a[1] <= b[1]) {
      if (py <= a[1] || py > b[1]) return false
      if (px >= a[0] && px >= b[0]) return false
      if (px < a[0] && px < b[0]) return true
      return (py - a[1]) / (px - a[0]) > (b[1] - a[1]) / (b[0] - a[0])
    }
    return intersects(b, a)
  }

判断三角形顶点环绕方向

设定三角形由点a、b、c构成。
向量B：a -> b，向量C：a -> c。
通过计算B x C得到的向量的Z值的正负可得三角形顶点环绕方向：正 -> 逆时针 负 -> 顺时针
由于点a、b、c的z值均可认为为0，所以B x C = (B.x * C.y - B.y * C.x) k
参见：https://baike.baidu.com/item/向量积/4601007 中坐标运算部分。

/**
 * 判断三角形顶点是否是CCW(counter-clock-wise)环绕
 * @param {[x: number, y: number]} a 
 * @param {[x: number, y: number]} b 
 * @param {[x: number, y: number]} c 
 */
function isCCW(a, b, c) {
    // vector: a -> b, a -> c
    var vb = [b[0]-a[0], b[1]-a[1]],
    vc = [c[0] - a[0], c[1] - a[1]]
    // a x b = (a.x * b.y - a.y * b.x)k
    return vb[0] * vc[1] - vb[1] * vc[0] > 0
}

测试

var p1 = [0,0], p2 = [1, 0], p3 = [0, 1]
console.log(isCCW(p2, p1, p3)) // false
console.log(isCCW(p3, p1, p2)) // true

/**
 * 将第三个点投影到由两个点所在的直线上
 * @param {Point} segP0
 * @param {Point} segP1
 * @param {Point} anotherPt
 */
function projectPointOnLine(segP0, segP1, anotherPt) {
  /**
   * x轴方向单位步进带来的y轴方向单位步进。
   * 当segP1[0] - segP0[0]为0时，JavaScript的计算结果为Infinity
   */
  var k = (segP1[1] - segP0[1]) / (segP1[0] - segP0[0])

  /**
   * 处理线段水平和竖直的情况
   */
  if (k === 0) {
    return [anotherPt[0], segP0[1]]
  } else if (!isFinite(k)) {
    return [segP0[0], anotherPt[1]]
  }

  /**
   * 垂线斜率
   */
  var prepK = -1 / k

  var n = anotherPt[1] - prepK * anotherPt[0]
  var m = segP0[1] - k * segP0[0]

  var ix = (n - m) / (k - prepK)
  var iy = k * ix + m

  return [ix, iy]
}

Mapping Irregular Quadrilateral to a Rectangle

求平面方程和直线交点

原坐标系下 三个两两垂直的向量构成新坐标系，求从原坐标系到新坐标系的空间变换矩阵

平面直角坐标系下情形：

坡面法向量在XOZ平面和YOZ平面上的投影线的坡度大小，可能实际用处比较小。

空间直角坐标系中点的平移、旋转、缩放后坐标计算

求经纬坐标系上某点沿经线、纬线方向前进单位距离所导致的经纬度变化量：

平面直角坐标系上点向左向右旋转90度以及以y=x的对称点：