經常使用空間分析函數

經常使用空間分析函數

包含 點、線、面 之間的相互關係

本文主要記錄一些本身開發中經常使用的一些分析函數，比較正宗和全面的能夠看一些空間分析庫，好比前端的 Turf 和 JSTS

一、點

1.一、點到點的距離

//點到點的距離
function dist2d(coord1, coord2) {
  let dx = coord1[0] - coord2[0];
  let dy = coord1[1] - coord2[1];
  return Math.sqrt(dx * dx + dy * dy)
}
dist2d([1, 1], [4, 5]) // 5

1.二、判斷兩個點集是否相等

//判斷兩個點集是否相等
function equals(coord1, coord2) {
  let equals = true;
  for (let i = coord1.length - 1; i >= 0; --i){
    if (coord1[i] != coord2[i]) {
      equals = false;
      break
    }
  }
  return equals
}
equals([1, 3], [1, 3]) // true
equals([1, 3, 1], [1, 3, 2]) // false

二、線

2.一、線段的長度

// 線段的長度
function formatLength(coords) {
  coords = coords || []
  let length = 0
  //經過遍歷座標計算兩點以前距離，進而獲得整條線的長度
  for (let i = 0, leng = coords.length - 1; i < leng; i++) {
    length += dist2d(coords[i], coords[i + 1]);
  }
  return length
}
formatLength([[1, 0], [2, 1], [3, 0], [4, 2]]) // 5.06449510224598

2.二、根據偏移量偏移一條線段

// 根據偏移量偏移一條線段
function lineOffset(coords, deltaX, deltaY){
  deltaX = deltaX || 0
  deltaX = isNaN(deltaX) ? 0 : deltaX
  deltaY = deltaY || 0
  deltaY = isNaN(deltaY) ? 0 : deltaY

  if(deltaX == 0 && deltaY == 0)return coords

  coords.forEach(coord => {
    coord[0] += deltaX;
    coord[1] += deltaY;
  })
  return coords
}
lineOffset([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]], 2, 3)
//         [[3, 4], [4, 6], [5, 6], [6, 5], [4, 3]]

2.三、線段上距離點P最近的一個點

