0

When searching for code to calculate the circle out of 3 points it leads me to this code:

def circleRadius(b, c, d):
  temp = c[0]**2 + c[1]**2
  bc = (b[0]**2 + b[1]**2 - temp) / 2
  cd = (temp - d[0]**2 - d[1]**2) / 2
  det = (b[0] - c[0]) * (c[1] - d[1]) - (c[0] - d[0]) * (b[1] - c[1])

  if abs(det) < 1.0e-10:
    return None

  # Center of circle
  cx = (bc*(c[1] - d[1]) - cd*(b[1] - c[1])) / det
  cy = ((b[0] - c[0]) * cd - (c[0] - d[0]) * bc) / det

  radius = ((cx - b[0])**2 + (cy - b[1])**2)**.5

  return radius

Based on Stackoverflow and Dr.Math. The code works perfectly, but I don't understand how the code fits to the explanation given at Dr.Math.

Can anyone help me to understand why the code is working and what substeps are implemented in the variables?

Glorfindel
  • 21,988
  • 13
  • 81
  • 109
Seppelandrio
  • 11
  • 1
  • 1
  • It computes the intersection of the perpendicular bisectors of the lines `BC` and `CD`. There's a fair amount of algebra involved but it is not too hard to prove – only requires grade-school geometry. – meowgoesthedog Oct 25 '18 at 13:31
  • But where does it calculate the perpendicular bisectors. BC and CD are values and no points and of what do I calculate the determinante (to determine the slope?) and how do they come together at cx,cy? I do understand the principle idea behind this approach, but can't map the steps in code to it... – Seppelandrio Oct 25 '18 at 14:54
  • If you "de-compress" this code into the relevant steps it boils down to the following: **1)** compute the midpoints of `BC` and `CD`; **2)** calculate the perpendicular directions to `BC` and `CD`; **3)** use the results of (1) and (2) to define the perpendicular bisector lines; **4)** solve a pair of simultaneous equations to obtain the intersection point coordinates – this is where the determinant is used. Work out these variables, do a bunch of tedious algebra to simplify, and you get the equations in the code above. – meowgoesthedog Oct 25 '18 at 15:01
  • What is the value of variable bc and dc and where does the slope come from? To calculate (1) I would suggest code like BC = ((B[0]+C[0])/2, (B[1],C[1])/2) and (2) with -1/(slope of BC)... but where does this calculation take place? – Seppelandrio Oct 25 '18 at 15:11
  • 1
    @Seppelandrio Forget about slopes. Rather simple math description for intersection point: https://en.wikipedia.org/wiki/Intersection_(Euclidean_geometry)#On_a_plane – MBo Oct 25 '18 at 15:33

2 Answers2

1

The code you see is a "simplified" and concise formula of the procedure described in the Dr. Math page.

Let us go over it step-by-step.

  1. For the sake of simplicity and adhering to mathematical notations, let points be the points on the triangle. [b is point 1, c is point 2, and d is point 3]

  2. For such triangle, the area is defined as:

    area

    This variable det in the function is equal to 2 * area of the triangle.

    The if abs(det) < 1.0e-10: is checking for collinearity. If the area is close to zero, the points given are collinear i.e. they are points on a single line.

  3. Find the slopes of the lines L1, L2 passing through points b, c & c, d

    slope

  4. Find the equations for lines L3, L4, which are the perpendicular bisectors of line sL1 and L2 respectively. bisectors1 bisectors2

  5. Find the intersection of lines L3 and L4, which is nothing but the center of the cirle.

    intersection

Do all the substitutions and you can see it all come together.

  1. Compute radius of the circle by finding the euclidean distance between the center and one of the three points.
pbskumar
  • 1,127
  • 12
  • 14
  • What are you going to do if `y1=y2 or y3=y2 or x2=x1` etc? – MBo Oct 26 '18 at 07:05
  • Thank you! So I understood the term! – Seppelandrio Oct 26 '18 at 09:06
  • @MBo, After all the variables, are substituted, the difference terms causing a divide by zero will be eliminated. That is evident in the formula used to calculate `cx` and `cy`. I think *divide by zero* is the reason this specific piece of code does not calculate slopes. By avoiding calculating slopes beforehand, we can avert that catastrophic divide/multiply by zero problems. – pbskumar Oct 26 '18 at 14:45
0

Not really an answer, I am showing here an alternative method.

Let p, q, r be the three points. We translate them so that p comes to the origin. Vectorially, q-= p, r-= p.

Now the equation of a point by the origin is

2.xc.x + 2.yc.y = x² + y²

where xc, yc are the coordinates of the center.

The plugging the coordinates of p, q, we get a 2x2 system

xc.xp + yc.yp = xp² + yp²
xc.xq + yc.yq = xq² + yq²

The code that implements this is given below

# Translate to the origin
xq-= xp
yq-= yp
q2= xp * xp + yp * yp
xr-= xq
yr-= yq
r2= xr * xr + yr * yr

# Solve for the center coordinates
d= 2 * (xp * yq - xq * yp)
xc= (p2 * yq - q2 * yp) / d
yc= (p2 * xp - q2 * xq) / d

# Radius
r= math.sqrt(xc * xc + yc * yc)

# Untranslate
xc+= xp
yc+= yp