如何在 Scala 中编写适合的圆

Question

我想用 Scala 中给定的二维点拟合一个圆。

Apache commons math 在 java 中有一个例子，我正在尝试将其翻译成 scala（没有成功，因为我对 Java 的了解几乎不存在）。

我从“http://commons.apache.org/proper/commons-math/userguide/leastsquares.html”中获取了示例代码（见页尾），我试图将其翻译成 scala：

  import org.apache.commons.math3.linear._
  import org.apache.commons.math3.fitting._
  import org.apache.commons.math3.fitting.leastsquares._
  import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer._
  import org.apache.commons.math3._
  import org.apache.commons.math3.geometry.euclidean.twod.Vector2D
  import org.apache.commons.math3.util.Pair
  import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum

  def circleFitting: Unit = {
    val radius: Double = 70.0

    val observedPoints = Array(new Vector2D(30.0D, 68.0D), new Vector2D(50.0D, -6.0D), new Vector2D(110.0D, -20.0D), new Vector2D(35.0D, 15.0D), new Vector2D(45.0D, 97.0D))

    // the model function components are the distances to current estimated center,
    // they should be as close as possible to the specified radius

    val distancesToCurrentCenter = new MultivariateJacobianFunction() {
      //def value(point: RealVector): (RealVector, RealMatrix) = {
      def value(point: RealVector): Pair[RealVector, RealMatrix] = {

        val center = new Vector2D(point.getEntry(0), point.getEntry(1))

        val value: RealVector = new ArrayRealVector(observedPoints.length)
        val jacobian: RealMatrix = new Array2DRowRealMatrix(observedPoints.length, 2)

        for (i <- 0 to observedPoints.length) {
          var o = observedPoints(i)
          var modelI: Double = Vector2D.distance(o, center)
          value.setEntry(i, modelI)
          // derivative with respect to p0 = x center
          jacobian.setEntry(i, 0, (center.getX() - o.getX()) / modelI)
          // derivative with respect to p1 = y center
          jacobian.setEntry(i, 1, (center.getX() - o.getX()) / modelI)
        }
        new Pair(value, jacobian)
      }
    }

    // the target is to have all points at the specified radius from the center
    val prescribedDistances = Array.fill[Double](observedPoints.length)(radius)
    // least squares problem to solve : modeled radius should be close to target radius
    
    val problem:LeastSquaresProblem = new LeastSquaresBuilder().start(Array(100.0D, 50.0D)).model(distancesToCurrentCenter).target(prescribedDistances).maxEvaluations(1000).maxIterations(1000).build()
    
    val optimum:Optimum = new LevenbergMarquardtOptimizer().optimize(problem) //LeastSquaresOptimizer.Optimum
    val fittedCenter: Vector2D = new Vector2D(optimum.getPoint().getEntry(0), optimum.getPoint().getEntry(1))
    println("circle fitting wurde aufgerufen!")
    println("CIRCLEFITTING: fitted center: " + fittedCenter.getX() + " " + fittedCenter.getY())
    println("CIRCLEFITTING: RMS: " + optimum.getRMS())
    println("CIRCLEFITTING: evaluations: " + optimum.getEvaluations())
    println("CIRCLEFITTING: iterations: " + optimum.getIterations())
    
  }

这不会产生编译错误，但会崩溃：

Exception in thread "main" java.lang.NullPointerException
    at org.apache.commons.math3.linear.EigenDecomposition.<init>(EigenDecomposition.java:119)
    at org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory.squareRoot(LeastSquaresFactory.java:245)
    at org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory.weightMatrix(LeastSquaresFactory.java:155)
    at org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory.create(LeastSquaresFactory.java:95)
    at org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder.build(LeastSquaresBuilder.java:59)
    at twoDhotScan.FittingFunctions$.circleFitting(FittingFunctions.scala:49)
    at twoDhotScan.Main$.delayedEndpoint$twoDhotScan$Main(hotScan.scala:14)
    at twoDhotScan.Main$delayedInit$body.apply(hotScan.scala:11)
    at scala.Function0.apply$mcV$sp(Function0.scala:34)
    at scala.Function0.apply$mcV$sp$(Function0.scala:34)
    at scala.runtime.AbstractFunction0.apply$mcV$sp(AbstractFunction0.scala:12)
    at scala.App.$anonfun$main$adapted(App.scala:76)
    at scala.collection.immutable.List.foreach(List.scala:389)
    at scala.App.main(App.scala:76)
    at scala.App.main$(App.scala:74)
    at twoDhotScan.Main$.main(hotScan.scala:11)
    at twoDhotScan.Main.main(hotScan.scala)