// 線段上距離點P最近的一個點
function getNearestCoord(point, lines) {
  let d, res = { dist: Infinity }
    if(!(Array.isArray(lines[0]) && Array.isArray(lines[0][0]))){
        lines = [lines]
    }
    lines.forEach(function(coords){
        for (let i = 0; i < coords.length; i++) {
            d = dist2d(point, coords[i]);
            if (d < res.dist) {
        res.dist = d;
        res.index = i
        res.point = coords[i]
            }
        }
    })
    
    return res.point ? res : null;
}
getNearestCoord([2.2, 3.1], [[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// {dist: 0.5099019513592785, index: 1, point: [2, 3]}

2.四、點P到線段AB的最短距離

/* 
* 點P到線段AB的最短距離
* 使用矢量算法，計算線AP在線段AB方向上的投影，當須要計算的數據量很大時，這種方式優點明顯
* 特殊狀況如點在線段上、點在端點、點在線段延長線上等等的狀況所有適用於此公式，只是做爲特殊狀況出現，無需另做討論。這也是矢量算法思想的優點所在。
* 函數返回值：point 投影座標  dist 點P到投影距離  type 垂足位置，不爲0表示垂足在線段外
*/
function pointToSegmentDist(point, point1, point2){
  let x = point[0], x1 = point1[0], x2 = point2[0]
  let y = point[1], y1 = point1[1], y2 = point2[1]

  //線段AB 爲一個點
  if(x1 == x2 && y1 == y2) return {
    type: 0,
    point: point1,
    dist: 0
  }

  let cross = (x2 - x1) * (x - x1) + (y2 - y1) * (y - y1);
  //let r = cross / d2
  //r < 0 點P的垂足在線段AB外，且點P距離線段AB最近的點爲A
  //r = 0 點P的垂足和點P距離線段AB最近的點爲A
  if (cross <= 0) return {
    type: 1,
    point: point1,
    dist: Math.sqrt((x - x1) * (x - x1) + (y - y1) * (y - y1))
  };
    
  let d2 = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
  //r > 1 點P的垂足在線段AB外，且點P距離線段AB最近的點爲B
  //r = 1 點P的垂足和點P距離線段AB最近的點爲B
  if (cross >= d2) return {
    type: 2,
    point: point2,
    dist: Math.sqrt((x - x2) * (x - x2) + (y - y2) * (y - y2))
  };
    
  let r = cross / d2;
  let px = x1 + (x2 - x1) * r;
  let py = y1 + (y2 - y1) * r;
  return {
    type: 0,
    point: [px, py],
    dist: Math.sqrt((x - px) * (x - px) + (py - y) * (py - y))
  };
}
pointToSegmentDist([1, 3], [0, 1], [3, 1]); 
//{ dist: 2, point: [1, 1], type: 0 }

2.五、點P到直線AB的垂足及最短距離

/* 
* pointToSegmentDist變種：點P到直線AB的垂足及最短距離
*/
function pointToFootDist(point, point1, point2){
  let x = point[0], x1 = point1[0], x2 = point2[0]
  let y = point[1], y1 = point1[1], y2 = point2[1]

  //線段AB 爲一個點
  if(x1 == x2 && y1 == y2) return {
    type: true,
    point: point1,
    dist: 0
  }

  let dx = x2 - x1;
  let dy = y2 - y1;
  //r < 0 點P的垂足在線段AB外，且點P距離線段AB最近的點爲A
  //r > 1 點P的垂足在線段AB外，且點P距離線段AB最近的點爲B
  //r = 0 點P的垂足和點P距離線段AB最近的點爲A
  //r = 1 點P的垂足和點P距離線段AB最近的點爲B
  let r = (dx * (x0 - x1) + dy * (y0 - y1)) / (dx * dx + dy * dy || 0);

  let px = x1 + dx * r;
  let py = y1 + dy * r;
  return {
    type: r >= 0 && r <= 1, //true  垂足在線段內   false 垂足在線段外
    point: [px, py], //垂足
    dist: Math.sqrt((x - px) * (x - px) + (py - y) * (py - y)) //點P到垂足距離
  };
}
pointToFootDist([3, 1], [1, 0], [3, 0])
// {dist: 1, point: [3, 0], type: true}

2.六、計算點P到直線AB的距離

/*
* 計算點P到直線AB的距離
* 使用三角形面積，計算線AP在直線AB方向上的投影
* Area   = |(1/2)(x1 * y2 + x2 * y3 + x3 * y1 - x2 * y1 - x3 * y2 - x1 * y3)|
* Bottom = Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
* Area = .5 * Bottom * H
* Height = Area / .5 / Bottom
*/
function verticalDistance(point, point1, point2){
  if(equals(point1, point2)){
    return dist2d(point, point1)
  }
  let area = Math.abs(0.5 * (point1[0] * point2[1] + point2[0] * point[1] + point[0] * point1[1] - point2[0] * point1[1] - point[0] * point2[1] - point1[0] * point[1]));
  let bottom = dist2d(point1, point2)
  let height = area / 0.5 / bottom;

  return height;
}
verticalDistance([4, 2], [1, 0], [3, 0]) // 2

2.七、計算點P到直線AB的距離

/* 
* 計算點P到直線AB的距離
* 根據公式，計算線AP在直線AB方向上的投影
* d = Math.abs((A * x + B * y + C) / Math.sqrt( A * A + B * B))
* 公式介紹：https://zhuanlan.zhihu.com/p/26307123
* 直線公式：y = kx + b
* 直線公式：kx - y + b = 0
* 假設 B = -1，則 A = k，C = b   
* 直線斜率：k = (y1 - y2) / (x1 - x2)
* 常數計算：b = y - kx
*/
function verticalDistance2(point, point1, point2){
    //若是point1[0] == point2[0] 說明是條豎着的線
    if(point1[0] == point2[0]){
        return Math.abs(point[0] - point1[0])
    }else{
    let k = (point1[1] - point2[1]) / (point1[0] - point2[0])
    let b = point1[1] - k * point1[0]
        return Math.abs((k * point[0] - point[1] + b) / Math.sqrt( k * k + 1))
    }
}
verticalDistance2([4, 2], [1, 0], [3, 0]) // 2

2.八、計算緯度和經度指定的兩點之間的距離(米)

/* 
* 計算緯度和經度指定的兩點之間的距離(米)
* 在計算涉及到經緯度的計算時(長度、面積等)，若是隻是粗略的計算，那麼咱們使用座標點計算，將計算結果按照每經緯度111km轉換便可。這個結果越往兩極失真越大
* 若是想要精細點，就須要用到地球長半軸參數去計算
* WGS84 座標系下的長半軸參數 radius = 6378137 m
* WGS84(World Geodetic System 1984) 座標系是爲GPS全球定位系統使用而創建的座標系統。是國際公認的座標系
* 北京54座標系下的長半軸參數 radius = 6378245 m
* 1954年北京座標系能夠認爲是前蘇聯1942年座標系的延伸。它的原點不在北京而是在前蘇聯的普爾科沃。因此偏差較大，缺點較多。
* 西安80座標系下的長半軸參數 radius = 6378140 ± 5 m
* 1980西安座標系所採用的IAG1975橢球，其長半軸要比WGS84橢球長半軸的值大3米左右，而這可能引發地表長度偏差達10倍左右
* CGCS2000(2000 國家大地座標系) 長半軸參數 radius = 6378137 m
* 2000國家大地座標系，是我國當前最新的國家大地座標系，是全球地心座標系在我國的具體體現，其原點爲包括海洋和大氣的整個地球的質量中心。
* 正弦曲線投影下的長半軸參數 radius = 6370997 m
* 正弦曲線投影是一種等面積的僞圓柱投影。規定緯線投影爲平行直線，經線投影爲對稱於中央經線的正弦曲線，同一緯線上經距相等，緯距向兩極縮小。主要用於小比例尺世界地圖
*/
function haversineDistance(c1, c2) {
  let radius = 6378137
  let lat1 = toRadians(c1[1]);
  let lat2 = toRadians(c2[1]);
  let deltaLatBy2 = (lat2 - lat1) / 2;
  let deltaLonBy2 = toRadians(c2[0] - c1[0]) / 2;
  let a = Math.sin(deltaLatBy2) * Math.sin(deltaLatBy2) + Math.sin(deltaLonBy2) * Math.sin(deltaLonBy2) * Math.cos(lat1) * Math.cos(lat2);
  return 2 * radius * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a))
}
function toRadians(_){
  return _ * Math.PI / 180
}
haversineDistance([117, 34], [118, 39])  //563732.8911197125米
haversineDistance([118, 39], [117, 34])  //563732.8911197125米

2.九、根據一個點切割線段

/* 
* 根據一個點切割線段
* 點沒必要落在線段上
* 經常使用於計算鼠標位置距離線段起點或終點的距離
* 簡單版的是使用 getNearestCoord，而不是 pointToSegmentDist
* getNearestCoord 尋找的是線上真實存在的一點
* pointToSegmentDist 尋找的點位於線上，但不必定真實存在與線上
*/
function getClosestPoint(coords, point){
  if(!coords || coords.length < 2) return [[], []]

  let squaredDistance = {dist: Infinity}, index
  for(let i = 1; i < coords.length; i++){
    let d = pointToSegmentDist(point, coords[i - 1], coords[i])
    if(d.dist < squaredDistance.dist){
      squaredDistance = d
      index = i
    }
  }

  if(index === undefined){
    return [coords, []]
  }

  let prearr = coords.slice(0, index)
  if(prearr.length && !equals(squaredDistance.point, prearr[prearr.length - 1])){
    prearr.push(squaredDistance.point)
  }

  let nextarr = coords.slice(index)
  if(nextarr.length && !equals(squaredDistance.point, nextarr[0])){
    nextarr.unshift(squaredDistance.point)
  }

  return [prearr, nextarr]
}
getClosestPoint([[0, 1], [2, 3], [5, 6], [3, 4], [4, 4]], [3, 5])
/*
[
    [[0,1],[2,3],[3.5,4.5]],
    [[3.5,4.5],[5,6],[3,4],[4,4]]
]
*/

2.十、線段與線段的交點

/* 
* 線段與線段的交點
* 解線性方程組, 求線段AB與線段CD的交點座標，若是沒有交點，返回null
*/
function intersects(coords1, coords2) {
  let x1 = coords1[0][0];
  let y1 = coords1[0][1];
  let x2 = coords1[1][0];
  let y2 = coords1[1][1];
  let x3 = coords2[0][0];
  let y3 = coords2[0][1];
  let x4 = coords2[1][0];
  let y4 = coords2[1][1];
  //斜率交叉相乘 k1 = (y4 - y3) / (x4 - x3)    k2 = (y2 - y1) / (x2 - x1)
  //k1 k2 同乘 (x4 - x3) * (x2 - x1) 並相減獲得denom
  let denom = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1)); 
  // 若是分母爲0 則平行或共線, 不相交 
  if (denom === 0) {
    return null;
  }

  let numeA = ((x4 - x3) * (y1 - y3)) - ((y4 - y3) * (x1 - x3));
  let numeB = ((x2 - x1) * (y1 - y3)) - ((y2 - y1) * (x1 - x3));
  let uA = numeA / denom;
  let uB = numeB / denom;

  // 交點在線段1上，且交點也在線段2上  
  if (uA >= 0 && uA <= 1 && uB >= 0 && uB <= 1) {
    let x = x1 + (uA * (x2 - x1));
    let y = y1 + (uA * (y2 - y1));
    return [x, y];
  }
  return null;
}
intersects([[0,0],[1,1]], [[3,0],[2,1]])  //null
intersects([[0,0],[1,1]], [[3,0],[0,1]])  //[0.75, 0.75]

