これまでのあらすじ

私はQ&Aサイトのteratailを覗いて、クイズのつもりで時折回答するのをささやかな趣味の一つにしております。
以前「unity 3D 切り取り」とのご質問を見かけまして、それに対して「モデルをTexture3Dを使ったボクセル的なデータとして表現し、それを加工することで削り取る」というのはどうかと提案いたしました。

ご質問者さんにはひとまず回答として受け入れていただけたようなのですが、解像度が低いと見た目がガタガタ、かといって解像度を上げるとエラーが発生する...といった問題点のご指摘を受けました。高解像度化の案として、別の回答者さんから八分木構造を使ってみてはとのコメントをいただき、ご質問者さんも新たに「Unity 3D Octree関連の実装」との質問を投稿されたご様子でした。
恥ずかしながら調べたことのない分野で、私には荷が重い...と感じ挑戦は控えていたのですが、興味深い話題で勉強になりそうだとの思いもありまして、実験用シーンを作ってチビチビいじっておりました。

現時点では速度が遅すぎて彫刻なんてできない状態ではありますが、Sparse Voxel Octree表現の一例として多少なりともご参考になれば...と思い、ひとまずここで区切りとして結果を報告させていただこうと考えました（さすがにそろそろモアイとにらめっこするのにもくたびれてしまったのです...ご容赦ください）。

本来であれば直接teratail上に投稿するべきなのでしょうが、コードも載せようと思うと回答が分割されすぎてしまいそうで、また最近teratailもリニューアルに伴っていまいち挙動がおかしい様子で長文を投稿するのははばかられ（過去の投稿文を見てみると追記分が消えてしまっていたりして、ご質問者さんとのコミュニケーションにも支障がありそうでした）、代わりに作成以来ずっと持ち腐れにしていたQiitaのアカウントを利用させていただこうと思った次第です。

本題

使用した各機能のバージョンは...

Unity：2021.2.8f1
Burst：1.6.4
Collections：1.1.0
Mathematics：1.2.5

となっています。なお、Collectionsはパッケージ一覧に出現しなかったため、「Add package from git URL...」メニューからインストールしました。

シチュエーションは「unity 3D 切り取り」の時と同じで、モアイ像を謎の棍棒で削るというものです。
モアイ像は以前と同じ大英博物館のHoa Hakananai'aですが、棍棒は先端の形をトゲトゲに変更してみました。削り取りの精密さを以前よりも向上させるのを主眼としましたので、効果を確認するにはある程度の鋭さを持たせた凶悪な形の方が適しているんじゃないかな...と思ってのことです。

各ファイルの概要

今回は、関連する各ファイルは一つのフォルダにまとめることにしました。

Carvee.cs：削られる側のスクリプト。モアイ像にアタッチする。
Carver.cs：削る側のスクリプト。棍棒にアタッチする。以前と同様に削り取り処理の中核はCarvee側が担っており、Carverは削り取りを実演する上での脇役。
Carving.asmdef：アセンブリ定義ファイル。
Filling.compute：モデル表面のボクセル化ののち、モデル内部を埋め立てる作業を担当するシェーダー。
MortonUtility.cs：3次元の整数座標にZ階数曲線に従った番号（モートンインデックス／モートンコード）を振る、あるいは逆にモートンインデックスから3次元座標を復元する機能を担当するユーティリティクラス。
TrianglePartitioner.cs：一度にボクセル化するのにちょうどいい大きさに空間を切り分ける際に使用するクラス。
Volume.shader：ツリー化されたモデルを画面上にレンダリングするためのマテリアルに用いるシェーダー。
Voxelization.compute：TrianglePartitionerによって切り分けられた小空間をボクセル化する作業を担当するシェーダー。
VoxelUtility.cs：ボクセル化やツリー構築の中核を担当するユーティリティクラス。

処理の流れの概要

まずCarveeのStartタイミングで、削られる側のモデルメッシュからツリー構造を作成しました。手順としては下記のような流れになっています。

ボクセル密度とモデルのサイズをもとに、モデルを包む一辺が2の冪乗の空間を定める（ルート空間）。以前のTexture3Dを使ったバージョンではモデルのバウンディングボックスに沿ったサイズだったが、今回はこの2の冪乗の縛りのため、いくらか余白のある大きめの空間になりがち。
ルート空間をさいの目状に切り、一度にボクセル化するのに手頃なサイズの小空間に分ける。
モデルメッシュの各ポリゴンを小空間に振り分ける。
各小空間についてボクセル化を行う。
1. 小空間を構成する各ボクセルについて、その小空間に振り分けられた各ポリゴンとの交差判定を行う。
2. 交差がある場合、そのボクセルの座標に基づくモートンインデックス、法線、Z表裏（ポリゴンのZ方向に関する向き）をひとまとめにしたボクセルデータを出力する。
ツリー構造を作る。
1. ツリーにルートノードを追加する。これが前述のルート空間に相当する。
2. 出力されたボクセルデータそれぞれについて、モートンインデックスを頼りにルートノードからボクセルの位置まで潜行させる。
3. 目的の位置まで到達できなかった場合は要分割ということになる。行き止まったノードのインデックスを分割リストに追加する。
  1. 分割リストが空なら（一つもノード分割要求がなければ）すべてのボクセルが目的位置まで到達できる状態になったと判断し、手順5-4に進む。
  2. 空でなければ、分割リストに載っているノードそれぞれに8個の子ノードを追加し、分割リストは空にする。その後5-2の手順に戻る。
4. 各ボクセルデータと対応する位置のノードに、そのボクセルの法線とZ表裏を埋め込む。
モデルの内部を埋め立てる。
1. Z=0平面上の各点を起点に、+Z方向にレイキャストを行う。ボクセルにヒットしたら、そのボクセルのZ表裏を調べる。
2. 向きが負なら（レイに対して表面を向けているなら）、その地点から先はモデルの内部だと判断する（入射点）。そこを起点にふたたび+Z方向へのレイキャストを行い、同じくヒットしたボクセルのZ表裏を調べる。
3. 向きが正なら（レイに対して裏面を向けているなら）、その地点から先はモデルの外部だと判断する（出射点）。
4. 入射点から出射点までをZ方向に走査し、間のノードに充填フラグを立てる。
5. 出射点からふたたび+Z方向へのレイキャストを行う。ルート空間から外にレイが出るまで、手順6-2に戻って埋め立て作業を繰り返す。

そしてCarveeのCarveを実行すると下記の手順で削り取りを行い、モデルのツリーを削られた形に更新します。

ルート空間とCarverが交差する領域を調べ、その領域を包含するノード（共通ノード）を求める。
共通ノードの空間をさいの目状に切り、一度にボクセル化するのに手頃なサイズの小空間に分ける。
モデルメッシュの各ポリゴンを小空間に振り分ける。
各小空間についてボクセル化を行う。
1. 小空間を構成する各ボクセルについて、その小空間に振り分けられた各ポリゴンとの交差判定を行う。
2. 交差がなく、かつその位置に既存のCarvee表面ボクセルが存在し、さらにその地点がCarverモデルの外部であるなら、そのCarvee表面ボクセルデータを出力する。
3. 交差があり、かつその位置が既存のCarveeモデルの内部であるなら、Carver表面ボクセルデータを生成・出力する（なお、削り取られたくぼみなので面の向きは反転させておく）。さらに、その位置がCarveeモデルの表面であるなら、データを既存のCarvee表面ボクセルデータと混合したうえで出力する。
既存のツリー構造から共通ノード以下を削除した上で、前述と同様の手順でツリー構造を作る。このとき、共通ノードの削除で空いた領域は再利用する。
前述と同様の手順でモデルの内部を埋め立てる。

ボクセル化時に出力されるボクセルデータは下図のような形にしてみました。

法線のx、y成分は14ビットの整数で表現し、-8191～8191が-1.0～1.0に対応します。z成分は符号だけを保存し、0なら正、1なら負を表します。
最上位の2ビットがポリゴンの表裏を表しており、左ビットが1なら-Z、右ビットが1なら+Zを向いています。Z軸に平行で表裏を決定できない場合はどちらも0、モデルの表面ではない領域はどちらも1とすることにしました。

シーンのセットアップ

モアイ像は以前と同様に身長2m程度に調整、今回作成した新しいCarveeスクリプトをアタッチしています。なお、もう一つアタッチされているMoaiControllerはCarving外のスクリプトで、単にモアイ像を回転させるだけのものです（以前のTexture3D版でも用いたものと同一です）。

インスペクター上の「Volume Shader」にはVolume.shader、「Voxelization Shader」にはVoxelization.compute、「Filling Shader」にはFilling.computeがそれぞれセットされています。「Voxels Per Unit」と「Color」は以前のCarveeと同じくボクセルの密度および表面の色を設定するためのものです。

「Cubic Appearance」以降はボクセルの配置やZ表裏、充填状況を視覚的に確認する目的で作成した機能をそのまま残したものです。これをオンにすると、Z表裏が-Zなら青、+Zなら赤、中立なら緑、非表面領域ならマゼンタ色の立方体として描画します。
オンの状態で拡大してみますと、モデルがキューブの集まりになっている様子をご覧いただけるかと思います。

陰影再現の程度

「unity 3D 切り取り」の時のTexture3D上の値は充填状態を表すだけで、法線はレンダリング時に周辺をサンプリングして充填状態の変化率から算出していました。一方、今回は各ボクセルに法線情報を埋め込みましたのでモデル本来の陰影付けに近い陰影になったかと思います。
本来のメッシュによる陰影が下図のようになるのに対して...

Voxels Per Unitが128だと下図のように、

Voxels Per Unitが2048だと下図のようになりました。

Texture3D版の時と同じ密度128でも、密度の低さゆえのガタガタ感は相変わらずですが、陰影の再現度はけっこう向上できたのではないでしょうか。

削り取りの様子

冒頭でも申し上げましたが、今のままではリアルタイム彫刻には使えないかと思います。私としてはSVO表現の入門のつもりで実験してみた...という認識でして、実用的な速度が出せるかどうかは二の次、Web上で検索した論文やらをつまみ食いして書き散らかした段階です...
処理実行時も完了を同期的に待機しており、つまりStartやCarveを実行すると終わるまでUnityがフリーズします。

Voxels Per Unitが128（ルート空間のサイズは512³）だと下図のような見た目となりました。

スピード面ではTexture3D版に劣り、削るたびに一拍引っかかりが生じています。リアルタイムの削り取りには足りず、Texture3D版はCarver中でキー入力をInput.GetKey(KeyCode.Space)で判定して連続的にCarveを実行していたところを、今回はInput.GetKeyDown(KeyCode.Space)に変更して散発的に削っています。
解像度をあきらめてもかまわなければ、削り取り方式はTexture3D版と同様にした上でTexture3D上に法線データを埋め込んで見た目の改善を狙った方がマシかもしれません。

また、Voxels Per Unitが2048（ルート空間のサイズは8192³）の場合はこうなりました。

ひと掘り40分（ボクセル化39分、ツリー更新1分、内部埋め立て30秒）という馬鹿げた状態ではありますが、八分木構造のおかげでサイズ8192³の空間...つまり以前のTexture3D版なら5500億個のボクセルを要するケースであってもノード総数は8000万個程度に抑えることができ、私の廉価ノートPCでもなんとか実行できました。

参考（と称してつまみ食いした）情報

Fast 3D Triangle-Box Overlap Testing
- ボクセルとポリゴンの交差判定にはこの方式を使わせていただきました。
Octree-Based Sparse Voxelization Using the GPU Hardware Rasterizer
- 八分木構造を作る手順の参考にさせていただきました。前半にはGPUのラスタライズ機能を利用したボクセル化の手法も解説されていたのですが、今回の例ではその方法ではなく、素朴に各ボクセルとポリゴンの交差判定をとるやり方になっています。
Efficient Ray Tracing of Sparse Voxel Octrees on an FPGA
- レンダリング用シェーダーでのツリー内ノードの探索には、この手順を利用させていただきました。また、埋め立てシェーダー内でもこの手順をアレンジした探索方法を使っています。
Fast Parallel Surface and Solid Voxelization on GPUs
- 中身の詰まったモデルを作ることについて言及があり、これは使える！と思ったのですが、読解力不足により自身の脳内でイメージを確立できず挫折しました。代わりの苦肉の策として、前述のようにZ表裏の情報を手がかりにちまちまと埋め立てを行うことにしたのです...

各ファイルの詳細

Carvee.cs

using System.Diagnostics;
using Unity.Burst;
using Unity.Collections;
using Unity.Collections.LowLevel.Unsafe;
using Unity.Jobs;
using Unity.Mathematics;
using UnityEngine;
using UnityEngine.Rendering;
using Debug = UnityEngine.Debug;

namespace Carving
{
    [RequireComponent(typeof(MeshFilter), typeof(MeshRenderer))]
    public class Carvee : MonoBehaviour
    {
        private const int RootResolutionMin = 8;
        private const int PartitionResolutionMax = 128;
        private static readonly int ClipMinProperty = Shader.PropertyToID("_ClipMin");
        private static readonly int ClipMaxProperty = Shader.PropertyToID("_ClipMax");
        private static readonly int NodePoolProperty = Shader.PropertyToID("_NodePool");
        private static readonly int VolumeMinProperty = Shader.PropertyToID("_VolumeMin");
        private static readonly int VolumeSizeProperty = Shader.PropertyToID("_VolumeSize");
        private static readonly int NormalizedLeafSizeProperty = Shader.PropertyToID("_NormalizedLeafSize");

        [SerializeField] private Shader volumeShader;
        [SerializeField] private ComputeShader voxelizationShader;
        [SerializeField] private ComputeShader fillingShader;
        [SerializeField][Min(0)] private int voxelsPerUnit = 128;
        [SerializeField] private Color color = Color.white;
        [SerializeField] private bool cubicAppearance;
        [SerializeField] private bool hideInternalNodes;
        [SerializeField][Range(0.0f, 1.0f)] private float clipXMin;
        [SerializeField][Range(0.0f, 1.0f)] private float clipXMax = 1.0f;
        [SerializeField][Range(0.0f, 1.0f)] private float clipYMin;
        [SerializeField][Range(0.0f, 1.0f)] private float clipYMax = 1.0f;
        [SerializeField][Range(0.0f, 1.0f)] private float clipZMin;
        [SerializeField][Range(0.0f, 1.0f)] private float clipZMax = 1.0f;

