intersectionLinePlane <- function(LM, PM) { # Bryan Hanson, DePauw University, June 2013 # Use as you see fit. # Function to find the intersection of a line with a plane # LM is a 2 x 3 matrix with columns x, y, z giving 2 pts on the line # PM is a similar 3 x 3 matrix, giving 3 points on the plane ### Helper function from RSEIS (3D vector cross product) xprod <- function(A1, A2) { x = A1[2] * A2[3] - A1[3] * A2[2] y = A1[3] * A2[1] - A1[1] * A2[3] z = A1[1] * A2[2] - A2[1] * A1[2] return(c(x, y, z)) } # Based upon the method described # paulbourke.net/geometry/pointlineplane/ La <- LM[1,] Lb <- LM[2,] P0 <- PM[1,] P1 <- PM[2,] P2 <- PM[3,] N <- xprod(P2-P0, P1-P0) # correct based upon graphical testing n <- sum(N[1]*(P2 - La)[1], N[2]*(P2 - La)[2], N[3]*(P2 - La)[3]) d <- sum(N[1]*(Lb - La)[1], N[2]*(Lb - La)[2], N[3]*(Lb - La)[3]) if (isTRUE(all.equal(d, 0.0))) stop("Line is parallel to or in the plane") u <- n/d if ((u >= 0.0) && (u <= 1.0)) message("Plane is between the two points") P <- La + u*(Lb-La) P # correct based upon graphical testing } # Some test runs ML <- matrix(rnorm(6), nrow = 2) MP <- matrix(rnorm(9), nrow = 3) chk <- intersectionLinePlane(ML, MP) ML <- matrix(rnorm(6), nrow = 2) MP <- matrix(rnorm(9), nrow = 3) chk <- intersectionLinePlane(ML, rbind(ML, MP[3,])) # check for d = 0 ML <- matrix(c(0,0,0,1,1,1), nrow = 2, byrow = TRUE) MP <- diag(nrow = 3) chk <- intersectionLinePlane(ML, MP) # check for betweenness # graphical checks library("rgl") open3d() segments3d(ML, col = "red", lwd = 3) segments3d(rbind(MP[1:2,], MP[2:3,], MP[c(3,1),]), col = "blue", lwd = 3) axes3d(box = TRUE) title3d(xlab = "x", ylab = "y", zlab = "z") points3d(chk[1], chk[2], chk[3], col = "green", size = 10)