Calculating the inverse of the radius of a circle in R

Question

I'm trying to replicate a function used in a study, but don't really have the mathematical background to fully appreciate how this ought to be done. The measure takes three points from a tongue contour and uses these three points to calculate the radius of a circle that would pass through them. I have looked here and found something that does this in python. I've tried to modify the code so it would work in R with my own data. (Posted at the bottom)

The problem is, based on the study I am reading, I then need to calculate the concavity of the circumference of the circle and find the inverse of the radius of the circle passing through the three points. I'm googling and googling but honestly this means nothing to me. The only thing I have found is that I seem to need to calculate the first and second derivatives of the tongue surface curve. I'm really hoping somebody might be able to help explore how I would do this in R. To be brutally honest, I am not overly interested in understanding the mathematics here, just how to actually implement it.

Edit: I thought below was the formula that I need to replicate. As MBo points out, this isn't the case.

I'll repeat something from another study that used a very, very similar method in case that helps.

'Any three points (A, B, C) can be conceived as lying on the circumference of a circle. The circle will have a radius, the inverse of which represents the curvature of the circle passing through those three points.' The set of three points 'yields a curvature numeber which is the inverse of the radius of the circle passing through them. Three points which lie along a straight line have a curvature of zero, since their concavity is zero and this becomes the numerator of of the curvature equation'. It's this that I need to do, but don't know where to begin operationalising it in R.

The code below is the python code I'm attempting to replicate for my purposes in R to obtain the radius from three points. I have no idea how to proceed after that.

def define_circle(p1, p2, p3):
    """
    Returns the center and radius of the circle passing the given 3 points.
    In case the 3 points form a line, returns (None, infinity).
    """
    temp = p2[0] * p2[0] + p2[1] * p2[1]
    bc = (p1[0] * p1[0] + p1[1] * p1[1] - temp) / 2
    cd = (temp - p3[0] * p3[0] - p3[1] * p3[1]) / 2
    det = (p1[0] - p2[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p2[1])

    if abs(det) < 1.0e-6:
        return (None, np.inf)

    # Center of circle
    cx = (bc*(p2[1] - p3[1]) - cd*(p1[1] - p2[1])) / det
    cy = ((p1[0] - p2[0]) * cd - (p2[0] - p3[0]) * bc) / det

    radius = np.sqrt((cx - p1[0])**2 + (cy - p1[1])**2)
    return ((cx, cy), radius)

Here's my R attempt. I haven't written the function yet, but I will be looking at three points along a curve, A, B and C. The function will extract x and y values for each of these three points (called x_value_a, y_value_a etc.). Once this is done. I will run the code that follows. It's after this that I am properly stumped.

temp = x_value_b ^ 2 + y_value_b ^ 2

bc = (x_value_a ^ 2 + y_value_a ^ 2 - temp) / 2

cd = (temp - x_value_c ^ 2 - y_value_c ^ 2) / 2

det = (x_value_a - x_value_b) * (y_value_b - y_value_c) - (x_value_b - x_value_c) * (y_value_a - y_value_b)

cx = (bc * (y_value_b - y_value_c) - cd * (y_value_a - y_value_b)) / det 

cy = ((x_value_a - x_value_b) * cd - (x_value_b - x_value_c) * bc) / det

radius = sqrt((cx - x_value_a)^2 + (cy - y_value_a)^2)

Any help would be greatly appreciated. I'm sorry for my mathematical ignorance.

Answer 1

If you only want the Python script translated into R, that's pretty straightforward (I don't quite understand why you split it up in the R code you added).

