I'm working on a program in R to calculate the Gabriel Graph for up to 1000 data points. I used a program I found online first (GabrielGraph based on Bhattacharya et al. 1981 lines 781-830).
Unfortunately it takes quite a bit of time to get the result so I tried reprogramming it using Rcpp. For this I wrote a couple of small programs and a big one called edges which is used to calculate the edges of the Gabriel Graph. I'm also new to programming with Rcpp so I probably did everything more complicated than necessary but I didn't know how to do it any better.
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double vecnorm(NumericVector x){
  //to calculate the vectornorm sqrt(sum of (vector entries)^2)
  double out;
  out = sqrt(sum(pow(x,2.0)));
  return out;
}
// [[Rcpp::export]]
NumericVector vektorzugriff(NumericMatrix  xy,int i){
  //to return a row of the Matrix xy
  int col = xy.ncol();
  NumericVector out(col);
  for(int j=0; j<=col; j++){
    out[j] = xy(i-1,j);
  }
  return out;
}
// [[Rcpp::export]]
IntegerVector vergl(NumericVector eins, NumericVector zwei){
  //to see if two Vectors have any identical entries
  IntegerVector out = match(eins, zwei);
  return out;
}
// [[Rcpp::export]]
IntegerVector verglInt(int eins, NumericVector zwei){
  NumericVector dummy =  NumericVector::create( eins ) ;
  IntegerVector out = match(dummy, zwei);
  return out;
}
// [[Rcpp::export]]
NumericVector toVec(NumericVector excluded, int k){
  //to append int k to a Vector excluded
  NumericVector dummy =  NumericVector::create( k ) ;
  int len = excluded.size();
  int len2 = dummy.size();
  int i=0;
  NumericVector out(len+len2);
  while(i<len+len2){
    if(i<len){
      out[i]=excluded[i];
      i++;
    }
    else{
      out[i]=dummy[i-len];
      i++;
    }
  }
  return out;
}
// [[Rcpp::export]]
LogicalVector isNA(IntegerVector x) {
  //to see which Vector Entries are NAs
  int n = x.size();
  LogicalVector out(n);  
  for (int i = 0; i < n; ++i) {
    out[i] = IntegerVector::is_na(x[i]);
  }
  return out;
}
// [[Rcpp::export]]
NumericMatrix Gab(NumericMatrix Gabriel, NumericVector edges1, NumericVector edges2, int anz){
  //to fill a Matrix with the Gabrieledges
  for(int i=0; i<anz; i++) {
    Gabriel(edges1[i]-1, edges2[i]-1)  = 1 ; 
    Gabriel(edges2[i]-1, edges1[i]-1)  = 1 ; 
  }
  return Gabriel;
}
// [[Rcpp::export]]
NumericVector edges(NumericMatrix  xy,NumericVector vertices,NumericVector excluded, int i){
  //actual function to calculate the edges of the GabrielGraph
  int npts = xy.nrow()+1;
  double d1; 
  double d2;
  double d3;
  for(int r=i+1; r<npts; r++) {
    // Skip vertices in excluded
    if(!is_true(any(isNA(verglInt(r,excluded))))){
      continue;}
    d1 = vecnorm(vektorzugriff(xy,i) - vektorzugriff(xy,r));
    for(int k=1; k<npts; k++) {
      if((k!=r) && (k!=i)){
        d2 = vecnorm(vektorzugriff(xy,i) - vektorzugriff(xy,k));
        d3 = vecnorm(vektorzugriff(xy,r) - vektorzugriff(xy,k));
        //Betrachte vertices, die noch nicht excluded sind
        if(!is_true(any(isNA(verglInt(k,vertices[isNA(vergl(vertices,excluded))]))))){
          //Wenn d(x,z)^2 > d(x,y)^2+d(y,z)^2 -> Kante gehoert nicht zum GG
          if( pow(d2,2.0) > pow(d1,2.0) + pow(d3,2.0) ) {
            excluded = toVec(excluded,k);
          }
        }
        if( pow(d1,2.0) > pow(d2,2.0) + pow(d3,2.0) ){
          excluded = toVec(excluded,r);
          break;
        }
      }
    }
  }
  return excluded;
}
I used these Rcpp programs in this R program:
GabrielGraphMatrix <- function(X,Y,PlotIt=FALSE){
# Heuristic rejection Algorithm for Gabriel Graph Construction (Bhattacharya et al. 1981)
# Algorithm is ~ O(d n^2)
  #loading Rcpp functions
  library(Rcpp)
  sourceCpp("... .cpp")
  XY <- cbind(X,Y)
  ndim <- ncol(XY)
  npts <- nrow(XY)
  edges1<- c()
  edges2<- c()
  for( i in 1:(npts-1) ) {
    # Candidate set of Gabriel neighbors
    vertices <- (i+1):npts
    # Initialize list of vertices to be excluded from Ni
    excluded <- edges(XY,vertices,vector(),i);
    adj <- vertices[which(!match(vertices,excluded,nomatch=F)>0)]
    if(length(adj) > 0) {
      edges1=c(edges1,rep(i,length(adj)))
      edges2=c(edges2,adj)  
    }
  }
  anz <- length(edges1)
  Gabriel <- Gab(matrix(0, npts, npts),edges1,edges2,anz)
  return(list(Gabriel=Gabriel,edges=cbind(edges1,edges2)))
}
For a sample data of ten data points it worked fine, for example:
z <- 10
X <- runif(z)*100
Y <- runif(z)*100
GabrielGraphMatrix(X,Y)
returns
> GabrielGraphMatrix(X,Y)
$Gabriel
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    1    0    0    0    0    0    0    0     0
 [2,]    1    0    0    1    0    0    1    0    0     0
 [3,]    0    0    0    1    1    0    0    0    0     1
 [4,]    0    1    1    0    0    0    0    0    0     0
 [5,]    0    0    1    0    0    0    0    0    0     0
 [6,]    0    0    0    0    0    0    0    1    0     0
 [7,]    0    1    0    0    0    0    0    0    0     0
 [8,]    0    0    0    0    0    1    0    0    1     1
 [9,]    0    0    0    0    0    0    0    1    0     1
[10,]    0    0    1    0    0    0    0    1    1     0
$edges
      edges1 edges2
 [1,]      1      2
 [2,]      2      4
 [3,]      2      7
 [4,]      3      4
 [5,]      3      5
 [6,]      3     10
 [7,]      6      8
 [8,]      8      9
 [9,]      8     10
[10,]      9     10
But if I try to put in bigger data sets I get this error message:
Error: Value of SET_STRING_ELT() must be a 'CHARSXP' not a 'builtin'
I would be amazingly grateful if anybody had at least an idea of what I did wrong.
 
    