        private Bounds carverBox;
        private TrianglePartitioner carverPartitioner;
        private LocalKeyword cubicAppearanceKeyword;
        private NativeQueue<int> deadNodes;
        private LocalKeyword hideInternalNodesKeyword;
        private Material initialMaterial;
        private Vector3 initialScale;
        private MeshFilter meshFilter;
        private MeshRenderer meshRenderer;
        private NativeList<VoxelUtility.TreeNode> nodePool;
        private ComputeBuffer nodePoolBuffer;
        private Bounds previousClip = new Bounds(new Vector3(float.NaN, float.NaN, float.NaN), new Vector3(float.NaN, float.NaN, float.NaN));
        private Color previousColor;
        private bool previousCubicAppearance;
        private bool previousHideInternalNodes;
        private int rootResolution;
        private Mesh sourceMesh;
        private Bounds volumeBounds;
        private Material volumeMaterial;
        private Mesh volumeMesh;
        private bool isInitialized;

        public Vector3 Center => this.transform.TransformPoint(this.volumeBounds.center);

        private void Start()
        {
            // ボクセル化作業用コンピュートシェーダーを準備する
            if (this.voxelizationShader != null)
            {
                VoxelUtility.EnableComputeShader(this.voxelizationShader, this.fillingShader);
            }

            // MeshFilterから元となるメッシュを取得する
            // なお、メッシュは「Read/Write」をオンにしておく必要がある
            this.meshFilter = this.GetComponent<MeshFilter>();
            this.sourceMesh = this.meshFilter!.sharedMesh;
            if ((this.sourceMesh == null) || !this.sourceMesh.isReadable)
            {
                Debug.LogError("Source mesh is null or not readable!");
                this.enabled = false;
                return;
            }

            // ボリュームの形は立方体にする
            // 辺の長さは2の冪乗にしておく
            this.initialScale = this.transform.localScale;
            var sourceBounds = this.sourceMesh.bounds;
            var scaledSourceBounds = new Bounds(Vector3.Scale(sourceBounds.center, this.initialScale), Vector3.Scale(sourceBounds.size, this.initialScale));
            var rootResolution3 = scaledSourceBounds.size * this.voxelsPerUnit;
            var rootResolution3Int = new Vector3Int(
                Mathf.Max(Mathf.CeilToInt(rootResolution3.x), RootResolutionMin),
                Mathf.Max(Mathf.CeilToInt(rootResolution3.y), RootResolutionMin),
                Mathf.Max(Mathf.CeilToInt(rootResolution3.z), RootResolutionMin));
            this.rootResolution = Mathf.NextPowerOfTwo(
                Mathf.Max(Mathf.Max(rootResolution3Int.x, rootResolution3Int.y), rootResolution3Int.z));
            var rootSize = (float)this.rootResolution / this.voxelsPerUnit;
            Debug.Log($"Root resolution: {this.rootResolution}^3");
            if (this.rootResolution > MortonUtility.MaxResolution)
            {
                Debug.LogError($"Root resolution must not exceed {MortonUtility.MaxResolution}^3!");
                this.enabled = false;
                return;
            }

            // 本来のメッシュを包むようにボリューム描画用メッシュを作る
            this.volumeBounds = new Bounds(scaledSourceBounds.center, new Vector3(rootSize, rootSize, rootSize));
            var vMin = this.volumeBounds.min;
            var vMax = this.volumeBounds.max;
            this.volumeMesh = new Mesh
            {
                name = this.sourceMesh.name + " Volume",
                vertices = new[]
                {
                    vMin,
                    new Vector3(vMax.x, vMin.y, vMin.z),
                    new Vector3(vMin.x, vMax.y, vMin.z),
                    new Vector3(vMax.x, vMax.y, vMin.z),
                    new Vector3(vMin.x, vMin.y, vMax.z),
                    new Vector3(vMax.x, vMin.y, vMax.z),
                    new Vector3(vMin.x, vMax.y, vMax.z),
                    vMax
                },
                triangles = new[]
                {
                    0, 1, 2, 3, 2, 1,
                    0, 2, 4, 6, 4, 2,
                    0, 4, 1, 5, 1, 4,
                    7, 5, 6, 4, 6, 5,
                    7, 3, 5, 1, 5, 3,
                    7, 6, 3, 2, 3, 6
                }
            };
            this.volumeMesh.RecalculateBounds();
            this.meshFilter.sharedMesh = this.volumeMesh;
            this.transform.localScale = Vector3.one;

            // ボリューム全体を小区画に分ける
            var partitionResolution = Mathf.Min(this.rootResolution, PartitionResolutionMax);
            var partitionCountPerAxis = Mathf.Max(this.rootResolution / PartitionResolutionMax, 1);
            var partitionSize = rootSize / partitionCountPerAxis;

            // 元のメッシュの頂点の位置・法線と、面のインデックスを取得し
            // 各ポリゴンを小区画に振り分ける
            var stopwatch = new Stopwatch();
            using var sourceMeshDataArray = Mesh.AcquireReadOnlyMeshData(this.sourceMesh);
            var sourceMeshData = sourceMeshDataArray[0];
            var sourceVerticesVector3 = new NativeArray<Vector3>(sourceMeshData.vertexCount, Allocator.TempJob);
            var sourceNormalsVector3 = new NativeArray<Vector3>(sourceMeshData.vertexCount, Allocator.TempJob);
            var sourceIndices = new NativeArray<int>(sourceMeshData.GetSubMesh(0).indexCount, Allocator.TempJob);
            sourceMeshData.GetVertices(sourceVerticesVector3);
            sourceMeshData.GetNormals(sourceNormalsVector3);
            sourceMeshData.GetIndices(sourceIndices, 0);
            var sourceVertices = sourceVerticesVector3.Reinterpret<float3>();
            var sourceNormals = sourceNormalsVector3.Reinterpret<float3>();
            stopwatch.Restart();
            new SourceScalingJob
            {
                SourceVertices = sourceVertices,
                SourceNormals = sourceNormals,
                SourceScale = this.initialScale,
                NormalScale = math.float3(1.0f) / this.initialScale
            }.Schedule(sourceVertices.Length, 64).Complete();
            stopwatch.Stop();
            Debug.Log($"Scaling finished in {stopwatch.ElapsedMilliseconds} ms.");
            var sourceTriangles = sourceIndices.Reinterpret<int3>(sizeof(int));
            using var partitioner = new TrianglePartitioner(
                sourceVertices,
                sourceTriangles,
                this.volumeBounds,
                partitionCountPerAxis,
                int3.zero,
                Allocator.TempJob);
            Debug.Log($"Split volume into {partitionCountPerAxis}^3 partitions. Each partition size: {partitionSize}, Resolution: {partitionResolution}^3");
            stopwatch.Restart();
            partitioner.ExecutePartitioning().Complete();
            stopwatch.Stop();
            Debug.Log($"Partitioning finished in {stopwatch.ElapsedMilliseconds} ms.");
            var partitions = partitioner.Result;

            // 個々の小区画を最高解像度でボクセル化する
            stopwatch.Restart();
            var needsUpdateSourceBuffers = true;
            var voxelDataQueue = new NativeQueue<VoxelUtility.VoxelData>(Allocator.TempJob);
            for (var i = 0; i < partitions.Length; i++)
            {
                var partition = partitions[i];
                if (partition.Triangles.IsEmpty)
                {
                    continue;
                }

                MortonUtility.Decode32((uint)i, out var partitionCoord);
                VoxelUtility.Voxelize(
                    sourceVertices,
                    sourceNormals,
                    sourceTriangles,
                    partition.Triangles,
                    partitionCoord,
                    partition.Min,
                    partitionSize,
                    partitionResolution,
                    this.rootResolution,
                    voxelDataQueue,
                    needsUpdateSourceBuffers);
                needsUpdateSourceBuffers = false;
            }
            var voxelDataPool = voxelDataQueue.ToArray(Allocator.TempJob);
            voxelDataQueue.Dispose();
            stopwatch.Stop();
            Debug.Log($"Voxelization finished in {stopwatch.ElapsedMilliseconds} ms. Data count: {voxelDataPool.Length}");

            // 元のメッシュのデータは不要になったので廃棄する
            sourceVerticesVector3.Dispose();
            sourceNormalsVector3.Dispose();
            sourceIndices.Dispose();

            // 八分木を構築する
            stopwatch.Restart();
            this.deadNodes = new NativeQueue<int>(Allocator.Persistent);
            this.nodePool = new NativeList<VoxelUtility.TreeNode>(Allocator.Persistent);
            this.nodePool.Add(VoxelUtility.TreeNode.Create());
            var rootLevel = math.floorlog2(this.rootResolution) - 1;
            while (VoxelUtility.SubdivideNodesOneLevel(rootLevel, voxelDataPool, this.nodePool))
            {
                Debug.Log($"Pool size: {this.nodePool.Length}/{this.nodePool.Capacity}");
            }
            VoxelUtility.EmbedVoxelData(rootLevel, voxelDataPool, this.nodePool);
            voxelDataPool.Dispose();
            stopwatch.Stop();
            Debug.Log($"Octree building finished in {stopwatch.ElapsedMilliseconds} ms. Node count: {this.nodePool.Length}");

            // まだモデルの内部はがらんどうなので、Z方向に走査してモデル内ノードに充填フラグを立てる
            stopwatch.Restart();
            VoxelUtility.Fill(this.nodePool, this.rootResolution);
            stopwatch.Stop();
            Debug.Log($"Model filling finished in {stopwatch.ElapsedMilliseconds} ms.");

            // ボリューム描画用マテリアルを作り、レンダラーにセットする
            this.UpdateNodePoolBuffer();
            this.meshRenderer = this.GetComponent<MeshRenderer>();
            this.initialMaterial = this.meshRenderer!.sharedMaterial;
            this.volumeMaterial = new Material(this.volumeShader);
            this.meshRenderer.sharedMaterial = this.volumeMaterial;
            this.volumeMaterial.color = this.color;
            this.previousColor = this.color;
            this.cubicAppearanceKeyword = new LocalKeyword(this.volumeShader, "_CUBIC_APPEARANCE");
            this.hideInternalNodesKeyword = new LocalKeyword(this.volumeShader, "_HIDE_INTERNAL_NODES");
            this.volumeMaterial.SetKeyword(this.cubicAppearanceKeyword, this.cubicAppearance);
            this.previousCubicAppearance = this.cubicAppearance;
            this.volumeMaterial.SetKeyword(this.hideInternalNodesKeyword, this.hideInternalNodes);
            this.previousHideInternalNodes = this.hideInternalNodes;
            this.volumeMaterial.SetVector(VolumeMinProperty, this.volumeBounds.min);
            this.volumeMaterial.SetFloat(VolumeSizeProperty, rootSize);
            this.volumeMaterial.SetFloat(NormalizedLeafSizeProperty, 2.0f / this.rootResolution);
            this.volumeMaterial.SetBuffer(NodePoolProperty, this.nodePoolBuffer);
            this.isInitialized = true;
        }

        private void Update()
        {
            if (this.color != this.previousColor)
            {
                this.volumeMaterial!.color = this.color;
                this.previousColor = this.color;
            }
            if (this.cubicAppearance != this.previousCubicAppearance)
            {
                this.volumeMaterial!.SetKeyword(this.cubicAppearanceKeyword, this.cubicAppearance);
                this.previousCubicAppearance = this.cubicAppearance;
            }
            if (this.hideInternalNodes != this.previousHideInternalNodes)
            {
                this.volumeMaterial!.SetKeyword(this.hideInternalNodesKeyword, this.hideInternalNodes);
                this.previousHideInternalNodes = this.hideInternalNodes;
            }
            var clipMin = new Vector3(this.clipXMin, this.clipYMin, this.clipZMin);
            var clipMax = new Vector3(this.clipXMax, this.clipYMax, this.clipZMax);
            clipMin = Vector3.Min(clipMin, clipMax);
            this.clipXMin = clipMin.x;
            this.clipYMin = clipMin.y;
            this.clipZMin = clipMin.z;
            var clip = new Bounds
            {
                min = clipMin,
                max = clipMax
            };
            if (clip != this.previousClip)
            {
                this.volumeMaterial!.SetVector(ClipMinProperty, clip.min);
                this.volumeMaterial.SetVector(ClipMaxProperty, clip.max);
                this.previousClip = clip;
            }
        }

        private void OnDestroy()
        {
            if (this.isInitialized)
            {
                VoxelUtility.DisableComputeShader();
                if (this.nodePool.IsCreated)
                {
                    this.nodePool.Dispose();
                }

                if (this.deadNodes.IsCreated)
                {
                    this.deadNodes.Dispose();
                }

                this.nodePoolBuffer?.Dispose();
                this.meshRenderer!.sharedMaterial = this.initialMaterial;
                this.meshFilter!.sharedMesh = this.sourceMesh;
                Destroy(this.volumeMaterial);
                Destroy(this.volumeMesh);
                this.carverPartitioner?.Dispose();
            }
        }

