Search code examples
rpolygoncomputational-geometrycurve

Calculate curvature of a closed object in R

I have a polygon that consists of 1,000 points. Is it possible to calculate the curvature at each point? The polygon originally contains only 13 points:

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

Then I use the smooth function in smoothr package for interpolation now that the polygon has 1,000 points and looks like: enter image description here And now I want to calculate curvature at each point. But since this is a closed object, how to actually perform the calculation?

enter image description here

EDIT I finally found a cell with protrusions to test the robustness. The cell looks like: enter image description here

And the corresponding K values are: enter image description here

Indeed, this plot captures two protrusions but can the curvature value be that high? I read a paper and seems like their values are all within 1: enter image description here paper link: https://www.biorxiv.org/content/10.1101/623793v1.full

Solution

  • Your example is not fully reproducible on its own, though it can be made so with reference to 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()

    enter image description here

    Now we have all the X and Y co-ordinates of the smoothed polygon in a data frame called smooth_df. We can calculate the x and y components of the curvature vectors like this:

    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

    These last two are the x and y components of the curvature vectors at each point on the periphery of the polygon. Since this polygon is smooth, the curvatures are small:

    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

    Adding these vectors to a plot would just give the inside of the polygon some "fur", since there are so many of them and they are so small, so instead we can show that the directions are correct by plotting a subset of them, magnified by 10, with arrowheads. The arrows should start on the periphery and point directly in the direction of the concavity of the polygon at that point. We should also see longer arrows where the curves are tight, and shorter arrows where the polygon is flat.

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

    enter image description here

    If you want the curvature as a single dimensional number 𝜿, you can do:

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

    Which then allows you to plot the curvature as a colour. This will also give negative values when the curve is concave outwards, though I have here just plotted the convex hull again. The red indicates areas with high curvature, the blue areas are flatter.

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

    enter image description here