define_circle = function(p1, p2, p3) {

  # Returns the center and radius of the circle passing the given 3 points.
  # In case the 3 points form a line, returns warning.
  
  temp = p2[1] * p2[1] + p2[2] * p2[2]
  bc = (p1[1] * p1[1] + p1[2] * p1[2] - temp) / 2
  cd = (temp - p3[1] * p3[1] - p3[2] * p3[2]) / 2
  det = (p1[1] - p2[1]) * (p2[2] - p3[2]) - (p2[1] - p3[1]) * (p1[2] - p2[2])
  
  if (abs(det) < 1.0e-6) {
    
    return(c("Three points form a line"))
    
  } else {
    
    # Center of circle
    cx = (bc*(p2[2] - p3[2]) - cd*(p1[2] - p2[2])) / det
    cy = ((p1[1] - p2[1]) * cd - (p2[1] - p3[1]) * bc) / det
    
    radius = sqrt((cx - p1[1])**2 + (cy - p1[2])**2)
    
    return(list("center" = c(cx, cy), "radius" = radius))
    
  }

}

Note that p1-3 represents a vector containing an x- and y-coordinate. I have to trust the original Python code here but a quick check using desmos.com seems to indicate it works:

> define_circle(c(0,1), c(2,2), c(0.5,5))
$center
[1] 0.25 3.00

$radius
[1] 2.015564

Example circle plot

By leaving the function intact, you can calculate the inverse radius for any set of points you want. I agree that the inverse radius simply means 1/radius.

Answer 2

Here's a geometric approach. Suppose I have three random points in a data frame:

set.seed(1)

df <- setNames(as.data.frame(matrix(rnorm(6), nrow = 3)), c("x", "y"))
df
#>            x          y
#> 1 -0.6264538  1.5952808
#> 2  0.1836433  0.3295078
#> 3 -0.8356286 -0.8204684

plot(df$x, df$y, xlim = c(-3, 2), ylim = c(-2, 2))

Now, I can draw lines between these points and find the mid-point arithmetically:

lines(df$x, df$y)

mid_df <- data.frame(x = diff(df$x)/2 + df$x[-3],
                     y = diff(df$y)/2 + df$y[-3],
                     slope = -diff(df$x)/diff(df$y))
mid_df$intercept <- mid_df$y - mid_df$slope * mid_df$x

points(mid_df$x, mid_df$y)

If I draw lines perpendicular to these lines through the midpoint, then the resulting point should be equidistant from my three starting points:

abline(a = mid_df$intercept[1], b = mid_df$slope[1], col = "red", lty = 2)
abline(a = mid_df$intercept[2], b = mid_df$slope[2], col = "red", lty = 2)

center_x <- (mid_df$intercept[2] - mid_df$intercept[1]) /
            (mid_df$slope[1] - mid_df$slope[2])

center_y <- mid_df$slope[1] * center_x + mid_df$intercept[1]

points(center_x, center_y)

As is indeed the case:

distances <- sqrt((center_x - df$x)^2 + (center_y - df$y)^2)

distances
#> [1] 1.136489 1.136489 1.136489

So, the radius of the circle is given by distances[1], and its center is at center_x, center_y. The curvature that is your end result is given by 1/distances[1]

To prove this, let's draw the circle this describes:

xvals <- seq(center_x - distances[1], center_x + distances[1], length.out = 100)
yvals <- center_y + sqrt(distances[1]^2 - (xvals - center_x)^2)
yvals <- c(yvals, center_y - sqrt(distances[1]^2 - (xvals - center_x)^2))
xvals <- c(xvals, rev(xvals))
lines(xvals, yvals)

Answer 3

My favorite resolution:

• subtract the coordinates of one point from the two others;

• now your circle is through the origin and has the simplified equation

  2 Xc X + 2 Yc Y = X² + Y²
  

• you have a standard and easy system of two equations in two unknowns.

  X1 Xc + Y1 Yc = (X1² + Y1²) / 2 = Z1
  X2 Xc + Y2 Yc = (X2² + Y2²) / 2 = Z2
  

• when you have computed Xc and Yc, the radius is √Xc²+Yc².

Answer 4

Using complex numbers:

We map the points Z1, Z2 to -1 and 1 by the transformation Z = (2Z - Z1 - Z2) / (Z2 - Z1). Now the center of the circle is on the imaginary axis, let iH. We express that the center is equidistant to 1 and to the third point (2 Z3 - Z0 - Z1) / (Z1 - Z0) = X + iY,

H² + 1 = X² + (Y - H)²

or

H = (X² + Y² - 1) / 2Y

and

R = √H²+1.