计算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

然后，现在多边形具有1,000个点，看起来像：现在，我想计算每个点的曲率。但是由于这是一个封闭的对象，如何实际执行计算？

Answer 1

尽管您的示例可以完全复制，但是可以参考your previous question进行复制：

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

现在，在名为X的数据框中，具有平滑多边形的所有Y和smooth_df坐标。我们可以像这样计算曲率向量的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

将这些向量添加到绘图中只会给多边形内部带来一些“毛发”，因为它们太多了，而且它们是如此之小，所以我们可以通过绘制它们的一个子集来表明方向是正确的，被箭头放大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()