[Matplotlib] 平面と法線ベクトル

平面の方程式

　この記事では三次元空間内の 平面 を表示するコードを掲載しています。

平面と法線ベクトル

　固定点 $A(x_0,y_0,z_0)$ と平面に垂直なベクトル $\vec{n}=(p,q,r)$ が与えられたとします。
　このようなベクトルを平面の 法線ベクトル とよびます。
　平面上の点を $P(x,y,z)$ で表すと、ベクトル $\vec{AP}$ は法線ベクトルに垂直です。

　つまり $\vec{AP}$ と $\vec{n}$ の内積は $0$ です：
　
\[\vec{n}\cdot\vec{AP}=0\tag{1}\]
　具体的に成分で書き表すと
　
\[\begin{bmatrix}p\\q\\r\end{bmatrix}\cdot
\begin{bmatrix}x-x_0\\y-y_0\\z-z_0\end{bmatrix}=0\tag{2}\]
となります。よって、平面の方程式は
　
\[p(x-x_0)+q(y-y_0)+r(z-z_0)=0\tag{3}\]
で与えられることになります。

　Matplotlib を使って平面を描いてみましょう。以下に定義する nvector_plane() は法線ベクトルと任意の点を与えて平面をプロットする関数です。

# PLANE_01-1

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# 法線ベクトルを指定して平面をプロットする関数
# axes：サブプロット
# vector：法線ベクトル
# point：平面上の点
# xrange,yrange,zrange：x軸,y軸,z軸の範囲
# loc：法線ベクトルの始点(デフォルトは原点)
# vcolor：法線ベクトルの色
# pcolor,alpha：平面の色,透明度

def nvector_plane(axes, vector, point,
                  xrange, yrange, zrange,
                  loc = [0, 0, 0],
                  vcolor="red", pcolor="blue", alpha=0.5):

    # 軸ラベルの設定
    axes.set_xlabel("x", fontsize = 16)
    axes.set_ylabel("y", fontsize = 16)
    axes.set_zlabel("z", fontsize = 16)

    # 軸範囲の設定
    axes.set_xlim(xrange[0], xrange[1])
    axes.set_ylim(yrange[0], yrange[1])
    axes.set_zlim(zrange[0], zrange[1])

    # 格子点の作成
    x = np.arange(xrange[0], xrange[1], 0.2)
    y = np.arange(yrange[0], yrange[1], 0.2)
    xx, yy = np.meshgrid(x, y)

    # 平面の方程式
    zz = point[2] - (vector[0]*(xx-point[0])+vector[1]*(yy-point[1])) / vector[2]

    # 平面をプロット
    ax.plot_surface(xx, yy, zz, color=pcolor, alpha=alpha)

    # 法線ベクトルをプロット
    axes.quiver(loc[0], loc[1], loc[2],
                vector[0], vector[1], vector[2],
                color = vcolor, length = 1, arrow_length_ratio = 0.2)

　点 $(0,0,0)$ を通り、法線ベクトル $\vec{n}=(2,5,3)$ に垂直な平面を表示させてみましょう。

# PLANE_01-2

# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection='3d')

# 平面上の点
p = (0, 0, 0)

# 法線ベクトル
v = (2, 5, 3)

# 法線ベクトルの始点
loc = (0, 0, 0)

# x,y,z軸の範囲
xr = [-5, 5]
yr = [-5, 5]
zr = [-5, 5]

# 法線ベクトルと平面をプロット
nvector_plane(ax, loc, v, p, xr, yr, zr)

# 画像の保存
plt.savefig("normal_vector_plane.png", bbox_inches = "tight")

　$3$ 点 $A,\ B,\ C$ が指定されると 1 つの平面が定まります。
　この平面の法線ベクトルは $\vec{AB}$ と $\vec{AC}$ に垂直なので、
　
\[\vec{n}=\vec{AB}\times \vec{AC}\tag{4}\]
で与えられます。nvector_plane() を再利用して、$3$ 点を通る平面を描画する関数を定義できます。

# PLANE_01-3

# 3点を通る平面をプロットする関数
# axes：サブプロット
# point：平面上の点
# xrange,yrange,zrange：x軸,y軸,z軸の範囲
# loc：法線ベクトルの始点(デフォルトは原点)
# vcolor：法線ベクトルの色
# pcolor,alpha：平面の色,透明度

def point_plane(axes, p0, p1, p2,
                xrange, yrange, zrange, 
                loc = [0, 0, 0],
                vcolor="red", pcolor="blue", alpha=0.5):

    u = p1 - p0
    v = p2 - p0
    w = np.cross(u, v)
    w = 0.5 * np.abs(xrange[0]-yrange[1]) * w / np.sqrt(np.sum(w**2))
    
    nvector_plane(axes, w, p0,
                  xrange, yrange, zrange,
                  loc=loc, vcolor=vcolor, pcolor=pcolor, alpha=alpha)

　
　点 $(3,1,1),\ (0,2,0),\ (-1, -1, 3)$ の $3$ 点を通る平面をプロットします。

# PLANE_01-4

# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection='3d')

# x,y,z軸の範囲
xr = [-5, 5]
yr = [-5, 5]
zr = [-5, 5]

# 平面内の3点を指定
p0 = np.array([3, 1, 1])
p1 = np.array([0, 2, 0])
p2 = np.array([-1, -1, 3])

# 3点の重心
g = (p0 + p1 + p2) / 3

