e.blog

主にUnity/UE周りのことについてまとめていきます

セーブデータを暗号化して保存する

概要

ゲームには状態の保存など、いわゆる「セーブデータ」が必要なケースが多いです。
今回はそんな「セーブデータ」をシリアライズしたものをバイナリ化し、さらに暗号化して保存する方法を書いておきたいと思います。

なお、今回の実装にあたっては以下の記事を参考にさせていただきました。

qiita.com

developer.wonderpla.net

loumo.jp

使うクラス

セーブデータを暗号化して保存するために、以下のクラス群を利用します。 (もちろん、暗号化には様々なアルゴリズムが存在し、今回紹介する以外の方法でももちろん暗号化を行うことが可能です)

  • System.Runtime.Serialization.Formatters.Binary.BinaryFormatter
  • System.IO.MemoryStream
  • System.IO.FileStream
  • System.Security.Cryptography.MD5CryptoServiceProvider
  • System.Security.Cryptography.RijndaelManaged
  • System.Security.Cryptography.Rfc2898DeriveBytes
  • System.Security.Cryptography.ICryptoTransform

セーブデータをシリアライズしてバイナリデータとしてファイルに保存する

まずは暗号化の話をする前に、セーブデータをシリアライズしてバイナリデータとしてファイルに保存する方法を解説します。

SerializableAttributeでシリアライズ可能なことを明示する

System.SerializableAttributeをclassに指定することでそのクラスがシリアル化可能なことを明示することができます。

ドキュメント↓
docs.microsoft.com

バイナリ化する

データをシリアル化できるよう明示したら、次はそのオブジェクトをバイナリ化します。 バイナリ化にはMemoryStreamBinaryFormatterクラスを使います。

コード断片で示すと以下のようになります。

using (MemoryStream stream = new MemoryStream())
{
    BinaryFormatter formatter = new BinaryFormatter();
    formatter.Serialize(stream, data);

    byte[] source = stream.ToArray();

    using (FileStream fileStream = new FileStream(SavePath, FileMode.Create, FileAccess.Write))
    {
        fileStream.Write(source, 0, source.Length);
    }
}

コード量はそんなに多くないのでぱっと見でなんとなく分かるかと思います。 MemoryStreamを生成し、BinaryFormatterSerializeメソッドを利用してオブジェクトをシリアライズします。 シリアライズしたbyte配列はMemoryStreamに書き込まれます。 結果のbyte配列を取得するにはToArrayメソッドを使います。

そして最後にFileStreamを生成し、byte配列をファイルに書き込みます。 バイナリ化に関してはBinaryFormatterが行ってくれるので、IO周りがしっかり把握できていればさしてむずかしい処理ではないと思います。

データをAESで暗号化する

バイナリデータの保存が分かったところで、次は暗号化についてです。 今回取り上げるのは「AES暗号化」です。

AESとは

AESとは以下の記事から引用させていただくと、

www.atmarkit.co.jp

「AES(Advanced Encryption Standard)」は、DESの後継として米国の国立標準技術研究所(NIST:National Institute of Standards and Technology)によって制定された新しい暗号化規格である。

とのこと。

そしてさらに以下のように続いています。

そして最終的に2001年に「Rijndael(ラインダール)」という暗号化方式が選ばれた。開発者はベルギーの暗号学者、「Joan Daemen(ホァン・ダーメン)」と「Vincent Rijmen(フィンセント・ライメン)」であり、Rijndaelという名称は2人の名前から取られた(とされている)。

実際にコードを見てもらうと分かりますが、実装にはRijndaelManagedというクラスが利用されており、これがまさに上の暗号化方式の名前となっていますね。

暗号化についてのアルゴリズムなどについては上記記事を読んでみるとなんとなく雰囲気は分かるかと思います。
そして完全に余談ですが、以下のiOSアプリが、色々なアルゴリズム(暗号化や証明書などなど)についてアニメーション付きで分かりやすく解説してくれているのでよかったらダウンロードしてみてください。

アルゴリズム図鑑

アルゴリズム図鑑

  • Moriteru Ishida
  • 教育
  • 無料

実装

暗号化に関しては前段で解説した際に行ったバイナリ化したデータ(byte配列)に対して操作を行い、暗号化します。 なので暗号化している部分だけを抜粋してコードを紹介します。 (引数のbyte[] dataは、前段で生成したbyte配列です)

static SaveDataManager()
{
    _rijindeal = new RijndaelManaged();
    _rijindeal.KeySize = 128;
    _rijindeal.BlockSize = 128;

    byte[] bsalt = Encoding.UTF8.GetBytes(_salt);
    Rfc2898DeriveBytes deriveBytes = new Rfc2898DeriveBytes(_password, bsalt);
    deriveBytes.IterationCount = 1000;

    _rijindeal.Key = deriveBytes.GetBytes(_rijindeal.KeySize / 8);
    _rijindeal.IV = deriveBytes.GetBytes(_rijindeal.BlockSize / 8);
}

static private byte[] Encrypt(byte[] data)
{
    ICryptoTransform encryptor = _rijindeal.CreateEncryptor();
    byte[] encrypted = encryptor.TransformFinalBlock(data, 0, data.Length);

    encryptor.Dispose();

    // Console.WriteLine(string.Join(" ", encrypted));

    return encrypted;
}

static private byte[] Dencrypt(byte[] data)
{
    ICryptoTransform decryptor = _rijindeal.CreateDecryptor();
    byte[] plain = decryptor.TransformFinalBlock(data, 0, data.Length);

    // Console.WriteLine(string.Join(" ", plain));

    return plain;
}

コードの冒頭では静的コンストラクタによってRijndaelManagedオブジェクトを生成しています。

そしてEncryptDencryptメソッドでデータの暗号化、復号化を行っています。 ここでやっていること自体はとてもシンプルですね。 静的コンストラクタで生成したRijndaelManagedオブジェクトから、CreateEncryptorCreateDecryptorをそれぞれ生成し、ICryptorTransformインターフェースのTransformFinalBlockメソッドを実行しているだけです。

