一、球面插值

球面插值的原理大概就如下圖所示，大致理解就是計算球面角度的佔比，計算公式不是太複雜，如下所示：

當角度無限接近於0的時候，這個時候球面插值就演變為線性插值

下面用scala對球面插值進行一個簡單的實現：

class Slerp4scala[T <: Double](start: Vector[T], end: Vector[T], t: Double, omiga: Double) {
  def this() = this(Vector(), Vector(), 0, 0)
  def slerp(): Unit = {
    if (start.size < 1) throw new Exception("Please fill your Array.")
    if (t > 1 || t < 0) throw new Exception("Interpolator point error, its range is in " +
      "[0, 1].")
    val divisor = math.sin(omiga)
    val p1 = math.sin((1 - t) * omiga) / divisor
    val p2 = math.sin(t * omiga) / divisor
    start.map(_ * p1).zip(end.map(_ * p2)).map(k => k._1 + k._2)
  }
}

專案本意是要實現地球座標的球面插值，因為做過資料稀疏化的資料的距離都很大，不能近似的認為是線性插值（誤差較大），但是在求地球的omiga的時候，感覺實現起來還是有點麻煩，於是就找到了第三方的這個庫geotools，之前用的地理圖形操作都是採用的jts的包，這次對geotools又有了一個新的認識。下面是對Geotools的認識。

二、Geotools

高興的是在這個工具集中有很多的工具類，失望的我好像沒有發現計算omiga的介面，但是發現了一個計算路徑中平分點的api，雖然不能有效的解決我的問題，但也算是失之東隅收之桑榆吧，下面將測試程式碼貼出來：

"Geotools " should "be success " in{
    val geoToolsTest = new GeodeticCalculator()
    geoToolsTest.setStartingGeographicPoint(0, 0)
    geoToolsTest.setDestinationGeographicPoint(60, 0)
    assert(geoToolsTest.getAzimuth > 0)
    // 取路徑上的幾個平均分段點
    println(geoToolsTest.getGeodeticPath(5))
    // 計算航向（與正北方向的夾角）
    println(geoToolsTest.getAzimuth)
  }

執行結果：

[Point2D.Double[0.0, 0.0], Point2D.Double[10.000000000000002, 0.0], Point2D.Double[20.000000000000004, 0.0], Point2D.Double[29.999999999999996, 0.0], Point2D.Double[40.00000000000001, 0.0], Point2D.Double[50.00000000000001, 0.0], Point2D.Double[59.99999999999999, 0.0]]
90.0

由於我用的是本地模式（本地代理上網sbt實在是有點慢），先下載geotools的java包，在SourceForge上面，最新版本好像已經是20了，這裡廢話就不多說了。

最後

回到最開始的問題，其實計算經緯度的球面插值，最為核心的就是將經緯度轉化為球面座標，然後根據球面插值的公式進行插值計算，計算出的插值點也為球面座標，然後再將球面座標轉換為通用的經緯度即可。