annulus_demo

library(ashapesampler)
library(alphahull)
library(ggplot2)
library(doParallel)
#> Loading required package: foreach
#> Loading required package: iterators
#> Loading required package: parallel
library(parallel)
cores <- min(2L, detectCores())

In this document, we demonstrate the α-shape sampler pipeline by simulating the process of learning a set of two-dimensional shapes (in this case, annuli) and simulating a new shape from that. This vignette requires the packages alphahull, ggplot2, parallel, and doParallel in addition to ashapesampler.

We begin by setting the parameters for our simulation. We will fix α = 0.25 and n = 100, and draw 20 shapes for our data set. Our true underlying manifold will be the annulus with outer radius 0.75 and inner radius 0.25.

set.seed(201723)
my_alpha = 0.15
n = 1000
N= 10
r_maj <- 0.75
r_min <- 0.25

Next we will draw the shapes themselves.

ann_list <- list()
complex_list <- list()
tau_vec <- vector("numeric", N)
for (k in 1:N){
  ann_pts <- runif_annulus(n, r_maj, r_min)
  ann_list[[k]] <- ashape(ann_pts, alpha = my_alpha)
  complex_list[[k]] <- get_alpha_complex(ann_pts, my_alpha)
  tau_vec[k] <- tau_bound(ann_list[[k]]$x, complex_list[[k]])
}

Now that we have the shapes generated and imported, we want to sample point clouds to combine. We’ll choose 2.

choose_2 <- sample(N,2)
point_cloud = rbind(ann_list[[choose_2[1]]]$x, ann_list[[choose_2[[2]]]]$x)

Then we will have our τ bound be a summary statistic of the τ found for each input shape. Here, we will use mean, but one can tweak this to see different results. Note that if τ is too small, then the random walk won’t be able to execute around the point cloud, but if τ is too big, then we risk losing geometric and topological information in the reconstruction.

tau_vec2 = c(tau_vec[choose_2[1]], tau_vec[choose_2[2]])

Now we can take the parameters and generate a new shape and plot it. Here, we assume k_min=2 as we are in two dimensions.

new_annulus <- generate_ashape2d(point_cloud, J=2, tau=min(tau_vec2), 
                                 cores=cores)
#> [1] "Acceptance Rate is 0.9371"


tri_keep = new_annulus$delvor.obj$tri.obj$trlist[which(new_annulus$delvor.obj$tri.obj$cclist[,3]<new_annulus$alpha), 1:3]
dim_tri = dim(tri_keep)[1]
tri_keep = as.vector(t(tri_keep))
triangles = data.frame("id"=sort(rep(1:dim_tri, 3)), 
                       "x"=new_annulus$x[tri_keep, 1], 
                       "y"=new_annulus$x[tri_keep,2])
extremes = as.data.frame(new_annulus$x[new_annulus$alpha.extremes,])

edges = as.data.frame(new_annulus$edges[,3:6])

ggplot(data.frame(new_annulus$x), aes(x=X1, y=X2)) +
  geom_polygon(data=triangles, aes(x=x, y=y, group=id), fill="gray") +
  geom_segment(data=edges, aes(x=x1, y=y1, xend=x2, yend=y2), color="blue")+
  geom_point(data=extremes, aes(x=V1, y=V2), size=1.5)+
  theme_classic()