hierarchy in voronoi tessellations

Question

I am working with voronoi tessellations. I have different polygons representing regions in the tessellations.

The points below are used to draw the tessellation in the figure.

tessdata
    [,1]       [,2]
1  -0.4960583 -0.3529047
2  -2.4986929  0.8897895
3   3.6514561 -1.3533369
4  -1.7263101 -5.5341202
5   2.2140143  0.3883696
6  -2.5208933 -1.4881461
7  -3.2556913  4.4535629
8   0.6423109 -2.8350062
9  -0.4160715  1.2676151
10  4.4059361  4.5641771

Using tessdata as input to draw the tessellation as below:

library(deldir)
dd<-deldir(tessdata[,1], tessdata[,2])
plot(dd,wlines="tess")

Sammon coordinates are below.

       [,1]        [,2]
1   3.14162704 -1.45728604
2   2.35422623  2.46437927
3  -0.85051049  2.71503294
4   1.94310458 -0.45936958
5   0.08737757  3.74324701
6   1.23007799  1.34443842
7   0.01571924  2.19322032
8   1.43320754  2.64818631
9  -0.05463431  0.66980876
10   1.51344967  5.03351176

I want to construct the tessellations for which the sammon coordinate points are input. The tessellation using these points should be within one of the regions in the figure shown and for that, the above points should be scaled or we can restrict the plot of the tessellation within one of the regions in the above figure.

Hope i have covered all the necessary data.

P.S:

sammon's projection comes in "MASS" package. voronoi tessellations from "deldir" package.

dirsgs argument of the deldir function output will give the coordinates of the points forming the lines in the tessellations.

segments function of package graphics can be used to join the 2 points whose coordinates are extracted from dirsgs.

Answer 1

If you want to restrict the second set of points to one of the tiles of the tessellation, you can use tile.list to have a description of each tile, and then check which points are in this tile (there are many functions to do so: in the following example, I use secr::pointsInPolygon).

# Sample data
x <- matrix( rnorm(20), nc = 2 )
y <- matrix( rnorm(1000), nc=2 )

# Tessellation
library(deldir)
d <- deldir(x[,1], x[,2])
plot(d, wlines="tess")

# Pick a cell at random 
cell <- sample( tile.list(d), 1 )[[1]]
points( cell$pt[1], cell$pt[2], pch=16 )
polygon( cell$x, cell$y, lwd=3 )

# Select the points inside that cell
library(secr)
i <- pointsInPolygon(
  y, 
  cbind( 
    c(cell$x,cell$x[1]), 
    c(cell$y,cell$y[1])
  )
)
points(y[!i,], pch=".")
points(y[i,], pch="+")

# Compute a tessellation of those points
dd <- deldir(y[i,1], y[i,2])
plot(dd, wlines="tess", add=TRUE)

If, instead, you want to translate and rescale the points to fit them into the tile, that is trickier.

We need to somehow estimate how far away from the tile the points are: to this end, let us define a few auxilliary functions to compute, first the distance from a point to a segment, then the distance from a point to a polygon.

distance_to_segment <- function(M, A, B) {
  norm <- function(u) sqrt(sum(u^2))
  lambda <- sum( (B-A) * (M-A) ) / norm(B-A)^2
  if( lambda <= 0 ) {
    norm(M-A)
  } else if( lambda >= 1 ) {
    norm(M-B)
  } else {
    N <- A + lambda * (B-A)
    norm(M-N)
  }
}
A <- c(-.5,0)
B <- c(.5,.5)
x <- seq(-1,1,length=100)
y <- seq(-1,1,length=100)
z <- apply(
  expand.grid(x,y), 
  1, 
  function(u) distance_to_segment( u, A, B )
)
par(las=1)
image(x, y, matrix(z,nr=length(x)))
box()
segments(A[1],A[2],B[1],B[2],lwd=3)

library(secr)
distance_to_polygon <- function(x, poly) {
  closed_polygon <- rbind(poly, poly[1,])
  if( pointsInPolygon( t(x), closed_polygon ) )
    return(0)
  d <- rep(Inf, nrow(poly))
  for(i in 1:nrow(poly)) {
    A <- closed_polygon[i,]
    B <- closed_polygon[i+1,]
    d[i] <- distance_to_segment(x,A,B)
  }
  min(d)
}
x <- matrix(rnorm(20),nc=2)
poly <- x[chull(x),]
x <- seq(-5,5,length=100)
y <- seq(-5,5,length=100)
z <- apply(
  expand.grid(x,y), 
  1, 
  function(u) distance_to_polygon( u, poly )
)
par(las=1)
image(x, y, matrix(z,nr=length(x)))
box()
polygon(poly, lwd=3)

We can now look for a transformation of the form

x --> lambda * x + a
y --> lambda * y + b

that minimizes the (sum of the squared) distances to the polygon. That is actually not sufficient: we are likely to end up with scaling factor lambda equal to (or close to) zero. To avoid this, we can add a penalty if lambda is small.

# Sample data 
x <- matrix(rnorm(20),nc=2)
x <- x[chull(x),]
y <- matrix( c(1,2) + 5*rnorm(20), nc=2 )
plot(y, axes=FALSE, xlab="", ylab="")
polygon(x)

# Function to minimize:
# either the sum of the squares of the distances to the polygon, 
# if at least one point is outside, 
# or minus the square of the scaling factor.
# It is not continuous, but (surprisingly) that does not seem to be a problem.
f <- function( p ) {
  lambda <- log( 1 + exp(p[1]) )
  a <- p[2:3]
  y0 <- colMeans(y)
  transformed_points <- t( lambda * (t(y)-y0) + a )
  distances <- apply(
    transformed_points, 
    1, 
    function(u) distance_to_polygon(u, x)
  )
  if( all(distances == 0) ) - lambda^2
  else                      sum( distances^2 )
}
# Minimize this function
p <- optim(c(1,0,0), f)$par
# Compute the optimal parameters
lambda <- log( 1 + exp(p[1]) )
a <- p[2:3]
y0 <- colMeans(y)
# Compute the new coordinates
transformed_points <- t( lambda * (t(y)-y0) + a )
# Plot them
segments( y[,1], y[,2], transformed_points[,1], transformed_points[,2], lty=3 )
points( transformed_points, pch=3 )
library(deldir)
plot( 
  deldir( transformed_points[,1], transformed_points[,2] ), 
  wlines="tess", add=TRUE 
)