Anti-clustering II

We'll need the libraries again.

library(deldir) #For plotting the voronoi regions
library(randtoolbox) #For generating quasi-random points from low discrepancy sequences
library(MASS) #For multidimensional scaling (MDS)
library(geometry) #For 3D delaunay triangulation
library(rgl) #For 3D plot of delaunay triangulation

We'll generate the quasi-random points again and fix the colors. This time we will use various other quasi-random sequences to generate the points, such as the Hammersley and Sobol sequences.

#Generate some quasi-random points using the Sobol Sequence
qpoints <- sobol(100,dim=40,scrambling=3);dd<-deldir(qpoints[,1], qpoints[,2])

## 
##      PLEASE NOTE:  The components "delsgs" and "summary" of the 
##      object returned by deldir() are now DATA FRAMES rather than 
##      matrices (as they were prior to release 0.0-18). 
##      See help("deldir").
##  
##      PLEASE NOTE: The process that deldir() uses for determining
##      duplicated points has changed from that used in version
##      0.0-9 of this package (and previously). See help("deldir").

#Get 4 colors and repeat them to show the space filling properties of the points
ccc <- topo.colors(4)
ccc <- c(rep(ccc[1],25),rep(ccc[2],25),rep(ccc[3],25),rep(ccc[4],25))

Here we'll generate files that can be used offline for the creation of 3D voronoi regions.

qpointMat <- dist(qpoints)
qpointMDS <- isoMDS(qpointMat,k=3) #Create the 2D projection (the default)

## initial  value 30.168743 
## iter   5 value 25.481613
## final  value 25.111030 
## converged

#Scale the projection to be within the unit square. Add a fudge factor, epsilon, to make sure points don't lie on the [0,1] interval
epsilon <- 0.01
qmdspX <- (qpointMDS$points[,1] - min(qpointMDS$points[,1]) + epsilon)/(max(qpointMDS$points[,1]) - min(qpointMDS$points[,1]) + epsilon)

qmdspY <- (qpointMDS$points[,2] - min(qpointMDS$points[,2]) + epsilon)/(max(qpointMDS$points[,2]) - min(qpointMDS$points[,2]) + epsilon)

qmdspY <- (qpointMDS$points[,2] - min(qpointMDS$points[,2]) + epsilon)/(max(qpointMDS$points[,2]) - min(qpointMDS$points[,2]) + epsilon)

qmdspZ <- (qpointMDS$points[,3] - min(qpointMDS$points[,3]) + epsilon)/(max(qpointMDS$points[,3]) - min(qpointMDS$points[,3]) + epsilon)
#write 3D quasi-points 

write.table(cbind(qmdspX,qmdspY,qmdspZ),"PointsFile.txt",quote=FALSE,sep=" ",row.names=FALSE,col.names=FALSE)

Write out Sobol sequence quasi points.

write.table(qpoints[,1:3],"PointsFileQuasi.txt",quote=FALSE,sep=" ",row.names=FALSE,col.names=FALSE)

Here we show what a projection of 40D random uniform points looks like.

rpoints40D <- matrix(0,nrow=100,ncol=40)
for(i in 1:40){
  rpoints40D[,i] <- runif(100)
}
z <- deldir(rpoints40D[,1],rpoints40D[,2],rw=c(0,1,0,1))
w <- tile.list(z)
plot(w,fillcol=ccc,close=FALSE)

rpoints2DMat <- dist(rpoints40D[,1:2])

rpoints40DMat <- dist(rpoints40D)
rpoints40DMDS <- isoMDS(rpoints40DMat) #Create the 2D projection (the default)

## initial  value 41.097771 
## iter   5 value 36.237114
## final  value 35.283016 
## converged

#Scale the projection to be within the unit square. Add a fudge factor, epsilon, to make sure points don't lie on the [0,1] interval
epsilon <- 0.001
rmdspX <- (rpoints40DMDS$points[,1] - min(rpoints40DMDS$points[,1]) + epsilon)/(max(rpoints40DMDS$points[,1]) - min(rpoints40DMDS$points[,1]) + epsilon)

rmdspY <- (rpoints40DMDS$points[,2] - min(rpoints40DMDS$points[,2]) + epsilon)/(max(rpoints40DMDS$points[,2]) - min(rpoints40DMDS$points[,2]) + epsilon)

z <- deldir(rmdspX,rmdspY,rw=c(0,1,0,1))
w <- tile.list(z)
plot(w,fillcol=ccc,close=FALSE)

Lo! It seems to preserve more of the quasi-random feel, rather than the weak clustering of several dimensions of random uniform sampled points. But there are some anomalies. 1. Note how the corners and edges of the image have fewer points, and 2. some of the points lie nearly on the edge, though not in the corners. 3. The colors show that there is some ordering, but not a whole lot.

par(mfrow=c(1,2))
hist(rpoints2DMat,freq=FALSE)
lines(density(rpoints2DMat))

hist(rpoints40DMat,freq=FALSE)
lines(density(rpoints40DMat))

Hammersley sequence uses the same van der Corput-like sequence for each dimension with different bases, but the first dimension is uniformly partitioned.

hamqpoints <- halton(100,dim=40); 
hamqpoints[,2:40] <- hamqpoints[,1:39]
hamqpoints[,1] <- seq(from =0.01,to= 1.0, by =0.01)

z <- deldir(hamqpoints[,1],hamqpoints[,2],rw=c(0,1,0,1))
w <- tile.list(z)
plot(w,fillcol=ccc,close=FALSE)

hamqpoints2DMat <- dist(hamqpoints)

hamqpointsMat <- dist(hamqpoints)
hamqpointsMDS <- isoMDS(hamqpointsMat) 

## initial  value 20.327606 
## iter   5 value 16.618611
## final  value 16.109053 
## converged

epsilon <- 0.001
hamqmdspX <- (hamqpointsMDS$points[,1] - min(hamqpointsMDS$points[,1]) + epsilon)/(max(hamqpointsMDS$points[,1]) - min(hamqpointsMDS$points[,1]) + epsilon)

hamqmdspY <- (hamqpointsMDS$points[,2] - min(hamqpointsMDS$points[,2]) + epsilon)/(max(hamqpointsMDS$points[,2]) - min(hamqpointsMDS$points[,2]) + epsilon)

z <- deldir(hamqmdspX,hamqmdspY,rw=c(0,1,0,1))
w <- tile.list(z)
plot(w,fillcol=ccc,close=FALSE)

Write out Hammersley sequence quasi points.

write.table(hamqpoints[,1:3],"PointsFileHamQuasi.txt",quote=FALSE,sep=" ",row.names=FALSE,col.names=FALSE)

Plot Delaunay Tesselation with hamqpoints (this is the dual of the voronoi).

pc <- hamqpoints[,1:3]
tc <- delaunayn(pc)
rgl.viewpoint(60)
rgl.light(120,60)
tetramesh(tc,pc, alpha=0.4)

In principle, one should be able to generate the 3D voronoi regions with deldir, but I have not been successfull in doing so. Instead, I used the printed out 3D points above and slightly modified an example program (random_points.cc) in Voro++ (http://math.lbl.gov/voro++/) to show the 3D voronoi regions in gnuplot. random_points generated a random point set, but I substitue the read in quasi-random points, for the random points (call the new program, say, quasi_points) such that it generates two files quasi_points_v.gnu quasi_points_p.gnu giving the necessary geometry to produce the voronoi regions. The gnuplot command then to plot the 3D interactive plot is gnuplot> splot “random_points_p.gnu” u 2:3:4, “random_points_v.gnu” with lines