視錐台(Frustum)を利用し、ある高さにおいてオブジェクトが描画されている範囲を求める方法

#はじめに
Qiita初投稿になります。
間違っている点、不明な点などあれば優しく沢山マサカリを投げていただけると嬉しいです。

Unityを触っていて、<b>「オブジェクトが現在カメラに写っているか」</b>を知りたいということがあるかと思います。
その場合、最も簡単な方法はObjectにアタッチされている`Renderer`コンポーネントの`isVisible()`メソッドを呼び出すことです。

```csharp
public bool IsRendered()
{
    var renderer = this.GetComponent<Renderer>();
    return renderer.isVisible();
}
```

この方法は監視対象のオブジェクトが少ない場合は非常に便利です。

ただ、多くのオブジェクトや特定の座標を監視したい場合、毎回このメソッドを呼ぶためにRendererコンポーネントを取得するのは不便ですし、計算負荷もそれなりにかかります。

今回、この記事では<b>ある基準となる高さ（y座標）における、オブジェクトが描画されている範囲を求める方法</b>を解説します。
この方法を用いると複数オブジェクトの描画範囲の計算がまとめて行えるので、床面の高さが変わらないボードゲームやアクションゲームで特に効果を発揮すると思われます。

まあここまでめんどくさい処理をしないといけない機会もなかなかなさそうですし、探せばもっと楽な方法がありそうな気もしますが、個人的な備忘録程度にまとめておきます。

