python - 计算R中圆半径的倒数

Question

我正在尝试复制一项研究中使用的函数，但我并没有真正的数学背景来完全理解应该如何完成这项工作。该测量从舌头轮廓上取三个点，并使用这三个点来计算通过它们的圆的半径。我看过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)

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

Answer 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/radius。