1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:

import numpy as np
import matplotlib.pyplot as plt

p1 = np.array([3, 5])
p2 = np.array([5, 6])
p3 = np.array([8, 2])
pts = np.dstack([p1, p2, p3, p1])[0]

x1, y1 = p1[0], p1[1]
x2, y2 = p2[0], p2[1]
x3, y3 = p3[0], p3[1]

# Calculate distances between points
c = (x1 - x2)**2 + (y1 - y2)**2
a = (x2 - x3)**2 + (y2 - y3)**2
b = (x3 - x1)**2 + (y3 - y1)**2

# Calculate the cross-product
s = 2 * (a * b + b * c + c * a) - (a**2 + b**2 + c**2)

# Calculate the central coordinate of the circle
mx = (a * (b + c - a) * x1 + b * (c + a - b) * x2 + c * (a + b - c) * x3) / s
my = (a * (b + c - a) * y1 + b * (c + a - b) * y2 + c * (a + b - c) * y3) / s

# Calculate the radius of the circle
ar = a**0.5
br = b**0.5
cr = c**0.5
r = ar * br * cr / ((ar + br + cr) * (-ar + br + cr) * (ar - br + cr) * (ar + br - cr))**0.5


print("Central coordinate (x, y) and radius of the circle:")
print("cx={:>.3f}".format(mx), "cy={:>.3f}".format(my), "cr={:>.3f}".format(r)    )

# cx=5.227 cy=3.045 cr=2.963

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:

# draw
fig, axes = plt.subplots()
# draw points
plt.scatter(  pts[0], pts[1], color = 'red')

# connect the points
plt.plot(pts[0], pts[1])
plt.plot(pts[0], pts[1])

#draw centerpoint
plt.scatter(mx, my, color='green')

# connect centerpoint with outer points
x1 = np.array([mx, p1[0], mx, p2[0],mx,p3[0]])
y1 = np.array([my, p1[1], my, p2[1],my, p3[1]])
plt.plot(x1, y1, color='green')

# draw the circle
theta = np.linspace( 0 , 2 * np.pi , 150 )
a = r * np.cos( theta )+ mx
b = r * np.sin( theta )+ my

#draw Location vector:
x1 = np.array([0, mx ])
y1 = np.array([0, my ])
plt.plot(x1, y1, color='black')

# axes.plot( a, b )
axes.set_aspect( 1 )
plt.plot( a, b )
plt.xlim(-1, 10)
plt.ylim(-1, 10)
plt.grid()
plt.show()

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:

aq = (x2 - x3)**2 + (y2 - y3)**2
bq = (x3 - x1)**2 + (y3 - y1)**2
cq = (x1 - x2)**2 + (y1 - y2)**2

# Calculate the barycentric coordinates of the circle center
bar_a = aq * (bq + cq - aq)
bar_b = bq * (cq + aq - bq)
bar_c = cq * (aq + bq - cq)
bar_sum = bar_a + bar_b + bar_c

# Calculate the central coordinate of the circle
mx = bar_a * x1 + bar_b * x2 + bar_c * x3) / bar_sum
my = bar_a * y1 + bar_b * y2 + bar_c * y3) / bar_sum

# Calculate the radius of the circle
r = (aq * bq * cq / bar_sum)**0.5

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:

#      3   5   1 |  34
#      5   6   1 |  61
#      8   2   1 |  68
# 

import numpy as np
from numpy import linalg as LA

x = np.array([ 3, 5, 8])
y = np.array([ 5, 6, 2])
z = np.array([ 1, 1, 1])
t = np.array([ 34, 61, 68])

# Cramersche Regel
M = np.stack( [x, y, z], axis=1)
Mx = np.stack( [t, y, z], axis=1)
My = np.stack( [x, t, z], axis=1)
Mz = np.stack( [x, y, t], axis=1)

D  = LA.det(M)      #  -11.000000000000002 
Dx = LA.det(Mx)     # -115.00000000000001
Dy = LA.det(My)     #  -66.99999999999997 
Dz = LA.det(Mz)     #  306.0 

x = Dx/D      #  10.454545454545455 
y = Dy/D      #   6.090909090909087 
z = Dz/D      # -27.818181818181813

print( x/2, y/2, z)
# 5.2272727272727275 3.0454545454545436 -27.818181818181813

print((np.sqrt(abs(z)/np.pi)))       # 2.9757019824548037


# other methods result, out:
# r=  2.9632738205466587 center_x=  5.2272727272727275 center_y= 3.0454545454545454
# cx=5.227 cy=3.045 cr=2.963

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:

import numpy as np
from numpy import linalg as LA

x = np.array([ 3, 5, 8])
y = np.array([ 5, 6, 2])
z = np.array([ 1, 1, 1])
t = np.array([ x[0]**2+y[0]**2, x[1]**2+y[1]**2, x[2]**2+y[2]**2])

# Cramersche Regel
M = np.stack( [x, y, z], axis=1)
Mx = np.stack( [t, y, z], axis=1)
My = np.stack( [x, t, z], axis=1)
Mz = np.stack( [x, y, t], axis=1)

D  = LA.det(M)      #  -11.000000000000002 
Dx = LA.det(Mx)     # -115.00000000000001
Dy = LA.det(My)     #  -66.99999999999997 
Dz = LA.det(Mz)     #  306.0 

x = Dx/(2*D)
y = Dy/(2*D)
z = Dz/D      # -27.818181818181813
r = np.sqrt(x**2 + y**2 + z)

print( x, y, r)
# 5.2272727272727275 3.0454545454545436 2.9632738205466578