方向ベクトルと平面の内分点 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2016-06-03

とある3D オブジェクトであるうさぎを貫くベクトルがあるとする。
ベクトルは極座標表示なので、中心からある程度伸ばして、そしてxyz 座標に変換する。
うさぎ自体はmesh3d オブジェクトで、もともとは頂点オブジェクト bun.v, エッジオブジェクト bun.f からできたものとする。

library(cwhmisc)
pol <- c(1, 1.475076, 1.970379) # 極座標
xyz <- t(mapply(toXyz, r=seq(0, 0.9, length=300), theta=pol[2], phi=pol[3])) # 直交座標

mesh.tri <- tmesh3d(t(bun.v), t(bun.f), homogeneous=FALSE)
plot3d(bun.v, type="n", axes=FALSE, xlab="", ylab="", zlab="")
shade3d(mesh.tri)
points3d(xyz, col=4)

　
いま、ベクトルが貫いている三角形を取り出したい。というのも、ベクトルと三角形の交点である、三角形表面上の点の動きを追いたいと思っているから。
はじめに、極座標とxyz 座標の変換は
$\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}r\sin{\theta}\cos{\phi}\\r\sin{\theta}\sin{\phi}\\r\cos{\theta}\end{pmatrix}$
$\theta$ と $\phi$ は決まっているので、 $r$ が動いて三角形を貫く。これを決めたい。
　
一方、ベクトルが三角形を点 $P$ で貫くとするならば、三角形を構成するベクトル $(v_1,v_2,v_3)$ により
$P=av_1+bv_2+cv_3$ とかける。ただし、 $P$ は三角形の内分点なので、 $a+b+c=1,\hspace{3}a>0,b>0,c>0$ である。
また、空間上の平面の方程式は
$\alpha X+\beta Y+\gamma Z+\delta=0$
と書ける。3点だと変数4つは解けないが、ある文字で固めれば
$\alpha X+\beta Y+\gamma Z+1=0$
というような形に持っていける。これで、3点で決まる平面上にベクトルのさきっちょが乗るようには解ける。
　
次に、 $P=av_1+bv_2+cv_3$ と書けるように $(a,b,c)$ の組を決める。これも、 $P$ が三角形に乗っていて $(v_1,v_2,v_3)$ で表せられることと、極座標ベクトルで表せられることを方程式で解けば出る。
このとき、 $(v_1,v_2,v_3)$ で表現しているから、内分点の前提を満たしている $a>0,b>0,c>0$ の組だけ採用する。
スクリプトではなぜか $r<0$ も採用されていたので、 $(v_1,v_2,v_3)$ に対してすべて同方向(なす角が180度以下)を拾った。
　

norm_sph <- c(sin(pol[2])*cos(pol[3]), sin(pol[2])*sin(pol[3]), cos(pol[2])) # 極座標ベクトル
abc <- mapply(function(i){
    hoge_tri <- bun.v[bun.f[i,],] # bun.f は三角形を構成するvertex のID が入っている。
                                  # つまり、ある三角形をなす3つの頂点の座標を得る。
    Rad <- 1/sum(base::solve(hoge_tri, rep(-1, 3)) * norm_sph) # 平面の方程式
    abc <- -base::solve(t(hoge_tri), Rad*norm_sph)             # abc を求める
    if( all(abc > 0) & colSums(abc*hoge_tri)%*%norm_sph > 0) colSums(abc*hoge_tri)
  }, seq(nrow(bun.f)))
triangle_points <- do.call(rbind.data.frame, abc)
colnames(triangle_points) <- NULL
triangle_points

points3d(triangle_points, col=4)
# ほしい三角形だけかく
m2 <- tmesh3d(t(bun.v), t(bun.f[sapply(abc, length)>0,]), homogeneous=FALSE)
wire3d(m2, col=4)

　
総当りでなんとかできるけれども、三角形の数が多くなれば計算オーダーが増えるので、 $a+b+c=1,\hspace{3}a>0,>0,c>0$ の条件をうまくいれて線形計画法的に解きたかったけどゴリ押しした。
　
追記
疎行列の対角成分に3x3 行列、つまり三角形を構成するベクトルを納めて、疎行列のまま線形代数的に解けば早い。
ただし、これは3つのベクトルの組で表した場合のパラメータ表示が求まるだけなので、三角形平面に乗っている（和が1）なのは別途、もう一回疎3つ組ベクトルだけ取り出して計算する。

v <- c(1, 1, 1)
library(Matrix)
n <- ncol(g$it)
i0 <- c(mapply(function(z) rep((3*z-2):(z*3), 3), 1:n))
j0 <- rep(1:(n*3), each=3)
sm <- sparseMatrix(i=i0, j=j0, x=c(g$vb[-4, g$it]))
res <- matrix(solve(sm, rep(v, n)), nr=3)