Merge pull request #18 from lcrawlab/new_functions

New functions

R/2D.R

@@ -184,7 +184,7 @@ euclid_dists_point_cloud_2D <- function(point, point_cloud){
   }
   dist_vec = vector("numeric",m)
   for (j in 1:m){
-    sqr_dist = (point[1]-point_cloud[j,1])^2+(point[2]-point_cloud[j,2])^2
+    sqr_dist = as.numeric((point[1]-point_cloud[j,1])^2+(point[2]-point_cloud[j,2])^2)
     dist_vec[j] = sqrt(sqr_dist)
   }
   return(dist_vec)

R/mcmc.R

@@ -68,7 +68,9 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
     .export = c("runif_ball_3D", "euclid_dists_point_cloud_3D")
   ) %dopar% {
   #for (i in 1:n_vert){
-    new_points = runif_ball_3D(m, tau/8)+rep(point_cloud[i,], each=m)
+    new_points = runif_ball_3D(m, tau/8)+ cbind(rep(point_cloud[i,1],m),
+                                                rep(point_cloud[i,2],m),
+                                                rep(point_cloud[i,3],m))
     keep_pts = matrix(NA, nrow=0, ncol=3)
     for (j in 1:m){
       dist_list = euclid_dists_point_cloud_3D(new_points[j,], point_cloud)
@@ -85,13 +87,13 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
     }
     keep_pts
   }
-
+  my_points = unique(my_points) #keeps error free if necessary.
   if(dim(my_points)[1]<5){
     stop("Not enough points accepted in MCMC walk to make a shape. Need at least 5.")
   }
   rr = dim(my_points)[1]/(m*dim(point_cloud)[1])
   print(paste0("Acceptance Rate is ", rr))
-  new_ashape <- alphashape3d::ashape3d(my_points, alpha=my_alpha)
+  new_ashape <- alphashape3d::ashape3d(my_points, alpha=tau-eps)
   return(new_ashape)
 }
 
@@ -163,7 +165,7 @@ generate_ashape2d <- function(point_cloud, J, tau, delta=0.05,
           .export = c("runif_disk", "euclid_dists_point_cloud_2D"))%dopar%{
 
   #for(i in 1:n_vert){
-    new_points = runif_disk(m, tau/8)+rep(point_cloud[i,], each=m)
+    new_points = runif_disk(m, tau/8)+ cbind(rep(point_cloud[i,1],m), rep(point_cloud[i,2],m))
     keep_pts = matrix(NA, nrow=0, ncol=2)
      for (j in 1:m){
        dist_list = euclid_dists_point_cloud_2D(new_points[j,], point_cloud)
@@ -180,11 +182,12 @@ generate_ashape2d <- function(point_cloud, J, tau, delta=0.05,
      }
     keep_pts
           }
+  my_points = unique(my_points) #keeps error free if necessary.
   if(dim(my_points)[1]<3){
     stop("Not enough points accepted in MCMC walk to make a shape. Need at least 3.")
   }
   rr = dim(my_points)[1]/(m*dim(point_cloud)[1])
   print(paste0("Acceptance Rate is ", rr))
-  new_ashape <- alphahull::ashape(my_points, alpha=my_alpha)
+  new_ashape <- alphahull::ashape(my_points, alpha=tau-eps)
   return(new_ashape)
 }

R/tau_bound.R

@@ -68,16 +68,16 @@ tau_bound <- function(v_list, complex, extremes=NULL, cores = 1){
     m=length(extremes)
   }
   dist_matrix = as.matrix(stats::dist(v_list))
-  e_list = extract_complex_edges(complex,m)
+  e_list = extract_complex_edges(complex,n)
   if(is.null(e_list)){
     return(min(dist_matrix[dist_matrix>0]))
   }
-  f_list = extract_complex_faces(complex,m)
+  f_list = extract_complex_faces(complex,n)
   f_circ = circumcenter_face(v_list, f_list)
   t_list = NULL
   t_circ = NULL
   if(dimension>2){
-    t_list = extract_complex_tet(complex,m)
+    t_list = extract_complex_tet(complex,n)
     t_circ = circumcenter_tet(v_list, t_list)
   }
   tau_vec=vector("numeric", m)