我猜问题出在函数 distancesToCurrentCenter 的定义中。我什至不知道这个 MultivariateJacobianFunction 应该是一个真正的函数还是一个对象或者什么。

Answer 1

经过长时间修改代码，我搞定了运行

在我的 build.sbt 文件中将 apache-commons-math3 从版本 3.3 更新到版本 3.6.1 后，NullPointerException 消失了。不知道我是否忘记了它是否是错误的参数。在 apache-commons-math 网站上的示例中也有 2 个错误：他们有两次 .getX 运算符应该是 .getY。

所以这是一个运行已知半径的圆拟合示例：

import org.apache.commons.math3.analysis.{ MultivariateVectorFunction, MultivariateMatrixFunction }
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum
import org.apache.commons.math3.fitting.leastsquares.{ MultivariateJacobianFunction, LeastSquaresProblem, LeastSquaresBuilder, LevenbergMarquardtOptimizer }
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D
import org.apache.commons.math3.linear.{ Array2DRowRealMatrix, RealMatrix, RealVector, ArrayRealVector }

object Main extends App {
  val radius: Double = 20.0
  val pointsList: List[(Double, Double)] = List(
    (18.36921795, 10.71416674),
    (0.21196357, -22.46528791),
    (-4.153845171, -14.75588526),
    (3.784114125, -25.55910336),
    (31.32998899, 2.546924253),
    (34.61542186, -12.90323269),
    (19.30193011, -28.53185596),
    (16.05620863, 10.97209111),
    (31.67011956, -20.05020878),
    (19.91175561, -28.38748712))
/*******************************************************************************
 ***** Random values on a circle with centerX=15, centerY=-9 and radius 20 *****
 *******************************************************************************/

  val observedPoints: Array[Vector2D] = (pointsList map { case (x, y) => new Vector2D(x, y) }).toArray

  val vectorFunktion: MultivariateVectorFunction = new MultivariateVectorFunction {
    def value(variables: Array[Double]): Array[Double] = {
      val center = new Vector2D(variables(0), variables(1))
      observedPoints map { p: Vector2D => Vector2D.distance(p, center) }
    }
  }

  val matrixFunction = new MultivariateMatrixFunction {
    def value(variables: Array[Double]): Array[Array[Double]] = {
      val center = new Vector2D(variables(0), variables(1))
      (observedPoints map { p: Vector2D => Array((center.getX - p.getX) / Vector2D.distance(p, center), (center.getY - p.getY) / Vector2D.distance(p, center)) })
    }
  }

  // the target is to have all points at the specified radius from the center
  val prescribedDistances = Array.fill[Double](observedPoints.length)(radius)
  // least squares problem to solve : modeled radius should be close to target radius
  val problem = new LeastSquaresBuilder().start(Array(100.0D, 50.0D)).model(vectorFunktion, matrixFunction).target(prescribedDistances).maxEvaluations(25).maxIterations(25).build
  val optimum: Optimum = new LevenbergMarquardtOptimizer().optimize(problem)
  val fittedCenter: Vector2D = new Vector2D(optimum.getPoint.getEntry(0), optimum.getPoint.getEntry(1))

  println("Ergebnisse des LeastSquareBuilder:")
  println("CIRCLEFITTING: fitted center: " + fittedCenter.getX + " " + fittedCenter.getY)
  println("CIRCLEFITTING: RMS: " + optimum.getRMS)
  println("CIRCLEFITTING: evaluations: " + optimum.getEvaluations)
  println("CIRCLEFITTING: iterations: " + optimum.getIterations + "\n")
}

在 Scala 版本 2.12.6 上测试，使用 sbt 版本 1.2.8 编译

有人知道如何在没有固定半径的情况下做到这一点吗？

Answer 2