2.十一、線與線的全部交點

/* 
* 線與線的全部交點
* 接收兩條線的點集合，求取line2 和 line1的交點
* 若是規定count參數，表明返回前多少個交點，不然所有返回
* 返回交點集合，集合中的元素屬性包括 index1：交點在line1的位置，index2：交點在line2的位置，coords：交點座標
* 代碼參考truf.lineIntersect：http://turfjs.org/docs/#lineIntersect
* 若是交點爲線的點，則會重複返回，包括返回的數據結構，都是爲下面切割面函數splitPolygon服務的
*/
function lineIntersect(line1, line2, count){
  let result = []
    for(let i = 1; i < line1.length; i++){
    let coords1 = [line1[i-1], line1[i]]
    //求取數據的邊界範圍，函數放在了下面的多邊形中
    let bbox1 = bbox(coords1)
        for(let j = 1; j < line2.length; j++){
      let coords2 = [line2[j-1], line2[j]]
      let bbox2 = bbox(coords2)
      //判斷兩個邊界範圍的關係： bbox1 是否包含 bbox2，函數放在了下面的多邊形中
      if(isIntersects(bbox1, bbox2)){
        let p = intersects(coords1, coords2)
        if(p){
          result.push({index1: i, index2: j, coords: p})
          if(count && (result.length >= count)){
            return result
          }
        }
      }
        }
  }
    return result
}
lineIntersect([[0, 0], [0, 1], [1, 3], [3, 2], [5, 0]], [[-1, 0], [4, 3]])
//[{"index1":1,"index2":1,"coords":[0,0.6]},{"index1":3,"index2":1,"coords":[2.6363636363636367,2.1818181818181817]}]