        public void Carve(
            NativeArray<float3> vertices,
            NativeArray<float3> normals,
            NativeArray<int3> triangles,
            Transform carverTransform)
        {
            // Carverの頂点をボリューム空間に移し、バウンディングボックスを求める
            const int threadCount = 16;
            if (!vertices.IsCreated || !normals.IsCreated || !triangles.IsCreated || (carverTransform == null))
            {
                return;
            }
            var vertexCount = vertices.Length;
            using var carverVertices = new NativeArray<float3>(vertexCount, Allocator.TempJob, NativeArrayOptions.UninitializedMemory);
            using var carverNormals = new NativeArray<float3>(vertexCount, Allocator.TempJob, NativeArrayOptions.UninitializedMemory);
            using var carverMinMax = new NativeArray<float3x2>(threadCount, Allocator.TempJob, NativeArrayOptions.UninitializedMemory);
            var carverToCarvee = (float4x4)(this.transform.worldToLocalMatrix * carverTransform.localToWorldMatrix);
            var carverToCarveeNormal = math.transpose(math.inverse(math.float3x3(carverToCarvee)));
            new CarverTransformJob
            {
                SourceVertices = vertices,
                SourceNormals = normals,
                Vertices = carverVertices,
                Normals = carverNormals,
                MinMax = carverMinMax,
                VertexCountPerThread = (int)math.ceil((double)vertexCount / threadCount),
                CarverToCarvee = carverToCarvee,
                CarverToCarveeNormal = carverToCarveeNormal
            }.Schedule(threadCount, 1).Complete();
            var carverMin = math.float3(float.PositiveInfinity);
            var carverMax = math.float3(float.NegativeInfinity);
            foreach (var minMax in carverMinMax)
            {
                carverMin = math.min(carverMin, minMax.c0);
                carverMax = math.max(carverMax, minMax.c1);
            }

            // CarverとCarveeのバウンディングボックスの交差がなければ、処理不要なので中断
            var volumeMin = (float3)this.volumeBounds.min;
            var volumeSize = this.volumeBounds.size.x;
            var clampedMin = math.clamp(carverMin, volumeMin, volumeMin + volumeSize);
            var clampedMax = math.clamp(carverMax, volumeMin, volumeMin + volumeSize);
            var clampedSize = clampedMax - clampedMin;
            if ((clampedSize.x * clampedSize.y * clampedSize.z) <= 0.0f)
            {
                return;
            }

            // Carverのバウンディングボックスの最大・最小がどのボクセルに相当するか求め、モートンコードの排他論理和を得る
            // それの上の桁のゼロの数を数え、ルートから数えて何レベルまでが共通空間かを調べる
            var coordScale = this.rootResolution / volumeSize;
            var rootCoordMax = math.int3(this.rootResolution - 1);
            var carverCoordMin = math.clamp(math.int3((carverMin - volumeMin) * coordScale), int3.zero, rootCoordMax);
            var carverCoordMax = math.clamp(math.int3((carverMax - volumeMin) * coordScale), int3.zero, rootCoordMax);
            var rootCodeMax = ((ulong)this.rootResolution * (ulong)this.rootResolution * (ulong)this.rootResolution) - 1;
            var carverCodeMin = MortonUtility.Encode(carverCoordMin);
            var carverCodeMax = MortonUtility.Encode(carverCoordMax);
            var bitOffset = math.lzcnt(rootCodeMax);
            var sharedLevel = (math.lzcnt(carverCodeMin ^ carverCodeMax) - bitOffset) / 3;

            // 八分木を潜行し、共通空間ノードを得る
            // 途中で空ノードで行き止まったらCarveeとCarverの交差はない
            var nodeIndex = 0;
            var node = this.nodePool[nodeIndex];
            var nodeMin = volumeMin;
            var nodeSize = volumeSize;
            var nodeResolution = this.rootResolution;
            var nodeCoord = int3.zero;
            for (var i = 0; i < sharedLevel; i++)
            {
                if (node.IsEmpty)
                {
                    return;
                }
                if (node.IsLeaf)
                {
                    break;
                }
                var childIndex = (int)((carverCodeMin >> (64 - bitOffset - (3 * (i + 1)))) & 7);
                nodeIndex = node.FirstChild + childIndex;
                node = this.nodePool[nodeIndex];
                nodeSize *= 0.5f;
                nodeResolution >>= 1;
                var childCoord = math.int3(childIndex & 1, (childIndex & 2) >> 1, (childIndex & 4) >> 2);
                nodeMin += nodeSize * (float3)childCoord;
                nodeCoord = (nodeCoord * 2) + childCoord;
            }
            this.carverBox = new Bounds
            {
                min = nodeMin,
                max = nodeMin + nodeSize
            };

            // ノードを小区画に分ける
            var partitionResolution = Mathf.Min(nodeResolution, PartitionResolutionMax);
            var partitionCountPerAxis = Mathf.Max(nodeResolution / PartitionResolutionMax, 1);
            var partitionSize = nodeSize / partitionCountPerAxis;

            // 各ポリゴンを小区画に振り分ける
            var stopwatch = new Stopwatch();
            var partitionCoordOffset = nodeCoord * partitionCountPerAxis;
            this.carverPartitioner ??= new TrianglePartitioner(
                carverVertices,
                triangles,
                this.carverBox,
                partitionCountPerAxis,
                partitionCoordOffset
            );
            this.carverPartitioner.Resize(
                carverVertices,
                triangles,
                this.carverBox,
                partitionCountPerAxis,
                partitionCoordOffset);
            Debug.Log($"Split volume into {partitionCountPerAxis}^3 partitions. Each partition size: {partitionSize}, Resolution: {partitionResolution}^3");
            stopwatch.Restart();
            this.carverPartitioner.ExecutePartitioning().Complete();
            stopwatch.Stop();
            Debug.Log($"Partitioning finished in {stopwatch.ElapsedMilliseconds} ms.");
            var partitions = this.carverPartitioner.Result;

            // 個々の小区画を最高解像度でボクセル化する
            // ここでは面反転オプションをオンにし（faceSignを-1.0f）、くぼみとして扱う
            stopwatch.Restart();
            var needsUpdateSourceBuffers = true;
            var voxelDataQueue = new NativeQueue<VoxelUtility.VoxelData>(Allocator.TempJob);
            VoxelUtility.UploadNodePool(this.nodePool);
            for (var i = 0; i < partitions.Length; i++)
            {
                var partition = partitions[i];
                VoxelUtility.Voxelize(
                    carverVertices,
                    carverNormals,
                    triangles,
                    partition.Triangles,
                    partition.Coord,
                    partition.Min,
                    partitionSize,
                    partitionResolution,
                    this.rootResolution,
                    voxelDataQueue,
                    needsUpdateSourceBuffers,
                    nodeIndex,
                    nodeResolution,
                    -1.0f);
                needsUpdateSourceBuffers = false;
            }
            var voxelDataPool = voxelDataQueue.ToArray(Allocator.TempJob);
            voxelDataQueue.Dispose();
            stopwatch.Stop();
            Debug.Log($"Voxelization finished in {stopwatch.ElapsedMilliseconds} ms. Data count: {voxelDataPool.Length}");

            // 八分木から共通空間ノードを削除したのち、八分木を再構築する
            stopwatch.Restart();
            VoxelUtility.KillNode(this.nodePool, nodeIndex, this.deadNodes);
            var rootLevel = math.floorlog2(this.rootResolution) - 1;
            while (VoxelUtility.SubdivideNodesOneLevel(rootLevel, voxelDataPool, this.nodePool, this.deadNodes))
            {
                Debug.Log($"Pool size: {this.nodePool.Length}/{this.nodePool.Capacity}");
            }
            VoxelUtility.EmbedVoxelData(rootLevel, voxelDataPool, this.nodePool);
            voxelDataPool.Dispose();
            stopwatch.Stop();
            Debug.Log($"Octree building finished in {stopwatch.ElapsedMilliseconds} ms. Node count: {this.nodePool.Length}");

            // モデル内ノードに充填フラグを立てる
            stopwatch.Restart();
            VoxelUtility.Fill(this.nodePool, this.rootResolution, nodeResolution, nodeCoord.xy);
            stopwatch.Stop();
            Debug.Log($"Model filling finished in {stopwatch.ElapsedMilliseconds} ms.");

            // 描画用のノードプールへデータを転送する
            this.UpdateNodePoolBuffer();
            this.volumeMaterial!.SetBuffer(NodePoolProperty, this.nodePoolBuffer);
        }

        private void UpdateNodePoolBuffer()
        {
            this.nodePoolBuffer ??= new ComputeBuffer(256, UnsafeUtility.SizeOf<VoxelUtility.TreeNode>());
            VoxelUtility.ExtendBufferIfNeeded(ref this.nodePoolBuffer, this.nodePool.Length);
            this.nodePoolBuffer!.SetData((NativeArray<VoxelUtility.TreeNode>)this.nodePool, 0, 0, this.nodePool.Length);
        }

        [BurstCompile]
        private struct SourceScalingJob : IJobParallelFor
        {
            public NativeArray<float3> SourceVertices;
            public NativeArray<float3> SourceNormals;
            public float3 SourceScale;
            public float3 NormalScale;

            public void Execute(int i)
            {
                this.SourceVertices[i] *= this.SourceScale;
                this.SourceNormals[i] = math.normalizesafe(this.SourceNormals[i] * this.NormalScale);
            }
        }

        [BurstCompile]
        private struct CarverTransformJob : IJobParallelFor
        {
            [ReadOnly] public NativeArray<float3> SourceVertices;
            [ReadOnly] public NativeArray<float3> SourceNormals;
            [NativeDisableParallelForRestriction][WriteOnly] public NativeArray<float3> Vertices;
            [NativeDisableParallelForRestriction][WriteOnly] public NativeArray<float3> Normals;
            [WriteOnly] public NativeArray<float3x2> MinMax;
            public int VertexCountPerThread;
            public float4x4 CarverToCarvee;
            public float3x3 CarverToCarveeNormal;

            public void Execute(int threadIndex)
            {
                var min = math.float3(float.PositiveInfinity);
                var max = math.float3(float.NegativeInfinity);
                var fromIndex = math.min(threadIndex * this.VertexCountPerThread, this.Vertices.Length);
                var toIndex = math.min(fromIndex + this.VertexCountPerThread, this.Vertices.Length);
                for (var i = fromIndex; i < toIndex; i++)
                {
                    var vertex = math.mul(this.CarverToCarvee, math.float4(this.SourceVertices[i], 1.0f)).xyz;
                    this.Vertices[i] = vertex;
                    this.Normals[i] = math.normalizesafe(math.mul(this.CarverToCarveeNormal, this.SourceNormals[i]));
                    min = math.min(min, vertex);
                    max = math.max(max, vertex);
                }
                this.MinMax[threadIndex] = math.float3x2(min, max);
            }
        }
    }
}

Carver.cs

using System.Collections;
using Unity.Collections;
using Unity.Mathematics;
using UnityEngine;

namespace Carving
{
    [RequireComponent(typeof(MeshFilter), typeof(MeshRenderer))]
    public class Carver : MonoBehaviour
    {
        [SerializeField] private Carvee target;
        [SerializeField][Range(0.0f, 50.0f)] private float rotationSpeed = 10.0f;
        [SerializeField][Range(0.0f, 5.0f)] private float pushPullSpeed = 1.0f;
        [SerializeField] private Color inactiveColor = Color.green;
        [SerializeField] private Color activeColor = Color.red;

        private Material carverMaterial;
        private NativeArray<float3> carverNormals;
        private NativeArray<int3> carverTriangles;
        private NativeArray<float3> carverVertices;
        private float distance;
        private float horizontalAngle;
        private float verticalAngle;

        private IEnumerator Start()
        {
            if (this.target == null)
            {
                Debug.LogError("Target is not set!");
                Destroy(this);
                yield break;
            }
            yield return null;
            Cursor.lockState = CursorLockMode.Locked;
            var meshFilter = this.GetComponent<MeshFilter>();
            var carverMesh = meshFilter!.sharedMesh;
            if ((carverMesh == null) || !carverMesh.isReadable)
            {
                // Carverのメッシュ情報も必要なので、Carvee同様こちらも「Read/Write」をオンにしておく必要がある
                Debug.LogError("Carver mesh is null or not readable!");
                Destroy(this);
                yield break;
            }
            using (var meshData = Mesh.AcquireReadOnlyMeshData(carverMesh))
            {
                this.carverVertices = new NativeArray<float3>(carverMesh.vertexCount, Allocator.Persistent);
                this.carverNormals = new NativeArray<float3>(carverMesh.vertexCount, Allocator.Persistent);
                this.carverTriangles = new NativeArray<int3>(
                    (int)(carverMesh.GetIndexCount(0) / 3),
                    Allocator.Persistent);
                meshData[0].GetVertices(this.carverVertices.Reinterpret<Vector3>());
                meshData[0].GetNormals(this.carverNormals.Reinterpret<Vector3>());
                meshData[0].GetIndices(this.carverTriangles.Reinterpret<int>(sizeof(int) * 3), 0);
            }
            var meshRenderer = this.GetComponent<MeshRenderer>();
            this.carverMaterial = meshRenderer!.material;
            var relativePosition = this.transform.position - this.target!.Center;
            this.distance = relativePosition.magnitude;
            var direction = relativePosition.normalized;
            this.verticalAngle = Mathf.Asin(direction.y) * Mathf.Rad2Deg;
            this.horizontalAngle = Mathf.Atan2(direction.z, direction.x) * Mathf.Rad2Deg;
        }

        private void Update()
        {
            if (!this.carverVertices.IsCreated)
            {
                return;
            }

            // マウス左右で水平に、上下で垂直に回転し
            // 左ボタン・右ボタンで前後に動かす
            var mouseX = Input.GetAxis("Mouse X");
            var mouseY = Input.GetAxis("Mouse Y");
            var push = Input.GetMouseButton(0) ? this.pushPullSpeed * Time.deltaTime : 0.0f;
            var pull = Input.GetMouseButton(1) ? this.pushPullSpeed * Time.deltaTime : 0.0f;
            this.distance = Mathf.Max(0.0f, (this.distance - push) + pull);
            this.horizontalAngle = (this.horizontalAngle + (mouseX * this.rotationSpeed)) % 360.0f;
            this.verticalAngle = Mathf.Clamp(this.verticalAngle + (mouseY * this.rotationSpeed), -90.0f, 90.0f);
            var h = this.horizontalAngle * Mathf.Deg2Rad;
            var v = this.verticalAngle * Mathf.Deg2Rad;
            var sinH = Mathf.Sin(h);
            var cosH = Mathf.Cos(h);
            var sinV = Mathf.Sin(v);
            var cosV = Mathf.Cos(v);
            var direction = new Vector3(cosH * cosV, sinV, sinH * cosV);
            var center = this.target!.Center;
            this.transform.position = center + (direction * this.distance);
            this.transform.LookAt(center);

            // スペースキーを押し下げたらCarveを行う
            if (Input.GetKeyDown(KeyCode.Space))
            {
                this.carverMaterial!.color = this.activeColor;
                this.target.Carve(this.carverVertices, this.carverNormals, this.carverTriangles, this.transform);
            }
            else
            {
                this.carverMaterial!.color = this.inactiveColor;
            }
        }

