计算R中圆的半径的倒数答案

【问题标题】：Calculating the inverse of the radius of a circle in R计算R中圆的半径的倒数
【发布时间】：2020-09-17 06:07:03
【问题描述】：

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

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

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

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

'任何三个点（A、B、C）都可以被认为是位于一个圆的圆周上。圆将有一个半径，其倒数表示通过这三个点的圆的曲率。这三个点的集合产生一个曲率数，它是通过它们的圆的半径的倒数。位于一条直线上的三个点的曲率为零，因为它们的凹度为零，这成为曲率方程的分子。这是我需要做的，但不知道从哪里开始在 R 中操作它。

下面的代码是我试图在 R 中复制的 python 代码，以获取三个点的半径。我不知道之后该怎么做。

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)

任何帮助将不胜感激。对不起，我对数学的无知。

【问题讨论】：

C 公式（分析曲线曲率）与您的问题无关。
我的第一反应是“逆”的意思是“倒数”——所以你只需要1 / radius。这将与凹度为 0 的直线相匹配，因为半径是无限的。
顺便说一句-您似乎在 R 代码中缺少 - temp 来计算 bc
@MBo，感谢您指出这一点。在我正在阅读的一篇论文中，曲率测量似乎是必不可少的。我想使用的度量是基于那个度量的，所以我认为计算在这里也是必不可少的。很高兴犯错，不必追求那条路。谢谢！
@Hobo，感谢您的两位 cmets。我已经更新了代码以包含缺失值（好地方！）。我希望你对 1/半径是正确的。

标签： python r geometry curve

【解决方案1】：

如果您只想将 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 坐标的向量。我必须相信这里的原始 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/半径。

【讨论】：

感谢您的回复，伯特。您对代码的翻译以匹配原始代码是有道理的。我只是手动将所有 x 和 y 值分开，所以它非常难看，可能需要更多解释。我也没有 if 和 else 但我可能应该包括它，以防万一，所以我肯定会从你的代码中借用。不幸的是，主要问题不是 python 代码到 R 代码的翻译，所以我不能投票给这个作为答案，但我非常感谢你的回复。

【解决方案2】：

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

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

现在，我可以在这些点之间画线并用算术方法找到中点：

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)

【讨论】：

感谢您如此彻底的回复。看到它以这种方式布置真的很美。
可爱的解释

【解决方案3】：

我最喜欢的解决方案：

从另外两个中减去一个点的坐标；
现在你的圆通过原点并且有简化的方程
```
2 Xc X + 2 Yc Y = X² + Y²
```

你有一个标准且简单的系统，包含两个未知数的两个方程。

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

当您计算出Xc 和Yc 时，半径为√Xc²+Yc²。

【讨论】：

【解决方案4】：

使用复数：

我们通过变换Z = (2Z - Z1 - Z2) / (Z2 - Z1) 将点Z1、Z2 映射到-1 和1。现在圆的中心在虚轴上，设iH。我们表示中心到1和第三点(2 Z3 - Z0 - Z1) / (Z1 - Z0) = X + iY等距离，

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

或

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

和

R = √H²+1.

【讨论】：