寫在前面

嗨！很高興看到你點(diǎn)進(jìn)來閱讀這篇文章，請別介意，標(biāo)題有點(diǎn)長有點(diǎn)啰嗦（完全是為了seo考慮），但也算是概括了這篇文章的內(nèi)容。如果你是要開發(fā)如下圖所示的場景，但又苦于沒什么好的思路，那么這篇文章一定會(huì)幫助到你！

往復(fù)式運(yùn)動(dòng)航線

基于不規(guī)則凸多邊形地塊的往復(fù)式航線規(guī)劃

哦，對了，本文的實(shí)現(xiàn)是基于web平臺(tái)的地圖，使用javascript。如果你也是在web平臺(tái)上開發(fā)，而且任務(wù)時(shí)間非常緊急，沒有時(shí)間閱讀完全文的話。。。我已經(jīng)將本文的思路封裝成一個(gè)庫了，你可以猛戳下面的鏈接，開箱即用：
https://github.com/Char-Ten/cpRPA
兼容各大地圖平臺(tái)api（其實(shí)不同平臺(tái)api差異的影響很低）哦，不信的話戳demo：
百度地圖demo
高德地圖demo
leaflet地圖demo
覺得好用的話記得給個(gè)star～原創(chuàng)不易，謝謝支持

正文！

其實(shí)也是套公式

其實(shí)這種問題，實(shí)際上是數(shù)學(xué)幾何應(yīng)用題，既然是數(shù)學(xué)題啦，那按照考試的套路第一步肯定是套公式啊，這種場景，核心的公式不多，就兩條：

• 一次函數(shù)兩點(diǎn)表達(dá)式
• 繞（tx,ty）點(diǎn)旋轉(zhuǎn)n度之后縮放SxSy倍的變換矩陣

第一條沒什么可以說的，初二數(shù)學(xué)就開始教一次函數(shù)的知識(shí)，這一條是用來計(jì)算航線與地塊邊界的交點(diǎn)的。
第二條就是很經(jīng)典的復(fù)合變換矩陣了，分別是位移矩陣叉乘旋轉(zhuǎn)矩陣叉乘位移縮放矩陣，我們設(shè)其叉乘結(jié)果為A，那么我們就可以列出下面的等式：

計(jì)算過程就是。。。橫乘豎橫乘豎橫乘豎橫乘豎橫乘豎橫乘豎橫乘豎橫乘豎。。。。最后化為用于程序的代數(shù)式就是：

  function transform(x, y, tx, ty, deg, sx, sy) {
        var deg = deg * Math.PI / 180;
        if (!sy) sy = 1;
        if (!sx) sx = 1;
        return [
            sx * ((x - tx) * Math.cos(deg) - (y - ty) * Math.sin(deg)) + tx,
            sy * ((x - tx) * Math.sin(deg) + (y - ty) * Math.cos(deg)) + ty
        ]
  }

  function calcPointInLineWithY(p1,p2,y){
        var s = p1[1] - p2[1];
        var x;
        if (s) {
            x = (y - p1[1]) * (p1[0] - p2[0]) / s + p1[0]
        } else {
            return false
        }
        
        /**判斷x是否在p1,p2在x軸的投影里，不是的話返回false*/
        if (x > p1[0] && x > p2[0]) {
            return false
        }
        if (x < p1[0] && x < p2[0]) {
            return false
        }
        return [x, y]
  }

看到這里，恭喜你，你已經(jīng)完成了50%的工作量！如果是在考試，你把這兩條公式列出來，不寫答案也有一半的分?jǐn)?shù)（

先從最簡單的場景開始

一個(gè)矩形地塊，航線水平于x軸:

1. 兩條橫線的間隔是20
2. 每一條橫線都可以表示為 y=N，N為常數(shù)，表示某個(gè)緯度值
3. 每一條橫線段都是y=N與矩形相交產(chǎn)生，也就是每一條橫線段都是該矩形地塊與維度相交的結(jié)果

那么，現(xiàn)在矩形的四個(gè)頂點(diǎn)的經(jīng)緯度是已知的，無人機(jī)飛行的間隔也是已知的，這個(gè)矩形需要與多少條緯度線相交是未知的，每一條橫線的N是未知的，每一條橫線段左右兩個(gè)點(diǎn)的緯度是未知的。根據(jù)已知求未知，你的目標(biāo)已經(jīng)很明確了，一道很簡單的幾何題:

• 該矩形需要與多少條緯度線相交：

