计算R中圆的半径的倒数

问题描述

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

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

解决方法

如果只想将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 /半径。

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

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,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],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
```
计算了Xc和Yc时，半径为√Xc²+Yc²。

使用复数：

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

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

或

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

和

R = √H²+1.

curve geometry python r r