@@ -108,16 +108,16 @@ tau_bound <- function(v_list, complex, extremes=NULL, cores = 1){
       dist_vec = dist_vec_point[edge_list_zoom]
       dist_vec_b = c()
       if (dimension == 2){
-        if(!is.null(face_list_zoom)){
+        if(length(face_list_zoom)>0){
           points = matrix(f_circ[face_list_zoom,], ncol=2)
           dist_vec_b = c(dist_vec_b, 2*euclid_dists_point_cloud_2D(v_list[i,],
                                                               points ))
         }
       } else {
-        if(!is.null(face_list_zoom)){
+        if(length(face_list_zoom)>0){
           points = matrix(f_circ[face_list_zoom,], ncol=3)
           dist_vec_b = 2*euclid_dists_point_cloud_3D(v_list[i,],points)
-          if(!is.null(tet_list_zoom)){
+          if(length(tet_list_zoom)>0){
           points = matrix(t_circ[tet_list_zoom,], ncol=3)
           dist_vec_b = c(dist_vec_b, 2*euclid_dists_point_cloud_3D(v_list[i,], points))
           }
@@ -171,6 +171,7 @@ extreme_pts <- function(complex, n_vert, dimension){
     edge_face = data.frame(edge_face)
     colnames(edge_face)=c("ed1", "ed2")
     int_edge = edge_face[which(duplicated(edge_face)),]
+    int_edge = unique(int_edge)
     bd_edge = setdiff(edge_list, int_edge)
     bd_vert = unique(c(bd_edge$ed1, bd_edge$ed2))
     return(c(iso_vert, bd_vert))

README.md

@@ -1,2 +1,75 @@
 # ashapesampler
-R package for sampling alpha shapes
+
+R package for generating alpha shapes via either sampling from a known distribution or sampling based on an existing data set.
+
+## R Packages for ashapesampler and Tutorials
+
+The ashapesampler software requires the installation of the following R libraries:
+
+[alphahull](https://cran.r-project.org/web/packages/alphahull/index.html)
+
+[alphashape3d](https://cran.r-project.org/web/packages/alphashape3d/index.html)
+
+[doParallel](https://cran.r-project.org/web/packages/doParallel/index.html)
+
+[dplyr](https://cran.r-project.org/web/packages/dplyr/index.html)
+
+[foreach](https://cran.r-project.org/web/packages/foreach/index.html)
+
+[ggplot2](https://cran.r-project.org/web/packages/ggplot2/index.html)
+
+[Rvcg](https://cran.r-project.org/web/packages/Rvcg/index.html)
+
+[TDA](https://cran.r-project.org/web/packages/TDA/index.html)
+
+[truncnorm](https://cran.r-project.org/web/packages/truncnorm/index.html)
+
+Note that the latest BAKR and RATE functions are also included in the `Software` directory for this repo. Unless stated otherwise, the easiest method to install many of these packages is with the following example command entered in an R shell:
+
+    install.packages("alphahull", dependecies = TRUE)
+
+Alternatively, one can also [install R packages from the command line](http://cran.r-project.org/doc/manuals/r-release/R-admin.html#Installing-packages).
+
+## C++ Packages for ashapesampler and Tutorials
+
+The code in this repository assumes that basic C++ functions and applications are already set up on the running personal computer or cluster. If not, the functions and necessary TDA and alphashape3d packages to build alpha complexes and alpha shapes in three dimensions will not work properly. A simple option is to use [gcc](https://gcc.gnu.org/). macOS users may use this collection by installing the [Homebrew package manager](http://brew.sh/index.html) and then typing the following into the terminal:
+
+    brew install gcc
+
+For macOS users, the Xcode Command Line Tools include a GCC compiler. Instructions on how to install Xcode may be found [here](http://railsapps.github.io/xcode-command-line-tools.html). Additional installs for macOS users are automake, curl, glfw3, glew, xquartz, and qpdf. For extra tips on how to run C++ on macOS, please visit [here](http://seananderson.ca/2013/11/18/rcpp-mavericks.html). For tips on how to avoid errors dealing with "-lgfortran" or "-lquadmath", please visit [here](http://thecoatlessprofessor.com/programming/rcpp-rcpparmadillo-and-os-x-mavericks-lgfortran-and-lquadmath-error/).
+
+## R Package Installation
+
+To install the package, we will use the remotes package and run the command:
+
+	remotes::install_github('lcrawlab/ashapesampler')
+
+To load the package, use the command
+
+	library(ashapesampler)
+
+Other common installation procedures may apply.
+
+# Code Usage
+
+## Vignettes
+
+The `vignettes` folder contains the following demonstrations for running and analyzing results in the ashapesampler: 
+
+Generating new annuli from a simulated set of annuli
+
+Generating new tori from a simulated set of tori
+
+Converting binary masks to simplicial complexes for input into the alpha shape sampler
+
+Analyzing landmarks on teeth via procrustes analysis
+
+## Relevant Citations
+
+E.T. Winn-Nuñez, H. Witt, D. Bhaskar, R. Huang, I.Y. Wong, J. Reichner, and L. Crawford. A probabilistic method for sampling alpha-shapes.
+
+## Questions and Feedback
+
+Please send any questions or feedback to the corresponding authors [Emily Winn-Nuñez](mailto:[email protected]) or [Lorin Crawford](mailto:[email protected]).
+
+We appreciate any feedback you may have with our repository and instructions.