How to build the alpha shape complex for a fixed radius with the TDA package?

Question

The R package TDA (https://cran.r-project.org/package=TDA) is used to do topological data analysis.

I'm interested on the reconstruction of alpha shapes complexes given a fixed radius (estimated from a persistence barcode)

With this in mind I have tried this reconstructions with this toy dataset

library(TDA)

n <- 100
theta <- runif(n, 0, 2 * pi)
r1 <- sqrt(runif(4 * n)) + 1
r2 <- sqrt(runif(n)) * 0.5 + 0.5

X1c1 <- r1 * cos(theta)
Yc1 <- r1 * sin(theta)
X1c2 <- r2 * cos(theta) + 2
Yc2 <- r2 * sin(theta) + 2
X1 <- c(X1c1, X1c2)

X <- cbind(c(Yc1, Yc2), c(X1c1, X1c2))

plot(X)

The figure is formed of two circles of different diameters.

Then, I create the alpha complexes with TDA

alphaDiag <-   alphaComplexDiag(
  X = X,
  library = c("GUDHI", "Dionysus"),
  location = TRUE,
  maxdimension = 1,
  printProgress = TRUE
)

alphaFiltration <- alphaComplexFiltration(
  X = X, 
  printProgress = TRUE)

Plotting the persistence diagram I get

plot(alphaDiag$diagram)

Persistence:

And the complex reconstruction

lim <- rep(c(-2, 3), 2)
plot(
  NULL,
  type = "n",
  xlim = lim[1:2],
  ylim = lim[3:4],
  main = "Alpha Complex Filtration Plot", 
  asp=1
)
for (idx in seq(along = alphaFiltration[["cmplx"]])) {
  polygon(
    alphaFiltration[["coordinates"]][alphaFiltration[["cmplx"]][[idx]], , drop = FALSE],
    col = "pink",
    border = NA,
    xlim = lim[1:2],
    ylim = lim[3:4]
  )
}
for (idx in seq(along = alphaFiltration[["cmplx"]])) {
  polygon(alphaFiltration[["coordinates"]][alphaFiltration[["cmplx"]][[idx]], , drop = FALSE],
          col = NULL,
          xlim = lim[1:2],
          ylim = lim[3:4])
}
points(alphaFiltration[["coordinates"]], pch = 16)

Complex:

As you can see, the alpha Complex is not capturing well the two holes. I would like to adjust the radius (similarly to a Vietoris o Cech Rips) to create a better figure.

Do you have any idea or pointers for this problem?

Answer 1

Perhaps something like this?

library("ggplot2")

filtration <- alphaComplexFiltration(X)
alpha <- 0.3

simplices_idx <- which(filtration$values <= alpha)
edges_idx <- which(sapply(filtration$cmplx[simplices_idx], function(s) length(s) == 2))
triangles_idx <- which(sapply(filtration$cmplx[simplices_idx], function(s) length(s) == 3))

edges <- matrix(do.call("rbind", filtration$cmplx[simplices_idx[edges_idx]]), ncol=2)
edges <- data.frame(cbind(X[edges[,1],], X[edges[,2],]))
colnames(edges) <- c("x1", "y1", "x2", "y2")

triangles_vertex_idx <- unlist(filtration$cmplx[simplices_idx[triangles_idx]])
triangles_grouping <- rep(1:length(triangles_idx), each=3)
triangles <- data.frame(cbind(triangles_grouping, X[triangles_vertex_idx,]))
colnames(triangles) <- c("id", "x", "y")

ggplot(data.frame(X), aes(x=X1, y=X2)) +
  geom_polygon(data=triangles, aes(x=x, y=y, group=id), fill="brown") +
  geom_segment(data=edges, aes(x=x1, y=y1, xend=x2, yend=y2), color="black") +
  geom_point(size=1.5)