0

我想在 Scala 中为给定的 2D 点拟合一个圆。

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$1(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$1$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 是否应该是一个真正的函数或对象或任何东西。

4

2 回答 2

0

在对代码进行了长时间的摆弄之后，我让它运行起来了

在我将 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 编译

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

于 2019-02-08T09:20:40.973 回答
0

在对圆拟合进行了一些研究之后，我在论文中发现了一个很棒的算法：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......这个值的特征向量产生了正确的结果，但是因为它被整理出来了负特征值。所以这个并不总是稳定的（与许多其他圆拟合算法一样）

但是多么好的算法啊！！！在我看来有点像黑魔法。

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

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

于 2019-02-14T08:53:29.287 回答