/**@method 獲取兩個(gè)經(jīng)緯度點(diǎn)的距離
 * @param {Object} p1,p2 - 兩個(gè)經(jīng)緯度點(diǎn)
 * @return {Number} 單位 米
*/
function distance(p1,p2){
  /**leaflet提供的方法*/
  return L.latLng(p1.lat,p1.lng).distanceTo(L.latLng(p2.lat,p2.lng))
}

/** nw到sw的距離*/
var dist = distance(nw,sw);

/** 得出答案*/
var lines = parseInt(dist / 20)

• 求每一條橫線的N:

var N=[];
var stepLat=(nw.lat-sw.lat)/lines;
for(var i=0;i<lines;i++){
  N.push(nw.lat - i * stepLat)
}

• 因?yàn)榫匦蔚膬蓷l邊是垂直的，所以，橫線段左右兩個(gè)點(diǎn)的經(jīng)度分別為nw.lng，ne.lng。這樣我們就可以繪制出來了：

for(var i=0;i<N.length;i++){
  L.polyline([
    {
      lat:N[i],
      lng:nw.lng
    },{
      lat:N[i],
      lng:ne.lng
    }]).addTo(map)
}

場景開始變形！

鏘鏘，我們把矩形上面的邊往東挪50米，得到一個(gè)平行四邊形：

聰明的你一定發(fā)現(xiàn)了，平行四邊形在Y軸上的投影根本沒有發(fā)生變化嘛，即使變了之后，穿過地塊的緯度線數(shù)目還是不變嘛，只不過，這次因?yàn)閮蓷l邊不是垂直的，所以，我們需要計(jì)算斜邊與緯度線的交點(diǎn)。等等，你這時(shí)候想起了，最開始50%工作量里面所封裝的那個(gè)calcPointInLineWithY函數(shù)！
你已經(jīng)知道斜邊兩個(gè)點(diǎn)的坐標(biāo)，然后你又知道y=N，那你通過一次函數(shù)的兩點(diǎn)表達(dá)式，完全就可以知道x，也就是經(jīng)度是多少啦：

/**這里注意 緯度lat對應(yīng)y，經(jīng)度lng對應(yīng)x*/
for(var i=0;i<N.length;i++){
  var westPoint=calcPointInLineWithY(
    [nw.lng,nw.lat],
    [sw.lng,sw.lat],
    N[i]
  );
  var eastPoint=calcPointInLineWithY(
    [ne.lng,ne.lat],
    [se.lng,se.lat],
    N[i]
  );
  if(eastPoint&&westPoint){
    L.polyline([westPoint,eastPoint]).addTo(map)
  }
}

那你可以再變一變，讓y軸上的投影也發(fā)生變化，就像這樣：

好了，這下你觀察到，每條邊都跟緯度線相交了，也就是說，這次你要遍歷一下這個(gè)平行四邊形四個(gè)頂點(diǎn)。等等，你似乎忘記了一個(gè)問題，這個(gè)四邊形在y軸上的投影發(fā)生了變化，相交緯度線數(shù)目也跟著發(fā)生變化了。這時(shí)候你想到，要不給這個(gè)多邊形做個(gè)外接矩形？就像這樣：

這樣是不是又回歸了最開始的場景？只是把calcPointInLineWithY函數(shù)加上去之后，你可以得到任意凸多邊形與緯度線相交的模型。

• 首先是創(chuàng)建多邊形外接的矩形，將多邊形頂點(diǎn)全部遍歷一遍，取到最大和最小的經(jīng)緯度值。經(jīng)度最大為東，經(jīng)度最小為西。緯度最大為北，緯度最小為南。將這幾個(gè)組合一下，獲得西北點(diǎn)，東北點(diǎn)，東南點(diǎn)和西南點(diǎn)，就可以弄出一個(gè)矩形出來

  /**@method 創(chuàng)建多邊形外接矩形
 * @param {Array} latlngs - 格式為[{lat,lng}]的經(jīng)緯度數(shù)組，多邊形的頂點(diǎn)集合
 * @return {Array} .rect - 返回外接矩形四頂點(diǎn)
 * @return {Object} .center - 返回外接矩形的中心的 
*/
function createPolygonBounds(latlngs){
  var lats=[];
  var lngs=[];
  for(var i=0;i<latlngs.length;i++){
    lats.push(latlngs[i].lat);
    lngs.push(latlngs[i].lng);
  }
  var maxLat=Math.max.apply(Math,lats);
  var maxLng=Math.max.apply(Math,lngs);
  var minLat=Math.min.apply(Math,lats);
  var minLng=Math.min.apply(Math,lngs);
  return {
    center:{
      lat:(maxLat+minLat)/2,
      lng:(maxLng+minLng)/2
    },
    latlngs:[{
      lat:maxLat,
      lng:minLng//西北
    },{
      lat:maxLat,
      lng:maxLng//東北
    },{
      lat:minLat,
      lng:maxLng//東南
    },{
      lat:minLat,
      lng:minLng//西南
    }]
  }
}

