I have sequences of points forming arcs(?). I would like to deduce the best fit circular (or even ellipsoid) curve for these points. Each arc has a relatively consistent change of angle across its length (that is partly how I isolate them).
An algorithm would be great (especially if explained so I can understand it), but since this is probably known territory, pointers to answers would also be good. If all else fails some better search terms would be appreciated.
EDIT: I am not confident in using matrix transformations, and certainly not confident in coding solutions described in them..... a solution without matrices would be much appreciated.
(I am fluent in c/java/kotlin languages so solutions in something like that would be easiest for me, but I am happy to work from any language or pseudocode).
Thanks in advance
Thanks to Renat for his solution - it works a treat. I am a little miffed that I do not understand the algorithm, but there is a reference and I can live with the ignorance for now.
For anyone needing the code translated to kotlin:
class Circle(val x: Int, val y: Int, val radius: Int) {
companion object {
/**
* Find the best fit circle for a collection of points.
* With many thanks to Renat:
* https://stackoverflow.com/questions/74493481/deduce-center-of-circle-from-portion-of-circumference
*/
fun fit(points: List<Point>): Circle {
val sx = points.sumOf { it.x.toDouble() }
val sy = points.sumOf { it.y.toDouble() }
val sx2 = points.sumOf { it.x.toDouble() * it.x }
val sy2 = points.sumOf { it.y.toDouble() * it.y }
val sxy = points.sumOf { it.x.toDouble() * it.y }
val sx3 = points.sumOf { it.x.toDouble() * it.x * it.x }
val sy3 = points.sumOf { it.y.toDouble() * it.y * it.y }
val sx2y = points.sumOf { it.x.toDouble() * it.x * it.y }
val sxy2 = points.sumOf { it.x.toDouble() * it.y * it.y }
val n = points.size;
val s11 = n * sxy - sx * sy
val s20 = n * sx2 - sx * sx
val s02 = n * sy2 - sy * sy
val s30 = n * sx3 - sx2 * sx
val s03 = n * sy3 - sy * sy2
val s21 = n * sx2y - sx2 * sy
val s12 = n * sxy2 - sx * sy2
val a = ((s30 + s12)*s02 - (s03 + s21)*s11) / (2*( s20 * s02 - s11 * s11));
val b = ((s03 + s21)*s20 - (s30 + s12)*s11) / (2*( s20 * s02 - s11 * s11));
val c = (sx2 + sy2 - 2*a*sx - 2*b*sy) / n;
val radius = Math.sqrt(c + a*a + b*b);
return Circle(a.roundToInt(), b.roundToInt(), radius.roundToInt())
}
}
}
Ok, so it's a Fitting a circle (or ellipse) problem. Which is solvable in linear time.
For code below I used answer from this question from Math.stackexchange, which references Régressions coniques, quadriques. Régressions linéaires et apparentées, circulaire, sphérique, pp 12-13, which uses a method of least squares to minimize the sum
(press 'Run code snipppet', and draw with mouse)
var canvas = document.createElement('canvas');
document.body.appendChild(canvas);
document.body.style.margin = 0;
canvas.style.position = 'fixed';
// get canvas 2D context and set him correct size
var ctx = canvas.getContext('2d');
resize();
let points = [];
let unclicked = false;
window.addEventListener('resize', resize);
document.addEventListener('mousemove', draw);
document.addEventListener('mousedown', setPosition);
function getApproxCircle(p) {
// Régressions coniques, quadriques, régressions linéaires et apparentées, circulaire, sphérique
// pp 12-13
let Sx = 0
let Sy = 0
let Sx2 = 0
let Sy2 = 0
let Sxy = 0
let Sx3 = 0
let Sy3 = 0
let Sx2y = 0
let Sxy2 = 0
for(point of p ){
const [x,y] = [point.x,point.y]
const [x2,y2] = [x*x,y*y]
const [x3,y3] = [x2*x,y2*y]
const xy = x*y
const [x2y,xy2] = [xy*x,xy*y]
Sx += x
Sy += y
Sx2 += x2
Sy2 += y2
Sxy += xy
Sx3 += x3
Sy3 += y3
Sx2y += x2y
Sxy2 += xy2
}
const n = points.length;
const s11 = n * Sxy - Sx * Sy
const s20 = n * Sx2 - Sx * Sx
const s02 = n * Sy2 - Sy * Sy
const s30 = n * Sx3 - Sx2 * Sx
const s03 = n * Sy3 - Sy * Sy2
const s21 = n * Sx2y - Sx2 * Sy
const s12 = n * Sxy2 - Sx * Sy2
const a = ((s30 + s12)*s02 - (s03 + s21)*s11)
/ (2*( s20 * s02 - s11 * s11));
const b = ((s03 + s21)*s20 - (s30 + s12)*s11)
/ (2*( s20 * s02 - s11 * s11));
const c = (Sx2 + Sy2 - 2*a*Sx - 2*b*Sy) / n;
const radius = Math.sqrt(c + a*a + b*b);
return {
center: {
x: a,
y: b
},
radius: radius
};
}
function setPosition(e) {
if (unclicked && e.buttons === 1) {
unclicked = false;
points = [];
}
points.push({ x: e.clientX, y: e.clientY });
}
function resize() {
ctx.canvas.width = window.innerWidth;
ctx.canvas.height = window.innerHeight;
}
function clear() {
ctx.clearRect(0, 0, window.innerWidth, window.innerHeight);
}
function draw(e) {
if (e.buttons !== 1) {
unclicked = true;
return;
}
points.push({ x: e.clientX, y: e.clientY });
clear();
ctx.lineWidth = 5;
ctx.lineCap = 'round';
ctx.strokeStyle = '#c0392b';
//point of drawn arc
ctx.beginPath();
ctx.moveTo(points[0].x, points[0].y);
for(i=1; i<points.length; i++) {
ctx.lineTo(points[i].x, points[i].y);
}
ctx.stroke();
let approxCircle = getApproxCircle(points);
//circle
ctx.lineWidth = 2;
ctx.strokeStyle = '#3039cb77';
ctx.beginPath();
ctx.arc(approxCircle.center.x, approxCircle.center.y, approxCircle.radius, 0, 2 * Math.PI, false);
ctx.stroke();
// center of the circle
ctx.beginPath();
ctx.arc(approxCircle.center.x, approxCircle.center.y, 1, 0, 2 * Math.PI, false);
ctx.fillStyle = 'green';
ctx.fill();
ctx.lineWidth = 5;
ctx.strokeStyle = '#003300';
ctx.stroke();
}
And if it's needed to draw a circle by 3 points, there is a formula from Wikipedia:
in Kotlin (link to sandbox):
data class Point(
val x: Float,
val y: Float
) {
fun lenSquare() = x*x + y*y
}
fun main() {
val A = Point(0f, 0f)
val B = Point(0f, 2f)
val C = Point(2f, 0f)
val D = 2 * (A.x * (B.y - C.y)
+ B.x * (C.y - A.y)
+ C.x * (A.y - B.y))
var lenASqr = A.lenSquare()
var lenBSqr = B.lenSquare()
var lenCSqr = C.lenSquare()
val Ux = (lenASqr * (B.y - C.y)
+ lenBSqr * (C.y - A.y)
+ lenCSqr * (A.y - B.y)
) / D;
val Uy = (lenASqr * (C.x - B.x)
+ lenBSqr * (A.x - C.x)
+ lenCSqr * (B.x - A.x)
) / D;
val U = Point(Ux, Uy)
println("A: $A")
println("B: $B")
println("C: $C")
println("Circumcircle center: $U")
}