        private void OnDestroy()
        {
            if (this.carverVertices.IsCreated)
            {
                this.carverVertices.Dispose();
            }
            if (this.carverNormals.IsCreated)
            {
                this.carverNormals.Dispose();
            }
            if (this.carverTriangles.IsCreated)
            {
                this.carverTriangles.Dispose();
            }
            Cursor.lockState = CursorLockMode.None;
        }
    }
}

Carving.asmdef

下図のように「Allow 'unsafe' Code」をオンにしており、「Assembly Definition References」にはUnity.Burst、Unity.Collections、Unity.Mathematicsを追加しました。

Filling.compute

#pragma kernel Fill

#define TREE_TRAVERSAL_STACK_CAPACITY 8
#define FILLING_LOOP_COUNT 128

struct TreeNode
{
    uint flags;
    uint data;
};

inline bool IsEmpty(TreeNode node)
{
    return (node.flags & 1) > 0;
}

inline bool IsLeaf(TreeNode node)
{
    return (node.flags & 2) > 0;
}

inline void FillNode(inout TreeNode node)
{
    node.flags |= 4;
}

bool ZcastCore(
    uint rootResolution,
    int3 origin,
    RWStructuredBuffer<TreeNode> nodePool,
    out int enter,
    out int nodeIndex)
{
    enter = 0;
    nodeIndex = 0;
    bool hit = false;
    int t0Stack[TREE_TRAVERSAL_STACK_CAPACITY];
    int t1Stack[TREE_TRAVERSAL_STACK_CAPACITY];
    int nodeIndexStack[TREE_TRAVERSAL_STACK_CAPACITY];
    int childIndexStack[TREE_TRAVERSAL_STACK_CAPACITY];
    int2 centerStack[TREE_TRAVERSAL_STACK_CAPACITY];
    int childIndex = 0;
    int tStart = origin.z;
    while (true)
    {
        int t0 = 0;
        int t1 = (int)rootResolution;
        int2 center = (int2)(rootResolution >> 1);
        nodeIndex = 0;
        bool isFirstNode = true;
        int level = 0;
        uint stackIndex = 0;
        uint stackLength = 0;
        TreeNode node = nodePool[nodeIndex];
        bool exitRoot = false;
        while (true)
        {
            int tm = ((t1 - t0) >> 1) + t0;
            bool exitNode = false;
            bool stackUnderflow = stackLength == 0;
            childIndex = isFirstNode
                ? ((origin.x >= center.x ? 1 : 0) | (origin.y >= center.y ? 2 : 0))
                : childIndex + 4;
            isFirstNode = false;
            exitNode = (childIndex >> 3) > 0;
            bool positionZ = (childIndex & 4) > 0;
            int childT0 = positionZ ? tm : t0;
            int childT1 = positionZ ? t1 : tm;
            bool passNode = childT1 <= tStart;
            if (exitNode)
            {
                exitRoot = (level == 0);
                if (stackUnderflow || exitRoot)
                {
                    tStart = max(t1, tStart);
                    break;
                }
                level--;
                stackIndex = ((stackIndex + TREE_TRAVERSAL_STACK_CAPACITY) - 1) % TREE_TRAVERSAL_STACK_CAPACITY;
                stackLength--;
                t0 = t0Stack[stackIndex];
                t1 = t1Stack[stackIndex];
                childIndex = childIndexStack[stackIndex];
                center = centerStack[stackIndex];
                nodeIndex = nodeIndexStack[stackIndex];
                node = nodePool[nodeIndex];
            }
            else if (!passNode && !IsEmpty(node))
            {
                if (IsLeaf(node))
                {
                    enter = childT0;
                    hit = true;
                    exitRoot = true;
                    break;
                }
                stackLength = min(stackLength + 1, TREE_TRAVERSAL_STACK_CAPACITY);
                t0Stack[stackIndex] = t0;
                t1Stack[stackIndex] = t1;
                childIndexStack[stackIndex] = childIndex;
                centerStack[stackIndex] = center;
                nodeIndexStack[stackIndex] = nodeIndex;
                stackIndex = (stackIndex + 1) % TREE_TRAVERSAL_STACK_CAPACITY;
                level++;
                center += int2(((childIndex & 1) << 1) - 1, (childIndex & 2) - 1) * (int)(rootResolution >> (level + 1));
                isFirstNode = true;
                t0 = childT0;
                t1 = childT1;
                nodeIndex = (int)node.data + childIndex;
                node = nodePool[nodeIndex];
            }
        }
        if (exitRoot)
        {
            break;
        }
    }
    return hit;
}

bool Zcast(
    RWStructuredBuffer<TreeNode> nodePool,
    uint rootResolution,
    int3 origin,
    out int enter,
    out int nodeIndex)
{
    bool result;
    enter = 0;
    nodeIndex = 0;
    result = ZcastCore(rootResolution, origin, nodePool, enter, nodeIndex);
    return result;
}

inline uint2 ShiftRight(uint2 code, int digit)
{
    return digit <= 0
        ? code
        : digit >= 32
            ? uint2(code.y >> (digit - 32), 0)
            : uint2((code.x >> digit) | (code.y << (32 - digit)), code.y >> digit);
}

inline uint2 ExpandSplit(int component)
{
    uint c = (uint)(component & 0x001FFFFF);
    uint2 x = uint2(c & 0x0000FFFFu, c & 0x001F0000u);
    x = (x | uint2(x.x << 0x10, (x.y << 0x10) | (x.x >> 0x10))) & uint2(0xFF0000FFu, 0x001F0000u);
    x = (x | uint2(x.x << 0x08, (x.y << 0x08) | (x.x >> 0x18))) & uint2(0x0F00F00Fu, 0x100F00F0u);
    x = (x | uint2(x.x << 0x04, (x.y << 0x04) | (x.x >> 0x1C))) & uint2(0xC30C30C3u, 0x10C30C30u);
    x = (x | uint2(x.x << 0x02, (x.y << 0x02) | (x.x >> 0x1E))) & uint2(0x49249249u, 0x12492492u);
    return x;
}

inline uint2 EncodeCoord(int3 coord)
{
    uint2 x = ExpandSplit(coord.x);
    uint2 y = ExpandSplit(coord.y);
    uint2 z = ExpandSplit(coord.z);
    y.y = (y.y << 1) | (y.x >> 31);
    y.x <<= 1;
    z.y = (z.y << 2) | (z.x >> 30);
    z.x <<= 2;
    return x | y | z;
}

int GetNode(RWStructuredBuffer<TreeNode> nodePool, uint2 code, int maxLevel, out int level)
{
    int nodeIndex = 0;
    level = 0;
    TreeNode node = nodePool[nodeIndex];
    while (!IsLeaf(node))
    {
        int childIndex = (int)ShiftRight(code, (maxLevel - level - 1) * 3).x & 7;
        nodeIndex = (int)node.data + childIndex;
        node = nodePool[nodeIndex];
        level++;
    }
    return nodeIndex;
}

RWStructuredBuffer<TreeNode> _NodePool;
RWStructuredBuffer<int> _Head;
RWStructuredBuffer<int> _TempHead;
RWStructuredBuffer<int> _Tail;
RWStructuredBuffer<bool> _Finished;
int _MaxLevel;
uint _Resolution;
uint _RootResolution;
uint2 _VoxelCoordOffset;

[numthreads(8, 8, 1)]
void Fill(uint3 localVoxelCoord : SV_DispatchThreadID)
{
    uint columnIndex = (localVoxelCoord.y * _Resolution) + localVoxelCoord.x;
    uint2 voxelCoord = localVoxelCoord.xy + _VoxelCoordOffset;
    int head = _Head[columnIndex];
    int tempHead = _TempHead[columnIndex];
    int tail = _Tail[columnIndex];
    bool continueLoop = true;
    for (int i = 0; continueLoop && (i < FILLING_LOOP_COUNT); i++)
    {
        switch (tail)
        {
            case -1:
                // モデル内部への入射点を探す
                int headEnter = 0;
                int headNodeIndex = 0;
                if (Zcast(
                    _NodePool,
                    _RootResolution,
                    int3((int2)voxelCoord, head),
                    headEnter,
                    headNodeIndex))
                {
                    TreeNode headNode = _NodePool[headNodeIndex];
                    int headCandidateFound = (int)(headNode.data >> 31) & 1;
                    head = headEnter + 1;
                    tail -= headCandidateFound;
                }
                else
                {
                    head = (int)_RootResolution;
                    continueLoop = false;
                }
                tempHead = head;
                break;
            case -2:
                // モデル外部への出射点を探す
                int tailEnter = 0;
                int tailNodeIndex = 0;
                if (Zcast(
                    _NodePool,
                    _RootResolution,
                    int3((int2)voxelCoord, tempHead),
                    tailEnter,
                    tailNodeIndex))
                {
                    TreeNode tailNode = _NodePool[tailNodeIndex];
                    switch ((tailNode.data >> 30) & 3)
                    {
                        case 1:
                            tail = tailEnter;
                            break;
                        case 2:
                            head = tailEnter + 1;
                            tempHead = head;
                            break;
                        default:
                            tempHead = tailEnter + 1;
                            break;
                    }
                }
                else
                {
                    head = (int)_RootResolution;
                    tempHead = head;
                    continueLoop = false;
                }
                break;
            default:
                // 入射点と出射点の間のノードに充填フラグを立てる
                if (head < tail)
                {
                    int level = 0;
                    int nodeIndex = GetNode(
                        _NodePool,
                        EncodeCoord(int3((int2)voxelCoord, head)),
                        _MaxLevel,
                        level);
                    TreeNode node = _NodePool[nodeIndex];
                    FillNode(node);
                    _NodePool[nodeIndex] = node;
                    head += 1 << (_MaxLevel - level);
                }
                else
                {
                    head = tail + 1;
                    tail = -1;
                }
                tempHead = head;
                break;
        }
    }
    _Head[columnIndex] = head;
    _TempHead[columnIndex] = tempHead;
    _Tail[columnIndex] = tail;
    if (head < (int)_RootResolution)
    {
        _Finished[0] = false;
    }
}

MortonUtility.cs

using Unity.Burst;
using Unity.Mathematics;

namespace Carving
{
    // Jeroen Baert「Libmorton」をもとにしたビットマスクを使ったモートン階数のエンコード・デコード
    // https://github.com/Forceflow/libmorton/blob/main/include/libmorton/morton3D.h
    // 各レベル3ビット占有し、64ビットエンコードの場合は63ビットまで使用できるため
    // 21レベル...つまり各軸解像度0x200000(2097152)までエンコードできる
    // Libmortonのライセンス情報:
    /*
    MIT License

    Copyright (c) 2016 Jeroen Baert

    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:

    The above copyright notice and this permission notice shall be included in all
    copies or substantial portions of the Software.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
    SOFTWARE.
    */
    [BurstCompile]
    public static class MortonUtility
    {
        public const int MaxResolution = 0x200000;

        [BurstCompile]
        public static ulong Encode(in int3 coord)
        {
            return Expand(coord.x) | (Expand(coord.y) << 1) | (Expand(coord.z) << 2);
        }
        
        [BurstCompile]
        public static void Decode32(uint code, out int3 coord)
        {
            coord = new int3(Collapse32(code), Collapse32(code >> 1), Collapse32(code >> 2));
        }
        
        [BurstCompile]
        private static ulong Expand(int component)
        {
            var x = (ulong)(component & 0x001FFFFF);
            x = (x | (x << 0x20)) & 0x001F00000000FFFFul;
            x = (x | (x << 0x10)) & 0x001F0000FF0000FFul;
            x = (x | (x << 0x08)) & 0x100F00F00F00F00Ful;
            x = (x | (x << 0x04)) & 0x10C30C30C30C30C3ul;
            x = (x | (x << 0x02)) & 0x1249249249249249ul;
            return x;
        }
        
        [BurstCompile]
        private static int Collapse32(uint code)
        {
            var x = code & 0x09249249u;
            x = (x ^ (x >> 0x02)) & 0x030C30C3u;
            x = (x ^ (x >> 0x04)) & 0x0300F00Fu;
            x = (x ^ (x >> 0x08)) & 0x030000FFu;
            x = (x ^ (x >> 0x10)) & 0x000003FFu;
            return unchecked((int)x);
        }
    }
}

TrianglePartitioner.cs

using System;
using System.Runtime.InteropServices;
using Unity.Burst;
using Unity.Collections;
using Unity.Collections.LowLevel.Unsafe;
using Unity.Jobs;
using Unity.Mathematics;
using UnityEngine;

namespace Carving
{
    public class TrianglePartitioner : IDisposable
    {
        private NativeArray<float3> vertices;
        private NativeArray<int3> triangles;
        private NativeList<Partition> result;
        private float partitionSize;
        private readonly Allocator allocator;

        public TrianglePartitioner(
            NativeArray<float3> sourceVertices,
            NativeArray<int3> sourceTriangles,
            Bounds volumeBounds,
            int partitionCountPerAxis,
            int3 coordOffset = default,
            Allocator resultAllocator = Allocator.Persistent)
        {
            this.vertices = sourceVertices;
            this.triangles = sourceTriangles;
            var partitionCount = partitionCountPerAxis * partitionCountPerAxis * partitionCountPerAxis;
            this.result = new NativeList<Partition>(
                partitionCount,
                resultAllocator);
            var volumeMin = (float3)volumeBounds.min;
            var volumeSize = volumeBounds.size.x;
            this.partitionSize = volumeSize / partitionCountPerAxis;
            for (var i = 0; i < partitionCount; i++)
            {
                MortonUtility.Decode32((uint)i, out var coord);
                this.result.Add(
                    Partition.Create(
                        resultAllocator,
                        volumeMin + ((float3)coord * this.partitionSize),
                        coord + coordOffset));
            }
            this.allocator = resultAllocator;
        }