• 然后是計(jì)算這個(gè)外接矩形穿過了多少條緯度線，跟之前那個(gè)場景是一樣的。rect是上面那個(gè)函數(shù)創(chuàng)建的數(shù)組，可以看到西北點(diǎn)的索引是0，西南點(diǎn)的索引是3，所以計(jì)算西北到西南的點(diǎn)，也就是這個(gè)外接矩形的高度。這個(gè)方法返回有l(wèi)en條緯度線穿過，而且穿過的緯度相差lat。

/**@method 計(jì)算有多少條緯度線穿過
 * @param {Array} rect - 外接矩形
 * @param {Number} space - 間隔
*/
function calcLatsInPolygon(rect,space){
  var lines=parseInt(distance(rect[0],rect[3])/space/2);
  var lat=(rect[0].lat-rect[3].lat)/lines;
  return {
    len:lines,
    lat:lat
  }
}

• 最后就是結(jié)合上面幾個(gè)場景，改寫一下最終的生成函數(shù)，這里你是直接將生成的橫線畫了出來，如果需要做折線連接順序處理，還需要聲明一個(gè)數(shù)組去存儲(chǔ)生成的點(diǎn)，比如當(dāng)緯度線的索引是奇數(shù)時(shí)，橫線從西往東畫，也就是先放西邊的點(diǎn)再放東邊的點(diǎn)；如果緯度線索引是偶數(shù)時(shí)，橫線從東往西畫，也就是反過來。這里就留給你自己處理了

/**任意多邊形頂點(diǎn)集*/
var polygon=[{lat:23.13184566463993,lng:113.25901150703432},/**...其他點(diǎn)*/];

/**飛行間隔*/
var space=10;

/**外接矩形*/
var outRect=createPolygonBounds(polygon);

/**緯度線*/
var latLines=calcLatsInPolygon(outRect.latlngs,space);

var line=[];
/**遍歷每一條緯度線*/
for(var i=0;i<latLines.len;i++){
  line=[];
  /**遍歷每一個(gè)多邊形頂點(diǎn)*/
  for(var j=0;j<polygon.length;j++){
    var point=calcPointInLineWithY([
       polygon[i].lng,
       polygon[i].lat,
    ],[
      polygon[si(i+1,polygon.length)].lng
      polygon[si(i+1,polygon.length)].lat
    ],outRect.latlngs[0].lat - i*latLines.lat)
    if(point){
      line.push(point)
    }
  }
  
  /**去掉只有一個(gè)交點(diǎn)的緯度線*/
  if(line.length<2){
    continue
  }
  
  /**去掉兩個(gè)交點(diǎn)重合的緯度線*/
  if(line[0][0] === line[1][0]){
    continue
  }
  
  /**leaflet 繪制*/
  L.polyline([
    {lat:line[0][1],lng:line[0][0]},
    {lat:line[1][1],lng:line[1][0]}
  ]).addTo(map)
}

/**防止索引溢出*/
function si(i,l){
  if (i > l - 1) {
    return i - l;
  }
  if (i < 0) {
    return l + i;
  }
  return i;
}

到這里，基本上已經(jīng)完成了90%，剩下10%那部分最簡單了。簡單么？你歪著頭問，到現(xiàn)在為止只是多邊形與緯度線相交啊，現(xiàn)實(shí)中無人機(jī)又不可能都是沿著緯度在地圖上橫著飛！它可以在地圖上沿任意方向飛行的，可上面的做法，只能讓無人機(jī)橫著走??！
別急，你這時(shí)候要淡定，你要想想開頭不是還留了一個(gè)transform函數(shù)么？
這個(gè)transform的作用是讓坐標(biāo)(x,y)繞著(tx,ty)旋轉(zhuǎn)deg度后在縮放SySx倍得到一個(gè)的新坐標(biāo)(_x,_y)。沒錯(cuò)，這次我們要拿這些坐標(biāo)進(jìn)行旋轉(zhuǎn)運(yùn)動(dòng)。
而且牛頓告訴過你：

運(yùn)動(dòng)都是相對的

先放一張gif圖，看看如何繪制一個(gè)斜45度角的航線：

GIF.gif

