计算 R 中封闭对象的曲率

Question

我有一个由 1,000 个点组成的多边形。是否可以计算每个点的曲率？多边形原本只包含 13 个点：

43748.72 40714.19
43743.99 40716.16
43741.36 40720.19
43740.95 40726.46
43742.67 40729.28
43745.52 40730.97
43748.72 40731.14
43752.86 40729.43
43756.77 40723.24
43757.19 40719.73
43755.27 40716.68
43752.23 40714.76
43748.72 40714.19

然后我使用 smoothr 包中的 smooth 函数进行插值，因为多边形有 1,000 个点，看起来像：现在我想计算每个点的曲率。但是由于这是一个封闭的对象，实际如何进行计算呢？

编辑我终于找到了一个带有突起的电池来测试坚固性。细胞看起来像：

而对应的K值为：

确实，这个图捕捉到了两个突起，但是曲率值可以那么高吗？我读了一篇论文，似乎它们的值都在 1 以内：论文 link: https://www.biorxiv.org/content/10.1101/623793v1.full

Answer 1

您的示例本身无法完全重现，但可以参考 :

library(sf)
library(smoothr)
library(ggplot2)

data <- structure(list(x = c(43740.95, 43741.36, 43742.67, 43743.99, 
           43745.52, 43748.72, 43748.72, 43748.72, 43752.23, 43752.86, 43755.27, 
           43756.77, 43757.19), y = c(40726.46, 40720.19, 40729.28, 40716.16, 
           40730.97, 40714.19, 40731.14, 40714.19, 40714.76, 40729.43, 40716.68, 
           40723.24, 40719.73)), class = "data.frame", row.names = c(NA,  -13L))

smooth_poly <- data %>% 
  st_as_sf(coords=c("x", "y")) %>% 
  st_union() %>% 
  st_convex_hull() %>% 
  smooth(method='spline', n=1000)

smooth_df <- as.data.frame(sf::st_coordinates(smooth_poly))

ggplot(smooth_df, aes(X, Y)) + 
  geom_polygon(alpha = 0, colour = "black", size = 1) +
  coord_equal()

现在我们在名为 smooth_df 的数据框中获得了平滑多边形的所有 X 和 Y 坐标。我们可以这样计算曲率向量的 x 和 y 分量：

dx <- diff(c(smooth_df$X, smooth_df$X[1])) # Distance between x coords with wrap-around
dy <- diff(c(smooth_df$Y, smooth_df$Y[1])) # Distance between y coords with wrap-around
ds <- sqrt(dx^2 + dy^2)                    # Segment size between points
ddx <- dx/ds                               # Ratio of x distance to segment size
ddy <- dy/ds                               # Ratio of y distance to segment size
ds2 <- (ds + c(ds[-1], ds[1]))/2           # Mean segment length either side per point
smooth_df$Cx <- diff(c(ddx, ddx[1]))/ds2   # Change in ddx per unit length
smooth_df$Cy <- diff(c(ddy, ddy[1]))/ds2   # Change in ddy per unit length

这最后两个是多边形外围每个点的曲率向量的 x 和 y 分量。由于这个多边形是光滑的，所以曲率很小：

head(smooth_df)
#>          X        Y L1 L2         Cx        Cy
#> 1 43748.72 40714.19  1  1 0.02288753 0.1419567
#> 2 43748.67 40714.20  1  1 0.02324771 0.1375075
#> 3 43748.61 40714.21  1  1 0.02356064 0.1332985
#> 4 43748.56 40714.22  1  1 0.02383216 0.1293156
#> 5 43748.51 40714.23  1  1 0.02406747 0.1255458
#> 6 43748.45 40714.24  1  1 0.02427127 0.1219768

将这些矢量添加到绘图中只会给多边形的内部一些 "fur"，因为它们太多而且太小，所以我们可以通过以下方式证明方向是正确的用箭头绘制其中的一个子集，放大 10 倍。箭头应从外围开始并直接指向多边形在该点的凹面方向。我们还应该在曲线较紧的地方看到较长的箭头，在多边形平坦的地方看到较短的箭头。

smooth_df$Cx_plot <- 10 * smooth_df$Cx + smooth_df$X
smooth_df$Cy_plot <- 10 * smooth_df$Cy + smooth_df$Y

ggplot(smooth_df, aes(X, Y)) + 
  geom_polygon(alpha = 0, colour = "black", size = 1) +
  geom_segment(data = smooth_df[seq(1, nrow(smooth_df), 50),],
               mapping = aes(xend = Cx_plot, yend = Cy_plot), 
               arrow =  arrow(length = unit(0.3, "cm"))) +
  coord_equal()

如果你想把曲率作为一维数，你可以这样做：

smooth_df$K <- (ddy * smooth_df$Cx - ddx * smooth_df$Cy)/
               ((ddx^2 + ddy^2)^(3/2))

然后您可以将曲率绘制为一种颜色。当曲线向外凹时，这也会给出负值，尽管我在这里刚刚再次绘制了凸包。红色表示高曲率区域，蓝色区域较平坦。

ggplot(smooth_df, aes(X, Y)) + 
  geom_point(aes(colour = K)) +
  coord_equal() + scale_colour_gradient(low = "skyblue", high = "red")

计算 R 中封闭对象的曲率

Calculate curvature of a closed object in R

r

polygon

computational-geometry

curve