计算R中圆半径的倒数

Calculating the inverse of the radius of a circle in R

我正在尝试复制一项研究中使用的函数,但我并没有真正的数学背景来完全理解应该如何完成这项工作。该测量从舌头轮廓上取三个点,并使用这三个点来计算通过它们的圆的半径。我查看了 并在 python 中找到了这样做的东西。我试图修改代码,以便它可以在 R 中使用我自己的数据。 (贴在底部)

问题是,根据我正在阅读的研究,我需要计算圆的圆周的凹度,并求出通过这三个点的圆的半径的倒数。我在谷歌上搜索,但老实说,这对我来说毫无意义。我唯一发现的是,我似乎需要计算舌面曲线的一阶和二阶导数。我真的希望有人能够帮助探索我如何在 R 中做到这一点。老实说,我对理解这里的数学并不太感兴趣,只是对如何实际实现它感兴趣。

编辑:我认为下面是我需要复制的公式。正如 MBo 指出的那样,情况并非如此。

我将重复另一项研究中的内容,该研究使用了非常非常相似的方法以防有帮助。

'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.'三点集'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'。这是我需要做的,但不知道从哪里开始在 R 中操作它。

下面的代码是 python 我试图在 R 中复制的代码,用于从三个点获取半径。之后我不知道如何进行。

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)

这是我的 R 尝试。 我还没有编写该函数,但我将查看沿曲线的三个点 A、B 和 C。该函数将为这三个点中的每一个提取 x 和 y 值(称为 x_value_a,y_value_a 等)。一旦完成。我将 运行 下面的代码。在此之后,我被难倒了。

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)

如有任何帮助,我们将不胜感激。我为我的数学无知感到抱歉。

如果您只想将 Python 脚本翻译成 R,那非常简单(我不太明白为什么要在添加的 R 代码中拆分它)。

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

}

请注意,p1-3 表示包含 x- 和 y-coordinate 的向量。我必须相信这里的原始 Python 代码,但使用 desmos.com 的快速检查似乎表明它有效:

> 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

通过保持函数不变,您可以计算所需的任何一组点的反半径。我同意倒数半径只是表示 1/radius.

这是一个几何方法。假设我在一个数据框中有三个随机点:

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

现在,我可以在这些点之间画线并通过算术找到 mid-point:

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)

如果我通过中点绘制垂直于这些线的线,那么所得点应该与我的三个起点等距:

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)

确实如此:

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

distances
#> [1] 1.136489 1.136489 1.136489

所以,圆的半径由distances[1]给出,圆心在center_x, center_y。最终结果的曲率由 1/distances[1]

给出

为了证明这一点,让我们画出这个描述的圆圈:

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)

我最喜欢的分辨率:

  • 从另外两个点减去一个点的坐标;

  • 现在你的圆通过原点并且有简化的方程

    2 Xc X + 2 Yc Y = X² + Y²
    
  • 你有一个标准且简单的系统,包含两个未知数的两个方程。

    X1 Xc + Y1 Yc = (X1² + Y1²) / 2 = Z1
    X2 Xc + Y2 Yc = (X2² + Y2²) / 2 = Z2
    
  • 当您计算出 XcYc 时,半径为 √Xc²+Yc².

使用复数:

我们通过变换Z = (2Z - Z1 - Z2) / (Z2 - Z1)将点Z1Z2映射到-11。现在圆心在虚轴上,令iH。我们表示中心与1和第三点(2 Z3 - Z0 - Z1) / (Z1 - Z0) = X + iY,

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

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

R = √H²+1.