# 3点を通る平面をプロット
point_plane(ax, p0, p1, p2, xr, yr, zr,loc=g)

一般的な平面方程式

　(3) で表された方程式
　
\[p(x-x_0)+q(y-y_0)+r(z-z_0)=0\tag{3}\]
の定数部分を $s$ にまとめてしまうと、一般的な 平面方程式 が得られます。
　
\[px+qy+rz+s=0\tag{5}\]
　式変形によって点 $A(x_0,y_0,z_0)$ の情報は失われましたが、係数 $p,q,r$ は残っているので法線ベクトルの情報は保持されています。さらに、$z=f(x,y)$ の形に変形すると
　
\[z=-\frac{px+qy+s}{r}\tag{6}\]
となります。(5) を使って平面を描く関数を定義してみます。

# PLANE_02-1

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# px+qy+rz+s=0をプロットする関数
# axes：サブプロット
# param：p,q,r,sのリストまたはタプルなど
# xrange,yrange,zrange：x軸,y軸,z軸の範囲
# pcolor,alpha：平面の色,透明度

def plot_plane(axes, param, xrange, yrange, zrange,
               pcolor="blue", alpha=0.5):

    # 軸ラベルの設定
    axes.set_xlabel("x", fontsize = 16)
    axes.set_ylabel("y", fontsize = 16)
    axes.set_zlabel("z", fontsize = 16)

    # 軸範囲の設定
    axes.set_xlim(xrange[0], xrange[1])
    axes.set_ylim(yrange[0], yrange[1])
    axes.set_zlim(zrange[0], zrange[1])

    # 格子点の作成
    x = np.arange(xrange[0], xrange[1], 0.2)
    y = np.arange(yrange[0], yrange[1], 0.2)
    xx, yy = np.meshgrid(x, y)

    # 平面の方程式
    zz = -(param[0]*xx + param[1]*yy + param[3]) / param[2]

    # 平面をプロット
    ax.plot_surface(xx, yy, zz, color=pcolor, alpha=alpha)

　plot_plane() を使って、$2x-y+3z=0$ で表される平面をプロットします。

# PLANE_02-2

# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection='3d')

# x,y,z軸の範囲
xr = [-5, 5]
yr = [-5, 5]
zr = [-5, 5]

# p,q,r,sを設定
param = [2, -1, 3, 0]

# 平面2x-y+3z=0をプロット
plot_plane(ax, param, xr, yr, zr)

# 画像の保存
plt.savefig("plot_plane.png", bbox_inches = "tight")

基底ベクトルで表される平面

　三次元空間における平面は 2 つの独立な三次元ベクトル (基底ベクトル) の線型結合によって表すことができます。
　
\[\begin{bmatrix}x\\y\\z\end{bmatrix}
=s\begin{bmatrix}a\\b\\c\end{bmatrix}
+t\begin{bmatrix}d\\e\\f\end{bmatrix}\tag{7}\]
　パラメータ (媒介変数) で表される方程式なので、Python に実装するときにはパラメータの格子点データをつくります。

# PLANE_03-1

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# 基底ベクトルを指定して平面をプロットする関数
# axes：サブプロット
# v,w：基底ベクトル
# xrange,yrange,zrange：x軸,y軸,z軸の範囲
# vcolor,wcolor：基底ベクトルの色
# pcolor,alpha：平面の色,透明度

def bvector_plane(axes, v, w, xrange, yrange, zrange,
                  srange, trange,
                  vcolor="red", wcolor="blue",
                  pcolor="green", alpha=0.2):

    # 軸ラベルの設定
    axes.set_xlabel("x", fontsize = 16)
    axes.set_ylabel("y", fontsize = 16)
    axes.set_zlabel("z", fontsize = 16)

    # 軸範囲の設定
    axes.set_xlim(xrange[0], xrange[1])
    axes.set_ylim(yrange[0], yrange[1])
    axes.set_zlim(zrange[0], zrange[1])

    # 格子点の作成
    s = np.arange(srange[0], srange[1], 0.2)
    t = np.arange(trange[0], trange[1], 0.2)
    ss, tt = np.meshgrid(s, t)

    # 基底ベクトルの線型結合
    xx = v[0]*ss + w[0]*tt
    yy = v[1]*ss + w[1]*tt
    zz = v[2]*ss + w[2]*tt

    # 平面をプロット
    ax.plot_surface(xx, yy, zz, color=pcolor, alpha=alpha)

    # 基底ベクトルをプロット
    axes.quiver(0, 0, 0, v[0], v[1], v[2],
                color=vcolor, length=1, arrow_length_ratio=0.2)
    
    axes.quiver(0, 0, 0, w[0], w[1], w[2],
                color=wcolor, length=1, arrow_length_ratio=0.2)

　bvector_plane() を使って基底ベクトル $\vec{v}=(1, 1, 1)$ と $\vec{w}=(2, -1, 4)$ が張る平面をプロットしてみます。

# PLANE_03-2

# FigureとAxes
fig = plt.figure(figsize = (8, 6.5))
ax = fig.add_subplot(111, projection='3d')

# 基底ベクトル
v = [1,  1, 1]
w = [2, -1, 4]

# パラメータの範囲
srange = [-2, 2]
trange = [-2, 2]

# 軸範囲の設定
xrange = [-5, 5]
yrange = [-5, 5]
zrange = [-5, 5]

# 基底ベクトルと平面をプロット
bvector_plane(ax, v, w, xrange, yrange, zrange, srange, trange)