        public NativeArray<Partition> Result => this.result;

        public void Dispose()
        {
            this.ReleaseUnmanagedResources();
            GC.SuppressFinalize(this);
        }

        public void Resize(
            NativeArray<float3> sourceVertices,
            NativeArray<int3> sourceTriangles,
            Bounds volumeBounds,
            int partitionCountPerAxis,
            int3 coordOffset = default)
        {
            this.vertices = sourceVertices;
            this.triangles = sourceTriangles;
            var volumeMin = (float3)volumeBounds.min;
            var volumeSize = volumeBounds.size.x;
            this.partitionSize = volumeSize / partitionCountPerAxis;
            var previousPartitionCount = this.result.Length;
            var partitionCount = partitionCountPerAxis * partitionCountPerAxis * partitionCountPerAxis;
            for (var i = partitionCount; i < previousPartitionCount; i++)
            {
                this.result[i].Dispose();
            }
            this.result.Resize(partitionCount, NativeArrayOptions.UninitializedMemory);
            for (var i = previousPartitionCount; i < partitionCount; i++)
            {
                MortonUtility.Decode32((uint)i, out var coord);
                this.result[i] = Partition.Create(
                    this.allocator,
                    volumeMin + ((float3)coord * this.partitionSize),
                    coord + coordOffset);
            }
            var partitionCountToClear = math.min(partitionCount, previousPartitionCount);
            for (var i = 0; i < partitionCountToClear; i++)
            {
                var partition = this.result[i];
                MortonUtility.Decode32((uint)i, out var coord);
                partition.Min = volumeMin + ((float3)coord * this.partitionSize);
                partition.Coord = coord + coordOffset;
                partition.Triangles.Clear();
                this.result[i] = partition;
            }
        }

        public unsafe JobHandle ExecutePartitioning()
        {
            return new TrianglePartitioningJob
            {
                Vertices = this.vertices,
                Triangles = this.triangles,
                PartitionSize = this.partitionSize,
                Result = (Partition*)this.result.GetUnsafeReadOnlyPtr()
            }.Schedule(this.result.Length, 1);
        }

        private void ReleaseUnmanagedResources()
        {
            if (this.result.IsCreated)
            {
                for (var i = 0; i < this.result.Length; i++)
                {
                    this.result[i].Dispose();
                }
                this.result.Dispose();
            }
        }

        ~TrianglePartitioner()
        {
            this.ReleaseUnmanagedResources();
        }

        // ボクセル化の単位となる小区画
        [StructLayout(LayoutKind.Sequential)]
        public struct Partition : IDisposable
        {
            public float3 Min;
            public int3 Coord;
            public UnsafeList<int> Triangles;

            public static Partition Create(Allocator allocator, float3 min, int3 coord)
            {
                return new Partition
                {
                    Min = min,
                    Coord = coord,
                    Triangles = new UnsafeList<int>(4, allocator)
                };
            }

            public void Dispose()
            {
                this.Triangles.Dispose();
            }
        }

        [BurstCompile]
        private unsafe struct TrianglePartitioningJob : IJobParallelFor
        {
            [ReadOnly] public NativeArray<float3> Vertices;
            [ReadOnly] public NativeArray<int3> Triangles;
            public float PartitionSize;
            [NativeDisableUnsafePtrRestriction] public Partition* Result;

            public void Execute(int partitionIndex)
            {
                var partition = this.Result + partitionIndex;
                for (var j = 0; j < this.Triangles.Length; j++)
                {
                    var t = this.Triangles[j];
                    if (VoxelUtility.Intersects(
                            partition->Min,
                            this.PartitionSize,
                            this.Vertices[t.x],
                            this.Vertices[t.y],
                            this.Vertices[t.z]))
                    {
                        partition->Triangles.Add(j);
                    }
                }
            }
        }
    }
}

Volume.shader

// Audun Wilhelmsen「Efficient Ray Tracing of Sparse Voxel Octrees on an FPGA」のノード探索方法をアレンジしたもの
// https://ntnuopen.ntnu.no/ntnu-xmlui/bitstream/handle/11250/2370620/567036_FULLTEXT01.pdf
Shader "Carving/Volume"
{
    Properties
    {
        _Color ("Color", Color) = (1.0, 1.0, 1.0, 1.0)
        [Toggle(_CUBIC_APPEARANCE)] _CubicAppearance ("Cubic Appearance", Float) = 0.0
        [Toggle(_HIDE_INTERNAL_NODES)] _HideInternalNodes ("Hide internal Nodes", Float) = 0.0
        _ClipMin ("Clip Min", Vector) = (0.0, 0.0, 0.0, 0.0)
        _ClipMax ("Clip Max", Vector) = (1.0, 1.0, 1.0, 1.0)
    }
    SubShader
    {
        Tags { "Queue" = "AlphaTest" "RenderType" = "Opaque" }

        Pass
        {
            Tags { "LightMode" = "ForwardBase" }

            Cull Back
            ZTest Always

            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #pragma require compute
            #pragma exclude_renderers gles
            #pragma multi_compile_local _ _CUBIC_APPEARANCE
            #pragma multi_compile_local _ _HIDE_INTERNAL_NODES
            #include "UnityCG.cginc"
            #include "Lighting.cginc"

            #define TREE_TRAVERSAL_STACK_CAPACITY 8
            #define DIRECTIONAL_EPSILON 0.0009765625
            #define EDGE_TINT_EXPONENT 8.0

            struct TreeNode
            {
                uint flags;
                uint data;
            };

            fixed4 _Color;
            float3 _ClipMin;
            float3 _ClipMax;
            float3 _VolumeMin;
            float _VolumeSize;
            float _NormalizedLeafSize;
            sampler2D _CameraDepthTexture;
            StructuredBuffer<TreeNode> _NodePool;

            bool isEmpty(TreeNode node) {return (node.flags & 1) > 0;}
            bool isLeaf(TreeNode node) {return (node.flags & 2) > 0;}
            bool isFilled(TreeNode node) {return (node.flags & 4) > 0;}
            float3 worldToObjectPos(float3 worldPos) {return mul(unity_WorldToObject, float4(worldPos, 1.0)).xyz;}
            float3 objectToWorldPos(float3 localPos) {return mul(unity_ObjectToWorld, float4(localPos, 1.0)).xyz;}
            float3 worldToObjectVec(float3 worldVec) {return mul((float3x3)unity_WorldToObject, worldVec);}
            float3 objectToWorldVec(float3 localVec) {return mul((float3x3)unity_ObjectToWorld, localVec);}
            float3 getWorldCameraPos() {return UNITY_MATRIX_I_V._14_24_34;}
            float3 getWorldCameraDir() {return -UNITY_MATRIX_V._31_32_33;}
            float3 getLocalCameraPos() {return worldToObjectPos(getWorldCameraPos());}
            float3 getLocalCameraDir() {return normalize(worldToObjectVec(getWorldCameraDir()));}
            bool isPerspective() {return any(UNITY_MATRIX_P._41_42_43);}
            float3 getLocalCameraOrigin(float3 localPos)
            {
                float3 localCameraPos = getLocalCameraPos();
                float3 localCameraDir = getLocalCameraDir();
                return isPerspective()
                    ? localCameraPos
                    : localPos + dot(localCameraPos - localPos, localCameraDir) * localCameraDir;
            }
            bool3 bitsToBool3(int childIndex) {return bool3((childIndex & 1) > 0, (childIndex & 2) > 0, (childIndex & 4) > 0);}
            float getMaxComponent(float3 v) {return max(max(v.x, v.y), v.z);}
            float getMinComponent(float3 v) {return min(min(v.x, v.y), v.z);}
            int getFirstChildIndex(float3 t0, float3 tm) {return (int)dot(lerp(0.0, float3(1.0, 2.0, 4.0), tm < getMaxComponent(t0)), 1.0);}
            int getNextChildIndex(float3 tm, float3 t1, int childIndex)
            {
                bool3 position = bitsToBool3(childIndex);
                t1 = lerp(tm, t1, position);
                return (t1.x < t1.y) && (t1.x < t1.z)
                    ? position.x ? 8 : childIndex ^ 1
                    : t1.y < t1.z
                        ? position.y ? 8 : childIndex ^ 2
                        : position.z ? 8 : childIndex ^ 4;
            }
            
            bool raycastCore(
                float3 rootT0,
                float3 rootT1,
                float2 enterExitRange,
                float3 clippedT0,
                int directionMask,
                StructuredBuffer<TreeNode> nodePool,
                int rootNodeIndex,
                out float enter,
                out TreeNode foundNode
                #ifdef _CUBIC_APPEARANCE
                , out float3 normal
                , out float edgeTint
                #endif
                )
            {
                enter = 0.0;
                foundNode = (TreeNode)0;
                #ifdef _CUBIC_APPEARANCE
                normal = 0.0;
                edgeTint = 1.0;
                #endif
                bool hit = false;
                float3 t0Stack[TREE_TRAVERSAL_STACK_CAPACITY];
                float3 t1Stack[TREE_TRAVERSAL_STACK_CAPACITY];
                int nodeIndexStack[TREE_TRAVERSAL_STACK_CAPACITY];
                int childIndexStack[TREE_TRAVERSAL_STACK_CAPACITY];
                int childIndex = 0;
                float tStart = enterExitRange.x;
                while (true)
                {
                    float3 t0 = rootT0;
                    float3 t1 = rootT1;
                    int nodeIndex = rootNodeIndex;
                    bool isFirstNode = true;
                    int level = 0;
                    uint stackIndex = 0;
                    uint stackLength = 0;
                    TreeNode node = nodePool[nodeIndex];
                    bool exitRoot = false;
                    while (true)
                    {
                        float3 tm = (t0 + t1) * 0.5;
                        bool exitNode = false;
                        bool stackUnderflow = stackLength == 0;
                        childIndex = isFirstNode ? getFirstChildIndex(t0, tm) : getNextChildIndex(tm, t1, childIndex);
                        isFirstNode = false;
                        exitNode = childIndex == 8;
                        bool3 position = bitsToBool3(childIndex);
                        int flippedChildIndex = childIndex ^ directionMask;
                        float3 childT0 = lerp(t0, tm, position);
                        float3 childT1 = lerp(tm, t1, position);
                        float childTEnter = getMaxComponent(childT0);
                        float childTExit = getMinComponent(childT1);
                        bool passNode = (childTExit <= tStart) || (childTEnter >= enterExitRange.y);
                        if (exitNode)
                        {
                            exitRoot = exitRoot || level == 0;
                            if (stackUnderflow || exitRoot)
                            {
                                tStart = getMinComponent(t1);
                                break;
                            }
                            level--;
                            stackIndex = ((stackIndex + TREE_TRAVERSAL_STACK_CAPACITY) - 1) % TREE_TRAVERSAL_STACK_CAPACITY;
                            stackLength--;
                            t0 = t0Stack[stackIndex];
                            t1 = t1Stack[stackIndex];
                            childIndex = childIndexStack[stackIndex];
                            nodeIndex = nodeIndexStack[stackIndex];
                            node = nodePool[nodeIndex];
                        }
                        else if (!passNode)
                        {
                            if (isFilled(node) && isLeaf(node)
                            #ifdef _HIDE_INTERNAL_NODES
                            && ((node.data >> 30) != 3)
                            #endif
                            )
                            {
                                enter = childTEnter;
                                foundNode = node;
                                #ifdef _CUBIC_APPEARANCE
                                bool clipped = enterExitRange.x > enter;
                                float3 plane = clipped ? ((float3)(enterExitRange.x == clippedT0)) : ((float3)(enter == childT0));
                                normal = plane * (2.0 * (float3)bitsToBool3(directionMask) - 1.0);
                                t0 = max(t0, clippedT0);
                                float3x4 d = min(float3x4(
                                    abs(float4(t0.yz, t1.yz) - t0.xxxx),
                                    abs(float4(t0.zx, t1.zx) - t0.yyyy),
                                    abs(float4(t0.xy, t1.xy) - t0.zzzz)) / _NormalizedLeafSize, 1.0);
                                float3 edgeTint3 = min(min(min(d._11_21_31, d._12_22_32), d._13_23_33), d._14_24_34) * plane;
                                edgeTint = 1.0 - pow(1.0 - dot(edgeTint3, 1.0), EDGE_TINT_EXPONENT);
                                #endif
                                hit = true;
                                exitRoot = true;
                                break;
                            }
                            if (!isEmpty(node))
                            {
                                stackLength = min(stackLength + 1, TREE_TRAVERSAL_STACK_CAPACITY);
                                t0Stack[stackIndex] = t0;
                                t1Stack[stackIndex] = t1;
                                childIndexStack[stackIndex] = childIndex;
                                nodeIndexStack[stackIndex] = nodeIndex;
                                stackIndex = (stackIndex + 1) % TREE_TRAVERSAL_STACK_CAPACITY;
                                level++;
                                isFirstNode = true;
                                t0 = childT0;
                                t1 = childT1;
                                nodeIndex = node.data + flippedChildIndex;
                                node = nodePool[nodeIndex];
                            }
                        }
                    }
                    if (exitRoot)
                    {
                        break;
                    }
                }
                return hit;
            }

            bool raycast(
                StructuredBuffer<TreeNode> nodePool,
                float3 volumeMin,
                float volumeSize,
                float3 localOrigin,
                float3 localDirection,
                float2x3 volumeRange,
                out float enter,
                out TreeNode foundNode
                #ifdef _CUBIC_APPEARANCE
                , out float3 normal
                , out float edgeTint
                #endif
                )
            {
                float3 volumeCenter = volumeMin + (volumeSize * 0.5);
                float3 rayOrigin = ((localOrigin - volumeCenter) * 2.0) / volumeSize;
                float3 rayDirection = localDirection;
                rayDirection += lerp(0.0, DIRECTIONAL_EPSILON, rayDirection == 0.0);
                float3 directionSign = sign(rayDirection);
                rayDirection *= directionSign;
                rayOrigin *= directionSign;
                int directionMask = (int)dot(max(-directionSign, 0.0) * float3(1.0, 2.0, 4.0), 1.0);
                float3 t0 = (-1.0 - rayOrigin) / rayDirection;
                float3 t1 = (1.0 - rayOrigin) / rayDirection;
                float tEnter = getMaxComponent(t0);
                float tExit = getMinComponent(t1);
                float3 flipped = 0.5 * (1.0 - directionSign);
                float2x3 invertedVolumeRange = float2x3(
                    directionSign * volumeRange[0] + flipped,
                    directionSign * volumeRange[1] + flipped);
                float2x3 flippedVolumeRange = float2x3(
                    lerp(invertedVolumeRange[0], invertedVolumeRange[1], flipped),
                    lerp(invertedVolumeRange[1], invertedVolumeRange[0], flipped));
                flippedVolumeRange *= 2.0;
                flippedVolumeRange -= 1.0;
                float3 clippedT0 = (flippedVolumeRange[0] - rayOrigin) / rayDirection;
                float3 clippedT1 = (flippedVolumeRange[1] - rayOrigin) / rayDirection;
                float2 enterExitRange = float2(getMaxComponent(clippedT0), getMinComponent(clippedT1));
                bool result;
                enter = 0.0f;
                foundNode = (TreeNode)0;
                #ifdef _CUBIC_APPEARANCE
                normal = 0.0;
                edgeTint = 1.0;
                #endif
                result = (tEnter < tExit) && raycastCore(t0, t1, enterExitRange, clippedT0, directionMask, nodePool, 0, enter, foundNode
                #ifdef _CUBIC_APPEARANCE
                , normal, edgeTint
                #endif
                );
                enter *= 0.5f * volumeSize;
                return result;
            }

            float sampleOpaqueZ(float4 screenPos)
            {
                float rawDepth = SAMPLE_DEPTH_TEXTURE_PROJ(_CameraDepthTexture, UNITY_PROJ_COORD(screenPos));
                return isPerspective()
                    ? LinearEyeDepth(rawDepth)
                    : -dot(unity_CameraInvProjection._33_34, float2(_ProjectionParams.x * (rawDepth * 2.0 - 1.0), 1.0));
            }

            float3 decodeNormal(uint data)
            {
                float2 xy = float2((int)((data & 0x3FFF) << 18), (int)(((data >> 14) & 0x3FFF) << 18)) / 0x7FFFFFFF;
                float z = (1 - (int)(((data >> 28) & 1) << 1)) * sqrt(max(1.0 - dot(xy, xy), 0.0));
                return float3(xy, z);
            }

            struct v2f
            {
                float4 vertex : SV_POSITION;
                float3 localPos : TEXCOORD0;
                float4 screenPos : TEXCOORD1;
            };

            v2f vert(float4 vertex : POSITION)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(vertex);
                o.localPos = vertex.xyz;
                o.screenPos = ComputeScreenPos(o.vertex);
                return o;
            }

            half4 frag(v2f i) : SV_Target
            {
                float3 localOrigin = getLocalCameraOrigin(i.localPos);
                float3 localDir = normalize(i.localPos - localOrigin);
                float enter = 0.0;
                TreeNode foundNode = (TreeNode)0;
                #ifdef _CUBIC_APPEARANCE
                float3 normal = 0.0;
                float edgeTint = 1.0;
                #endif
                clip((float)raycast(_NodePool, _VolumeMin, _VolumeSize, localOrigin, localDir, float2x3(_ClipMin, _ClipMax), enter, foundNode
                #ifdef _CUBIC_APPEARANCE
                , normal, edgeTint
                #endif
                ) - 0.5);
                float3 worldDir = objectToWorldVec(localDir);
                float opaqueZ = sampleOpaqueZ(i.screenPos);
                clip(opaqueZ - length(worldDir * enter));
                float3 worldNormal = UnityObjectToWorldNormal(
                    #ifdef _CUBIC_APPEARANCE
                    normal
                    #else
                    decodeNormal(foundNode.data)
                    #endif
                );
                float3 voxelColor;
                #ifdef _CUBIC_APPEARANCE
                switch ((foundNode.data >> 30) & 3)
                {
                    case 0:
                        voxelColor = float3(0.0, 1.0, 0.0);
                        break;
                    case 1:
                        voxelColor = float3(1.0, 0.0, 0.0);
                        break;
                    case 2:
                        voxelColor = float3(0.0, 0.0, 1.0);
                        break;
                    case 3:
                        voxelColor = float3(1.0, 0.0, 1.0);
                        break;
                }
                #else
                voxelColor = _Color.rgb;
                #endif
                half3 color = max(dot(_WorldSpaceLightPos0.xyz, worldNormal), 0.0) * _LightColor0 * voxelColor
                #ifdef _CUBIC_APPEARANCE
                * edgeTint
                #endif
                + voxelColor * ShadeSH9(half4(worldNormal, 1.0));
                return half4(color, 1.0);
            }
            ENDCG
        }
    }
}

