После некоторых исследований по подгонке окружности я нашел замечательный алгоритм в статье: "Анализ ошибок для алгоритмов подгонки окружности" Х. Аль-Шарадки и Н. Чернова (доступно здесь: http://people.cas.uab.edu/~mosya/cl/) Iреализовал это в 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-алгоритм) и попытался ограничить операции Matrix.
Это все еще не то, что мне нужно, потому что у алгоритма такого типа (алгебраическое соответствие) есть известные проблемы с подгонкой частично окружностей (дуг) и, возможно, больших окружностей.совершенно неверные результаты, и я обнаружил, что у меня есть собственное значение -0,1 ... Собственный вектор этого значения дал правильный результат, но он был отсортирован из-за отрицательного собственного значения.Так что этот не всегда стабилен (как и многие другие алгоритмы подбора окружности)
Но что за хороший алгоритм !!!Для меня это выглядит как темная магия.
Если кому-то не нужна большая точность и большая скорость (и у него есть данные от полного круга, а не большого), это был бы мой выбор.
Следующее, что я постараюсь, это реализовать алгоритм Левенберга Марквардта на той же странице, о которой я упоминал выше.