在对圆拟合进行一些研究后，我在论文中发现了一个很棒的算法："Error alalysis for circle fitting algorithms"，作者 H. Al-Sharadqah 和 N. Chernov（可在此处获取：http://people.cas.uab.edu/~mosya/cl/）我在 scala 中实现了它：

import org.apache.commons.math3.linear.{ Array2DRowRealMatrix, RealMatrix, RealVector, LUDecomposition, EigenDecomposition }

object circleFitFunction {
  def circleFit(dataXY: List[(Double, Double)]) = {

    def square(x: Double): Double = x * x
    def multiply(pair: (Double, Double)): Double = pair._1 * pair._2

    val n: Int = dataXY.length
    val (xi, yi) = dataXY.unzip
    //val S: Double = math.sqrt(((xi map square) ++ yi map square).sum / n)
    val zi: List[Double] = dataXY map { case (x, y) => x * x + y * y }
    val x: Double = xi.sum / n
    val y: Double = yi.sum / n
    val z: Double = ((xi map square) ++ (yi map square)).sum / n
    val zz: Double = (zi map square).sum / n
    val xx: Double = (xi map square).sum / n
    val yy: Double = (yi map square).sum / n
    val xy: Double = ((xi zip yi) map multiply).sum / n
    val zx: Double = ((zi zip xi) map multiply).sum / n
    val zy: Double = ((zi zip yi) map multiply).sum / n

    val N: RealMatrix = new Array2DRowRealMatrix(Array(
      Array(8 * z, 4 * x, 4 * y, 2),
      Array(4 * x, 1, 0, 0),
      Array(4 * y, 0, 1, 0),
      Array(2.0D, 0, 0, 0)))

    val M: RealMatrix = new Array2DRowRealMatrix(Array(
      Array(zz, zx, zy, z),
      Array(zx, xx, xy, x),
      Array(zy, xy, yy, y),
      Array(z, x, y, 1.0D)))

    val Ninverse = new LUDecomposition(N).getSolver().getInverse()
    val eigenValueProblem = new EigenDecomposition(Ninverse.multiply(M))
    // Get all eigenvalues
    // As we need only the smallest positive eigenvalue, all negative eigenvalues are replaced by Double.MaxValue
    val eigenvalues: Array[Double] = eigenValueProblem.getRealEigenvalues() map (lambda => if (lambda < 0) Double.MaxValue else lambda)

    // Now get the index of the smallest positive eigenvalue, to get the associated eigenvector
    val i: Int = eigenvalues.zipWithIndex.min._2
    val eigenvector: RealVector = eigenValueProblem.getEigenvector(3)

    val A = eigenvector.getEntry(0)
    val B = eigenvector.getEntry(1)
    val C = eigenvector.getEntry(2)
    val D = eigenvector.getEntry(3)

    val centerX: Double = -B / (2 * A)
    val centerY: Double = -C / (2 * A)
    val Radius: Double = math.sqrt((B * B + C * C - 4 * A * D) / (4 * A * A))
    val RMS: Double = (dataXY map { case (x, y) => (Radius - math.sqrt((x - centerX) * (x - centerX) + (y - centerY) * (y - centerY))) } map square).sum / n
    (centerX, centerY, Radius, RMS)
  }
}

我保留了论文中的所有名称（请参阅第 4 章和第 8 章并寻找 Hyperfit-Algorithm）并尝试限制矩阵运算。

它仍然不是我需要的，因为这种算法（代数拟合）在拟合部分圆（弧）和大圆方面存在已知问题。

根据我的数据，我曾经遇到过完全错误的结果，我发现我的特征值为 -0.1... 这个Value的Eigenvector产生了正确的结果，但是因为Eigenvalue为负而被整理出来了。所以这个并不总是稳定的（就像许多其他圆拟合算法一样）

但是多么好的算法！！！我觉得有点像黑魔法。

如果有人不需要太高的精度和速度（并且从一个完整的圆圈得到的数据不是很大），这将是我的选择。

接下来我将尝试从我上面提到的同一页面实施 Levenberg Marquardt 算法。

如何在 Scala 中编写适合的圆

How to program a circle fit in scala

math

geometry

scala

function-fitting