Voxelization.compute

#pragma kernel Voxelize

struct VoxelData
{
    double mortonCode;
    uint data;
};

struct TreeNode
{
    uint flags;
    uint data;
};

inline bool IsEmpty(TreeNode node)
{
    return (node.flags & 1) > 0;
}

inline bool IsLeaf(TreeNode node)
{
    return (node.flags & 2) > 0;
}

inline bool IsFilled(TreeNode node)
{
    return (node.flags & 4) > 0;
}

// VoxelUtility.Intersectsを移植したもの
// AABBと三角形の交差があればtrue、なければfalse
inline bool Intersects(float3 voxelMin, float voxelSize, float3 u0, float3 u1, float3 u2)
{
    float h = voxelSize * 0.5;
    float3 c = voxelMin + h;
    float3x3 v = float3x3(u0 - c, u1 - c, u2 - c);
    float3 vMin = min(min(v[0], v[1]), v[2]);
    float3 vMax = max(max(v[0], v[1]), v[2]);
    bool result = true;
    if (any((vMin > h) || (vMax < -h)))
    {
        result = false;
    }
    else
    {
        float3x3 f = float3x3(v[1] - v[0], v[2] - v[1], v[0] - v[2]);
        float3x3 identity = float3x3(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
        for (int j = 0; j < 3; j++)
        {
            for (int i = 0; i < 3; i++)
            {
                float3 a = cross(identity[i], f[j]);
                float r = dot(abs(a), 1.0) * h;
                float3 p = mul(v, a);
                if ((min(min(p.x, p.y), p.z) > r) || (max(max(p.x, p.y), p.z) < -r))
                {
                    result = false;
                    break;
                }
            }
            if (!result)
            {
                break;
            }
        }
        if (result)
        {
            float3 n = cross(f[0], f[1]);
            float3 d = sign(n) * h;
            result = (dot(n, -d - v[0]) > 0.0) ? false : (dot(n, d - v[0]) >= 0.0);
        }
    }
    return result;
}

inline uint2 ExpandSplit(int component)
{
    uint c = (uint)(component & 0x001FFFFF);
    uint2 x = uint2(c & 0x0000FFFFu, c & 0x001F0000u);
    x = (x | uint2(x.x << 0x10, (x.y << 0x10) | (x.x >> 0x10))) & uint2(0xFF0000FFu, 0x001F0000u);
    x = (x | uint2(x.x << 0x08, (x.y << 0x08) | (x.x >> 0x18))) & uint2(0x0F00F00Fu, 0x100F00F0u);
    x = (x | uint2(x.x << 0x04, (x.y << 0x04) | (x.x >> 0x1C))) & uint2(0xC30C30C3u, 0x10C30C30u);
    x = (x | uint2(x.x << 0x02, (x.y << 0x02) | (x.x >> 0x1E))) & uint2(0x49249249u, 0x12492492u);
    return x;
}

inline uint2 EncodeCoord(int3 coord)
{
    uint2 x = ExpandSplit(coord.x);
    uint2 y = ExpandSplit(coord.y);
    uint2 z = ExpandSplit(coord.z);
    y.y = (y.y << 1) | (y.x >> 31);
    y.x <<= 1;
    z.y = (z.y << 2) | (z.x >> 30);
    z.x <<= 2;
    uint2 code = x | y | z;
    return code;
}

inline uint EncodeNormal(float3 normal)
{
    int2 xy = (int2)((normal.xy * 0x1FFF) + 0.5);
    uint d = (uint)(xy.x & 0x3FFF);
    d |= (uint)(xy.y & 0x3FFF) << 14;
    d |= normal.z < 0.0 ? 1 << 28 : 0;
    return d;
}

inline float3 DecodeNormal(uint data)
{
    float2 xy = float2((int)((data & 0x3FFF) << 18), (int)(((data >> 14) & 0x3FFF) << 18)) / 0x7FFFFFFF;
    float z = (1 - (int)(((data >> 28) & 1) << 1)) * sqrt(max(1.0 - dot(xy, xy), 0.0));
    return float3(xy, z);
}

inline uint2 ShiftRight(uint2 code, int digit)
{
    return digit <= 0
        ? code
        : digit >= 32
            ? uint2(code.y >> (digit - 32), 0)
            : uint2((code.x >> digit) | (code.y << (32 - digit)), code.y >> digit);
}

// 特定のモートンコードの位置に存在するノードをツリー内から探す
int GetNode(StructuredBuffer<TreeNode> nodePool, int subRoot, int subRootLevel, uint2 code, int maxLevel)
{
    int nodeIndex = subRoot;
    int level = subRootLevel;
    TreeNode node = nodePool[nodeIndex];
    while (!IsLeaf(node))
    {
        int childIndex = (int)ShiftRight(code, (maxLevel - level - 1) * 3).x & 7;
        nodeIndex = (int)node.data + childIndex;
        node = nodePool[nodeIndex];
        level++;
    }
    return nodeIndex;
}

// ある点がポリゴンメッシュの内部に存在するかを判定する
// メッシュは多様体であることを前提とし、その点から+Z方向にレイキャストしたときの
// ポリゴンとの交差回数の偶奇から内外を判定する
// 全ポリゴンを走査する仕様のため、メッシュはなるべくローポリゴンであることが望ましい
bool Contains(float3 p, StructuredBuffer<float3> vertices, StructuredBuffer<int3> triangles, uint triangleCount)
{
    int result = 0;
    float3 d = float3(0.0f, 0.0f, 1.0f);
    for (uint i = 0; i < triangleCount; i++)
    {
        int3 t = triangles[i];
        float3 v0 = vertices[t.x];
        float3 v1 = vertices[t.y];
        float3 v2 = vertices[t.z];
        float3 r = p - v0;
        float3 e1 = v1 - v0;
        float3 e2 = v2 - v0;
        float3 c1 = cross(r, e1);
        float3 c2 = cross(d, e2);
        float3 tuv = float3(dot(c1, e2), dot(c2, r), dot(c1, d)) / dot(c2, e1);
        result += (all(tuv >= 0.0) && (dot(tuv.yz, 1.0) <= 1.0)) ? 1 : 0;
    }
    return (result & 1) > 0;
}

uint CombineData(uint data0, uint data1)
{
    float3 n = DecodeNormal(data0) + DecodeNormal(data1);
    float nLengthSq = dot(n, n);
    if (nLengthSq > 0.0)
    {
        n /= sqrt(nLengthSq);
    }
    uint facingBits = (data0 >> 30) | (data1 >> 30);
    facingBits ^= (facingBits < 3) ? 0 : 3;
    return EncodeNormal(n) | (facingBits << 30);
}

int _Resolution;
float _PartitionSize;
float3 _PartitionMin;
int3 _PartitionCoord;
int2 _TriangleIndexToVoxelize;
int _SubRoot;
int _SubRootLevel;
float _FaceSign;
int _MaxLevel;
uint _TriangleCount;
StructuredBuffer<float3> _SourceVertices;
StructuredBuffer<float3> _SourceNormals;
StructuredBuffer<int3> _Triangles;
StructuredBuffer<int> _TrianglesToVoxelize;
AppendStructuredBuffer<VoxelData> _VoxelDataStorage;
StructuredBuffer<TreeNode> _NodePool;

[numthreads(8, 8, 8)]
void Voxelize(uint3 voxelCoord : SV_DispatchThreadID)
{
    uint3 globalVoxelCoord = (_PartitionCoord * _Resolution) + voxelCoord;
    uint2 globalCodeUint = EncodeCoord(globalVoxelCoord);
    double globalCode = asdouble(globalCodeUint.x, globalCodeUint.y);
    bool subtractive = _SubRoot >= 0;
    TreeNode existingNode = (TreeNode)0;
    if (subtractive)
    {
        existingNode = _NodePool[GetNode(_NodePool, _SubRoot, _SubRootLevel, globalCodeUint, _MaxLevel)];
    }
    bool voxelFound = false;
    VoxelData foundVoxel = (VoxelData)0;
    foundVoxel.mortonCode = globalCode;
    if (subtractive && !IsFilled(existingNode))
    {
        // 削り取り時は空の領域には新規ボクセルを配置する必要はないので、ここで処理を終える
        return;
    }
    float voxelSize = _PartitionSize / _Resolution;
    float voxelExtents = voxelSize * 0.5;
    float3 voxelMin = _PartitionMin + (voxelSize * (float3)voxelCoord);
    float3 normal = 0.0;
    uint facingBits = 0;
    for (int i = _TriangleIndexToVoxelize.x; i < _TriangleIndexToVoxelize.y; i++)
    {
        int3 t = _Triangles[_TrianglesToVoxelize[i]];
        float3x3 vertices = float3x3(
        _SourceVertices[t.x],
        _SourceVertices[t.y],
        _SourceVertices[t.z]);
        if (!Intersects(voxelMin, voxelSize, vertices[0], vertices[1], vertices[2]))
        {
            continue;
        }
        float3 fx = vertices[1] - vertices[0];
        float3 fy = vertices[2] - vertices[0];
        float3 fz = cross(fx, fy);
        float det = dot(fz, fz);
        if (det <= 0.0)
        {
            continue;
        }   
        // ボクセルが面積を持つポリゴンと交差している場合、ボクセルの座標をポリゴンの重心座標に変換して
        // その比率で法線を合成する
        // また、ポリゴンがZ方向に対して表を向いているか裏を向いているかを調べてfacingBitsに積算していく
        // こちらの情報は内部充填時に内部入射点・外部出射点を判定するのに使う
        // facingBitsは0でZ符号0、1で+Z、2で-Zを表し、0よりも+Zを優先、+Zよりも-Zを優先するルールにする
        float3x3 localToFace = float3x3(cross(fy, fz), cross(fz, fx), cross(fx, fy)) / det;
        float2 xy = mul(localToFace, voxelMin + voxelExtents - vertices[0]).xy;
        float xySum = dot(xy, 1.0);
        float z = 1.0 - xySum;
        float3x3 normals = float3x3(
        _SourceNormals[t.x],
        _SourceNormals[t.y],
        _SourceNormals[t.z]);
        normal += mul(float3(z, xy), normals);
        facingBits |= (uint)(int)(sign(fz.z) * _FaceSign);
        voxelFound = true;
    }
    if (voxelFound)
    {
        float normalLengthSq = dot(normal, normal);
        normal = normalLengthSq > 0.0 ? normal / sqrt(normalLengthSq) : 0.0;
        facingBits ^= (facingBits < 3) ? 0 : 1;
        foundVoxel.data = EncodeNormal(_FaceSign * normal) | (facingBits << 30);
    }
    if (subtractive && !IsEmpty(existingNode) && IsLeaf(existingNode))
    {
        if (voxelFound)
        {
            // 削り取り時、既存のボクセルと新規ボクセルが競合した場合は
            // 頂点法線は両者の平均とするが、Z符号については-Z側を覆い被せるようにする
            // +Zボクセルが-Z側に露出するのを防ぐことで、内部充填時の充填ミス...本来はモデル外部であるはずの
            // 領域を埋め立ててしまう現象を回避したつもりだが、はたしてこれでいいのかは十分検証していない
            uint sourceBits = foundVoxel.data & 0xC0000000;
            uint destinationBits = existingNode.data & 0xC0000000;
            foundVoxel.data = CombineData(foundVoxel.data, existingNode.data);
            uint combinedNormal = foundVoxel.data & 0x3FFFFFFF;
            sourceBits |= combinedNormal;
            destinationBits |= combinedNormal;
            foundVoxel.data = ((sourceBits >> 31) > 0) ? sourceBits : ((((sourceBits >> 30) & 1) > 0) ? destinationBits : foundVoxel.data);
        }
        else 
        {
            // 新規ボクセルなし・既存ボクセルありの場合は、既存ボクセル中心がCarverモデルの内部にあるかどうかを判定し
            // 内部にない場合だけ既存ボクセルを有効なボクセルとして出力させる
            if (!Contains(voxelMin + voxelExtents, _SourceVertices, _Triangles, _TriangleCount))
            {
                foundVoxel.data = existingNode.data;
                voxelFound = true;
            }
        }
    }
    if (voxelFound)
    {
        _VoxelDataStorage.Append(foundVoxel);
    }
}

VoxelUtility.cs

using System;
using System.Runtime.InteropServices;
using Unity.Burst;
using Unity.Collections;
using Unity.Collections.LowLevel.Unsafe;
using Unity.Jobs;
using Unity.Mathematics;
using UnityEngine;
using UnityEngine.Rendering;

namespace Carving
{
    [BurstCompile]
    public static class VoxelUtility
    {
        private const int TriangleCountPerVoxelizationDispatchMax = 256;
        private const int ResolutionPerFillingDispatchMax = 1024;
        private const string VoxelizationKernel = "Voxelize";
        private const string FillingKernel = "Fill";

        private static int voxelizationKernelIndex = -1;
        private static int fillingKernelIndex = -1;
        private static bool computeShaderIsEnabled;
        private static ComputeBuffer sourceVerticesBuffer;
        private static ComputeBuffer sourceNormalsBuffer;
        private static ComputeBuffer sourceTrianglesBuffer;
        private static ComputeBuffer trianglesToVoxelizeBuffer;
        private static ComputeBuffer voxelDataStorageBuffer;
        private static ComputeBuffer resultCountBuffer;
        private static ComputeBuffer nodePoolBuffer;
        private static ComputeBuffer headBuffer;
        private static ComputeBuffer tempHeadBuffer;
        private static ComputeBuffer tailBuffer;
        private static ComputeBuffer finishedBuffer;
        private static ComputeShader voxelizationShader;
        private static ComputeShader fillingShader;
        private static uint3 voxelizationThreadsPerGroup;
        private static uint3 fillingThreadsPerGroup;
        private static VoxelData[] voxelDataStorageArray = new VoxelData[8];

        private static readonly int ResolutionProperty = Shader.PropertyToID("_Resolution");
        private static readonly int RootResolutionProperty = Shader.PropertyToID("_RootResolution");
        private static readonly int PartitionSizeProperty = Shader.PropertyToID("_PartitionSize");
        private static readonly int PartitionMinProperty = Shader.PropertyToID("_PartitionMin");
        private static readonly int PartitionCoordProperty = Shader.PropertyToID("_PartitionCoord");
        private static readonly int TriangleIndexToVoxelizeProperty = Shader.PropertyToID("_TriangleIndexToVoxelize");
        private static readonly int SourceVerticesProperty = Shader.PropertyToID("_SourceVertices");
        private static readonly int SourceNormalsProperty = Shader.PropertyToID("_SourceNormals");
        private static readonly int TrianglesProperty = Shader.PropertyToID("_Triangles");
        private static readonly int TrianglesToVoxelizeProperty = Shader.PropertyToID("_TrianglesToVoxelize");
        private static readonly int VoxelDataStorageProperty = Shader.PropertyToID("_VoxelDataStorage");
        private static readonly int SubRootProperty = Shader.PropertyToID("_SubRoot");
        private static readonly int SubRootLevelProperty = Shader.PropertyToID("_SubRootLevel");
        private static readonly int FaceSignProperty = Shader.PropertyToID("_FaceSign");
        private static readonly int NodePoolProperty = Shader.PropertyToID("_NodePool");
        private static readonly int HeadProperty = Shader.PropertyToID("_Head");
        private static readonly int TempHeadProperty = Shader.PropertyToID("_TempHead");
        private static readonly int TailProperty = Shader.PropertyToID("_Tail");
        private static readonly int FinishedProperty = Shader.PropertyToID("_Finished");
        private static readonly int MaxLevelProperty = Shader.PropertyToID("_MaxLevel");
        private static readonly int TriangleCountProperty = Shader.PropertyToID("_TriangleCount");
        private static readonly int VoxelCoordOffsetProperty = Shader.PropertyToID("_VoxelCoordOffset");
        private static readonly int[] Int3Args = new int[3];
        private static readonly int[] ResultCountArray = new int[3];

        // Tomas Akenine-Möller「Fast 3D Triangle-Box Overlap Testing」による
        // AABB対三角形のヒット判定を立方体専用に変更したもの
        [BurstCompile]
        public static bool Intersects(in float3 min, in float size, in float3 u0, in float3 u1, in float3 u2)
        {
            var h = size * 0.5f;
            var c = min + h;
            var v = new float3x3(u0 - c, u1 - c, u2 - c);
            var vMin = math.min(math.min(v.c0, v.c1), v.c2);
            var vMax = math.max(math.max(v.c0, v.c1), v.c2);
            if (math.any((vMin > h) | (vMax < -h)))
            {
                return false;
            }
            var f = new float3x3(v.c1 - v.c0, v.c2 - v.c1, v.c0 - v.c2);
            for (var j = 0; j < 3; j++)
            {
                for (var i = 0; i < 3; i++)
                {
                    var a = math.cross(float3x3.identity[i], f[j]);
                    var r = math.csum(math.abs(a)) * h;
                    var p = math.mul(a, v);
                    if ((math.cmin(p) > r) || (math.cmax(p) < -r))
                    {
                        return false;
                    }
                }
            }
            var n = math.cross(f.c0, f.c1);
            var d = math.sign(n) * h;
            if (math.dot(n, -d - v.c0) > 0.0f)
            {
                return false;
            }
            return math.dot(n, d - v.c0) >= 0.0f;
        }

        public static void DisableComputeShader()
        {
            if (!computeShaderIsEnabled)
            {
                return;
            }
            sourceVerticesBuffer!.Dispose();
            sourceNormalsBuffer!.Dispose();
            sourceTrianglesBuffer!.Dispose();
            trianglesToVoxelizeBuffer!.Dispose();
            voxelDataStorageBuffer!.Dispose();
            resultCountBuffer!.Dispose();
            nodePoolBuffer!.Dispose();
            headBuffer!.Dispose();
            tempHeadBuffer!.Dispose();
            tailBuffer!.Dispose();
            finishedBuffer!.Dispose();
            sourceVerticesBuffer = null;
            sourceNormalsBuffer = null;
            sourceTrianglesBuffer = null;
            trianglesToVoxelizeBuffer = null;
            voxelDataStorageBuffer = null;
            nodePoolBuffer = null;
            headBuffer = null;
            tailBuffer = null;
            finishedBuffer = null;
            voxelizationShader = null;
            fillingShader = null;
            computeShaderIsEnabled = false;
        }

        public static void EnableComputeShader(ComputeShader shaderForVoxelization, ComputeShader shaderForFilling)
        {
            if (computeShaderIsEnabled || (shaderForVoxelization == null) || (shaderForFilling == null))
            {
                return;
            }
            voxelizationShader = shaderForVoxelization;
            voxelizationKernelIndex = voxelizationShader.FindKernel(VoxelizationKernel);
            if (!voxelizationShader.IsSupported(voxelizationKernelIndex))
            {
                Debug.LogError($"Kernel {voxelizationKernelIndex} of {voxelizationShader.name} is not supported.");
                return;
            }
            voxelizationShader.GetKernelThreadGroupSizes(
                voxelizationKernelIndex,
                out voxelizationThreadsPerGroup.x,
                out voxelizationThreadsPerGroup.y,
                out voxelizationThreadsPerGroup.z);
            fillingShader = shaderForFilling;
            fillingKernelIndex = fillingShader.FindKernel(FillingKernel);
            if (!fillingShader.IsSupported(fillingKernelIndex))
            {
                Debug.LogError($"Kernel {fillingKernelIndex} of {fillingShader.name} is not supported.");
                return;
            }
            fillingShader.GetKernelThreadGroupSizes(
                fillingKernelIndex,
                out fillingThreadsPerGroup.x,
                out fillingThreadsPerGroup.y,
                out fillingThreadsPerGroup.z);
            sourceVerticesBuffer = new ComputeBuffer(64, UnsafeUtility.SizeOf<float3>(), ComputeBufferType.Default);
            sourceNormalsBuffer = new ComputeBuffer(64, UnsafeUtility.SizeOf<float3>(), ComputeBufferType.Default);
            sourceTrianglesBuffer = new ComputeBuffer(64, UnsafeUtility.SizeOf<int3>(), ComputeBufferType.Default);
            trianglesToVoxelizeBuffer = new ComputeBuffer(64, sizeof(int), ComputeBufferType.Default);
            voxelDataStorageBuffer = new ComputeBuffer(64, UnsafeUtility.SizeOf<VoxelData>(), ComputeBufferType.Append);
            resultCountBuffer = new ComputeBuffer(3, sizeof(int), ComputeBufferType.IndirectArguments); // IndirectArgumentsは12バイト以上必要らしい
            resultCountBuffer.SetData(ResultCountArray);
            nodePoolBuffer = new ComputeBuffer(64, UnsafeUtility.SizeOf<TreeNode>(), ComputeBufferType.Default);
            headBuffer = new ComputeBuffer(64, sizeof(int), ComputeBufferType.Default);
            tempHeadBuffer = new ComputeBuffer(64, sizeof(int), ComputeBufferType.Default);
            tailBuffer = new ComputeBuffer(64, sizeof(int), ComputeBufferType.Default);
            finishedBuffer = new ComputeBuffer(1, sizeof(int), ComputeBufferType.Default);
            computeShaderIsEnabled = true;
        }

        public static void Voxelize(
            NativeArray<float3> sourceVertices,
            NativeArray<float3> sourceNormals,
            NativeArray<int3> sourceTriangles,
            UnsafeList<int> trianglesToVoxelize,
            int3 partitionCoord,
            float3 partitionMin,
            float partitionSize,
            int resolution,
            int rootResolution,
            NativeQueue<VoxelData> voxelDataStorage,
            bool needsUpdateSourceBuffers = true,
            int existingSubRoot = -1,
            int subRootResolution = 1,
            float faceSign = 1.0f)
        {
            if (!computeShaderIsEnabled)
            {
                Debug.LogError("Voxelization shader is not enabled.");
                return;
            }
            if (needsUpdateSourceBuffers)
            {
                ExtendBufferIfNeeded(ref sourceVerticesBuffer, sourceVertices.Length);
                sourceVerticesBuffer!.SetData(sourceVertices);
                ExtendBufferIfNeeded(ref sourceNormalsBuffer, sourceNormals.Length);
                sourceNormalsBuffer!.SetData(sourceNormals);
                ExtendBufferIfNeeded(ref sourceTrianglesBuffer, sourceTriangles.Length);
                sourceTrianglesBuffer!.SetData(sourceTriangles);
            }
            ExtendBufferIfNeeded(ref trianglesToVoxelizeBuffer, trianglesToVoxelize.Length);
            unsafe
            {
                var trianglesToVoxelizeArray = NativeArrayUnsafeUtility.ConvertExistingDataToNativeArray<int>(
                    trianglesToVoxelize.Ptr,
                    trianglesToVoxelize.Length,
                    Allocator.None);
                #if ENABLE_UNITY_COLLECTIONS_CHECKS
                var atomicSafetyHandle = AtomicSafetyHandle.Create();
                NativeArrayUnsafeUtility.SetAtomicSafetyHandle(ref trianglesToVoxelizeArray, atomicSafetyHandle);
                #endif
                trianglesToVoxelizeBuffer!.SetData(trianglesToVoxelizeArray);
                #if ENABLE_UNITY_COLLECTIONS_CHECKS
                AtomicSafetyHandle.Release(atomicSafetyHandle);
                #endif
            }
            ExtendBufferIfNeeded(ref voxelDataStorageBuffer, resolution * resolution * resolution, ComputeBufferType.Append);
            var triangleCountToVoxelize = trianglesToVoxelize.Length;
            var dispatchCount = math.max((int)math.ceil((double)triangleCountToVoxelize / TriangleCountPerVoxelizationDispatchMax), 1);
            voxelizationShader!.SetInt(ResolutionProperty, resolution);
            voxelizationShader.SetFloat(PartitionSizeProperty, partitionSize);
            voxelizationShader.SetVector(PartitionMinProperty, (Vector3)partitionMin);
            Int3Args![0] = partitionCoord.x;
            Int3Args[1] = partitionCoord.y;
            Int3Args[2] = partitionCoord.z;
            voxelizationShader.SetInts(PartitionCoordProperty, Int3Args);
            voxelizationShader.SetBuffer(voxelizationKernelIndex, SourceVerticesProperty, sourceVerticesBuffer);
            voxelizationShader.SetBuffer(voxelizationKernelIndex, SourceNormalsProperty, sourceNormalsBuffer);
            voxelizationShader.SetBuffer(voxelizationKernelIndex, TrianglesProperty, sourceTrianglesBuffer);
            voxelizationShader.SetBuffer(voxelizationKernelIndex, TrianglesToVoxelizeProperty, trianglesToVoxelizeBuffer);
            voxelizationShader.SetBuffer(voxelizationKernelIndex, VoxelDataStorageProperty, voxelDataStorageBuffer);
            voxelizationShader.SetBuffer(voxelizationKernelIndex, NodePoolProperty, nodePoolBuffer);
            voxelizationShader.SetInt(SubRootProperty, existingSubRoot);
            voxelizationShader.SetInt(SubRootLevelProperty, math.floorlog2(rootResolution / subRootResolution));
            voxelizationShader.SetFloat(FaceSignProperty, faceSign);
            voxelizationShader.SetInt(MaxLevelProperty, 32 - math.lzcnt(rootResolution - 1));
            voxelizationShader.SetInt(TriangleCountProperty, sourceTriangles.Length);
            var groupCount = (int3)resolution / (int3)voxelizationThreadsPerGroup;

            // 1回のディスパッチで処理するポリゴンが多すぎると
            // クラッシュしたりするので、手頃な分量に分割する
            for (var i = 0; i < dispatchCount; i++)
            {
                voxelDataStorageBuffer!.SetCounterValue(0);
                Int3Args[0] = i * TriangleCountPerVoxelizationDispatchMax;
                Int3Args[1] = math.min(Int3Args[0] + TriangleCountPerVoxelizationDispatchMax, triangleCountToVoxelize);
                if (dispatchCount > 1)
                {
                    Debug.Log($"Dispatch: {i + 1}/{dispatchCount}, Triangles: {Int3Args[1]}/{triangleCountToVoxelize}");
                }
                voxelizationShader.SetInts(TriangleIndexToVoxelizeProperty, Int3Args);
                voxelizationShader.Dispatch(voxelizationKernelIndex, groupCount.x, groupCount.y, groupCount.z);
                ComputeBuffer.CopyCount(voxelDataStorageBuffer, resultCountBuffer, 0);
                resultCountBuffer!.GetData(ResultCountArray);
                var resultCount = ResultCountArray![0];
                if (voxelDataStorageArray!.Length < resultCount)
                {
                    var newLength = voxelDataStorageArray.Length << 1;
                    while (newLength < resultCount)
                    {
                        newLength <<= 1;
                    }
                    Array.Resize(ref voxelDataStorageArray, newLength);
                }
                voxelDataStorageBuffer.GetData(voxelDataStorageArray, 0, 0, resultCount);
                for (var j = 0; j < resultCount; j++)
                {
                    voxelDataStorage.Enqueue(voxelDataStorageArray![j]);
                }
            }
        }

        public static void Fill(
            NativeList<TreeNode> nodePool,
            int rootResolution,
            int nodeResolution = default,
            int2 nodeCoord = default)
        {
            if (computeShaderIsEnabled)
            {
                UploadNodePool(nodePool);
                var resolution = nodeResolution > 0 ? nodeResolution : rootResolution;
                var dispatchCountPerAxis = (int)math.ceil(
                    (double)resolution / ResolutionPerFillingDispatchMax);
                var patchResolution = math.min(resolution, ResolutionPerFillingDispatchMax);
                var patchArea = patchResolution * patchResolution;
                var groupCounts = math.int2(patchResolution) / (int2)fillingThreadsPerGroup.xy;
                ExtendBufferIfNeeded(ref headBuffer, patchArea);
                ExtendBufferIfNeeded(ref tempHeadBuffer, patchArea);
                ExtendBufferIfNeeded(ref tailBuffer, patchArea);
                var finishedValue = new int[1];
                var initialFinishedValue = new[] { 1 };
                var initialHeadValues = new int[patchArea];
                var initialTailValues = new int[patchArea];
                Array.Fill(initialTailValues, -1);
                fillingShader!.SetBuffer(fillingKernelIndex, NodePoolProperty, nodePoolBuffer);
                fillingShader.SetBuffer(fillingKernelIndex, HeadProperty, headBuffer);
                fillingShader.SetBuffer(fillingKernelIndex, TempHeadProperty, tempHeadBuffer);
                fillingShader.SetBuffer(fillingKernelIndex, TailProperty, tailBuffer);
                fillingShader.SetBuffer(fillingKernelIndex, FinishedProperty, finishedBuffer);
                fillingShader.SetInt(MaxLevelProperty, 32 - math.lzcnt(rootResolution - 1));
                fillingShader.SetInt(ResolutionProperty, patchResolution);
                fillingShader.SetInt(RootResolutionProperty, rootResolution);
                for (var j = 0; j < dispatchCountPerAxis; j++)
                {
                    for (var i = 0; i < dispatchCountPerAxis; i++)
                    {
                        headBuffer!.SetData(initialHeadValues);
                        tempHeadBuffer!.SetData(initialHeadValues);
                        tailBuffer!.SetData(initialTailValues);
                        Int3Args![0] = (resolution * nodeCoord.x) + (patchResolution * i);
                        Int3Args[1] = (resolution * nodeCoord.y) + (patchResolution * j);
                        fillingShader.SetInts(VoxelCoordOffsetProperty, Int3Args);
                        var k = 0;
                        do
                        {
                            k++;
                            Debug.Log($"Dispatch filling for ({i}, {j}), trial: {k}");
                            finishedBuffer!.SetData(initialFinishedValue);
                            fillingShader.Dispatch(fillingKernelIndex, groupCounts.x, groupCounts.y, 1);
                            finishedBuffer.GetData(finishedValue);
                            if (k > 10000)
                            {
                                Debug.LogWarning("Too many dispatch loops.");
                                break;
                            }
                        } while (finishedValue[0] == 0);
                    }
                }
                var nodePoolArray = nodePool.AsArray();
                AsyncGPUReadback.RequestIntoNativeArray(
                    ref nodePoolArray,
                    nodePoolBuffer,
                    nodePool.Length * UnsafeUtility.SizeOf<TreeNode>(),
                    0).WaitForCompletion();
            }
            else
            {
                Debug.LogError("Filling shader is not enabled.");
            }
        }

        public static void UploadNodePool(NativeList<TreeNode> nodePool)
        {
            ExtendBufferIfNeeded(ref nodePoolBuffer, nodePool.Length);
            nodePoolBuffer!.SetData(nodePool.AsArray());
        }

        public static void KillNode(NativeList<TreeNode> nodePool, int nodeIndex, NativeQueue<int> deadNodes)
        {
            var node = nodePool[nodeIndex];
            nodePool[nodeIndex] = TreeNode.Create();
            if (!node.IsLeaf)
            {
                deadNodes.Enqueue(node.FirstChild);
                for (var i = 0; i < 8; i++)
                {
                    KillNode(nodePool, node.FirstChild + i, deadNodes);
                }
            }
        }

        public static unsafe bool SubdivideNodesOneLevel(
            int rootLevel,
            NativeArray<VoxelData> fragments,
            NativeList<TreeNode> nodePool,
            NativeQueue<int>? deadNodes = null)
        {
            using var requests = new NativeHashSet<int>(fragments.Length, Allocator.TempJob);
            new NodeMarkingJob
            {
                RootLevel = rootLevel,
                Fragments = fragments,
                NodePool = nodePool.GetUnsafeList()->Ptr,
                SubdivisionRequests = requests.AsParallelWriter()
            }.Schedule(fragments.Length, 8).Complete();
            var subdivided = !requests.IsEmpty;
            if (subdivided)
            {
                foreach (var i in requests)
                {
                    var node = nodePool[i];
                    node.Flags &= ~3u;
                    if (deadNodes.HasValue && deadNodes.Value.TryDequeue(out var deadNode))
                    {
                        node.FirstChild = deadNode;
                    }
                    else
                    {
                        node.FirstChild = nodePool.Length;
                        for (var j = 0; j < 8; j++)
                        {
                            nodePool.Add(TreeNode.Create());
                        }
                    }

                    nodePool[i] = node;
                }
            }
            return subdivided;
        }

        public static unsafe void EmbedVoxelData(int rootLevel, NativeArray<VoxelData> fragments, NativeList<TreeNode> nodePool)
        {
            new DataEmbeddingJob
            {
                RootLevel = rootLevel,
                Fragments = fragments,
                NodePool = nodePool.GetUnsafeList()->Ptr
            }.Schedule(fragments.Length, 8).Complete();
        }

        public static void ExtendBufferIfNeeded(ref ComputeBuffer buffer, int request, ComputeBufferType type = ComputeBufferType.Default)
        {
            var currentCapacity = buffer!.count;
            if (request <= currentCapacity)
            {
                return;
            }
            var newCapacity = currentCapacity * 2;
            while (newCapacity < request)
            {
                newCapacity *= 2;
            }
            var stride = buffer.stride;
            buffer.Dispose();
            buffer = new ComputeBuffer(newCapacity, stride, type);
        }

        private static uint CombineData(uint data0, uint data1)
        {
            var n = math.normalizesafe(DecodeNormal(data0) + DecodeNormal(data1));
            var facingBits = (data0 >> 30) | (data1 >> 30);
            facingBits ^= facingBits < 3 ? 0u : 3u;
            return EncodeNormal(n) | (facingBits << 30);
        }

        private static uint EncodeNormal(in float3 normal)
        {
            var xy = (int2)((normal.xy * 0x1FFF) + 0.5f);
            var d = (uint)(xy.x & 0x3FFF);
            d |= (uint)(xy.y & 0x3FFF) << 14;
            d |= normal.z < 0.0f ? 1u << 28 : 0u;
            return d;
        }

        private static float3 DecodeNormal(uint data)
        {
            var xy = math.float2((int)((data & 0x3FFF) << 18), (int)(((data >> 14) & 0x3FFF) << 18)) / 0x7FFFFFFF;
            var z = (1 - (int)(((data >> 28) & 1) << 1)) * math.sqrt(math.max(1.0f - math.dot(xy, xy), 0.0f));
            return math.float3(xy, z);
        }

        // 個々のボクセルに埋め込むデータ
        [StructLayout(LayoutKind.Sequential, Pack = 4)]
        public readonly struct VoxelData
        {
            public readonly ulong MortonCode;
            public readonly uint Data;
        }

        // 八分木のノード
        [StructLayout(LayoutKind.Explicit)]
        public struct TreeNode
        {
            private const uint InitialFlags = 3;

            [FieldOffset(0)] public uint Flags;
            [FieldOffset(4)] public int FirstChild;
            [FieldOffset(4)] public uint Data;

            public bool IsEmpty
            {
                get => (this.Flags & 1) > 0;
                set
                {
                    if (value)
                    {
                        this.Flags |= 1;
                    }
                    else
                    {
                        this.Flags &= ~1u;
                    }
                }
            }

            public bool IsLeaf => (this.Flags & 2) > 0;

            public void SetFilled(bool value)
            {
                if (value)
                {
                    this.Flags |= 4;
                }
                else
                {
                    this.Flags &= ~4u;
                }
            }

            public static TreeNode Create()
            {
                return new TreeNode
                {
                    Flags = InitialFlags,
                    FirstChild = -1
                };
            }
        }

        [BurstCompile]
        private unsafe struct NodeMarkingJob : IJobParallelFor
        {
            public int RootLevel;
            [ReadOnly] public NativeArray<VoxelData> Fragments;
            [NativeDisableUnsafePtrRestriction] public TreeNode* NodePool;
            [WriteOnly] public NativeHashSet<int>.ParallelWriter SubdivisionRequests;

            public void Execute(int fragmentIndex)
            {
                var fragment = this.Fragments[fragmentIndex];
                var currentNodeIndex = 0;
                for (var level = this.RootLevel; level >= 0; level--)
                {
                    var currentNode = this.NodePool + currentNodeIndex;
                    var firstChildIndex = currentNode->FirstChild;
                    if (firstChildIndex < 0)
                    {
                        this.SubdivisionRequests.Add(currentNodeIndex);
                        break;
                    }
                    currentNodeIndex = firstChildIndex + (int)((fragment.MortonCode >> (level * 3)) & 7);
                }
            }
        }

        [BurstCompile]
        private unsafe struct DataEmbeddingJob : IJobParallelFor
        {
            public int RootLevel;
            [ReadOnly] public NativeArray<VoxelData> Fragments;
            [NativeDisableUnsafePtrRestriction] public TreeNode* NodePool;

            public void Execute(int fragmentIndex)
            {
                var fragment = this.Fragments[fragmentIndex];
                var currentNodeIndex = 0;
                for (var level = this.RootLevel; level >= 0; level--)
                {
                    currentNodeIndex = (this.NodePool + currentNodeIndex)->FirstChild + (int)((fragment.MortonCode >> (level * 3)) & 7);
                }
                var targetNode = this.NodePool + currentNodeIndex;
                targetNode->Data = targetNode->IsEmpty ? fragment.Data : CombineData(targetNode->Data, fragment.Data);
                targetNode->IsEmpty = false;
                targetNode->SetFilled(true);
            }
        }
    }
}