[R] Plot polygon in 3D with rgl

[Duncan Murdoch]

On 09/05/2010 6:41 PM, Remko Duursma wrote:
> Dear R-helpers, an rgl-ers in particular,
>
>
> what is the easiest way to plot a section of a plane in 3D, that is
> given by the xyz coordinates of the outline?
>
> Suppose I have a polygon - which I know for sure is a set of
> coordinates on the same plane. One method I found is to use surf.tri
> from the geometry package, and then plot the triangles with
> rgl.triangles. This method is not perfect though, because it makes the
> convex hull - no good for the polygon below.
>
> thanks for your help,
>
> Remko
>
> # Example polygon in 3D (matrix with tree columns x,y,z)
> m <- structure(c(-24.71, -36.45, -59.54, -83.97, -112.63, -126.66,
> -152.79, -171.04, -178.92, -183.71, -189.27, -191.56, -195.98,
> -203.09, -207.89, -212.12, -216.92, -221.73, -210.58, -199.43,
> -195.58, -175.38, -161.72, -152.09, -140.75, -123.63, -97.87,
> -79.04, -69.63, -61.75, -54.08, -38.34, -30.47, -24.71, -15.44,
> -13.33, -14.26, -17.39, -24.5, -31.92, -52.3, -70.45, -79.41,
> -89.16, -102.08, -111.61, -119.76, -125.55, -128.56, -132.55,
> -135.55, -138.55, -132.93, -127.31, -124.91, -125.78, -126.69,
> -127.42, -124.38, -117.74, -105.69, -88.52, -79.94, -70.97, -59.43,
> -34.75, -25.78, -15.44, 689.49, 686.54, 680.81, 674.79, 667.78,
> 664.42, 658.27, 654.05, 652.25, 651.22, 650.06, 649.66, 648.7,
> 647.04, 645.89, 644.91, 643.77, 642.63, 645.3, 647.98, 648.89,
> 653.93, 657.35, 659.75, 662.52, 666.66, 672.86, 677.25, 679.44,
> 681.24, 682.95, 686.44, 688.24, 689.49), .Dim = c(34L, 3L))
>
> # One (wrong) way
> library(geometry)
> tm <- t(surf.tri(m, delaunayn(m)))
> plot3d(m[,1],m[,2],m[,3], type='l', col="black", size=2)
> rgl.triangles(m[tm, 1], m[tm, 2], m[tm, 3],  col="green")
>
>   

There's a function triangulate() in the gpclib package that can 
triangulate a 2d polygon.  So you could pick two out of your three 
dimensions, triangulate those, then compute the 3rd coordinate from a 
fitted plane to the other two.  For example:

library(gpclib)
x <- m[,1]
y <- m[,2]
z <- m[,3]
triangles <- triangulate(as(cbind(x,y), "gpc.poly"))
zfit <- predict(lm(z ~ x + y), newdata=data.frame(x=triangles[,1], 
y=triangles[,2]))
triangles3d(cbind(triangles, zfit), col="green")

Duncan Murdoch