circlefitting.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
 
"""
Circular regression (Circle-fitting)
Calculates the best fitting circle by the given Points.
Based on http://www.block-net.de/messtechnik/p-messtechnik.html
"""
 
__author__      = "Kevin Gliewe aka KillerGoldFisch"
__copyright__   = "Copyright 2014, Kevin Gliewe"
__credits__     = ["Kevin Gliewe",]
__license__     = "MIT"
__version__     = "1.0.0"
__date__        = "2014-05-26"
__maintainer__  = "Kevin Gliewe"
__email__       = "kevingliewe@gmail,com"
__status__      = "Production"
 
"""
Changelog:
    no changes
"""

import math

#      ____  ____  _____   ________
#     / __ \/ __ \/  _/ | / /_  __/
#    / /_/ / / / // //  |/ / / /   
#   / ____/ /_/ // // /|  / / /    
#  /_/    \____/___/_/ |_/ /_/     
#                                  

class Point:

    def __init__(self, x, y):
        self.X = x
        self.Y = y

    def __str__(self):
        return "Point(X=%s,Y=%s)"%(self.X, self.Y) 

#     _____ __  _______           __   
#    / __(_) /_/ ____(_)_________/ /__ 
#   / /_/ / __/ /   / / ___/ ___/ / _ \
#  / __/ / /_/ /___/ / /  / /__/ /  __/
# /_/ /_/\__/\____/_/_/   \___/_/\___/ 

def fitCircle(points):
    xsum = sum([p.X for p in points])
    ysum = sum([p.Y for p in points])
    xsumsq = sum([(p.X**2) for p in points])
    ysumsq = sum([(p.Y**2) for p in points])
    n = len(points)

    #print "xsum=%d ysum=%d xsumsq=%d ysumsq=%d n=%d"%(xsum, ysum, xsumsq, ysumsq, n)

    def getA(points, xsum, ysum, n):
        a = 0.0
        for p in points:
            a += p.X * (2*n*p.X - 2*xsum)
        return a

    def getB(points, xsum, ysum, n):
        b = 0.0
        for p in points:
            b += p.Y * (2*n*p.X - 2*xsum)
        return b

    def getC(points, xsum, ysum, n):
        c = 0.0
        for p in points:
            c += p.X * (2*n*p.Y - 2*ysum)
        return c

    def getD(points, xsum, ysum, n):
        d = 0.0
        for p in points:
            d += p.Y * (2*n*p.Y - 2*ysum)
        return d

    def getE(points, xsum, ysum, xsumsq, ysumsq, n):
        e = 0.0
        for p in points:
            e += p.X * (n*p.X**2 + n*p.Y**2 - xsumsq - ysumsq)
        return e

    def getF(points, xsum, ysum, xsumsq, ysumsq, n):
        f = 0.0
        for p in points:
            f += p.Y * (n*p.X**2 + n*p.Y**2 - xsumsq - ysumsq)
        return f

    def getX0(a, b, c, d, e, f):
        return (d*e-c*f)/(a*d-b*c)

    def getY0(a, b, c, d, e, f):
        return (a*f-b*e)/(a*d-b*c)

    def getR(points, n, x0, y0):
        r = 0.0
        for p in points:
            r += ((p.X - x0)**2 + (p.Y - y0)**2)
        r = (r / n)
        r = math.sqrt(r)
        return r

    a = getA(points, xsum, ysum, n)
    b = getB(points, xsum, ysum, n)
    c = getC(points, xsum, ysum, n)
    d = getD(points, xsum, ysum, n)
    e = getE(points, xsum, ysum, xsumsq, ysumsq, n)
    f = getF(points, xsum, ysum, xsumsq, ysumsq,  n)
    x0 = getX0(a, b, c, d, e, f)
    y0 = getY0(a, b, c, d, e, f)
    r = getR(points, n, x0, y0)

    #print "a=%d b=%d c=%d d=%d e=%d f=%f"%(a,b,c,d,e,f)

    return (Point(x0, y0), r) #dict(x0=x0, y0=y0, r=r)




#      __  ______    _____   __
#     /  |/  /   |  /  _/ | / /
#    / /|_/ / /| |  / //  |/ / 
#   / /  / / ___ |_/ // /|  /  
#  /_/  /_/_/  |_/___/_/ |_/   
#                            

def main(argv):
    import random
    a45 = math.sin(math.pi/4) * 2
    x0=88
    y0=-880
    """
    points = [
        #Point(0.0, 1.0),
        Point(  a45 + x0, -a45 + y0),
        Point( -a45 + x0, -a45 + y0),
        Point( -a45 + x0,  a45 + y0),
        Point(  a45 + x0,  a45 + y0),
    ]
    """
    points = []
    for n in range(10000):
        r = 5.8 + random.random()*0.4
        angle = 2 * math.pi * random.random()
        x = x0 +  r * math.sin(angle)
        y = y0 + r * math.cos(angle)
        points.append(Point(x, y))
    #"""

    if len(points) < 101:
        print [str(p) for p in points]
    
    (center, radius) = fitCircle(points)
    print "Center: %s Radius:%s"%(center, radius)
 
 
if __name__ == "__main__":
    import sys
    main(sys.argv[1:])