scala实现球面插值（Slerp）

Posted 2022-12-14 通凡

一、球面插值

球面插值的原理大概就如下图所示，大致理解就是计算球面角度的占比，计算公式不是太复杂，如下所示：

当角度无限接近于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了，这里废话就不多说了。

最后

回到最开始的问题，其实计算经纬度的球面插值，最为核心的就是将经纬度转化为球面坐标，然后根据球面插值的公式进行插值计算，计算出的插值点也为球面坐标，然后再将球面坐标转换为通用的经纬度即可。