# persistence_script.R # Peter Bubenik, January 2017 # Load triangulated figures, filter, visualize, calculate persistence diagrams # and persistence landscapes # read persistence landscape files, assemble them, and save to R data file pl0_list <- list() pl1_list <- list() pose <- 1 # pose number between 1 and num_figures (148) for (pose in 1:num_figures){ # for each pose # read the data file <- paste(figures[pose,1],figures[pose,2],sep="") print(file) vertices <- as.matrix(read.table(paste(file.path(data_directory,file),".vert",sep=""))) triangles <- as.matrix(read.table(paste(file.path(data_directory,file),".tri",sep=""))) num_triangles <- dim(triangles)[1] mymesh <- tmesh3d(t(vertices), t(triangles), homogeneous = FALSE) if (show_individual_figures == TRUE){ open3d(); wire3d(mymesh) } # filter the triangulation # perseus will construct a filtered simplicial complex from values on the top dimensional faces # here we give the vertices values and use those to give the triangles values vertex_values <- knn.dist(vertices,k=10)[,10] # filter by k-nearest-neighbor triangle_values <- vector(length = num_triangles) for (i in 1:num_triangles) triangle_values[i] <- mean(vertex_values[triangles[i,]]) # assign a triangle the average value of its vertices if (show_individual_figures == TRUE) draw_filtered_triangulation(triangles,vertices,triangle_values) #construct animated gif - need imageMagick # dir.create(file.path(main_directory, animation_directory)) # comment this line after running once # setwd(animation_directory) # vertices <- -1 + 2*(vertices-min(vertices))/(max(vertices)-min(vertices)) # rescale to [-1,1]. If this is inside the function, then rgl doesn't center on the figure for some reason. # animate_filtered_triangulation(triangles,vertices,triangle_values,num_frames,file,animation_delay) # setwd("..") #use perseus to calculate persistent homology setwd(perseus_directory) ph <- pers_hom_filtered_triangulation(triangles,triangle_values,grid_size,file,operating_system) setwd("..") if (show_individual_figures == TRUE){ plot_persistence_diagram(ph[[1]],paste("PH in degree 0 of",file)) plot_persistence_diagram(ph[[2]],paste("PH in degree 1 of",file)) } #use Persistence Landscape Toolbox to calculate persistence landscape setwd(pl_directory) persistence_landscape(ph,file,max_filtration_value,max_depth,degrees=0:1) if (show_individual_figures == TRUE){ plot_persistence_landscape_from_files(file,degree=0,max_depth) plot_persistence_landscape_from_files(file,degree=1,max_depth) } pl0 <- persistence_landscape_vector(paste(file,"_0", sep = ""),pl_param_vals, max_depth) pl1 <- persistence_landscape_vector(paste(file,"_1", sep = ""),pl_param_vals, max_depth) # plot_persistence_landscape_from_vector(pl0,length(pl_param_vals),"PL in degree 0") # plot_persistence_landscape_from_vector(pl1,length(pl_param_vals),"PL in degree 1") pl0_list <- c(pl0_list,list(pl0)) pl1_list <- c(pl1_list,list(pl1)) setwd("..") } plm0 <- matrix_from_list_of_vectors(pl0_list) plm1 <- matrix_from_list_of_vectors(pl1_list) rm(pl0_list,pl1_list) plot_persistence_landscape_from_vector(colMeans(plm0),length(pl_param_vals),"Average PL in degree 0") plot_persistence_landscape_from_vector(colMeans(plm1),length(pl_param_vals),"Average PL in degree 1") plm_list <- list(plm0,plm1) rm(plm0,plm1) #save to R Data file save(figures,num_figures,plm_list,grid_size,max_depth,pl_param_vals,file="nonrigid3d_landscape.RData")