2.十二、判斷線段路徑是順時針仍是逆時針

/* 
* 判斷線段路徑是順時針仍是逆時針
* 使用格林公式計算
* 返回值 d < 0  '順時針' : d > 0  '逆時針'
* truf對應實現truf.booleanClockwise://turfjs.org/docs/#booleanClockwise
*/
function isClockwise(line){
  let length = line.length
  let d = 0

  if(length.length < 3)return d;

  //若是不是閉合線，則改成閉合線
  if(!equals(line[0], line[length - 1])){
    length = (line = line.slice()).push(line[0]);
  }

  //沿着多邊形的邊求曲線積分,若積分爲正,則是沿着邊界曲線正方向(逆時針),反之爲順時針
  //最後一個點是開始點，只參與末尾點(倒數第二個點)的計算，不單獨計算
  for(let i = 0; i < length - 1; i++){
    d += -0.5 * (line[i + 1][0] - line[i][0]) * (line[i + 1][1] + line[i][1])
  }

  return d
}
isClockwise([[0,0],[1,1],[1,0],[0,0]]) // -0.5 順時針
isClockwise([[0,0],[1,0],[1,1],[0,0]]) //  0.5 逆時針

2.1三、判斷線段路徑是順時針仍是逆時針

/* 
* 判斷線段路徑是順時針仍是逆時針
* 使用端點判斷
* 爲正時，p1-p2-p3   路徑的走向爲逆時針，  
* 爲負時，p1-p2-p3   走向爲順時針，  
* 爲零時，p1-p2-p3   所走的方向不變，亦即三點在一直線上
* 返回值 d < 0  '順時針' : d > 0  '逆時針'
*/
function isClockwise2(line){
  let length = line.length
  let d = 0

  if(length.length < 3)return d;

  //若是不是閉合線，則改成閉合線
  if(!equals(line[0], line[length - 1])){
    length = (line = line.slice()).push(line[0]);
  }

  //循環遍歷多邊形的座標選取X或者Y值中最大或者最小的點，這個點必然是凸點
  //而後取該點先後各一個點Pm-一、Pm+1，組成向量(Pm-1,Pm)、(Pm,Pm+1)。而後進行向量叉乘便可判斷出順時針或逆時針
  let max = line[0][0], maxIndex = 0;
  for(let i = 1; i < length - 1; i++){
    if(line[i][0] > max){
      maxIndex = i
      max = line[i][0]
    }
  }

  //找到端點 p2 ，而後取改點先後各一點 p1 p3
  let p1 = line[maxIndex ? (maxIndex - 1) : (length - 2)]
  let p2 = line[maxIndex]
  let p3 = line[maxIndex + 1]
  //而後根據三個點組成的兩個向量乘積判斷，進行向量叉乘便可判斷出順時針或逆時針
  //d = (x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1)
  d = (p2[0] - p1[0]) * (p3[1] - p2[1]) - (p3[0] - p2[0]) * (p2[1] - p1[1])

  return d
}
isClockwise2([[0,0],[1,1],[1,0],[0,0]]) // -1 順時針
isClockwise2([[0,0],[1,0],[1,1],[0,0]]) //  1 逆時針

