geometrybest-fit-curvebest-fit

Deduce center of circle from portion of circumference


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())
        }
    }
}

Solution

  • 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 enter image description here

    (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:

    enter image description here enter image description here

    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")
    }