#視錐台（Frustum）
シーンにおける3次元空間の中で、カメラが描画している範囲は<b>視錐台（Frustum）</b>と呼ばれ、以下の図のように三角柱の上部分を削り取ったような形をしています。この図ではFar clipping planeと呼ばれる面とNear clipping planeに挟まれる、台形の立体版のような範囲がこれにあたります。
![ViewFrustum.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/433211/00c2f2b0-f609-4ae2-1392-9275f1d1efa1.png)
この図形は6つの面から構成されており、それぞれカメラの「左」「右」「上」「下」「奥」「手前」を表します。
詳しくは[Unity公式マニュアル](https://docs.unity3d.com/ja/current/Manual/UnderstandingFrustum.html)をご覧ください。画像もこちらから引用させていただきました。

Frustumを構成する面は`GeometryUtility.CalculateFrustumPlanes()`メソッドに現在のCameraを引数として渡してあげることで取得できます。

```csharp
// メソッドの戻り値の配列は要素数6で、それぞれカメラ側から見て
// [0]: 左の面, [1]: 右の面, [2]: 下の面, [3]: 上の面, [4]: 手前の面, [5]: 奥の面
Plane[] planes = GeometryUtility.CalculateFrustumPlanes(Camera.main);
```

ここで得られる戻り値は`UnityEngine.Plane`クラスなるものの配列で、今回は要素として<b>平面の法線ベクトル</b>と<b>原点（Vector3.zero）からの距離</b>を持つ、という理解のみで大丈夫です。
詳しくは[UnityスクリプティングAPI](https://docs.unity3d.com/jp/460/ScriptReference/Plane.html)をご覧ください。

今回はこの得られたplaneを平面の方程式に変換し、（上と左）（上と右）といった接する2平面の交線を求めます。
そしてその交線上のある高さにおける座標を求めることで、その２つの面の接合部のある高さにおける座標を求めます。

#ある高さにおけるオブジェクトの描画判定
ここからの説明は算数の知識が少し必要となるので（私は半分理解をあきらめました）、結論まで飛ばしていただいて大丈夫です。

##平面の方程式

３次元空間上のある平面は

```math
ax+by+cz+d=0
```

という数式で表すことができます。
ここで (a, b, c) は平面の法線を表し、`unityEngine.Plane`クラスでは`normal`変数にVector3型で格納されています。

dは、平面上に存在する座標 (x<sub>on</sub>, y<sub>on</sub>, z<sub>on</sub>) を (x, y, z) にそれぞれ代入したときに求められるもので、つまり

```math
d = -(a(x-x_{on})+b(y-y_{on})+c(z-z_{on}))
```

となります。

が、今回はこの平面に関する情報として「原点からの距離」が分かっているので別の方法を用いてdを求めます。
詳細は省きますが、平面外のある座標 (x<sub>out</sub>, y<sub>out</sub>, z<sub>out</sub>) からの距離Dは以下のようにあらわされます。

```math
D = \frac{|ax_{out}+by_{out}+cz_{out}-d|}{\sqrt{a^2+b^2+c^2}}
```

今回の場合は座標 (x<sub>out</sub>, y<sub>out</sub>, z<sub>out</sub>) は原点、距離Dは`UnityEngine.Plane`クラスの`distance`変数に格納されているので、以下のようにdの式に変形できます。

```math
d = -D * \sqrt{a^2+b^2+c^2}
```

これで`Plane`クラスを平面の方程式に変換することができました。

C#で平面の方程式の各要素を求めるコードは以下のようになります。

```csharp
Plane plane;
float a = plane.normal.x;
float b = plane.normal.y;
float c = plane.normal.z;
float d = -plane.distance * Mathf.Pow(a*a+b*b+c*c, 0.5f);
```

※6/10追記　式が誤っていたので修正しました

##2平面の交線を求める
まだ数式のターンは続きます。
2つの平面があったとき、その平面は一つの交線で交わるか平行であるかのどちらかです。
交線は理系の人だったら多分一度は求めたことがあると思いますが<font color="gray">（2式を連立させて1変数を消して媒介変数表示して...）</font>、私は遠い過去の記憶に消えてしまっているので[ここのサイト様](https://www.hiramine.com/programming/graphics/3d_planeintersection.html)に全面的にお世話になりました。
一般的に、3次元空間上の直線は以下の式で表されます。

```math
\begin{bmatrix}
x \\
y \\
z 
\end{bmatrix} = A+te
```

ここでAは直線上に存在する点の座標、tは媒介変数、eは交線の方向ベクトルです。

ちょっと何を言っているかよくわからないですが、定数であるAとeを求めることでこの式は定められます。
eは2平面の法線ベクトルと共に垂直なベクトルであり、外積計算で求められるらしいです。

2つの平面の平面の方程式をそれぞれ

```math
a_1x+b_1y+c_1z+d_1=0 \\
a_2x+b_2y+c_2z+d_2=0
```

とおくと、eは以下の式で定められます。

```math
\begin{bmatrix}
e_x \\
e_y \\
e_z 
\end{bmatrix}
=
\begin{bmatrix}
b_1c_2 - c_1b_2 \\
c_1a_2 - a_1c_2 \\
a_1b_2 - b_1a_2
\end{bmatrix}
```

つぎにAの値を求めます。
交線上の点Aは、両方の平面上に存在する点、ということができます。

詳しい求め方に関する話は上記のサイト様を参考にしていただくとすると、Aの値は
e<sub>z</sub> != 0 の場合、

```math
\begin{bmatrix}
A_x \\
A_y \\
A_z 
\end{bmatrix}
=
\begin{bmatrix}
\frac{d_1b_2-d_2b_1}{e_z} \\
\frac{d_1a_2-d_2a_1}{-e_z} \\
0
\end{bmatrix}
```

e<sub>y</sub> != 0 の場合、

```math
\begin{bmatrix}
A_x \\
A_y \\
A_z 
\end{bmatrix}
=
\begin{bmatrix}
\frac{d_1c_2-d_2c_1}{-e_y} \\
0 \\
\frac{d_1a_2-d_2a_1}{e_y}
\end{bmatrix}
```

e<sub>x</sub> != 0 の場合、

```math
\begin{bmatrix}
A_x \\
A_y \\
A_z 
\end{bmatrix}
=
\begin{bmatrix}
0 \\
\frac{d_1c_2-d_2c_1}{e_x} \\
\frac{d_1b_2-d_2b_1}{-e_x}
\end{bmatrix}
```

このどれでもない場合、つまり e<sub>x</sub> = e<sub>y</sub> = e<sub>z</sub> = 0 となる場合、2つの平面は平行です。

お疲れ様です、これで2平面の交線の方程式を求めることができました！~~大変すぎて記事書くのやめようかと思った~~
平面`plane0`および平面`plane1`の交線の方程式を求めるコードは以下のようになります

```csharp
var e = new Vector3(plane0.normal.y * plane1.normal.z - plane0.normal.z * plane1.normal.y,
        plane0.normal.z * plane1.normal.x - plane0.normal.x * plane1.normal.z,
        plane0.normal.x * plane1.normal.y - plane0.normal.y * plane1.normal.x);
var d0 =
    -Mathf.Pow(
        plane0.normal.x * plane0.normal.x + plane0.normal.y * plane0.normal.y +
        plane0.normal.z * plane0.normal.z, 0.5f) * plane0.distance;
var d1 =
    -Mathf.Pow(
        plane1.normal.x * plane1.normal.x + plane1.normal.y * plane1.normal.y +
        plane1.normal.z * plane1.normal.z, 0.5f) * plane1.distance;
            
Vector3 A;
if (e.z != 0)
{
    A = new Vector3((d0 * plane1.normal.y - d1 * plane0.normal.y) / e.z,
        (d0 * plane1.normal.x - d1 * plane0.normal.x) / (-e.z), 0);
}
else if (e.y != 0)
{
    A = new Vector3((d0 * plane1.normal.z - d1 * plane0.normal.z) / (-e.y), 0,
        (d0 * plane1.normal.x - d1 * plane0.normal.x) / e.y);
}
else if (e.x != 0)
{
    A = new Vector3(0, (d0 * plane1.normal.z - d1 * plane0.normal.z) / e.x,
        (d0 * plane1.normal.y - d1 * plane0.normal.y) / (-e.x));
}
else
{
    A = Vector3.positiveInfinity;
}
```

##ある高さ（y座標）における描画範囲を計算する
ここまでくれば後は簡単です。
先ほど求めた直線の方程式の一つの次元（今回はy）に定数を代入し、tの値を求めることでxおよびz座標を求めます。

```math
\begin{bmatrix}
x \\
y_0 \\
z 
\end{bmatrix}
=
\begin{bmatrix}
A_x \\
A_y \\
A_z
\end{bmatrix}
+
t
\begin{bmatrix}
e_x \\
e_y \\
e_z
\end{bmatrix}
```

より、

```math
t = \frac{y_0-a_y}{e_y}
```

となるので、求める座標 (x<sub>0</sub>, z<sub>0</sub>) は

```math
\begin{bmatrix}
x_0 \\
z_0 
\end{bmatrix}
=
\begin{bmatrix}
A_x + \frac{y_0-a_y}{e_y}e_x \\
A_z + \frac{y_0-a_y}{e_y}e_z
\end{bmatrix}
```

と求めることができました。
ある高さ y0 におけるｘおよびｚ座標を求めるコードは以下のようになります。

```csharp
var t = (y0 - A.y) / e.y;
var x0 = A.x + t * e.x;
var z0 = A.z + t * e.z;
```

#サンプルコードと動作例
サンプルコードを以下に示します。
ここでは、描画範囲の端を表す4つGameObject(`leftUp`, `leftDown`, `rightUp`, `rightDown`)を基準となる高さ`baseHeight`上で動かしています。
視錐台での近い面及び遠い面の2つの面に関しては計算を省略し、上下左右の4面の交線から各座標を計算しています。

```csharp
using UnityEngine;

namespace Sample
{
    public class DrawFrustumCorner : MonoBehaviour
    {
        public float baseHeight;

        public GameObject leftUp;
        public GameObject leftDown;
        public GameObject rightUp;
        public GameObject rightDown;

        private void Update()
        {
            var planes = GeometryUtility.CalculateFrustumPlanes(Camera.main);
            
            var leftUpPosition = GetPlaneIntersention(planes[0], planes[3]);
            var rightUpPosition = GetPlaneIntersention(planes[1], planes[3]);
            var leftDownPosition = GetPlaneIntersention(planes[0], planes[2]);
            var rightDownPosition = GetPlaneIntersention(planes[1], planes[2]);

            leftUp.transform.position = new Vector3(leftUpPosition.x, baseHeight, leftUpPosition.y);
            rightUp.transform.position = new Vector3(rightUpPosition.x, baseHeight, rightUpPosition.y);
            leftDown.transform.position = new Vector3(leftDownPosition.x, baseHeight, leftDownPosition.y);
            rightDown.transform.position = new Vector3(rightDownPosition.x, baseHeight, rightDownPosition.y);
        }

        private Vector2 GetPlaneIntersention(Plane plane0, Plane plane1)
        {
            var e = new Vector3(plane0.normal.y * plane1.normal.z - plane0.normal.z * plane1.normal.y,
                plane0.normal.z * plane1.normal.x - plane0.normal.x * plane1.normal.z,
                plane0.normal.x * plane1.normal.y - plane0.normal.y * plane1.normal.x);
            var d0 =
                -Mathf.Pow(
                    plane0.normal.x * plane0.normal.x + plane0.normal.y * plane0.normal.y +
                    plane0.normal.z * plane0.normal.z, 0.5f) * plane0.distance;
            var d1 =
                -Mathf.Pow(
                    plane1.normal.x * plane1.normal.x + plane1.normal.y * plane1.normal.y +
                    plane1.normal.z * plane1.normal.z, 0.5f) * plane1.distance;
            
            Vector3 A;
            if (e.z != 0)
            {
                A = new Vector3((d0 * plane1.normal.y - d1 * plane0.normal.y) / e.z,
                    (d0 * plane1.normal.x - d1 * plane0.normal.x) / (-e.z), 0);
            }
            else if (e.y != 0)
            {
                A = new Vector3((d0 * plane1.normal.z - d1 * plane0.normal.z) / (-e.y), 0,
                    (d0 * plane1.normal.x - d1 * plane0.normal.x) / e.y);
            }
            else if (e.x != 0)
            {
                A = new Vector3(0, (d0 * plane1.normal.z - d1 * plane0.normal.z) / e.x,
                    (d0 * plane1.normal.y - d1 * plane0.normal.y) / (-e.x));
            }
            else
            {
                A = Vector3.positiveInfinity;
            }
            
            var t = (baseHeight - A.y) / e.y;
            
            return new Vector2(A.x + t * e.x, A.z + t * e.z);
        }
    }
}
```

このクラスを適当なGameObjectにアタッチして実行させると、以下のようになりました。
分かりやすいようにベースとなる高さにplaneを置き、カメラを四角柱、求められた四隅を球で表現しています
![result.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/433211/4429aa0d-38eb-ae15-ea0b-53bb4f3ce7b5.gif)

今回は一定の高さにおける描画範囲を可視化させただけですが、<b>この4点の内側判定は[この記事](http://www.sousakuba.com/Programming/gs_hittest_point_triangle.html)の三角形の内側判定を2回行う</b>ことで実現できます。
カメラ位置が変化しないのであれば、この手法を用いることで各オブジェクトに対してこれまでの計算を省略できるため、<b>オブジェクト数が多いプロジェクト</b>では特に効果を発揮するかと思います。

#結論
今回は、一定の高さにおけるオブジェクトの描画判定をFrustumの各面の法線ベクトル及び距離から行う方法について説明しました。
この手法は<b>オブジェクトが多く、かつ一定の高さに存在するようなプロジェクト</b>において、「カメラの中にいる間だけアニメーションを実行する」などといった状況で<b>高速に描画判定を行う</b>ことができます。

ここまで書いてきて思いましたが、Shaderを書くときに出てくるMVP行列などは確かこれと同じような処理をしてたような気がするので、もしかしたらさらに簡潔に書ける方法があるかもしれませんね、まあ一例として見ていただけたらと思います。

#補足：Planeを用いた別の描画判定方法
今回の主題とは異なる方法ですが、それぞれのPlaneオブジェクトは視錐台の内側を向いた法線を持つため、次の関数にオブジェクトの範囲（Bounds）を引数として与えることで内外判定（true/false）を行うこともできます。

```csharp
public bool IsRendered(Bounds bounds)
{
    var cam = Camera.main;
    var planes = GeometryUtility.CalculateFrustumPlanes(cam);
    return GeometryUtility.TestPlanesAABB(planes, bounds);
}
```

Boundsは[Unityスクリプトリファレンス](https://docs.unity3d.com/ja/current/ScriptReference/Bounds.html)に詳しい内容があるのでここでは割愛しますが、3次元空間のある範囲を表します。
`TestPlanesAABB()`メソッドは[ここに詳しく書いてありますが](https://docs.unity3d.com/jp/460/ScriptReference/GeometryUtility.TestPlanesAABB.html)、Boundsが引数として渡されたPlaneすべての法線方向に存在しているときにtrueを返す関数です。

ある範囲が描画されているか知りたい場合は、こちらの方法を使うのがよいでしょう。

#参考文献
[視錐台を理解する -Unity公式マニュアル](https://docs.unity3d.com/ja/current/Manual/UnderstandingFrustum.html)
[GeometryUtility.CalculateFrustumPlanes -UnityScriptingAPI](https://docs.unity3d.com/540/Documentation/ScriptReference/GeometryUtility.CalculateFrustumPlanes.html)
[GeometryUtility.TestPlanesAABB -UnityスクリプティングAPI](https://docs.unity3d.com/jp/460/ScriptReference/GeometryUtility.TestPlanesAABB.html)
[２平面の交線](https://www.hiramine.com/programming/graphics/3d_planeintersection.html)
[点と三角形の当たり判定( 内外判定 )](http://www.sousakuba.com/Programming/gs_hittest_point_triangle.html)