先讓多邊形繞著中心點(diǎn)旋轉(zhuǎn)想要的角度，將得到的新多邊形再與緯度線做相交操作，獲取到那些交點(diǎn)之后，再將那些交點(diǎn)旋轉(zhuǎn)回來。換句話說，變換前它是一個(gè)任意多邊形，變換后，它還是一個(gè)任意多邊形，都是滿足上面已經(jīng)預(yù)設(shè)好的場景的。這樣的好處顯而易見，你不需要修改上面的任何一個(gè)函數(shù)，也不需要去多寫一條兩個(gè)一次函數(shù)求交點(diǎn)的公式。把問題化到最簡單的場景去，只需要添加變換的代碼：

/**@method 創(chuàng)建一個(gè)旋轉(zhuǎn)后的多邊形
 * @param {Array}  latlngs - 經(jīng)緯度頂點(diǎn)集
 * @param {Object} bounds - 由上文  createPolygonBounds 函數(shù)創(chuàng)建的對象
 * @param {Number} rotate - 旋轉(zhuǎn)角度
 * @return {Array} - 變換后的經(jīng)緯度數(shù)組
 */
function createRotatePolygon(latlngs, bounds, rotate){
  var res=[];
  for(var i=0;i<latlngs.length;i++){
    var tr=transform(
      latlngs[i].lng,
      latlngs[i].lat,
      bounds.center.lng,
      bounds.center.lat,
      rotate
    );
    res.push({lat:tr[1],lng:tr[0]});
  }
  return res
}

這里你一定要牢記，lng是經(jīng)度，對應(yīng)x；lat是緯度，對應(yīng)y（被leaflet框架坑哭的我QuQ。。。
然后，將原來的代碼加上去

var line=[];

/**折線*/
var polyline=[];

/**旋轉(zhuǎn)45度*/
var rotate=45;

/**創(chuàng)建變換后的多邊形*/
var rPolygon=createRotatePolygon(polygon,outRect,-rotate)
/**遍歷每一條緯度線*/
for(var i=0;i<latLines.len;i++){
  line=[];
  /**遍歷每一個(gè)多邊形頂點(diǎn)*/
  for(var j=0;j<rPolygon.length;j++){
    var point=calcPointInLineWithY([
      rPolygon[i].lng,
      rPolygon[i].lat,
    ],[
      rPolygon[si(i+1,rPolygon.length)].lng
      rPolygon[si(i+1,rPolygon.length)].lat
    ],outRect.latlngs[0].lat - i*latLines.lat)
    if(point){
      line.push(point)
    }
  }
  
  /**去掉只有一個(gè)交點(diǎn)的緯度線*/
  if(line.length<2){
    continue
  }
  
  /**去掉兩個(gè)交點(diǎn)重合的緯度線*/
  if(line[0][0] === line[1][0]){
    continue
  }
  
  /**這里不繪制了，先塞進(jìn)polyline這個(gè)數(shù)組*/
  polyline.push(
    {lat:line[0][1],lng:line[0][0]},
    {lat:line[1][1],lng:line[1][0]}
  )
}

/**然后就可以直接轉(zhuǎn)回來繪制*/
L.polyline(
  createRotatePolygon(polyline,outRect,rotate)
).addTo(map)

小瑕疵

也許細(xì)心的你發(fā)現(xiàn)了，gif圖里面轉(zhuǎn)是轉(zhuǎn)了，斜45度的角也是畫出來了，但是旋轉(zhuǎn)后的圖形好像是被拉長了！旋轉(zhuǎn)回來時(shí)又會(huì)被壓肥回來。這個(gè)嘛。。。
這個(gè)問題其實(shí)我在寫驗(yàn)證的時(shí)候也發(fā)現(xiàn)了，粗略計(jì)算多邊形面積和航線掃過的面積的比值，總是發(fā)現(xiàn)0度的比值最接近于1，90度的比值最小，不管多邊形是什么形狀。
最開始我一直以為是面積算法的原因，直到寫這篇文章的時(shí)候去寫那個(gè)gif動(dòng)畫后才發(fā)現(xiàn)，多邊形變換后變形了，轉(zhuǎn)換后與緯度相交的數(shù)量變多了。而我猜想變形的原因可能是，地球并不是一個(gè)平面，它是彎的你知道吧，這些經(jīng)緯度雖然是投影到了平面上了，但實(shí)際上它們是在一個(gè)球上。直接拿經(jīng)緯度變換相當(dāng)于先在球上做了旋轉(zhuǎn)，然后在投影到地圖的平面上，這樣看起來就像是被拉長了。
所以，你不能直接變換經(jīng)緯度，而是要將經(jīng)緯度換算成地圖上的像素坐標(biāo)，變換完之后再轉(zhuǎn)回來，這樣圖像就不會(huì)被拉長了。
因此先來兩個(gè)像素系與經(jīng)緯度坐標(biāo)系轉(zhuǎn)換的方法壓壓驚：

