View Shiny app.
Get the code.

We start by downloading the data for each bone. (This command reads the space delimited file from a website into R, removes the header content, and accounts for the fact that the coordinate data, and the data to reconstruct the surface of the bone are of different dimensions.)

bone1 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/001.off",sep = " ", skip = 2,fill = TRUE)
bone2 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/002.off",sep = " ", skip = 2,fill = TRUE)
bone3 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/003.off",sep = " ", skip = 2,fill = TRUE)
bone4 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/004.off",sep = " ", skip = 2,fill = TRUE)
bone5 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/005.off",sep = " ", skip = 2,fill = TRUE)

Next, we extract the coordinate data of each bone from .

bone1 = bone1[1:5128,]
bone2 = bone2[1:5086,]
bone3 = bone3[1:5091,]
bone4 = bone4[1:5092,]
bone5 = bone5[1:5154,]

Now, we plot the each point cloud.

Finally, we compute persistent homology using the Vietoris-Rips complex.

for (i in 1:5) {
     cat(paste0("Persistent Homology for Bone ", i,"\n"))
     #~ Compute Homology
     eval(parse(text = paste0("boneRips.",i," = ripsDiag(X = bone",i,",maxdimension = 1, maxscale = 0.015)")))
     #~ Plot Persistence Diagrams and Barcodes
     eval(parse(text = paste0("par(mfrow = c(2,1));plot(boneRips.",i,"[['diagram']],main = 'Bone ",i," Persistence Diagram');plot(boneRips.",i,"[['diagram']],barcode = TRUE,main = 'Bone ",i," Barcode Diagram')")))
}
## Persistent Homology for Bone 1

## Persistent Homology for Bone 2

## Persistent Homology for Bone 3

## Persistent Homology for Bone 4

## Persistent Homology for Bone 5

The persistence diagrams show a large number of holes. From my understanding of the structure we are studying, many of these holes may not be significant and complicating our understanding of the bone structure. To filter some of these features out, we evaluate homology over a grid.

#~ Evaluate over a grid
# Some constants used for evaluation
by = 0.01
buffer = 0.05
seqlen = 10
band = 0.01

# Limit of x coordinates
Xlim = max(bone1[,1],bone2[,1],bone3[,1],bone4[,1],bone5[,1]) + buffer
Xlim = c(min(bone1[,1],bone2[,1],bone3[,1],bone4[,1],bone5[,1]) - buffer,Xlim)
Xgrid = seq(Xlim[1],Xlim[2],by = by)
# y
Ylim = max(bone1[,2],bone2[,2],bone3[,2],bone4[,2],bone5[,2]) + buffer
Ylim = c(min(bone1[,2],bone2[,2],bone3[,2],bone4[,2],bone5[,2]) - buffer,Ylim)
Ygrid = seq(Ylim[1],Ylim[2],by = by)
# Z
Zlim = max(bone1[,3],bone2[,3],bone3[,3],bone4[,3],bone5[,3]) + buffer
Zlim = c(min(bone1[,3],bone2[,3],bone3[,3],bone4[,3],bone5[,3]) - buffer,Zlim)
Zgrid = seq(Zlim[1],Zlim[2],by = by)

# construct grid/mesh over the generated sequences
grid = expand.grid(Xgrid,Ygrid,Zgrid)

#~ Compute homology over a grid
for (i in 1:5) {
cat(paste0("Persistent Homology for Bone ", i,"\n"))
#~ Compute bootstrap bandwidth
eval(parse(text = paste0("bootband = bootstrapBand(X = bone",i,", FUN = kde, Grid = grid, B = 10, parallel = TRUE, alpha = 0.05, h = ",band,")")))
#~ Compute Homology over a grid
eval(parse(text = paste0("boneGrid.",i," = gridDiag(X = bone",i,", FUN = kde, h = ",band,", lim = cbind(Xlim,Ylim,Zlim), by = ",by,", sublevel = FALSE, library = 'Dionysus')")))
#~ Plot Persistence Diagrams and Barcodes; save plots as a pdf
eval(parse(text = paste0("par(mfrow = c(2,1));plot(boneGrid.",i,"[['diagram']], band = ",2*bootband$width,",main = 'Bone ",i," Persistence Diagram');plot(boneGrid.",i,"[['diagram']],barcode = TRUE,main = 'Bone ",i," Barcode Diagram')")))
}
## Persistent Homology for Bone 1

## Persistent Homology for Bone 2

## Persistent Homology for Bone 3

## Persistent Homology for Bone 4

## Persistent Homology for Bone 5