2.1四、簡化線 與 平滑線

/* 
* 簡化線 與 平滑線
* 簡化線 與 平滑線 比較複雜，這裏就不展現代碼了
* 簡化線 又能夠叫作線的抽稀，通常使用的方法爲 道格拉斯-普克，這也是我使用的方法
* 道格拉斯-普克算法(Douglas–Peucker algorithm，亦稱爲拉默-道格拉斯-普克算法、迭代適應點算法、分裂與合併算法)是將曲線近似表示爲一系列點，並減小點的數量的一種算法。其思想主要是保留關鍵點。

* Douglas-Peucker算法描述：
* 一、在線首尾兩點A，B之間鏈接一條直線AB，該直線爲曲線的弦；
* 二、獲得曲線上離該直線段距離最大的點C，計算其與AB的距離d；
* 三、比較該距離與預先給定的閾值tolerance的大小，若是小於tolerance，則該直線段做爲曲線的近似，該段曲線處理完畢。
* 四、若是距離大於閾值tolerance，則用C將曲線分爲兩段AC和BC，並分別對兩段曲線進行1~3的處理。
* 五、當全部曲線都處理完畢時，依次鏈接各個分割點造成的折線，便可以做爲曲線的近似。

* 除了可使用道格拉斯-普克算法對線進行抽稀外，還可使用Wang-Müller算法(保留關鍵折彎)、Zhou-Jones算法(保留加權有效面積)、Visvalingam-Whyatt算法(保留有效面積)等。
* 在實際應用中發現，現實是殘酷的，若是閾值太小，簡化後的控制點會不少，用戶很難進行操做。而閾值變大，控制點變少了，可是通過平滑以後，線面變形嚴重，比較失真。真實和美觀之間很難平衡，難以兼得。


* 平滑線 通常使用的方法爲 貝賽爾曲線插值和B樣條曲線插值
* 在應用中發現，二者在某些時候都會有些問題，因此本身對貝賽爾曲線算法作了一點改變，具體的能夠看https://segmentfault.com/a/1190000031626358
*/

三、多邊形

3.一、獲取多邊形邊界範圍