戻り値は暗号化、復号化されたbyte配列となります。 あとはこれを、前段のファイル保存の処理で保存してやれば晴れて、セーブデータが暗号化されて保存されたことになります。 (結局のところ、最終的に保存されるのは01で表されるバイナリ表現のデータなので、それ自体が暗号化されているか否かに関わらず、ファイルの保存・読み込みは問題なく行えるというわけですね。(というか、FileStreamはそれを関知しない)

ソースコード

最後に、コンソールアプリとして実行するといくつかの項目を入力するとそれを保存、復元できるものを作ったのでソースコードを載せておきます。
もちろん、Unity上でも動作します。

using System;
using System.IO;
using System.Text;
using System.Runtime.Serialization.Formatters.Binary;
using System.Security.Cryptography;

[System.Serializable]
public class SaveData
{
    public float Number = 0.5f;
    public string Name = "Hoge";
    public int Count = 5;

    public override string ToString()
    {
        return string.Format("Name: {0}, Number: {1}, Count: {2}", Name, Number, Count);
    }
}

static public class SaveDataManager
{
    public const string SavePath = "./test.bytes";
    private const string _password = "passwordstring";
    private const string _salt = "saltstring";
    static private RijndaelManaged _rijindeal;

    static SaveDataManager()
    {
        _rijindeal = new RijndaelManaged();
        _rijindeal.KeySize = 128;
        _rijindeal.BlockSize = 128;

        byte[] bsalt = Encoding.UTF8.GetBytes(_salt);
        Rfc2898DeriveBytes deriveBytes = new Rfc2898DeriveBytes(_password, bsalt);
        deriveBytes.IterationCount = 1000;

        _rijindeal.Key = deriveBytes.GetBytes(_rijindeal.KeySize / 8);
        _rijindeal.IV = deriveBytes.GetBytes(_rijindeal.BlockSize / 8);
    }

    static public void Save(SaveData data)
    {
        using (MemoryStream stream = new MemoryStream())
        {
            BinaryFormatter formatter = new BinaryFormatter();
            formatter.Serialize(stream, data);

            byte[] source = stream.ToArray();

            source = AESlize(source);

            using (FileStream fileStream = new FileStream(SavePath, FileMode.Create, FileAccess.Write))
            {
                fileStream.Write(source, 0, source.Length);
            }

            Console.WriteLine("Done [" + data.ToString() + "]");
        }
    }

    static public SaveData Load(string name)
    {
        SaveData data = null;

        using (FileStream stream = new FileStream(name, FileMode.Open, FileAccess.Read))
        {
            using (MemoryStream memStream = new MemoryStream())
            {
                const int size = 4096;
                byte[] buffer = new byte[size];
                int numBytes;

                while ((numBytes = stream.Read(buffer, 0, size)) > 0)
                {
                    memStream.Write(buffer, 0, numBytes);
                }

                byte[] source = memStream.ToArray();
                source = DeAESlize(source);

                using (MemoryStream memStream2 = new MemoryStream(source))
                {
                    BinaryFormatter formatter = new BinaryFormatter();
                    data = formatter.Deserialize(memStream2) as SaveData;

                    Console.WriteLine("Loaded.");
                }
            }
        }

        return data;
    }

    static private byte[] AESlize(byte[] data)
    {
        ICryptoTransform encryptor = _rijindeal.CreateEncryptor();
        byte[] encrypted = encryptor.TransformFinalBlock(data, 0, data.Length);

        encryptor.Dispose();

        // Console.WriteLine(string.Join(" ", encrypted));

        return encrypted;
    }

    static private byte[] DeAESlize(byte[] data)
    {
        ICryptoTransform decryptor = _rijindeal.CreateDecryptor();
        byte[] plain = decryptor.TransformFinalBlock(data, 0, data.Length);

        // Console.WriteLine(string.Join(" ", plain));

        return plain;
    }
}

static public class EntryPoint
{
    static public void Main()
    {
        Console.WriteLine("Save? [y/n]");
        string cond = Console.ReadLine();

        if (cond == "y")
        {
            Save();
        }
        else
        {
            Load();
        }
    }

    static void Save()
    {
        Console.WriteLine("Name?");
        string name = Console.ReadLine();

        Console.WriteLine("Number?");
        float number;
        if (!float.TryParse(Console.ReadLine(), out number))
        {
            Console.WriteLine("Must input float value.");
            return;
        }

        Console.WriteLine("Count?");
        int count;
        if (!int.TryParse(Console.ReadLine(), out count))
        {
            Console.WriteLine("Must input int value.");
            return;
        }

        SaveData data = new SaveData
        {
            Name = name,
            Number = number,
            Count = count,
        };
        SaveDataManager.Save(data);
    }

    static void Load()
    {
        SaveData data = SaveDataManager.Load(SaveDataManager.SavePath);
        Console.WriteLine(data.ToString());
    }
}

UNETのLow Level APIを使ってシンプルな位置同期の仕組みを作る

概要

UNETとはUnityが提供しているネットワークゲームを作るためのシステムです。
UNETには「High Level API」と「Low Level API」があり、High Level APIについては以前、記事に簡単にまとめました。

edom18.hateblo.jp

今回はLow Level APIについて書きたいと思います。
サンプルとして、ごく簡単に位置を同期してみます。

実際に動かしてみたのがこちら↓

なお、今回の実装のサンプルはGithubにアップしてあるので興味がある人は見てみてください。

github.com

また、今回の実装にあたっては以下の記事を参考にさせていただきました。ありがとうございます。

tarukosu.hatenablog.com

概観する

まず、基本的なコネクションのフローなどをざっくりと見ていきたいと思います。

登場人物

ネットワーク対応にするにあたり、必要となるクラスなどを紹介します。

クラス

メインとなるクラスです。NetworkTransportsealedクラスになっていて、実質staticクラスになっています。

ネットワーク関連の初期化(IP AddressやPort番号の設定など)、データの送受信などネットワークに関連する処理は本クラスを用いて行います。

以下のHostTopologyクラスのコンストラクタ引数に指定するコンフィグ情報です。

トポロジーとは

こちらの記事から引用させてもらうと、

元々は位相幾何学の意味であるが、パソコン用語としては、LANの接続形態の総称。ネットワークトポロジーとも呼ばれる。LANには接続形態によって、リング型やバス型、スター型などがあり、接続する機器によって用途や目的に合った接続形態を使用する必要がある。

つまりは、どういう形態のネットワークなのか、を定義します。
コンストラクタにはコンフィグオブジェクトと最大接続数を指定します。

ConnectionConfig config = new ConnectionConfig();

// ... 中略 ...

int maxConnection = 10;
HostTopology topology = new HostTopology(config, maxConnection);

そして生成したホストトポロジーオブジェクトをホストとして登録します。

int hostId = NetworkTransport.AddHost(topology, _localPort);

登録すると、Host IDが返されるので保持しておきます。(後の通信時に利用します)

enum

QosType 説明
Unreliable There is no guarantee of delivery or ordering.
UnreliableFragmented There is no guarantee of delivery or ordering, but allowing fragmented messages with up to 32 fragments per message.
UnreliableSequenced There is no guarantee of delivery and all unordered messages will be dropped. Example: VoIP.
Reliable Each message is guaranteed to be delivered but not guaranteed to be in order.
ReliableFragmented Each message is guaranteed to be delivered, also allowing fragmented messages with up to 32 fragments per message.
ReliableSequenced Each message is guaranteed to be delivered and in order.
StateUpdate An unreliable message. Only the last message in the send buffer is sent. Only the most recent message in the receive buffer will be delivere |d.
ReliableStateUpdate A reliable message. Note: Only the last message in the send buffer is sent. Only the most recent message in the receive buffer will |be delivered.
AllCostDelivery A reliable message that will be re-sent with a high frequency until it is acknowledged.
UnreliableFragmentedSequenced There is garantee of ordering, no guarantee of delivery, but allowing fragmented messages with up to 32 fragments per mess |age.
ReliableFragmentedSequenced Each message is guaranteed to be delivered in order, also allowing fragmented messages with up to 32 fragments per message.

Quality of serviceの略です。到達の保証や順序の保証についてのクオリティを指定します。

NetworkEventType 説明
DataEvent Data event received. Indicating that data was received.
ConnectEvent Connection event received. Indicating that a new connection was established.
DisconnectEvent Disconnection event received.
Nothing No new event was received.
BroadcastEvent Broadcast discovery event received. To obtain sender connection info and possible complimentary message from them, call NetworkTransport.GetBroadcastConnectionInfo() and NetworkTransport.GetBroadcastConnectionMessage() functions.

データを受信した際のイベントタイプです。

フロー

簡単なフローは以下の通りです。

  1. NetworkTransport.Init();で初期化
  2. ConnectionConfigでコンフィグ情報を生成する
  3. 必要なQosType(Quality of service)をチャンネルに追加する
  4. HostTopolosyトポロジーオブジェクトを生成
  5. NetworkTransport.AddHost(...);でホスト(トポロジー)を登録
  6. NetworkTransport.Connect(...);でサーバへ接続(クライアントの場合)

セットアップ

さて、では実際のコードを参考にフローを見てみます。

ネットワークの初期化

サーバ/クライアントに関わらず、ネットワークを開始するためのセットアップを行います。

/// <summary>
/// ネットワークのセットアップ
/// </summary>
private void NetworkSetup()
{
    Debug.Log("Setup network.");

    _isStarted = true;

    // ネットワークの初期化処理
    // NetworkTransportクラスの他のメソッドを呼び出す前に初期化が必要
    NetworkTransport.Init();

    // コンフィグ情報の生成
    ConnectionConfig config = new ConnectionConfig();
    _channelId = config.AddChannel(_qosType);

    // トポロジーオブジェクトの生成
    int maxConnections = 10;
    HostTopology topology = new HostTopology(config, maxConnections);

    _hostId = NetworkTransport.AddHost(topology, _localPort);
}

コネクション、データ受信などのネットワークイベントをハンドリングする

ネットワークの初期化が終わったら、サーバへの接続およびデータの受信についてのハンドリングを行います。

前述の通りNetworkEventTypeというイベントタイプが定義されており、データを受信した際のタイプとして使用します。

データの受信にはNetworkTransport.Receive(...);メソッドを利用します。
詳細はコードを見てください。

/// <summary>
/// ネットワークイベントのハンドリング
/// </summary>
private void Receive()
{
    int receiveHostId;
    int connectionId;
    int channelId;
    byte[] receiveBuffer = new byte[_maxBufferSize];
    int bufferSize = _maxBufferSize;
    int dataSize;
    byte error;

    // 受信データは複数ある場合があるので、ループでデータが無くなるまで処理する
    do
    {
        NetworkEventType eventType = NetworkTransport.Receive(out receiveHostId, out connectionId, out channelId, receiveBuffer, bufferSize, out dataSize, out error);

        switch (eventType)
        {
            case NetworkEventType.Nothing:

                // 受信データなし

                return;
            case NetworkEventType.ConnectEvent:
                Debug.Log("Connected.");

                // Connection処理

                break;
            case NetworkEventType.DisconnectEvent:
                Debug.Log("Disconnected.");

                // Disconnection処理

                break;
            case NetworkEventType.DataEvent:
                
                // データ受信の処理

                break;
        }
    } while(true);
}

サーバへ接続する

サーバへ接続するにはNetworkTransport.Connect(...);メソッドを利用します。
以下のようにして、Host IDServer AddressServerPortExeption Connection IDを指定して接続を試みます。

あくまで試みるだけで、実際に接続が完了した場合はデータ受信(ハンドリング)のところで、ConnectEventが渡ってくるので、そこで適切に処理します。

/// <summary>
/// サーバへコネクションを張る
/// </summary>
private void Connect()
{
    byte error;
    _connectionId = NetworkTransport.Connect(_hostId, _serverAddress, _serverPort, 0, out error);
}

シリアライズしてデータを送受信する

接続が完了したら、次は必要なデータを送信します。
ただ注意が必要なのは、Vector3などの構造体はもちろん、floatなどのプリミティブな値もそのままでは送信することができません。

それらをbyte配列に変換し、シリアライズしてからデータを送信する必要があります。

/// <summary>
/// 位置を同期する
/// </summary>
private void SyncPosition()
{
    if (!_isServer && !_hasConnected)
    {
        return;
    }

    if (!_hasAuthorized)
    {
        return;
    }

    byte[] x = ConversionUtil.ToBytes(_target.position.x);
    byte[] y = ConversionUtil.ToBytes(_target.position.y);
    byte[] z = ConversionUtil.ToBytes(_target.position.z);

    byte[] pos = ConversionUtil.Serialize(x, y, z);

    if (_isServer)
    {
        for (int i = 0; i < _clientIds.Count; i++)
        {
            SendData(pos, _clientIds[i], pos.Length);
        }
    }
    else
    {
        SendData(pos, _connectionId, pos.Length);
    }
}
/// <summary>
/// データ送信
/// </summary>
/// <param name="data">送信するデータ配列</param>
public void SendData(byte[] data, int connectionId, int dataSize)
{
    byte error;
    NetworkTransport.Send(_hostId, connectionId, _channelId, data, dataSize, out error);
}

データのシリアライズ

シリアライズとは、こちらの記事から引用させてもらうと

シリアライズとは、複数の並列データを直列化して送信することである。

具体的には、メモリ上に存在する情報を、ファイルとして保存したり、ネットワークで送受信したりできるように変換することである。他方、既にファイルとして存在しているデータや、一旦シリアライズされたデータがネットワークから送られてきた際に、プログラムで扱えるようにする作業をデシリアライズと呼ぶ。

これが前述の「byte配列に変換して」の部分ですね。
intfloatなどは4バイトで表現されます。
つまり、これを1バイトごとに分解してシリアライズし、それをまとめて送る必要があります。

今回のサンプルでは、位置データ(Vector3)の各要素(x, y, z)をbyte配列に変換し、さらにそれをシリアル化してまとめて送信しています。

変換用のユーティリティクラスは以下のように実装しました。
(あくまで今回の簡易的な位置同期のための実装です)

[StructLayout(LayoutKind.Explicit)]
public struct UIntFloat
{
    [FieldOffset(0)]
    public uint intValue;

    [FieldOffset(0)]
    public float floatValue;

    [FieldOffset(0)]
    public double doubleValue;

    [FieldOffset(0)]
    public ulong longValue;
}

/// <summary>
/// 値をByte列に、Byte列を値に変換する
/// </summary>
public class ConversionUtil
{
    static public float ToFloat(uint value)
    {
        UIntFloat uif = new UIntFloat();
        uif.intValue = value;
        return uif.floatValue;
    }

    static public float ToFloat(byte[] values)
    {
        int x = 0;

        if (System.BitConverter.IsLittleEndian)
        {
            x |= (values[3] <<  0);
            x |= (values[2] <<  8);
            x |= (values[1] << 16);
            x |= (values[0] << 24);
        }
        else
        {
            x |= (values[0] <<  0);
            x |= (values[1] <<  8);
            x |= (values[2] << 16);
            x |= (values[3] << 24);
        }

        UIntFloat uif = new UIntFloat();
        uif.intValue = (uint)x;
        return uif.floatValue;
    }

    static public double ToDouble(ulong value)
    {
        UIntFloat uif = new UIntFloat();
        uif.longValue = value;
        return uif.doubleValue;
    }

    static public double ToDouble(byte[] values)
    {
        long x = 0;

        if (System.BitConverter.IsLittleEndian)
        {
            x |= ((long)values[7] <<  0);
            x |= ((long)values[6] <<  8);
            x |= ((long)values[5] << 16);
            x |= ((long)values[4] << 24);
            x |= ((long)values[3] << 32);
            x |= ((long)values[2] << 40);
            x |= ((long)values[1] << 48);
            x |= ((long)values[0] << 56);
        }
        else
        {
            x |= ((long)values[0] <<  0);
            x |= ((long)values[1] <<  8);
            x |= ((long)values[2] << 16);
            x |= ((long)values[3] << 24);
            x |= ((long)values[4] << 32);
            x |= ((long)values[5] << 40);
            x |= ((long)values[6] << 48);
            x |= ((long)values[7] << 56);
        }

        UIntFloat uif = new UIntFloat();
        uif.longValue = (ulong)x;
        return uif.doubleValue;
    }

    /// <summary>
    /// Floatの値をByte列に変換、ビッグエンディアンとして返す
    /// </summary>
    /// <param name="value">変換するfloat値</param>
    /// <returns>ビッグエンディアンのByte配列</returns>
    static public byte[] ToBytes(float value)
    {
        UIntFloat uif = new UIntFloat();
        uif.floatValue = value;

        uint x = uif.intValue;
        byte a = (byte)((x >>  0) & 0xff);
        byte b = (byte)((x >>  8) & 0xff);
        byte c = (byte)((x >> 16) & 0xff);
        byte d = (byte)((x >> 24) & 0xff);

        if (System.BitConverter.IsLittleEndian)
        {
            return new[] { d, c, b, a };
        }
        else
        {
            return new[] { a, b, c, d };
        }
    }

    /// <summary>
    /// Doubleの値をByte列に変換、ビッグエンディアンとして返す
    /// </summary>
    /// <param name="value">変換するdouble値</param>
    /// <returns>ビッグエンディアンのByte配列</returns>
    static public byte[] ToBytes(double value)
    {
        UIntFloat uif = new UIntFloat();
        uif.doubleValue = value;

        ulong x = uif.longValue;
        byte a = (byte)((x >>  0) & 0xff);
        byte b = (byte)((x >>  8) & 0xff);
        byte c = (byte)((x >> 16) & 0xff);
        byte d = (byte)((x >> 24) & 0xff);
        byte e = (byte)((x >> 32) & 0xff);
        byte f = (byte)((x >> 40) & 0xff);
        byte g = (byte)((x >> 48) & 0xff);
        byte h = (byte)((x >> 56) & 0xff);

        if (System.BitConverter.IsLittleEndian)
        {
            return new[] { h, g, f, e, d, c, b, a };
        }
        else
        {
            return new[] { a, b, c, d, e, f, g, h };
        }
    }

    static public byte[] Serialize(params byte[][] bytesArray)
    {
        int total = 0;
        for (int i = 0; i < bytesArray.Length; i++)
        {
            total += bytesArray[i].Length;
        }

        byte[] result = new byte[total];
        int index = 0;
        for (int i = 0; i < bytesArray.Length; i++)
        {
            System.Array.Copy(bytesArray[i], 0, result, index, bytesArray[i].Length);
            index += bytesArray[i].Length;
        }

        return result;
    }

    static public float Deserialize(byte[] bytes, int start, int end)
    {
        int num = end - start;
        byte[] data = new byte[num];

        for (int i = 0; i < num; i++)
        {
            data[i] = bytes[start + i];
        }

        return ToFloat(data);
    }
}

実際に利用に足るものにするには、以下のUnityのNetworkBufferのように、様々な値を簡単に取り扱えるようにする必要があります。
各値をbyte配列にしたり、などを実装しているので参考にしてみるといいと思います。

Unity-Technologies / Networking / source / Runtime / NetworkBuffer.cs — Bitbucket

ちなみに上で登場しているUIntFloat構造体ですが、メモリレイアウトを明示することによって効率的にintfloatなどを相互変換しています。

細かい点については以下の記事を見てください。

ufcpp.net

そこで示されているのは以下のようなUnionライクな構造体です。

using System;
using System.Runtime.InteropServices;

[StructLayout(LayoutKind.Explicit)]
struct Union
{
    [FieldOffset(0)]
    public byte A;
    
    [FieldOffset(1)]
    public byte B;
    
    [FieldOffset(2)]
    public byte C;
    
    [FieldOffset(3)]
    public byte D;
    
    [FieldOffset(0)]
    public int N;
}

public class Hello
{
    public static void Main()
    {
        Union x = new Union{ N = 0x12345678 };
        Console.WriteLine(x.A.ToString("x")); // => 78
        Console.WriteLine(x.B.ToString("x")); // => 56
        Console.WriteLine(x.C.ToString("x")); // => 34
        Console.WriteLine(x.D.ToString("x")); // => 12
    }
}

プロトコルについて考える

最後に少しだけプロトコルについて。

上記記事から引用させてもらうと、

コンピューター同士が通信をする際の手順や規約などの約束事。ネットワークでコンピューターが使う言語のようなもので、双方が理解できる同じプロトコルを使わないと通信は成立しない。

ということ。

Low Level APIではデータ送信のルールは実装者が決める必要があります。
どういう順番でデータを格納し、先頭から何バイトまでがヘッダを表すのか、などです。

コンピュータが認識できるのは01のビット信号のみです。
そしてそれを「バイト」という単位で扱い、それがどういう順番で並び、それがどういう意味を持っているのか、は人間が考えないとなりません。

今回のサンプルでは位置同期だけを行い、さらには同期するオブジェクトはひとつだけに限定していたため、送られてきたデータが位置を表すx, y, zのfloat値が3要素分、と固定して解釈していました。

しかし、実際に使えるものにするためには、「どのオブジェクトの位置なのか」や「誰から送られてきたのか」などの制御が必要になってきます。

そしてそれらは、ヘッダなど「データの構造を示すデータ」を有することで実現します。
このヘッダなどの「メタデータ」が何バイトから何バイトまでか、を規定しそれに基づいて処理を行う仕様が「プロトコル」というわけですね。

なので、Low Level APIを使って実際に使えるレベルの通信を行うためにはまだ実装しないとならないことがたくさんあります。
が、最低限の通信に関してはそこまで複雑なことはやっていません。

これを元に、簡単な通信の仕組みを作れれば色々と捗りそうですね。

頂点座標とUV座標から接ベクトルを求める

概要

法線マップを利用する際に使われる「接ベクトル空間」。
今回はその接ベクトルを頂点座標とUV座標から求めるくだりを書いてみます。

先に断っておくと、今回の記事の内容はマルペケさんのこちらの記事(その12 頂点座標とUV座標から接ベクトルを求めるちょっと眠い話)を参考に、Unityで実装したものになります。

なので、今回の記事は上記記事を参考にUnityで実装するにあたって自分の理解をメモするためのものです。

ちなみにUnityで動かすとこんな感じ↓

今回実装したやつはGithubにアップしてあります。

github.com

今回は接ベクトルについての話のみになります。
バンプマップについては以前、Qiitaの記事で書いているのでそちらをご覧ください。

qiita.com

余談ですが、前回の記事では(法線マップは使っていませんが)バンプマップを使って波紋エフェクトを表示する内容で記事を書いたのでよかったら読んでみてください。

edom18.hateblo.jp

接ベクトル空間

法線マップを利用する上で「接ベクトル空間」の話は外せません。
なぜなら、法線マップの情報はこの接ベクトル空間での値が格納されているからです。

ということで、今回は接ベクトル空間、接ベクトルを求める話です。

接ベクトルはその名の通り、「とある点(ピクセル)に接している」ベクトルのことです。
図にすると以下のような感じ。

f:id:edo_m18:20180827084802p:plain

見てもらうと分かる通り、とある点(ピクセル)の上に乗っている直交座標系と見ることができます。
ただ、「乗せる」と一口に言っても、「どう乗せるか」が問題です。

もちろん座標系をどう取ろうが本来は自由です。
しかし、今回は「法線マップ」の法線を適切に扱うことを目標としているので、自分勝手に定義してしまっては適切な法線を得ることはできません。

そしてこれが今回の主旨である「頂点情報から接ベクトルを求める」になりますが、座標系をどう取ったら適切な法線が取れるのか、を解説していきます。

接平面

接平面とは、今回の接ベクトル空間のUVベクトルが成す平面です。
つまり接ベクトルが成す平面ということですね。

接平面については以下の記事が分かりやすいでしょう。

接平面と接ベクトル - 物理とか

詳細は記事を読んでいただくとして、ざっくりと説明します。

記事を引用すると、

接平面を求めるには、曲面のある点において2つの接ベクトルを求めればいい

ということなのでまずは曲面を表す式を以下のように定義します。


p = p ( u_1, u_2 )

これはふたつのパラメータによって曲面が表されています。
方針としては曲面のある点において2つの接ベクトルを求めるので、この点を含むふたつの曲線を表現しそれの接線を求めることで接平面を求めることにします。

そしてできるだけ分かりやすい曲線を選ぶようにするために、以下のようにふたつの関数を定義してみます。


x_1 = p( u_1(t), u_2 )  \\\
x_2 = p( u_1, u_2(t) )  \\

これは、曲面の定義から媒介変数 tによってパラメータが変化する2つの曲線の関数を取り出したと見ることができます。

さて、この曲線の変化を見ることで接線を求めることができるわけですね。
曲線の変化は微分の出番なので、以下のようにして偏微分で求めます。


v_i = \frac{\partial p}{\partial u_i}

頂点にある5要素から接ベクトルを求める

接ベクトルを求めるのに「接平面」を考えることを書きました。
だいぶざっくり書いてしまったので、細かな点については紹介した記事を読んでください。

ここで言いたかったのは、接平面をどう定義し、どうやって必要なベクトルを求めるか、のヒントを得ることです。

さて、話は変わって。
3Dの世界では頂点はポリゴンを形成するひとつの点を表します
なにを当たり前のことを、と思うかもしれませんが、頂点を表す要素のうち、(x, y, z)要素が頂点の「位置」を決めます。
昔ながらの3Dであれば「頂点カラー」と呼ばれる、頂点ごとに設定された色情報を読み出して処理することもありました。(最近ではカラーとしてそのまま使うケースは稀でしょう)

このように、「頂点ごとに様々な情報」を付与してレンダリングを行うのが3Dです。
そしてテクスチャ空間をあらわす(u, v)値も頂点に付与される情報です。
これはもちろん、テクスチャの位置を決める値です。

つまり、位置を司る情報としては(x, y, z)に加え、(u, v)の合計5要素がある、というわけですね。

そして見出しの通り、この(u, v)値を利用して接ベクトルを求めます。

テクスチャ座標のUV値の変化を観察する

まず、以下の図を見てみてください。

f:id:edo_m18:20180830131534p:plain

図ではU軸の方向に移動すると(x, y)の値がどう変化するかを示したものです。
(本来は(z)の値も同様に考える必要がありますが、図的に分かりづらいので(x, y)に絞って説明しています)

よくよく見てみると、UVの値の変化の方向はすでに「接平面」に存在していることが分かります。
そしてu方向に移動したとき、xyの値がどう変化するかを図示したのが上の図です。

u方向に少しだけ移動したとき、xの値がどう変化するのか、あるいはyの値がどう変化するのか。
それぞれの要素に絞って変化を見ているわけです。
これって「偏微分」の考え方ですよね。

そうです。U軸方向への(x, y, z)それぞれの要素の変化量を見ることでU軸の方向ベクトルを得ることができるのです。
これは言ってみれば、UV値の勾配を求めることでそれぞれのベクトルが求まる、ということですね。

そして今回求めたいのは接ベクトル空間の接ベクトルおよび従法線ベクトルです。
これはまさにU軸、V軸ベクトルの方向ですね。

偏微分の記号を使ってU軸およびV軸を表してみると、


U軸 = \biggr(\frac{\partial x}{\partial u}, \frac{\partial y}{\partial u}, \frac{\partial z}{\partial u}\biggl)

V軸 = \biggr(\frac{\partial x}{\partial v}, \frac{\partial y}{\partial v}, \frac{\partial z}{\partial v}\biggl)

と表すことができます。

ローカル座標をUV座標で表す

前述のように位置に関する頂点の情報は全部で5つあります。


P_0 = (x_0, y_0, z_0, u_0, v_0)

ですね。

マルペケさんの記事を引用させていただくと、

ここでうまい事を考えます。5つの成分のうち、例えば(x, u, v)の3成分だけに注目し、3頂点から平面を作ってみます。平面は点の数が3つあればできるわけです。平面の方程式にすると、

\(A_0 x + B_0 u + C_0 v + D_0 = 0\)

です。両辺を\(D_0\)で割ると正規化されて、

\(A_0 x + B_0 u + C_0 v + 1 = 0\)

となります。(ABCの記号が同じ添え字なのは目を瞑ってください(^-^;)。このABCは未知の係数ですが、今頂点が3つあるので解く事ができます。連立方程式を立てても良いのですが、平面の方程式の(A,B,C)は平面の法線である事を利用すると、ポリゴンの法線から一発で求まります。

とのこと。

さて、上記の式を整理すると以下のようになります。


\begin{eqnarray}
A_0 x + B_0 u + C_0 v + 1 &=& 0 \\\
A_0 x &=& -B_0 u - C_0 v - 1 \\\
x &=& -\frac{B_0}{A_0}u - \frac{C_0}{A_0}v - \frac{1}{A_0}
\end{eqnarray}

 xについての式に変形したわけですね。
さて、これを uについて偏微分してみます。すると以下のようになります。


\frac{\partial x}{\partial u} = -\frac{B_0}{A_0}

お。これは先に書いた偏微分 x要素じゃないですか。
そう、この値を導き出すために上記のように平面の方程式を持ち出したのですね。

なお、引用した文章でも示されている通り、この ABCは平面の方程式の意味から「法線ベクトルの各要素( x, y, z)」です。
そして平面の法線は3頂点からベクトルを作り、その外積によって求めることができます。

ちなみに平面の方程式についてはこちらを参照↓

mathtrain.jp

そしてU軸およびV軸のベクトルの各要素は、上の例を(y, z)にも適用してやることで以下のように求まります。


U = \biggr(\frac{\partial x}{\partial u}, \frac{\partial y}{\partial u}, \frac{\partial z}{\partial u}\biggl) = \biggr(-\frac{B_0}{A_0}, -\frac{B_1}{A_1}, -\frac{B_2}{A_2}\biggl)

V = \biggr(\frac{\partial x}{\partial v}, \frac{\partial y}{\partial v}, \frac{\partial z}{\partial v}\biggl) = \biggr(-\frac{C_0}{A_0}, -\frac{C_1}{A_1}, -\frac{C_2}{A_2}\biggl)

さぁ、これを元に実際のコードに落とし込むと以下のようになります。

実装コード

今回はUnityで実装したのでC#コードです。

static public Vector3[] GetTangentSpaceVectors(Vector3[] p, Vector2[] uv)
{
    Vector3[] cp0 = new[]
    {
        new Vector3(p[0].x, uv[0].x, uv[0].y),
        new Vector3(p[0].y, uv[0].x, uv[0].y),
        new Vector3(p[0].z, uv[0].x, uv[0].y),
    };

    Vector3[] cp1 = new[]
    {
        new Vector3(p[1].x, uv[1].x, uv[1].y),
        new Vector3(p[1].y, uv[1].x, uv[1].y),
        new Vector3(p[1].z, uv[1].x, uv[1].y),
    };

    Vector3[] cp2 = new[]
    {
        new Vector3(p[2].x, uv[2].x, uv[2].y),
        new Vector3(p[2].y, uv[2].x, uv[2].y),
        new Vector3(p[2].z, uv[2].x, uv[2].y),
    };

    Vector3 u = Vector3.zero;
    Vector3 v = Vector3.zero;

    for (int i = 0; i < 3; i++)
    {
        Vector3 v1 = cp1[i] - cp0[i];
        Vector3 v2 = cp2[i] - cp0[i];
        Vector3 ABC = Vector3.Cross(v1, v2).normalized;

        if (ABC.x == 0)
        {
            Debug.LogWarning("ポリゴンかUV上のポリゴンが縮退しています");
            return new[] { Vector3.zero, Vector3.zero };
        }

        u[i] = -(ABC.y / ABC.x);
        v[i] = -(ABC.z / ABC.x);
    }

    u.Normalize();
    v.Normalize();

    return new[] { u, v };
}

やっていることは前述の文章の説明をそのままプログラムしただけです。
それぞれのベクトルの外積を取って ABCを求め、そこからUV軸のベクトル要素としているのが分かるかと思います。

以上で接ベクトルを求めることができました。
もしバンプマップなどの計算に用いたい場合は、この接ベクトル空間へライトなどの位置を変換してやり、そこでのシェーディングを計算することで法線マップを利用したライティングが可能になる、というわけです。

パースペクティブコレクト

さて、最後は少しだけ余談です。

今回の色々を実装する際に初めて聞いた単語。

以下の記事に詳しく書かれていました。

ラスタライザを作る人の古文書集 - ushiroad

どういうものかざっくり書くと。

通常、3DCGで画面にレンダリングを行う際は頂点情報を入力しそれをいくつかの座標変換を経て、最終的にスクリーンに表示します。
そしてこのとき、ポリゴンに指定された値(色とUV値など)はSlopeとSpanという処理を行い補間します。

冒頭の記事から引用させていただくと、

ポリゴンを描画するとき、頂点の属性(典型的には色、UV座標)を滑らかに変化させながら各ピクセルを描きます。このとき、頂点の属性を頂点間を結ぶ辺上で内挿したもの(属性値の"坂")はSlopeと呼ばれます。 さらに、2本のSlopeの間にSpan(橋)を架け、この両端の値を内挿することにより、三角形内部の全ての点で内挿値を得ることができます。

ということ。
この補間処理の際、UV値を色などと同じ方法で補間してしまうと問題が生じる、ということのようです。
紹介した記事にはどんな感じになるのかの画像があるので見てみてください。

さて、ではどういうふうに補間処理をしたらいいのかというと。
これも冒頭の記事に解説があります。引用させていただくと、

パースペクティブコレクトの具体的な方法ですが、テクスチャ座標(u, v)をそのまま補間するのではなく、透視変換で出てくる斉次のwを使い

(u/w, v/w, 1/w) という値を頂点毎に作り、これでSlope/Spanの処理を行います。補間結果を

{ (u/w)' , (v/w)' , (1/w)' } としたら、(u/w)' および (v/w)' を (1/w)' で割ります(つまり、補間前にwの逆数をとり、補間後にもう一度wをひっくり返します)

自分の理解を書くと、射影変換された空間で補間処理を行い、かつ1/wも同様に補間しておく。そして補間後の値((u/w)', (v/w)')を、補間した(1/w)'で元に戻すことで正常な値が得られる、ということのようです。
(射影空間上で補完処理をして戻す、ということかもしれません)

パースペクティブコレクトを用いるとき

なぜパースペクティブコレクトについて書いたかというと。
前回書いた記事↓

edom18.hateblo.jp

これを実装しようと考えたとき、法線の方向の計算がおかしくてバンプマップについて改めて調べていたときにInk Painter - Asset Storeを作られている方のブログを読んでいて知った単語でした。

そのときに読んだ記事↓

esprog.hatenablog.com

最終的にはこの話はまったく関係なかったのですが、調べていた過程で見つけたってことで備忘録的に書いてみました。

ちなみに、普通はパースペクティブコレクトされた状態でUV座標は補間されるので意識することはあまりないと思いますが、上記記事では「とある点」に近い位置のUV値を算出するようにしています。

その際にまさにこのパースペクティブコレクトが必要となります。
UV値をいじる必要が出た場合はこれを思い出すといいかと思います。

CustomRenderTextureを使って波紋エフェクトを作る

概要

何番煎じか分かりませんが、CustomRenderTextureを使ってシェーダでお絵かきができると色々表現の幅が増えるので、それの練習のために表題のサンプルを作ってみました。

↓実際に動かしたの動画

サンプルはCustomRenderTextureを使って波動方程式を解き波形を描画、さらにそれをレンダリング用シェーダの法線マップとして入力し、歪みとライティング(ランバート反射とフォン反射)を適用してみたものです。

ダウンロード

今回のサンプルもGithubにアップしてあります。

github.com

CustomRenderTextureとは

ドキュメントは以下です。
https://docs.unity3d.com/ja/current/Manual/CustomRenderTextures.htmldocs.unity3d.com

ドキュメントから引用させてもらうと、

カスタムレンダーテクスチャはレンダーテクスチャの拡張機能で、これを使うと簡単にシェーダー付きのテクスチャを作成できます。これは、コースティクス、雨の効果に使われるリップルシミュレーション、壁面へぶちまけられた液体など、あらゆる種類の複雑なシミュレーションを実装するのに便利です。また、カスタムレンダーテクスチャはスクリプトやシェーダーのフレームワークを提供し、部分的更新、または、マルチパスの更新、更新頻度の変更などのさらに複雑な設定をサポートします。

つまり、波形などの複雑なシミュレーションを必要とする演算をシェーダベースで行い、それをテクスチャとして利用できるようにしてくれる機能です。
この機能の面白いところは、updateや初期化などはある程度テクスチャ側の設定項目でまかなえるので、シェーダだけでも完結できる点にあります。

つまりC#スクリプトはいらないということです。
もちろん、細やかな制御をするためにスクリプトから制御することも可能です。

今回の波紋シミュレーションではスクリプトによる更新はしていません。
そしてそれを法線マップとして利用することで冒頭の動画のようなエフェクトを生成しています。

CustomRenderTexture用のcgincがある

さて、では実際に実装を行っていきます。
CustomRenderTextureは前述の通り、シェーダ(マテリアル)を設定するだけでその結果をRenderTextureに保持してくれる便利なものです。

そのため、RenderTextureに保存したい絵を描くためのシェーダを書く必要があります。
(ちなみに、CustomRenderTextureが登場する前は、C#も動員して自分で保存などの処理を書く必要がありました)

まずはCustomRenderTexture向けのシェーダで利用するものを紹介します。

// 専用のcgincファイル
#include "UnityCustomRenderTexture.cginc"

// 専用の定義済みvertexシェーダ関数
#pragma vertex CustomRenderTextureVertexShader

// 専用の構造体
struct v2f_customrendertexture { ... }

// テクスチャサイズを取得
float width = _CustomRenderTextureWidth;
float height = _CustomRenderTextureHeight;

// UV座標を取得
float2 uv = i.globalTexcoord;

// CustomRenderTextureからテクセルをフェッチ
float3 c = tex2D(_SelfTexture2D, uv);

Unityでシェーダを書いたことがある人であればどう使うのかイメージ付くかと思います。 CustomRenderTexture向けにシェーダを書く場合、上記のようにいくつか専用のものが定義されているのでそれを利用します。

これらを使ってCurstomRenderTextureにお絵かきするシェーダを書いていきます。

CustomRenderTexture用シェーダを書く

では実際のシェーダを書いていきます。

波動方程式を解く

CustomRenderTextureで利用するシェーダが分かったところで、実用的な使い方として波動方程式を解いた波形を描いてみたいと思います。

ちなみに、以前自分もQiitaの記事で波動方程式について書いているので、よかったらそちらも参考にしてみてください。(実装はWebGLですが)

qiita.com

上の記事は以下の動画を元に実装したものになります。
波動方程式について、高校数学でも分かるような感じで丁寧に説明してくれていてとても分かりやすいのでオススメです。

www.nicovideo.jp

なお今回の実装はこちらの凹みさんの記事を参考にさせていただきました。

tips.hecomi.com

波動方程式を解くシェーダ(for CustomRenderTexture)

下記シェーダは上記の凹みさんの記事で書かれているものを利用させていただきました。
それを少し整形し、自分がコメントを追記したものです。

Shader "Hidden/WaveShader"
{
    Properties
    {
        _S2("PhaseVelocity^2", Range(0.0, 0.5)) = 0.2
        _Atten("Attenuation", Range(0.0, 1.0)) = 0.999
        _DeltaUV("Delta UV", Float) = 3
    }

        SubShader
    {
        Cull Off
        ZWrite Off
        ZTest Always

        Pass
        {
            CGPROGRAM
            #pragma vertex CustomRenderTextureVertexShader
            #pragma fragment frag

            #include "UnityCustomRenderTexture.cginc"

            half _S2;
            half _Atten;
            float _DeltaUV;
            sampler2D _MainTex;

            float4 frag(v2f_customrendertexture  i) : SV_Target
            {
                float2 uv = i.globalTexcoord;

                // 1pxあたりの単位を計算する
                float du = 1.0 / _CustomRenderTextureWidth;
                float dv = 1.0 / _CustomRenderTextureHeight;
                float3 duv = float3(du, dv, 0) * _DeltaUV;

                // 現在の位置のテクセルをフェッチ
                float2 c = tex2D(_SelfTexture2D, uv);

                // ラプラシアンフィルタをかける
                // |0  1 0|
                // |1 -4 1|
                // |0  1 0|
                float k = (2.0 * c.r) - c.g;
                float p = (k + _S2 * (
                    tex2D(_SelfTexture2D, uv - duv.zy).r +
                    tex2D(_SelfTexture2D, uv + duv.zy).r +
                    tex2D(_SelfTexture2D, uv - duv.xz).r +
                    tex2D(_SelfTexture2D, uv + duv.xz).r - 4.0 * c.r
                )) * _Atten;

                // 現在の状態をテクスチャのR成分に、ひとつ前の(過去の)状態をG成分に書き込む。
                return float4(p, c.r, 0, 0);
            }
            ENDCG
        }
    }
}

これを適切にCustomRenderTextureにセットすると以下のようにシミュレーションが進んでいきます。

f:id:edo_m18:20180823084141g:plain

CustomRenderTextureのセットアップ

シェーダが書けたところで、実際に使用するためにセットアップを行っていきます。

CustomRenderTextureの生成

まずはCustomRenderTextureを生成します。
生成には以下のようにメニューから生成することができます。

f:id:edo_m18:20180823091039p:plain

CustomRenderTexutreの設定

生成したCustomRenderTextureのインスペクタを見ると以下のような設定項目があります。

f:id:edo_m18:20180823091057p:plain

設定項目抜粋

いくつかの設定項目について見ていきます。

Size

CustomRenderTextureのサイズです。

Color Format

今回は凹みさんの記事にならって「RG Float」としています。

これは、現在の波の高さをR要素に、ひとつ前の波の高さをG要素に格納して計算を行うためです。

Material

CustomRenderTextureで使用するシェーダを適用したマテリアルを設定します。
サブ項目に「Shader Pass」がありますが、これは複数パスを定義した場合に、どのパスを利用するか選択できるようになっています。

Initialize Mode

初期化モード。OnLoadとすることで、起動時に初期化されます。
スクリプトから制御する場合は「OnDemand」を選択します。

サブ項目として「Source」があります。
設定できる内容は「Texture and Color」か「Material」があります。
今回は初期化をテクスチャから指定しています。
その場合は初期化に利用するテクスチャを設定する必要があります。

Update Mode

テクスチャの更新に関する設定です。
サンプルでは「Realtime」を設定しています。
こうすることで、スクリプトなしに自動的に更新が行われるようになります。

なお、スクリプトから制御したい場合はこちらも「OnDemand」を選択します。

Double Buffered

ダブルバッファのオン/オフ。
これをオンにすることで、テクスチャの情報をピンポンするように更新することができるようになります。

バンプマップに関して

以上でCustomRenderTextureに関する話はおしまいです。

今回の例ではシミュレーションした波形を利用してバンプマップを行っています。
ただ、シミュレーションした結果はRG要素しかない上に、それぞれのチャンネルはHeightMapとなっているためそのままでは法線マップとして利用できません。

そのため、取得したHeightMapから法線を計算して利用する必要があります。

バンプマップとは

法線マップを使ってポリゴンに凹凸があるように見せる処理です。
以前、Qiitaの記事にも書いたのでよかったら見てみてください。
なお、今回は法線をハイトマップから計算で求めて適用するための方法について書いていきます。

qiita.com

法線をハイトマップから計算で求める

ということで、法線の計算についても書いておきたいと思います。
法線の計算にあたっては、以下のふたつの記事を参考にさせていただきました。

t-pot『動的法線マップ』

esprog.hatenablog.com

法線は接ベクトルに垂直

接ベクトルをWikipediaで調べると以下のようにあります。

数学において、接ベクトル(英: tangent vector)とは、曲線や曲面に接するようなベクトルのことである。

そして、参考にさせていただいた記事(t-pot『動的法線マップ』)から説明を引用させていただくと、

法線ベクトルの求め方ですが、接ベクトルが法線ベクトルに直行することを利用します。 今回の場合は、高さ情報だけしか変化しないので、X軸及びZ軸の方向に関する高さの偏微分が接ベクトルの方向になります。

ということ。
つまり、接ベクトルをX軸およびZ軸に対して求め、それの外積を取ることで法線を求める、ということです。

図にすると以下のようなイメージです。

f:id:edo_m18:20180823122548p:plain

数式にすると以下のようになります。


dy_x = \frac{(y\_{i+1j} - y\_{ij}) + (y\_{ij} - y\_{i-1j})}{2} = \frac{y\_{i+1j} - y\_{i-1j}}{2} \\\
dy_z = \frac{(y\_{ij+1} - y\_{ij}) + (y\_{ij} - y\_{ij-1})}{2} = \frac{y\_{ij+1} - y\_{ij-1}}{2}

よって、これを利用して以下のように法線を求めることができます。

float3 du = (1, dyx, 0)
float3 dv = (0, dyz, 1)
float3 n = normalize(cross(dv, du))

実装

コードにすると以下のようになります。(一部抜粋)

// _ParallaxMap_TexelSizeは、テクスチャサイズの逆数
// テクセルの「ひとつ隣(シフト)」分の値を計算する
float2 shiftX = float2(_ParallaxMap_TexelSize.x, 0);
float2 shiftZ = float2(0, _ParallaxMap_TexelSize.y);

// 現在計算中のテクセルの上下左右の隣のテクセルを取得
float3 texX = tex2D(_ParallaxMap, float4(i.uv.xy + shiftX, 0, 0)) * 2.0 - 1;
float3 texx = tex2D(_ParallaxMap, float4(i.uv.xy - shiftX, 0, 0)) * 2.0 - 1;
float3 texZ = tex2D(_ParallaxMap, float4(i.uv.xy + shiftZ, 0, 0)) * 2.0 - 1;
float3 texz = tex2D(_ParallaxMap, float4(i.uv.xy - shiftZ, 0, 0)) * 2.0 - 1;

// 偏微分により接ベクトルを求める
float3 du = float3(1, (texX.x - texx.x), 0);
float3 dv = float3(0, (texZ.x - texz.x), 1);

// 接ベクトルの外積によって「法線」を求める
float3 n = normalize(cross(dv, du));

実装は難しい点はありません。
上下の高さ差分と、左右の高さ差分からそれぞれ勾配(接ベクトル)を求め、その2つのベクトルから外積を求めるだけです。
法線は接ベクトルと垂直なので、これで無事、法線が計算できたことになりますね。

最後に、バンプマップを行ったシェーダの全文を載せておきます。

Shader "Unlit/HightMapNormal"
{
    Properties
    {
        _Color("Tint color", Color) = (1, 1, 1, 1)
        _MainTex("Texture", 2D) = "white" {}
        _ParallaxMap("Parallax Map", 2D) = "gray" {}
    }

        SubShader
    {
        Tags { "RenderType" = "Transparent" "Queue" = "Transparent" }
        LOD 100

        GrabPass { }

        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #include "UnityCG.cginc"
            #include "Lighting.cginc"

            struct appdata
            {
                float4 vertex : POSITION;
                float2 uv : TEXCOORD0;
            };

            struct v2f
            {
                float2 uv : TEXCOORD0;
                float4 vertex : SV_POSITION;
                float4 uvgrab : TEXCOORD1;
                float3 worldPos : TEXCOORD2;
            };

            sampler2D _MainTex;
            sampler2D _ParallaxMap;
            sampler2D _GrabTexture;
            float4 _MainTex_ST;
            fixed4 _Color;

            float2 _ParallaxMap_TexelSize;

            v2f vert(appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.worldPos = mul(UNITY_MATRIX_MV, v.vertex);
                o.uv = TRANSFORM_TEX(v.uv, _MainTex);

#if UNITY_UV_STARTS_AT_TOP
                float scale = -1.0;
#else
                float scale = 1.0;
#endif
                o.uvgrab.xy = (float2(o.vertex.x, o.vertex.y * scale) + o.vertex.w) * 0.5;
                o.uvgrab.zw = o.vertex.zw;

                return o;
            }

            fixed4 frag(v2f i) : SV_Target
            {
                float2 shiftX = float2(_ParallaxMap_TexelSize.x, 0);
                float2 shiftZ = float2(0, _ParallaxMap_TexelSize.y);

                float3 texX = tex2D(_ParallaxMap, float4(i.uv.xy + shiftX, 0, 0)) * 2.0 - 1;
                float3 texx = tex2D(_ParallaxMap, float4(i.uv.xy - shiftX, 0, 0)) * 2.0 - 1;
                float3 texZ = tex2D(_ParallaxMap, float4(i.uv.xy + shiftZ, 0, 0)) * 2.0 - 1;
                float3 texz = tex2D(_ParallaxMap, float4(i.uv.xy - shiftZ, 0, 0)) * 2.0 - 1;

                float3 du = float3(1, (texX.x - texx.x), 0);
                float3 dv = float3(0, (texZ.x - texz.x), 1);

                float3 n = normalize(cross(dv, du));

                i.uvgrab.xy += n * i.uvgrab.z;

                fixed4 col = tex2Dproj(_GrabTexture, UNITY_PROJ_COORD(i.uvgrab)) * _Color;

                float3 lightDir = normalize(_WorldSpaceLightPos0 - i.worldPos);
                float diff = max(0, dot(n, lightDir)) + 0.5;
                col *= diff;

                float3 viewDir = normalize(_WorldSpaceCameraPos - i.worldPos);
                float NdotL = dot(n, lightDir);
                float3 refDir = -lightDir + (2.0 * n * NdotL);
                float spec = pow(max(0, dot(viewDir, refDir)), 10.0);
                col += spec + unity_AmbientSky;

                return col;
            }
            ENDCG
        }
    }
}

このシェーダのマテリアルを平面に適用することで、冒頭の動画の効果を得ることができます。

参考記事

その他、参考になった記事を載せておきます。

esprog.hatenablog.com

esprog.hatenablog.com

カールノイズを使ったパーティクル表現

f:id:edo_m18:20180726190953p:plain

概要

以前、カールノイズについてふたつの記事を書きました。
カールノイズの「流体表現」についてと「衝突判定」についてです。

edom18.hateblo.jp

edom18.hateblo.jp

今回はこれを発展させて、上記のカールノイズを使った具体的なパーティクル表現について書きたいと思います。
実際に動いている様子はこんな感じです↓

移動量に応じて自動でエミットするタイプ↓

モデルの頂点位置からエミットするタイプ↓

上の例では移動距離に応じて発生させるパーティクルと、メッシュの頂点部分にパーティクルを発生させる2パターンを実装しました。

なお、今回のサンプルはGithubにアップしてあります。

github.com

ちなみに、後者の「頂点部分にパーティクル」というのは、ドラクエ11リレミトの表現をイメージしてみました↓


【ドラクエ11】目指せ最速クリア!ぶっ続けでプレイして世界を救う配信!【ネタバレ注意】

前回の実装ではパーティクルはすべて常にアップデートされていた

前の実装ではすべてのパーティクルが生きていて常に更新がかかる実装になっていました。
以前に投稿した動画はこんな感じです↓

見てもらうと分かりますがパーティクルが発生し続けていて「任意のタイミングで追加する」というケースがなかったわけですね。
つまりパーティクルは常に発生し続けライフタイムが0になったらまた新しく生成し直す、というサイクルで実装していました。

具体的には、ライフタイムがゼロになったらライフタイムを回復させ、位置を初期位置(ランダム性あり)として更新していました。

なので毎フレームパーティクルのデータを更新すればよく、休止中のパーティクルを起こす必要もなければ状態を管理する必要もなかったわけです。
しかし移動距離で発生させたりなど、「任意のタイミングで追加」する必要がある場合は「活動中」と「休止中」のパーティクルを管理し、追加の場合は休止中のパーティクルに対して処理を実行する必要が出てきます。

連続してパーティクルをエミットする

冒頭で紹介したパーティクル表現では(Unity標準のパーティクルシステムでも同様の機能がありますが)「移動した距離に応じてエミットする」という実装になっています。
つまり、前に発生させたパーティクルは維持しつつ、新しく追加でパーティクルを発生させる必要がある、というわけです。

パーティクルをプールして管理する

今回の実装ではパーティクルの状態を管理し、必要であれば追加でエミットする必要があることは書きました。
ではそれをどうやったら実現できるのでしょうか。

先に結論を書いてしまうと、休止しているパーティクルのIDを保持するバッファを用意し、休止になったパーティクルのIDはプールへ戻し、追加の必要が出た場合はそのプールからIDを取り出し利用する、という仕組みを作ります。

この実装にあたっては、凹みTipsの以下の記事を参考にさせていただきました。

tips.hecomi.com

まずはザッと今回新しく追加した処理を見てみます。

プールの状態を初期化するInitカーネル

Initカーネルで行っているのは、用意されたパーティクルバッファ分の初期化です。
処理はシンプルに、全パーティクルの非活性化およびそのIDのプールへの保存です。

////////////////////////////////////////////////////////////
///
/// 初期化処理のカーネル関数
///
[numthreads(8, 1, 1)]
void Init(uint id : SV_DispatchThreadID)
{
    _Particles[id].active = false;
    _DeadList.Append(id);
}

上記の_DeadList.Append(id);がプールへIDを保存している箇所ですね。
_DeadListの名前の通り、非活性化状態のパーティクルのIDを保存しているバッファです。

さて、このAppendを実行しているバッファはなんでしょうか。答えは以下です。

Append Buffer、Consume Bufferを利用してプールを管理する

Append BufferConsume Bufferはそれぞれ、追加・取り出し可能なLIFO(Last In First Out)コンテナです。(つまり、最後に入れた要素が最初に取り出されるタイプのコンテナです)

MSDNのドキュメントは以下です。

docs.microsoft.com

docs.microsoft.com

AppendStructuredBufferのドキュメントから引用すると、

Output buffer that appears as a stream the shader may append to. Only structured buffers can take T types that are structures.

シェーダから追加することができるストリームとしてのバッファ、とあります。
このバッファに対して、寿命が来たパーティクルをバッファに入れてあげることで死活管理ができるというわけです。

ちなみにシェーダ内では以下のように定義しています。

AppendStructuredBuffer<uint> _DeadList;
ConsumeStructuredBuffer<uint> _ParticlePool;

さて、ふたつのバッファを定義しているのでふたつのバッファを用意する必要があるのか、と思われるかもしれませんが、インターフェースが違うだけでバッファとしての実態は同じものを利用します。

CPU(C#)側での処理は以下のようになっています。

private ComputeBuffer _particlePoolBuffer;
// ---------------------
_particlePoolBuffer = new ComputeBuffer(_particleNumLimit, sizeof(int), ComputeBufferType.Append);
// ---------------------
_computeShader.SetBuffer(_curlnoiseKernel, _deadListId, _particlePoolBuffer);
// ---------------------
_computeShader.SetBuffer(_emitKernel, _particlePoolId, _particlePoolBuffer);

まず、ComputeBuffer型の_particlePoolBufferを定義し、タイプをComputeBufferType.Appendとしてバッファを生成しています。

その後SetBufferによってバッファをセットしていますが、見て分かる通り実際に渡しているバッファの参照はどちらも同じ_particlePoolBufferです。

つまり、追加する場合はAppendStructuredBufferとして参照を渡し、取り出すときはConsumeStructuredBufferとして参照を渡す、というわけですね。

これで、寿命が尽きたパーティクルはバッファ(リスト)に戻され、必要になったときに、非活性化しているパーティクルのIDを「取り出して」値を設定することで、今回のように追加でエミットできるようになる、というわけです。

パーティクルを追加でエミットする

さて、仕組みが分かったところで実際にエミットしているところを見てみましょう。

////////////////////////////////////////////////////////////
///
/// パーティクルをエミットさせるカーネル関数
///
[numthreads(8, 1, 1)]
void Emit()
{
    uint id = _ParticlePool.Consume();

    float2 seed = float2(id + 1, id + 2);
    float3 randomPosition = rand3(seed);

    Particle p = _Particles[id];

    p.active = true;
    p.position = _Position + (randomPosition * 0.05);
    p.velocity = _Velocity;
    p.color = _Color;
    p.scale = _Scale;
    p.baseScale = _BaseScale;
    p.time = 0;
    p.lifeTime = randRange(seed + 1, _MinLifeTime, _MaxLifeTime);
    p.delay = _Delay;

    _Particles[id] = p;
}

冒頭のuint id = _ParticlePool.Consume();がIDを取り出しているところですね。
こうして、非活性化しているパーティクルのIDを取り出し、取り出したパーティクルに対して新しくライフタイムや色、位置などを設定することで無事、該当のパーティクルが動き出す、というわけです。

パーティクルの死活管理

さて、パーティクル生成の最後に実際にパーティクルをアップデートしている箇所を見てみます。
前回の実装ではライフタイムがゼロになった場合はまた新しいライフタイムを即座に設定し、すぐに復活させていました。

が、今回は「非活性化状態リスト」に追加する必要があります。
該当の処理は以下になります。

////////////////////////////////////////////////////////////
///
/// カールノイズのカーネル関数
///
[numthreads(8, 1, 1)]
void CurlNoiseMain(uint id : SV_DispatchThreadID)
{
    // ... 前略

    if (p.active)
    {
        // ... 中略

        if (p.time >= p.lifeTime)
        {
            p.active = false;
            _DeadList.Append(id);
            p.scale = 0;
        }
    }

    // ... 後略
}

ライフタイムがゼロ以下になった場合は_DeadList.Append(id);として、非活性化リストに戻しているのが分かるかと思います。

これで晴れて、パーティクルの死活管理ができる、というわけです。

パーティクルの残り数を考慮する

以上でパーティクルを追加でエミットすることが可能になりました。
しかし今のままだと少し問題が残っています。

というのは、パーティクルの残数を管理する必要があることです。
パーティクルの計算にはComputeBufferによって渡されたパーティクル情報を元に計算を行います。
このとき、初期化のタイミングでパーティクル数を規定します。(バッファサイズを明示する必要があります)
サンプルでは30000という数字でバッファを生成しています。

つまり、(当然ですが)このサイズを超えるパーティクルは生成することができません。
しかし、今回の追加エミットの処理ではその上限を超えてパーティクルの追加を要求することができてしまいます。

当然、バッファサイズを超えた数をリクエストしてもプールには休止中のIDがなく、結果生成されることがありません。
自分が遭遇した問題としては、上限を超えてリクエストを行ってしまうとプールに戻るパーティクルがなくなり、結果的にパーティクルが再生されない、という問題がありました。
しかもその状態が発生してしまうと以後、まったくパーティクルが発生しなくなる、という問題に遭遇しました。

なので「上限を超えないように」エミットを調整しないとならない、というわけです。

現在のプールに残っているIDの数を取得する

ではどうするのか。
結論から先に書いてしまうと、プールに残っているIDの数を取得して調整を行います。

プールのカウントを取得するにはComputeBuffer.CopyCountを利用します。
取得、確認しているコードは以下です。

_particleArgsBuffer = new ComputeBuffer(1, sizeof(int), ComputeBufferType.IndirectArguments);
_particleArgs = new int[] { 0 };

// --------------------

_particleArgsBuffer.SetData(_particleArgs);
ComputeBuffer.CopyCount(_particlePoolBuffer, _particleArgsBuffer, 0);

_particleArgsBuffer.GetData(_particleArgs);

return _particleArgs[0];

最終的に_particleArgsint型の配列に結果が格納されます。
今回のサンプルでは30000のパーティクルを生成しているので、まったくの未使用状態だと_particleArgs[0] == 30000となります。

MSDNのドキュメント↓
ComputeBuffer.CopyCount - Unity スクリプトリファレンス

あとは、ここで取得した数と実際にエミットしたい数とを比較すればバッファ以上のパーティクルを生成してしまうのを防ぐことができるようになります。

ジオメトリシェーダ入門

概要

今回はジオメトリシェーダ入門という形で記事を書いていきます。
頂点シェーダ、フラグメントシェーダに比べるとあまり書く機会の少ないシェーダ。
その2シェーダに比べて違いが結構あるのでそのあたりについてのメモを書いていこうと思います。

実際に使った例はこんな感じ↓
f:id:edo_m18:20180711140346p:plain

ジオメトリシェーダとは

Wikipediaによると以下のように説明されています。

ジオメトリシェーダー
ジオメトリシェーダー(英: Geometry Shader, GS)はピクセルシェーダーに渡されるオブジェクト内の頂点の集合を加工するために使用される。ジオメトリシェーダーにより、実行時に頂点数を増減させたり、プリミティブの種類を変更したりすることが可能となる。OpenGLではプリミティブシェーダーとも呼ばれる。

ジオメトリシェーダーはポイント、ライン、トライアングルといった既存のプリミティブから新しいプリミティブを生成できる。

ジオメトリシェーダーは頂点シェーダーの後に実行され、プリミティブ全体または隣接したプリミティブの情報を持つプリミティブを入力する。例えばトライアングルを処理するとき、3つの頂点がジオメトリシェーダーの入力となる。ジオメトリシェーダーはラスタライズされるプリミティブを出力でき、そのフラグメントは最終的にピクセルシェーダーに渡される。またプリミティブを出力せずにキャンセルすることもできる。

ジオメトリシェーダーのよくある使い方としては、ポイントスプライトの生成、ジオメトリテセレーション、シャドウボリュームの切り出し、キューブマップあるいはテクスチャ配列へのシングルパスレンダリングなどがある。

ざっくりと要約すると、頂点単位ではなく、ポリゴン単位などで値を計算し、さらにはポリゴン枚数すら増やすことができるシェーダ、という感じでしょうか。

テッセレーションは分かりやすい使用例ですね。
その他の使い方としては、頂点郡からポリゴン(Quad)を生成しテクスチャを貼ることでパーティクルのような演出をする、なども使い道がありそうです。

まさにその発想で解説してくれている記事があるので紹介↓

wordpress.notargs.com

コードを見てみる

さて、ジオメトリシェーダのシンプルなコードを見てみましょう。
下の例ではジオメトリシェーダではある意味で「なにもしていません」。
つまり、ポリゴンは変形しません。

が、ジオメトリシェーダがなにをしてくれているのかを見るのにちょうどいいと思います。

appdata vert(appdata v)
{
    // ジオメトリシェーダで処理するため、頂点シェーダへの入力をそのまま返す
    return v;
}

// ジオメトリシェーダでは3頂点のinputと3頂点のoutputを行う=つまり普通の三角ポリゴン
// max vertex count 3が3頂点のoutputであることを伝えている
// また、triangle appdata input[3]が、「三角形」で3頂点のinputが必要であることを伝えている
[maxvertexcount(3)]
void geom(triangle appdata input[3], inout TriangleStream<vsout> outStream)
{
    [unroll]
    for (int i = 0; i < 3; i++)
    {
        // 頂点シェーダからもらった3頂点それぞれを射影変換して通常のレンダリングと同様にポリゴン位置を決める
        appdata v = input[i];
        vsout o;
        o.vertex = UnityObjectToClipPos(v.vertex);
        outStream.Append(o);
    }

    outStream.RestartStrip();
}

fixed4 frag(vsout i) : SV_Target
{
    return float4(1, 0, 0, 1);
}

頂点シェーダやフラグメントシェーダに比べるとやや特殊な感じになっています。
特に、関数に対してアトリビュートが指定されています。
細かな点についてはコメントで記載しました。

MSDNのドキュメント(ジオメトリ シェーダー オブジェクト (DirectX HLSL).aspx))を見ると以下のように書かれています。

ドキュメント

Syntax

[maxvertexcount(NumVerts)]
void ShaderName (
  PrimitiveType DataType Name [ NumElements ],
  inout StreamOutputObject
);

パラメーター

  • [maxvertexcount(NumVerts)]
    • [in] 作成する頂点の最大数の宣言です。
      • [maxvertexcount()] - 必須のキーワードです。正しい構文とするために、山かっこと丸かっこを必ず記述します。
    • NumVerts - 頂点の数を表す整数です。
  • ShaderName
    • [in] ジオメトリ シェーダー関数の一意の名前を含む ASCII 文字列。
  • PrimitiveType DataType Name [ NumElements ]
    • PrimitiveType - プリミティブ型。プリミティブ データの順序を決定します。
    • DataType - [in] 入力データ型。任意の HLSL データ型を指定できます。
    • Name - ASCII 文字列の引数名。引数が配列の場合は、NumElements で配列サイズを任意に指定できます。
プリミティブの種類 説明
point ポイント リスト
line ライン リストまたはライン ストリップ
triangle トライアングル リストまたはトライアングル ストリップ
lineadj 隣接性のあるライン リストまたは隣接性のあるライン ストリップ
triangleadj 隣接性のあるトライアングル リストまたは隣接性のあるトライアングル ストリップ

ストリーム出力オブジェクト

ここで、「ストリーム出力オブジェクト」は出力されるオブジェクトの種類に応じて指定するものが変わります。(例えばポリゴンとして出力するのか、ラインとして出力するのか、など)

ドキュメントから引用すると、以下の3つが定義されているようです。

Syntax
inout StreamOutputObject<DataType> Name;
パラメーター
  • StreamOutputObject Name
    • ストリーム出力オブジェクト (SO) 宣言
ストリーム出力オブジェクトの型 説明
PointStream ポイント プリミティブのシーケンス
LineStream ライン プリミティブのシーケンス
TriangleStream トライアングル プリミティブのシーケンス

違ったコードを見てみる

次に、ちょっとだけ違ったコードを見てみます。
これはポリゴン単位で法線を計算している例です。

通常の頂点シェーダでは頂点ごとの計算しかできないため、実際に生成されたポリゴンに対しての処理は書けません。
頂点だけではポリゴンを形成する「隣接する」他の頂点情報がないために計算が行えないためですね。

しかしジオメトリシェーダでは「ポリゴン単位」で処理が行えるため、こうしたことが可能になっています。

[maxvertexcount(3)]
void geom(triangle appdata input[3], inout TriangleStream<vsout> outStream)
{
    float3 edge1 = (input[1].vertex.xyz - input[0].vertex.xyz);
    float3 edge2 = (input[2].vertex.xyz - input[0].vertex.xyz);
    float3 normal = normalize(cross(edge1, edge2));

    [unroll]
    for (int i = 0; i < 3; i++)
    {
        appdata v = input[i];

        vsout o;

        o.pos = UnityObjectToClipPos(v.vertex);
        o.normal = normal;
        outStream.Append(o);
    }

    outStream.RestartStrip();
}

上記の例では、ポリゴンのエッジから法線を計算しています。
こんな感じで「ポリゴン単位」で計算が行えるのが特徴です。

また上記で紹介した記事のように、inputの数とouputの数を変えることができるので、1枚のポリゴンから数枚のポリゴン、というようにポリゴン数を増やすことも可能です。

実例を見てみる

さて最後に、実際に利用されているシェーダを覗いてみて終わりにします。
以下の例は、凹みさんが書かれている「凹みTips」の記事(HoloLens で使える Near Clip 表現について解説してみた)で紹介されているものです。
本人の承諾をいただいて解説させていただいています。

※ 以下のコードは凹みさんが公開しているものを若干変更したものになります。

凹みさんが作った実際のものを動かすと以下のような感じになります。
(凹みさんの画像を引用させていただきました)

ポリゴン分解の動作イメージ

ジオメトリシェーダの部分だけ抜粋しました↓
ポイントを絞ってコード内にコメントを記載しています。

// 前のサンプル同様、3頂点のinput / outputが指定されています。
[maxvertexcount(3)]
void geom(triangle appdata_t input[3], inout TriangleStream<g2f> stream)
{
    // ポリゴンの中心を計算。
    // ポリゴン単位で計算を行えるため、「ポリゴンの中心位置」も計算可能です。
    float3 center = (input[0].vertex + input[1].vertex + input[2].vertex) / 3;

    // ポリゴンの辺ベクトルを計算し、ポリゴンの法線を計算する。
    // 続いて、前のサンプルでもあった「ポリゴン法線」の計算です。
    float3 vec1 = input[1].vertex - input[0].vertex;
    float3 vec2 = input[2].vertex - input[0].vertex;
    float3 normal = normalize(cross(vec1, vec2));

    #ifdef _METHOD_PROPERTY
    // カメラ視点を利用しない場合はパラメータをそのまま利用する
    fixed destruction = _Destruction;
    #else
    // カメラ視点からの距離によって分解する場合は、中心位置をワールド座標に変換し、
    // カメラとの距離を計算して係数に利用する
    float4 worldPos = mul(Unity_ObjectToWorld, float4(center, 1.0));
    float3 dist = length(_WorldSpaceCameraPos - worldPos);
    fixed destruction = clamp((_StartDistance - dist) / (_StartDistance - _EndDistance), 0.0, 1.0);
    #endif

    // 省略していますが、独自で定義した「rand」関数を使って乱数を生成しています。
    // ここではポリゴン位置などをseedにして乱数を生成しています。
    fixed r = 2.0 * (rand(center.xy) - 0.5);
    fixed3 r3 = r.xxx;
    float3 up = float3(0, 1, 0);

    // コードをループじゃない状態に展開することを明示する(詳細は下の記事を参照)
    [unroll]
    for (int i = 0; i < 3; i++)
    {
        appdata_t v = input[i];
        g2f o;

        UNITY_SETUP_INSTANCE_ID(v);
        UNITY_INITIALIZE_VERTEX_OUTPUT_STEREO(o);

        // 以下では、各要素(位置、回転、スケール)に対して係数に応じて変化を与えます。

        // center位置を起点にスケールを変化させます。
        v.vertex.xyz = (v.vertex.xyz - center) * (1.0 - destruction * _ScaleFactor) + center + (up * destruction);

        // center位置を起点に、乱数を用いて回転を変化させます。
        v.vertex.xyz = rotate(v.vertex.xyz - center, r3 * destruction * _RotationFactor) + center;

        // 法線方向に位置を変化させます
        v.vertex.xyz += normal * destruction * _PositionFactor * r3;

        // 最後に、修正した頂点位置を射影変換しレンダリング用に変換します。
        o.vertex = UnityObjectToClipPos(v.vertex);

        #ifdef SOFTPARTICLES_ON
        o.projPos = ComputeScreenPos(o.vertex);
        COMPUTE_EYEDEPTH(o.projPos.z);
        #endif

        o.color = v.color;
        o.color.a *= 1.0 - destruction * _AlphaFactor;
        o.texcoord = TRANSFORM_TEX(v.texcoord, _MainTex);
        UNITY_TRANSFER_FOG(o, o.vertex);

        stream.Append(o);
    }

    stream.RestartStrip();
}

※ [unroll]については以下の記事がどうなるのか分かりやすいでしょう。

wlog.flatlib.jp

だいぶ駆け足ですが、ジオメトリシェーダの概要と実例を紹介しました。
最後の実例の詳細については凹みさんのブログ記事を参照してください。

tips.hecomi.com

ちなみにこうした動作を頂点シェーダで行ってしまうと、全頂点がむすばれてしまってトゲトゲした感じの見た目になるだけで、こうした動きは得られません。

凹みさんの記事には書いてありますが、ジオメトリシェーダに対応していないモバイルなどの環境では全頂点をポリゴン用に複製し、実際にポリゴンとして分解しておく必要があります。
こうした点でも、事前準備なしにこうした演出が行えるのはメリットですね。

実際に動いている動画

ちなみに以前投稿した以下の動画で動いているポリゴン分解は、凹みさんの記事を参考にさせてもらって実装したものです。SAO感が好きですw

OBB(Oriented Bounding Box)を計算する

概要

今回はOBB(Oriented Bounding Box)についてのまとめです。
今実装している中でボリュームのサイズをある程度正確に量りたいと思ったのが理由です。

今回の実装はこちらの記事を参考に、いくつかの記事を参照して行いました。

sssiii.seesaa.net

こちらの記事に書かれている手順を引用させていただくと以下のようになります。

  1. 頂点座標で行列を作る(xyz座標3データ×5000個のデータで3×5000の行列)
  2. xyz各行の平均を計算して引く
  3. 引いた行列の転置行列を作る(5000×3の行列)
  4. 掛けてNで割って共分散行列ができる((3×5000)×(5000×3)→3×3の行列)
  5. ヤコビ法などを用いて共分散行列の固有値固有ベクトルを計算する
  6. 固有値固有ベクトルが求まれば、固有値の大きい順に固有ベクトルを並べ替えます。
  7. 固有ベクトルを並べて3×3の行列を作っておきます。
  8. 固有ベクトルの行列と頂点の内積を計算し、各固有ベクトル毎に最大と最小の内積値を計算します。
  9. そのデータを使えば重心、OBBの座標が計算できます。(参考サイト参照)

ひとつひとつ見ていきましょう。

OBBとは

OBBとはOrientedと名前がついている通り、オブジェクトに対してできるだけ余白が出ないようにそれをすっぽりと覆うBoxを定義します。

それとは別にAABB(Axis Aligned Bounding Box)があります。Axis Alignedの名前の通り、座標空間のXYZ軸(Axis)に平行(Aligned)なBoxです。
つまり、オブジェクトが回転すると場合によっては無駄な余白ができてしまいます。

それぞれ図にすると以下のような感じの違いになります。

f:id:edo_m18:20180625140112p:plain
▲球体を半分にカットしてAABBを表示した状態

f:id:edo_m18:20180625140202p:plain
▲上記をOBBで表示した状態

いかがですか。オブジェクトの形状や回転状態によってだいぶ違いが出るのが分かるかと思います。

OBBの実装の概要についてはこちらの記事が参考になります。
記事から引用させていただくと、

簡単に言うと、あるポリゴン群に対してどんな方向を向いていてもかまわないので、それをすっぽりと囲む、できる限り小さな直方体を作りましょうというのがOBBの手法です。
この直方体をいかにして作るのかという所が腕の見せ所になります。
とりあえずここでは、上記論文で紹介されている主成分分析を使うことにしましょう。

と書かれています。

主成分分析とは

記事で書かれている主成分分析とはなんでしょうか。
Wikipediaによるとこう書かれています。

主成分分析(しゅせいぶんぶんせき、英: principal component analysis; PCA)とは、相関のある多数の変数から相関のない少数で全体のばらつきを最もよく表す主成分と呼ばれる変数を合成する多変量解析の一手法[1]。データの次元を削減するために用いられる。

主成分を与える変換は、第一主成分の分散を最大化し、続く主成分はそれまでに決定した主成分と直交するという拘束条件の下で分散を最大化するようにして選ばれる。主成分の分散を最大化することは、観測値の変化に対する説明能力を可能な限り主成分に持たせる目的で行われる。選ばれた主成分は互いに直交し、与えられた観測値のセットを線型結合として表すことができる。言い換えると、主成分は観測値のセットの直交基底となっている。主成分ベクトルの直交性は、主成分ベクトルが共分散行列(あるいは相関行列)の固有ベクトルになっており、共分散行列が実対称行列であることから導かれる。

なんだか分かったような分からないようなw
自分の理解としては、多数ある点群に対して方向性を見出し、それを計算によって求めた値を掛けることで分析しやすい形にする、といった感じでしょうか。

固有値固有ベクトルとは

主成分分析には「固有値固有ベクトル」を使って変換が行われるようです。
固有値固有ベクトルについては以下の記事がとても分かりやすかったです。

qiita.com

なんだかどんどん主題と離れていっている気がしなくもないですが、まぁせっかくなのでまとめておきましょうw

上で紹介した記事を見てもらうとアニメーションで示してくれているのでとても分かりやすいですが、自分の理解として言葉で説明しておくと。

とあるベクトルに対して、とある行列を掛けたとき、そのベクトルが回転せず定数倍になるベクトルのことを「固有ベクトル」、そしてその定数倍の定数が「固有値」、ということのようです。

数式で表すと以下になります。


Ax = \lambda x

※ ただし x \neq 0

 Aが行列で \lambdaが定数(スカラー)です。

行列での演算が定数(スカラー)の演算で表せる、というのが上記式の意味です。

(推測ですが)統計学や主成分分析に使える理由としては、傾向がある=回転しないベクトル(方向)、みたいなことから来ているんでしょうか。
ちょっと数学に疎いのでそのあたりの深い話は分かりません・・・。

閑話休題

さて、最初に紹介した記事を読み進めると以下のように続いています。

OBBではデータとはポリゴンの頂点群なので、x,y,zの3次元のデータが頂点の個数分存在します。 それを使って分散共分散行列を作って、その行列の固有ベクトルを計算すると、 OBBを表す3つの軸が得られます。

さぁ、また新しい単語が出てきましたw

分散共分散行列とは

ここでもWikipediaから引用させてもらうと、

統計学と確率論において分散共分散行列(ぶんさんきょうぶんさんぎょうれつ、英: variance-covariance matrix)や共分散行列(きょうぶんさんぎょうれつ、英: covariance matrix)とは、ベクトルの要素間の共分散の行列である。これは、スカラー値をとる確率変数における分散の概念を、多次元に自然に拡張したものである。

と書かれています。

「分散」と「共分散」があり、それを合成した行列のためこう呼ばれるようです。

ちなみに、「分散共分散行列」については以下の記事がとても分かりやすかったです。

zellij.hatenablog.com

分散

まずは分散について。
これまたWikipedia)から引用すると、

確率論および統計学において、分散(ぶんさん、英: variance)は、確率変数の2次の中心化モーメントのこと。これは確率変数の分布が期待値からどれだけ散らばっているかを示す非負の値である[1]。

つまり、確率変数がどれだけ「分散」しているか、ということかな。

また、上で紹介した記事を見てみると、以下のように説明されています。

具体的には、分散は「(各データの平均値からの距離)の2乗の平均」。
分散は2乗であることに注意。単位をそろえるために、分散の平方根を取ったものが標準偏差
標準偏差をσで表すと、分散はσ2で表される。

とのこと。なるほど。
ちなみに「標準偏差」については以下の記事が分かりやすかったです。

atarimae.biz

いわゆる「偏差値」とか、その類の話ですね。
二乗する辺りが「なるほど」と思わされます。
つまり、平均値からの距離を二乗することにより、大きく散らばっている(大きく分散している)場合ほど偏差が大きくなる、というわけです。(大きな値の二乗はより大きな数値になる)

そして完全に余談ですが、以前、「[数学] 最小二乗平面をプログラムで求める」という記事を書いていてそこでも似た話が出てくるので興味がある方は読んでみてください。

qiita.com

共分散

こちらもWikipediaから引用すると、

共分散(きょうぶんさん、covariance)は、2 組の対応するデータ間での、平均からの偏差の積の平均値である[1]。2 組の確率変数 X, Y の共分散 Cov(X, Y) は、E で期待値を表すことにして、 Cov(X, Y) = E[(X - E[X])(Y - E[Y])] で定義する。

上で、分かりやすいと紹介した記事を読んでの自分の理解を書くと、異なるデータに対して分散を計算し、どれくらい相関があるか、を計算したものが共分散かな、と思っています。

ちなみにその記事では以下のように例えています。非常に分かりやすい。

例えば、生徒の「数学の点数」と「英語の点数」がどのような関係にあるか知りたい。数学ができる生徒はやはり英語ができるのか?(正の相関)、それとも数学ができる生徒は英語が苦手なのか(負の相関)。

そこで、数学の点数(xの値)と英語の点数(yの値)という、2つのデータ群を考慮した分散を「共分散」と呼び、この共分散Sxyは次の式で表される。


S_{xy} = \frac{1}{n}(x - \overline{x})' (y - \overline{y})

 'は転置

つまり、XとYというふたつのデータがあった場合、上記式で求まるXとYの相関を見るのがXYの共分散となります。

さらにXとXの共分散が分散、ということのようです。
式にすると以下。


S_{xy} = \frac{1}{n}(x - \overline{x})' (x - \overline{x})

分散と共分散が分かりました。あとはこれを行列の形に合成すると「分散共分散行列」となります。
行列の形は以下のようになるようです。


S =
\begin{vmatrix}
S_{xx} &S_{xy} \\\
S_{xy} &S_{yy}
\end{vmatrix}

OBBの話に戻すと、「頂点郡」というバラけたデータから、XYZという3要素がどれくらい相関しているか、を主成分分析によって求める、みたいなことなのかなーと想像しています。(あくまで想像)

さて、では実際にどうこの分散共分散行列を求めるかと言うと、冒頭の記事の引用を再掲すると以下のように書かれています。

  1. 頂点座標で行列を作る(xyz座標3データ×5000個のデータで3×5000の行列)
  2. xyz各行の平均を計算して引く
  3. 引いた行列の転置行列を作る(5000×3の行列)
  4. 掛けてNで割って共分散行列ができる((3×5000)×(5000×3)→3×3の行列)

分散共分散行列の意味の解説を元に見てみるとやっていることはまさに同じことですね。

そして冒頭の記事が参考にしている記事には以下のようなコードが示されています。

Vertex配列に頂点データが入っていて、配列の長さがSizeになります。
mは全頂点の平均値です。

float C11 = 0, C22 = 0, C33 = 0, C12 = 0, C13 = 0, C23 = 0;
for( int i = 0; i < Size; ++i ) {
    C11 += ( Vertex[i].x - m.x ) * ( Vertex[i].x - m.x );
    C22 += ( Vertex[i].y - m.y ) * ( Vertex[i].y - m.y );
    C33 += ( Vertex[i].z - m.z ) * ( Vertex[i].z - m.z );
    C12 += ( Vertex[i].x - m.x ) * ( Vertex[i].y - m.y );
    C13 += ( Vertex[i].x - m.x ) * ( Vertex[i].z - m.z );
    C23 += ( Vertex[i].y - m.y ) * ( Vertex[i].z - m.z );
}

C11 /= Size;
C22 /= Size;
C33 /= Size;
C12 /= Size;
C13 /= Size;
C23 /= Size;

float Matrix[3][3] = {
    { C11, C12, C13 },
    { C12, C22, C23 },
    { C13, C23, C33 },
};

平均からの距離を求めて、それを行列にしている部分ですね。
ただ、いくつかの行列要素は同じ計算結果になるため必要な部分だけを計算して行列にしているのが分かります。

これで「分散共分散行列」を求めることができます。
次にここから「固有ベクトル」を求め、OBBの軸を求めます。

なお、コードを紹介してくれた記事では以下のように続いています。

こうして得られる分散共分散行列ことMatrix[3][3]は、
実対称行列なのでjacobi法で固有ベクトルを求めることができます。

固有ベクトルを求める

さて、だいぶ遠回りしてしまいましたが、ふんわりと、どういう処理を用いてOBBを計算するのかを見てきました。

最後に、求めたいのはOBBの軸と長さです。
3軸の向きと長さが分かればそこからOBBを求めることができますね。

ということで固有ベクトルを求める計算です。

jacobi法(ヤコビ法)で固有ベクトルを求める

前述の通り、分散共分散行列を求め、そこからjacobi法(ヤコビ法)という手法を用いて固有ベクトルを求めます。

jacobi法についてはこちらの記事(Jacobi法 - [物理のかぎしっぽ])を参考にさせていただきました。

該当記事の冒頭の解説を引用させてもらうと、

Jacobi法(ヤコビ法)を用いて行列の固有値固有ベクトルを求めてみたいと思います. Jacobi法では対称行列しか扱えないという制約があるものの,アルゴリズム自体も簡単な上に,解が必ず実数になる事もあり,大変シンプルです.

とあります。
そして分散共分散行列は対称行列のためこの方法が使える、というわけのようです。

さて、実際の求め方ですが、引用させていただくと、

 P^{-1} AP = B
となるとき,行列Aと行列Bの固有値は一致します.また,行列Aが対称行列である場合,直交行列Pを用いて

 P^{-1}AP = \Lambda
と,対角化することができ,行列Λの対角成分は行列Λの固有値となります.ここで直交行列Pの列ベクトルとそれに対応する行列Λの列の対角成分は固有値固有ベクトルの関係となります.

このことから行列Aの固有値固有ベクトルを求めるには,行列Aを対角化する直交行列Pを求めれば良いことになります.

Jacobi法ではこの直交行列Pを反復作業で求めています.つまり,行列Aの非対角成分の1つをゼロにする直交行列を用い,

 P^{-1}AP ⇒ A
とする事を繰り返します.

[tex: P^{-1}n … P^{-1}4 P^{-1}3 P^{-1}2 P^{-1}_1 AP_1 P_2 P_3 P_4 … P_n = \Lambda]
これを繰り返し全ての非対角成分をゼロにし,対角化をするのがJacobi法です.また,

 P_1 P_2 P_3 P_4 … P_n
この部分が行列Aを対角化する直交行列であり,列ベクトルは固有ベクトルとなります.

上記の説明の中で出てくる式や単語を整理しておきます。
直交行列」ですが、意味は逆行列と転置行列が等しくなる正方行列のことです。以下を満たす行列ですね。


M^T M = M M^T = E

また「対角化」を調べてみると以下の記事が分かりやすいです。

mathtrain.jp

記事から説明を引用すると、

与えられた正方行列 A に対して,正則行列 P をうまく取ってきて  P^{−1}AP を対角行列にする操作を対角化と言う。

ヤコビ法で語られていた以下の式が出てきました。


P^{-1} AP

ここで「正則行列Pをうまく取ってきて」とあります。
つまりこの「うまく取ってくる」というのがヤコビ法、ということになりますね。

さて実際のヤコビ法の計算方法ですが、詳細は紹介した記事を読んでみてください。
プログラムでの計算方法も掲載されています。

今回はそのプログラムコードを参考にさせていただき、Unity上で使えるようにC#で書き直しました。

C#コード

最後に、実際に書いたコードを掲載しておきます。よかったら参考にしてみてください。

以下のコードを貼り付けて実行すると、下図のようにOBBが表示されます。

f:id:edo_m18:20180627132122p:plain

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class OBBTest : MonoBehaviour
{
    private float _max = 0.0001f;

    private Vector3[] _edges = new Vector3[0];
    private Vector3 _origin;

    private void OnDrawGizmos()
    {
        if (_edges.Length == 0)
        {
            return;
        }

        Gizmos.color = Color.yellow;
        Gizmos.DrawWireSphere(_origin, 0.01f);

        Gizmos.color = Color.blue;

        Vector3 from, to;

        Vector3 ftr = _origin + (_edges[0] * 0.5f) + (_edges[1] * 0.5f) + (_edges[2] * 0.5f);
        Vector3 fbr = ftr - _edges[2];
        Vector3 ftl = ftr - _edges[0];
        Vector3 fbl = ftl - _edges[2];

        Vector3 btr = ftr - _edges[1];
        Vector3 bbr = btr - _edges[2];
        Vector3 btl = btr - _edges[0];
        Vector3 bbl = btl - _edges[2];

        Gizmos.DrawLine(ftr, fbr);
        Gizmos.DrawLine(ftr, ftl);
        Gizmos.DrawLine(ftl, fbl);
        Gizmos.DrawLine(fbl, fbr);

        Gizmos.DrawLine(btr, bbr);
        Gizmos.DrawLine(btr, btl);
        Gizmos.DrawLine(btl, bbl);
        Gizmos.DrawLine(bbl, bbr);

        Gizmos.DrawLine(ftr, btr);
        Gizmos.DrawLine(ftl, btl);
        Gizmos.DrawLine(fbr, bbr);
        Gizmos.DrawLine(fbl, bbl);
    }

    private void Awake()
    {
        CalcEdges();
    }

    private void Update()
    {
        if (Input.GetKeyDown(KeyCode.U))
        {
            CalcEdges();
        }
    }

    private void CalcEdges()
    {
        MeshFilter filter = GetComponent();

        Vector3[] vertices = filter.mesh.vertices;

        float[,] eigenVectors = GetEigenVectors(CollectMatrix(vertices));

        int rank = eigenVectors.Rank;

        Vector3 vec1, vec2, vec3;
        {
            float x = eigenVectors[0, 0];
            float y = eigenVectors[1, 0];
            float z = eigenVectors[2, 0];

            vec1 = new Vector3(x, y, z);
        }
        {
            float x = eigenVectors[0, 1];
            float y = eigenVectors[1, 1];
            float z = eigenVectors[2, 1];

            vec2 = new Vector3(x, y, z);
        }
        {
            float x = eigenVectors[0, 2];
            float y = eigenVectors[1, 2];
            float z = eigenVectors[2, 2];

            vec3 = new Vector3(x, y, z);
        }

        // 全頂点に対して内積を取り、最小値・最大値を計算する
        float min1 = float.MaxValue;
        float min2 = float.MaxValue;
        float min3 = float.MaxValue;
        float max1 = float.MinValue;
        float max2 = float.MinValue;
        float max3 = float.MinValue;

        for (int i = 0; i < vertices.Length; i++)
        {
            Vector3 pos = vertices[i];
            float dot1 = Vector3.Dot(vec1, pos);
            if (dot1 > max1)
            {
                max1 = dot1;
            }
            if (dot1 < min1)
            {
                min1 = dot1;
            }

            float dot2 = Vector3.Dot(vec2, pos);
            if (dot2 > max2)
            {
                max2 = dot2;
            }
            if (dot2 < min2)
            {
                min2 = dot2;
            }

            float dot3 = Vector3.Dot(vec3, pos);
            if (dot3 > max3)
            {
                max3 = dot3;
            }
            if (dot3 < min3)
            {
                min3 = dot3;
            }
        }

        float len1 = max1 - min1;
        float len2 = max2 - min2;
        float len3 = max3 - min3;

        Vector3 edge1 = transform.localToWorldMatrix.MultiplyVector(vec1 * len1);
        Vector3 edge2 = transform.localToWorldMatrix.MultiplyVector(vec2 * len2);
        Vector3 edge3 = transform.localToWorldMatrix.MultiplyVector(vec3 * len3);

        _edges = new[]
        {
            edge1, edge2, edge3,
        };

        Vector3 center1 = (vec1 * (max1 + min1)) * 0.5f;
        Vector3 center2 = (vec2 * (max2 + min2)) * 0.5f;
        Vector3 center3 = (vec3 * (max3 + min3)) * 0.5f;

        _origin = transform.localToWorldMatrix.MultiplyPoint(center1 + center2 + center3);
    }

    private float[,] CollectMatrix(Vector3[] vertices)
    {
        // 各成分の平均を計算
        Vector3 m = Vector3.zero;

        for (int i = 0; i < vertices.Length; i++)
        {
            m += vertices[i];
        }

        m /= vertices.Length;

        float c11 = 0; float c22 = 0; float c33 = 0;
        float c12 = 0; float c13 = 0; float c23 = 0;

        for (int i = 0; i < vertices.Length; i++)
        {
            c11 += (vertices[i].x - m.x) * (vertices[i].x - m.x);
            c22 += (vertices[i].y - m.y) * (vertices[i].y - m.y);
            c33 += (vertices[i].z - m.z) * (vertices[i].z - m.z);

            c12 += (vertices[i].x - m.x) * (vertices[i].y - m.y);
            c13 += (vertices[i].x - m.x) * (vertices[i].z - m.z);
            c23 += (vertices[i].y - m.y) * (vertices[i].z - m.z);
        }

        c11 /= vertices.Length;
        c22 /= vertices.Length;
        c33 /= vertices.Length;
        c12 /= vertices.Length;
        c13 /= vertices.Length;
        c23 /= vertices.Length;

        float[,] matrix = new float[,]
        {
            { c11, c12, c13 },
            { c12, c22, c23 },
            { c13, c23, c33 },
        };

        return matrix;
    }

    /// 
    /// 行列の中の絶対値の最大値とその位置を返す
    /// 
    /// 評価する行列
    /// 最大値の行位置
    /// 最大値の列位置
    /// 最大値
    private float GetMaxValue(float[,] matrix, out int p, out int q)
    {
        p = 0;
        q = 0;

        int rank = matrix.Rank;

        float max = float.MinValue;

        for (int i = 0; i < rank; i++)
        {
            int len = matrix.GetLength(i);

            for (int j = 0; j < len; j++)
            {
                // 対角成分は評価しない
                if (i == j)
                {
                    continue;
                }

                float absmax = Mathf.Abs(matrix[i, j]);
                if (max <= absmax)
                {
                    max = absmax;
                    p = i;
                    q = j;
                }
            }
        }

        if (p > q)
        {
            int temp = p;
            p = q;
            q = temp;
        }

        return max;
    }

    /// 
    /// 固有ベクトルを得る
    /// 
    /// 評価する行列
    private float[,] GetEigenVectors(float[,] matrix)
    {
        // 固有ベクトルのための行列を正規化
        float[,] eigenVectors = new float[3,3];

        for (int i = 0; i < 3; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                if (i == j)
                {
                    eigenVectors[i,j] = 1f;
                }
                else
                {
                    eigenVectors[i,j] = 0;
                }
            }
        }

        int limit = 100;
        int count = 0;
        int p, q;
        while (true)
        {
            count++;

            if (count >= limit)
            {
                Debug.Log("somothing was wrong.");
                break;
            }

            float max = GetMaxValue(matrix, out p, out q);
            if (max <= _max)
            {
                break;
            }

            float app = matrix[p, p];
            float apq = matrix[p, q];
            float aqq = matrix[q, q];

            float alpha = (app - aqq) / 2f;
            float beta = -apq;
            float gamma = Mathf.Abs(alpha) / Mathf.Sqrt(alpha * alpha + beta * beta);

            float sin = Mathf.Sqrt((1f - gamma) / 2f);
            float cos = Mathf.Sqrt((1f + gamma) / 2f);

            if (alpha * beta < 0)
            {
                sin = -sin;
            }

            for (int i = 0; i < 3; i++)
            {
                float temp = cos * matrix[p, i] - sin * matrix[q, i];
                matrix[q, i] = sin * matrix[p, i] + cos * matrix[q, i];
                matrix[p, i] = temp;
            }

            for (int i = 0; i < 3; i++)
            {
                matrix[i, p] = matrix[p, i];
                matrix[i, q] = matrix[q, i];
            }

            matrix[p, p] = cos * cos * app + sin * sin * aqq - 2 * sin * cos * apq;
            matrix[p, q] = sin * cos * (app - aqq) + (cos * cos - sin * sin) * apq;
            matrix[q, p] = sin * cos * (app - aqq) + (cos * cos - sin * sin) * apq;
            matrix[q, q] = sin * sin * app + cos * cos * aqq + 2 * sin * cos * apq;

            for (int i = 0; i < 3; i++)
            {
                float temp = cos * eigenVectors[i, p] - sin * eigenVectors[i, q];
                eigenVectors[i, + q] = sin * eigenVectors[i, p] + cos * eigenVectors[i, q];
                eigenVectors[i, + p] = temp;
            }
        }

        return eigenVectors;
    }
}