Die Ebene wird durch zwei (nicht parallele) Vektoren \(\vec{u}\) und \(\vec{v}\) definiert:

u = c(1, 0, 0) # basis vectors of the plane 
v = c(0, 1, 0)

Der Normalenvektor zu dieser Ebene ergibt sich aus dem Kreuzprodukt der Basisvektoren:

suppressWarnings(library(pracma))
n = cross(u,v)
n
## [1] 0 0 1
nx = n[1]
ny = n[2]
nz = n[3]  # a, b, and c define the normal vector

Wir wählen den Nullpunkt (Punkt auf der Ebene):

x0 = 0   # origin
y0 = 0
z0 = 0

Wir berechnen die Hilfsvariable \(d\):

d = -(nx*x0 + ny*y0 + nz*z0)  
d 
## [1] 0

Wir definieren den Punkt, der auf die Ebene projiziert werden soll:

x1 = 1  # we project this point
y1 = 0.5
z1 = 1.5

Eine weitere Hilfsvariable hatten wir zuvor k genannt:

k = (d + nx*x1 + ny*y1 + nz*z1) / (nx^2 + ny^2 + nz^2) 
k
## [1] 1.5

Nun können wir die Koordinaten des projizierten Punktes berechnen:

# calculate 3D-coordinates of the projected point:
xp = x1 - k*nx  
yp = y1 - k*ny  
zp = z1 - k*nz  
rp = c(xp, yp, zp)   # projected point
rp
## [1] 1.0 0.5 0.0

Wir plotten die Ebene, den zu projizierenden und den projizierten Punkt mit Hilfe des packages rgl:

suppressWarnings(library(rgl))
options(rgl.useNULL = TRUE) # don't open in a separate window
axes3d()
wire3d(scale3d(cube3d(color = "black", alpha = 0.1), 2, 2, 2))
title3d(xlab = "x", ylab = "y", zlab = "z") 
planes3d(nx, ny, nz, d, col = "green", alpha = 0.5, lwd = 2) 
plot3d(x = x1, y = y1, z = z1, col = "red", type = 's', radius = .1, add = TRUE) # the point to projec
plot3d(x = xp, y = yp, z = zp, col = "red", type = 's', radius = .1, add = TRUE) # the projected point
lines3d(c(x1, xp),  c(y1, yp), c(z1, zp), col = "yellow", lty = 2, lwd = 2)      # line connecting both points
rglwidget()

Wir stellen noch eine Funktion zur Berechnung der Koordinaten des projizierten Punktes bereit. Hier sind \(u\) und \(v\) die Basisvektoren der Ebene, \(r_0\) ist irgendein Punkt auf dieser Ebene, \(r_1\) ist der zu projizierende Punkt:

projectPoint <- function(u, v, r0, r1) {
  if(!is.vector(u) || !is.vector(v) || !is.vector(r0) || !is.vector(r1)) stop("All input must be vectors")
  suppressWarnings(library(pracma))  
  n = cross(u,v) 
  if(all(n == c(0,0,0))) stop("Basis vectors must not be parallel!")
  nx = n[1]; ny = n[2]; nz = n[3]      # nx, ny, and nz define the normal vector
  x0 = r0[1]; y0 = r0[2]; z0 = r0[3]   # some point on the plane
  x1 = r1[1]; y1 = r1[2]; z1 = r1[3]   # the point to be projected
  d = -(nx*x0 + ny*y0 + nz*z0)
  k = (d + nx*x1 + ny*y1 + nz*z1) / (nx^2 + ny^2 + nz^2)
  xp = x1 - k*nx  
  yp = y1 - k*ny  
  zp = z1 - k*nz  
  rp = c(xp, yp,zp)   # projected point  
  return(rp)
}

Wir testen die Funktion mit dem gleichen Input wie oben:

u = c(1, 0, 0)        # basis vectors of the plane 
v = c(0, 1, 0)
r0 = c(0, 0, 0)       # point on the plane
r1 = c(1, 0.5, 1.5)   # we project this point
rp <- projectPoint(u = u, v = v, r0 = r0, r1 = r1)
rp
## [1] 1.0 0.5 0.0

Die Funktion ist kann hier heruntergeladen werden.

Wir wollen nun mehrere Punkte simultan projizieren. Wir definieren zunächst die Basisvektoren \(u\) und \(v\) der Ebene, auf die projiziert werden soll, sowie deren Ursprung \(r_0\):

u = c(1, 0, 0)       # u,v : basis vectors of the plane 
v = c(0, 1, 0)
r0 = c(0, 0, 0)      # point on the plane

Wir definieren einige zu projizierende Punkte, welche wir in einer Matrix r1s zusammenfassen:

point1 = c(1, 1, 1)       # x,y,z-coordinates
point2 = c(1, 0.5, 1.5)   # same as above
point3 = c(3, 1, 1.5)
point4 = c(2, 2, 0.5)
all_points = c(point1, point2, point3, point4)
r1s = matrix(all_points, byrow = FALSE, nrow = 3) # four points, each column of the matrix represents a point to project
rownames(r1s) = c("x", "y", "z")
colnames(r1s) = paste0("p", 1:4)
r1s
##   p1  p2  p3  p4
## x  1 1.0 3.0 2.0
## y  1 0.5 1.0 2.0
## z  1 1.5 1.5 0.5

Um alle Punkte gleichzeitig zu projizieren, wenden wir die Funktion apply an:

rps = apply(r1s, 2, projectPoint, u = u, v = v, r0 = r0)
rps  # each column of rps represents a solution for a column of the input matrix
##   p1  p2 p3 p4
## x  1 1.0  3  2
## y  1 0.5  1  2
## z  0 0.0  0  0

Das Ergebnis für Punkt 2 stimmt mit dem oben erhaltenen Ergebnis überein.

Wir zeichnen noch das Ergebnis. (Die Festlegung der Achsenlimits ist leider etwas kompliziert)

xlim_low  = min(r1s[1,], rps[1,])*0.8
xlim_high = max(r1s[1,], rps[1,])*1.2 
xlimits = c(xlim_low, xlim_high)

ylim_low  = min(r1s[2,], rps[2,])*0.8
ylim_high = max(r1s[2,], rps[2,])*1.2 
ylimits = c(ylim_low, ylim_high)

zlim_low  = min(r1s[3,], rps[3,])*0.8
zlim_high = max(r1s[3,], rps[3,])*1.2 
zlimits = c(zlim_low, zlim_high)

plot3d(x = r1s[1,], y = r1s[2,], z = r1s[3,], col = "red", 
       type = 's', radius = .05, xlab = "x", ylab = "y", zlab = "z",
       xlim = xlimits, ylim = ylimits, zlim = zlimits)
plot3d(x = rps[1,], y = rps[2,], z = rps[3,], col = "red", 
       type = 's', radius = .05, add = TRUE)
planes3d(nx, ny, nz, d, col = "yellow", alpha = 1, lwd = 2)

x = c(r1s[1,1], rps[1,1])
for (i in 2:ncol(r1s)) x = c(x, r1s[1,i], rps[1,i])
y = c(r1s[2,1], rps[2,1])
for (i in 2:ncol(r1s)) y = c(y, r1s[2,i], rps[2,i])
z = c(r1s[3,1], rps[3,1])
for (i in 2:ncol(r1s)) z = c(z, r1s[3,i], rps[3,i])
segments3d(x, y, z, col = "yellow", lty = 3, lwd = 2) # lines connecting point pairs
rglwidget()


uwe.menzel@matstat.org