// 獲取多邊形邊界範圍
function bbox(coords) {
  // x/經度最小值 y/緯度最小值 x/經度最大值 y/緯度最大值
  let res = [Infinity, Infinity, -Infinity, -Infinity];
  coords.forEach(coord => {
    if (res[0] > coord[0]) res[0] = coord[0];
    if (res[2] < coord[0]) res[2] = coord[0];
    if (res[1] > coord[1]) res[1] = coord[1];
    if (res[3] < coord[1]) res[3] = coord[1];
  })
  return res;
}
bbox([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [1, 0, 4, 3]

3.二、判斷兩個邊界範圍的關係： a 是否包含 b

//判斷兩個邊界範圍的關係： a 是否包含 b
function isContains(a, b) {
  return a[0] <= b[0] &&
    a[1] <= b[1] &&
    b[2] <= a[2] &&
    b[3] <= a[3];
}
isContains([1, 0, 4, 3], [2, 2, 5, 5])  //false
isContains([1, 0, 4, 3], [2, 1, 3, 2])  //true

3.三、判斷兩個邊界範圍的關係： a 與 b 是否有交集

//判斷兩個邊界範圍的關係： a 與 b 是否有交集
function isIntersects(a, b) {
  return b[0] <= a[2] &&
    b[1] <= a[3] &&
    b[2] >= a[0] &&
    b[3] >= a[1];
}
isIntersects([1, 0, 4, 3], [2, 2, 5, 5])  //true
isIntersects([1, 0, 4, 3], [5, 2, 5, 5])  //false

3.四、獲取多邊形中心

// 獲取多邊形中心
function center(coords) {
  let ext = bbox(coords);
  let x = (ext[0] + ext[2]) / 2;
  let y = (ext[1] + ext[3]) / 2;
  return [x, y];
}
center([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [2.5, 1.5]

3.五、獲取多邊形重心

// 獲取多邊形重心
function centroid(coords) {
  let xSum = 0, ySum = 0, len = 0;
  coords.forEach(coord => {
    xSum += coord[0];
    ySum += coord[1];
    len++;
  })

  return [xSum / len, ySum / len];
}
centroid([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [2.4, 1.8]

3.六、獲取多邊形質心

// 獲取多邊形質心
function centerOfMass(coords) {
  let centr = centroid(coords);
  let neutralizedPoints = coords.map(function (point$$1) {
    return [
      point$$1[0] - centr[0],
      point$$1[1] - centr[1]
    ];
  });
  
  let sx = 0, sy = 0, sArea = 0;
  let x1, x2, y1, y2, a;
  for (let i = 0, len = coords.length - 1; i < len; i++) {
    x1 = neutralizedPoints[i][0]; 
    y1 = neutralizedPoints[i][1];
    x2 = neutralizedPoints[i + 1][0]; 
    y2 = neutralizedPoints[i + 1][1];

    // a 是計算有符號面積和最終座標的公因子
    a = x1 * y2 - x2 * y1;
    // sArea 用來計算有符號面積的和
    sArea += a;
    // sx 和 sy 是用於計算最終座標的總和
    sx += (x1 + x2) * a;
    sy += (y1 + y2) * a;
  }

  if (sArea === 0) {
    return centre;
  } else {
    // 計算有符號面積，並因式分解: x = 1 / 6A
    let area = sArea * 0.5;
    let areaFactor = 1 / (6 * area);

    return [
      centr[0] + areaFactor * sx,
      centr[1] + areaFactor * sy
    ]
  }
}
centerOfMass([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [2.569230769230769, 1.8735042735042735]

3.七、獲取多邊形面積

/* 
* 獲取多邊形面積
* 原理：△ABC的面積就是'向量AB'和'向量AC'兩個向量叉積的絕對值的一半。其正負表示三角形頂點是在右手系仍是左手系。
* 因此能夠把多邊形以多邊形內一點爲頂點，拆分出一個個三角形，再求取面積
* 設 n 邊形的頂點依次是 (x1，y1)(x2,y2)......(xn,yn)
* 那麼：s = (x1y2-x2y1)/2 + (x2y3-x3y2)/2 +......+ (xny1-x1yn)/2
* 獲取的 s 是有向面積，須要取絕對值
*/
function getArea(polygon){
  let area = 0;
  let len = polygon.length
  let x1 = polygon[len - 1][0]
  let y1 = polygon[len - 1][1]
  for(let i = 0; i < len; i++){
    let x2 = polygon[i][0];
    let y2 = polygon[i][1];
    area += y1 * x2 - x1 * y2;
    x1 = x2;
    y1 = y2
  }
  return Math.abs(area / 2)
}
getArea([[0,0],[0,10],[10,10],[10,0]])  //100
getArea([[0,0],[0,3],[4,0]])  //6

3.八、計算多邊形投影到地球上時的近似面積（平方米）

/* 
* 計算多邊形投影到地球上時的近似面積（平方米）
* 注意，若是環是順時針方向的，這個區域將是正的，不然它將是負的
* 注意，投影及其長半軸參數radius的介紹這裏不作敘述，能夠看上面haversineDistance的相關介紹
*/
function geodesicArea(coords) {
  let radius = 6378137
  let area = 0, len = coords.length;
  if(len <= 2){ return area }

  let x1 = coords[len - 1][0];
  let y1 = coords[len - 1][1];
  for (let i = 0; i < len; i++) {
    let x2 = coords[i][0], y2 = coords[i][1];
    area += toRadians(x2 - x1) * (2 + Math.sin(toRadians(y1)) + Math.sin(toRadians(y2)));
    x1 = x2;
    y1 = y2
  }
  area = area * radius * radius / 2
  return Math.abs(area)
}
geodesicArea([[118, 39], [117, 34], [117, 33], [116, 36], [117, 40]])  //64519860945.307144平方米

3.九、判斷點是否在多邊形內

/* 
* 判斷點是否在多邊形內
* 射線法（ray casting）或者奇偶規則法（even odd rule）
* 從這個點作一條射線，計算它跟多邊形邊界的交點個數，若是交點個數爲奇數，那麼點在多邊形內部，不然點在多邊形外。
* 對於射線法，須要排除如下幾種狀況：
* 一、點在多邊形的邊上
* 二、點和多邊形的頂點重合
* 三、射線通過多邊形頂點
* 四、射線恰好通過多邊形的一條邊
* truf對應實現truf.booleanPointInPolygon://turfjs.org/docs/#booleanPointInPolygon
*/
function pointInPolygon(point, polygon) {
  let px = point[0], py = point[1], flag = false

  let sx, sy, tx, ty
  for(let i = 0, l = polygon.length, j = l - 1; i < l; j = i, i++) {
    sx = polygon[i][0]
    sy = polygon[i][1]
    tx = polygon[j][0]
    ty = polygon[j][1]

    // 點與多邊形頂點重合       
    if((sx === px && sy === py) || (tx === px && ty === py)) {
      return true
    }

    // 判斷線段兩端點是否在射線兩側       
    if((sy < py && ty >= py) || (sy >= py && ty < py)) {
      // 線段上與射線 Y 座標相同的點的 X 座標         
      let x = sx + (py - sy) * (tx - sx) / (ty - sy)

      // 點在多邊形的邊上         
      if(x === px) {
        return true
      }

      // 射線穿過多邊形的邊界         
      if(x > px) {
        flag = !flag
      }
    }
  }

  // 射線穿過多邊形邊界的次數爲奇數時點在多邊形內     
  return flag
}

3.十、判斷點是否在多邊形內

/* 
* 判斷點是否在多邊形內
* 迴轉數法
* 用線段分別鏈接點和多邊形的所有頂點，計算全部點與相鄰頂點連線的夾角，計算全部夾角和，最後根據角度累加值計算迴轉數
* 注意每一個夾角都是有方向的，因此有多是負值。360°（2π）至關於一次迴轉。
* 當迴轉數爲 0 時，點在閉合曲線外部。
*/
function pointInPolygon2(point, polygon) {
  let px = point[0], py = point[1], sum = 0

  let sx, sy, tx, ty
  for(let i = 0, l = polygon.length, j = l - 1; i < l; j = i, i++) {
    sx = polygon[i][0]
    sy = polygon[i][1]
    tx = polygon[j][0]
    ty = polygon[j][1]

    // 點與多邊形頂點重合       
    if((sx === px && sy === py) || (tx === px && ty === py)) {
      return true
    }

    // 點與相鄰頂點連線的夾角       
    let angle = Math.atan2(sy - py, sx - px) - Math.atan2(ty - py, tx - px)

    // 確保夾角不超出取值範圍（-π 到 π）
    if(angle >= Math.PI) {
      angle = angle - Math.PI * 2
    } else if(angle <= -Math.PI) {
      angle = angle + Math.PI * 2
    }

    sum += angle
  }

  // 計算迴轉數並判斷點和多邊形的幾何關係     
  return Math.round(sum / Math.PI) !== 0
}

3.十一、判斷點是否在多邊形內

/* 
* 判斷點是否在多邊形內
* openlayers中的方法
*/
function pointInPolygon3(point, polygon) {
  let x = point[0], y = point[1]
  let wn = 0;
  let x1 = polygon[0][0];
  let y1 = polygon[0][1];
  for (let i = 0; i < polygon.length; i++) {
    let x2 = polygon[i][0];
    let y2 = polygon[i][1];
    if (y1 <= y) {
      if (y2 > y && (x2 - x1) * (y - y1) - (x - x1) * (y2 - y1) > 0){
        wn++;
      }
    } else if (y2 <= y && (x2 - x1) * (y - y1) - (x - x1) * (y2 - y1) < 0){
      wn--;
    }
    x1 = x2;
    y1 = y2
  }
  return wn !== 0
}

3.十二、兩個多邊形之間的關係：相交、包含、相離

/* 
* 兩個多邊形之間的關係：相交、包含、相離
* https://segmentfault.com/a/1190000020916225
* 一、遍歷A多邊形上的點，判斷是否有座標點在B多邊形內 --- 返回結果 a
* 二、遍歷B多邊形上的點，判斷是否有座標點在A多邊形內 --- 返回結果 b
* 若是a、b都爲true，則兩個多邊形相交
* 若是a爲true，b爲false，則多邊形B包含多邊形A
* 若是a爲false，b爲true，則多邊形A包含多邊形B
* 若是a、b都爲false，則兩個多邊形遠離
* truf對應實現 包含 truf.booleanContains://turfjs.org/docs/#booleanContains
* truf對應實現 交叉 truf.booleanCrosses://turfjs.org/docs/#booleanCrosses
* truf對應實現 相離 truf.booleanDisjoint://turfjs.org/docs/#booleanDisjoint
*/
function judge(coordsA, coordsB){
  let boola = coordsA.some(item => {
    return pointInPolygon3(item, coordsB)
  }) ? 1 : 0;
  let boolb = coordsB.some(item => {
    return pointInPolygon3(item, coordsA)
  }) ? 1 : 0;

  return ['相離', 'A包含B', 'B包含A', '相交'][boola * 2 + boolb]
}

3.1三、線切割多邊形(閉合線)

/* 
* 線切割多邊形(閉合線)
* 使用線將多邊形切割爲一個個小多邊形
* truf中相似的有：
* 差別 經過從第一個多邊形中剪切第二個多邊形來找到兩個多邊形之間的差別。
* truf.difference://turfjs.org/docs/#difference
* 相交 取兩個多邊形，並找到它們的交點。若是它們相交，則返回相交邊界。若是它們不相交，則返回undefined。
* truf.intersect://turfjs.org/docs/#intersect
*/
function splitPolygonByLine(polygon, line){
  let polygonItem = { coords: polygon }
  _splitPolygonByLine(polygonItem, line)
  return getChildrenPolygons([polygonItem])
}
function _splitPolygonByLine(polygonItem, lineCoords){
  if(polygonItem.children && polygonItem.children.length){
    return polygonItem.children.forEach(polygon => {
      _splitPolygonByLine(polygon, lineCoords)
    });
  }

  lineCoords = lineCoords.slice();
  let points = lineIntersect(lineCoords, polygonItem.coords, 2)
  //console.log('points', points)
  if(points.length >= 2){
    //判斷線的開始點在面內仍是面外，以判斷哪兩個點練成的線在面外，捨棄在外面的
    let startOut = pointInPolygon3(lineCoords[0], polygonItem.coords)
    let point1Out = pointInPolygon3(lineCoords[points[0].index1], polygonItem.coords)
    if(startOut && !point1Out){
      lineCoords = lineCoords.slice(points[0].index1);
      return _splitPolygonByLine(polygonItem, lineCoords)
    }

    let newPolygons = splitPolygon(polygonItem, lineCoords, points);
    //console.log('newPolygons', newPolygons)
    if(newPolygons.length >= 2){
      polygonItem.children = newPolygons

      lineCoords = lineCoords.slice(points[1].index1)
      newPolygons.forEach(polygon => {
        _splitPolygonByLine(polygon, lineCoords)
      })
    }
  }
}

function splitPolygon(polygonItem, lineCoords, points){
  let result = [];

  //lineCoords = lineCoords.slice()
  let polygonCoords = polygonItem.coords
  
  let startIndex = 0, endIndex = 0, lineIndex = 0

  points.sort((a,b)=>a.index1 - b.index1)
  
  points.reduce((prePoint, point, index) => {
    startIndex = prePoint.index1
    endIndex = point.index1

    let line;
    if(endIndex == startIndex){
      //return point
      line = [prePoint.coords, point.coords]
    }else{
      line = lineCoords.slice(startIndex, endIndex)
      line.unshift(prePoint.coords)
      line.push(point.coords)
    };

    lineIndex = endIndex

    let p1, p2 = polygonCoords.slice();
    if(prePoint.index2 > point.index2){
      startIndex = point.index2
      endIndex = prePoint.index2

      p1 = polygonCoords.slice(startIndex, endIndex)
      p1 = p1.concat(line)
      p1.push(p1[0])

      line.reverse()
      p2 = p2.slice(0, startIndex).concat(line).concat( p2.slice(endIndex) )
    }else if(prePoint.index2 < point.index2){
      startIndex = prePoint.index2
      endIndex = point.index2

      p2 = p2.slice(0, startIndex).concat(line).concat( p2.slice(endIndex) )

      line.reverse()
      p1 = polygonCoords.slice(startIndex, endIndex)
      p1 = p1.concat(line)
      p1.push(p1[0])
    }else{
      //startIndex = endIndex = prePoint.index1
      p1 = line.slice()

      let pc = polygonItem.coords[point.index2 - 1]
      if(dist2d(prePoint.coords, pc) > dist2d(point.coords, pc)){
        line.reverse()
      }

      p2.splice(point.index2, 0, ...line)
    }

    result.push({ coords: p1 }, { coords: p2 })

    return point
  })

  return result
}
function getChildrenPolygons(polygons){
  let result = []
  polygons.forEach(item => {
    if(item.children){
      result = result.concat( getChildrenPolygons(item.children) )
    }else if(getArea(item.coords) > 0){
      result.push(item.coords)
    }
  });
  
  return result
}
splitPolygonByLine([[0, 0], [0, 10], [10, 10], [10, 0], [0, 0]], [[-1, 0], [6, 12], [12, 2], [-1, 9]])
/* 
[
  [[0,1.7142857142857142],[0,8.461538461538462],[2.995121951219513,6.848780487804878],[0,1.7142857142857142]],
  [[0,10],[4.833333333333334,10],[2.995121951219513,6.848780487804878],[0,8.461538461538462],[0,10]],
  [[10,10],[10,5.333333333333334],[7.2,10],[10,10]],
  [[4.833333333333334,10],[7.2,10],[10,5.333333333333334],[10,3.076923076923077],[2.9951219512195113,6.848780487804879],[4.833333333333334,10]],
  [[0,0],[0,1.7142857142857142],[2.9951219512195113,6.848780487804879],[10,3.076923076923077],[10,0],[0,0]]
]
*/

github地址：https://github.com/hfhan/math