I want to fit a circle to given 2D points in Scala.
Apache commons math has an example for this in java, which I am trying to translate to scala (without success, because my knowledge of Java is almost non existent).
I took the example code from "http://commons.apache.org/proper/commons-math/userguide/leastsquares.html", (see end of page) which I tried to translate into 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())
}
This gives no compile errors, but crashes with:
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)
I guess the problem is somewhere in the definition of the function distancesToCurrentCenter. I don't even know if this MultivariateJacobianFunction is supposed to be a real function or an object or what ever.
After some long fiddeling with the code, I got it running
The NullPointerException was gone after I updated apache-commons-math3 from version 3.3 to version 3.6.1 in my build.sbt file. Don't know if I forgot a paramater of if it was a bug. There were also 2 bugs in the example on the apache-commons-math website: They had two times a .getX operator where should have been an .getY.
So here is a running example for a circle fit with known radius:
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")
}
Tested on Scala version 2.12.6, compiled with sbt version 1.2.8
Does anabody know how to do this without a fixed radius?