/**@method 設(shè)置經(jīng)緯度轉(zhuǎn)換成頁面像素坐標(biāo)的方法
 * @param {Object} latlng - {lat,lng}
*/
function latlng2px(latlng){
    /**百度，map為 new BMap.Map() 對象*/
    return map.pointToPixel(new BMap.Point(latlng.lng, latlng.lat))

    /**
    * 高德，map為 new AMap.Map() 對象
    * return map.lngLatToContainer(new AMap.LngLat(latlng.lng, latlng.lat))
    *
    * leaflet map 為 L.map對象
    * return map.latLngToLayerPoint(L.latLng(latlng.lat, latlng.lng)) 
    */
}

/**@method 設(shè)置像素坐標(biāo)轉(zhuǎn)換成經(jīng)緯度點(diǎn)的方法
 * @param {Array} px - [lng,lat] 
*/
function px2latlng(px){
    /**百度，map為 new BMap.Map() 對象*/
    return map.pixelToPoint(new BMap.Pixel(px[0], px[1]))

    /**
    * 高德，map為 new AMap.Map() 對象
    * return map.containerToLngLat(new AMap.Pixel(px[0], px[1]))
    * 
    * leaflet map 為 L.map對象
    * return map.layerPointToLatLng(L.point(px[0], px[1]))
    */
}

然后將原來的createRotatePolygon函數(shù)改為

function createRotatePolygon(latlngs,bounds,rotate){
  var res = [] , a , b;
  var c = latlng2Px(bounds.center);
  for (var i = 0; i < latlngs.length; i++) {
    a = latlng2Px(latlngs[i]);
    b = transform(a.x, a.y, c.x, c.y, rotate);
    res.push(px2Latlng(b));
  }
  return res;
}

這樣就解決了拉長的問題：

GIF.gif

寫在最后

終于寫完了此文，真是不容易。這個(gè)思路，從最開始思想混亂與Leaflet框架緊密耦合，到一步步解耦，到自己決定寫一個(gè)適用于各個(gè)地圖平臺(tái)的庫，到寫這篇文章，差不多已經(jīng)過去兩個(gè)星期了。我發(fā)現(xiàn)，很多問題是你調(diào)試過程中發(fā)現(xiàn)不了的，等到你調(diào)試好了，決定寫一篇裝逼的文章，在寫的過程中你就會(huì)發(fā)現(xiàn)各種調(diào)試中出現(xiàn)不了的問題，比如那個(gè)變換變形的問題。

在動(dòng)手前，關(guān)于此類的教程文章少之又少，唯一可以找到比較符合場景的竟然是百度文庫里面的一篇論文:
基于作業(yè)航向的不規(guī)則區(qū)域作業(yè)航線規(guī)劃算法研究
論文懂吧，長倒是不長，就是臭，里面堆砌著各種奇奇怪怪的術(shù)語。之后看了兩天后才明白他的實(shí)現(xiàn)思路，大同小異，只不過不知道他怎么搞的，新弄出一個(gè)坐標(biāo)系出來，感覺這樣是增加了思路的復(fù)雜度啊。所以在他的算法的基礎(chǔ)上我做了簡化，不做坐標(biāo)系偏移，取代的是變換地塊坐標(biāo)，這樣實(shí)現(xiàn)起來相對簡單些。

因此在這里，小小貢獻(xiàn)一下這個(gè)思路吧，也讓以后有人像我這樣被坑去寫這種無人機(jī)航線規(guī)劃的，能夠很快地實(shí)現(xiàn)，不用再去看那些奇奇怪怪的論文了。。。

最后最后最后，安利一下https://github.com/Char-Ten/cpRPA 這個(gè)庫（已經(jīng)無恥到這個(gè)地步了。。。），可以接入百度、高德、leaflet，然后算航線！什么？你不需要。。。那你需要計(jì)算地圖上多邊形地塊的面積不？我也提供算面積的方法，百度地圖、高德地圖都沒有這種算面積的方法！好評給個(gè)star！謝謝大家！（夠了喂！(╯‵□′)╯︵┻━┻）
如果你有更好的實(shí)現(xiàn)思路，或者新的使用場景，或者有使用問題，或者有bug，歡迎在下面評論。。。或者直接提issue:https://github.com/Char-Ten/cpRPA/issues