Unityでボリュームレンダリング ~CTデータビジュアライゼーション~

Unityでボリュームレンダリングをするのに必要な実装について簡易的に解説してみます。方式としてはダイレクトボリュームレンダリング・レイキャスティングです。

 

今回は以下2編です。

データ入力編

カスタムシェーダー編

コンピュートシェーダー編

 

データ入力編

ボリュームデータは、3次元画像の形で表せます。Unityでは、Texture3Dとして3次元画像がサポートされていますので、そちらを使います。

CTデータは、3次元×CT値です。CTの値域は-1000 ~ 1000といった値をとりますので、通常の8bitカラーでは不足します。16bit浮動小数点のTextureFormat.RHalfあたりが適当でしょう。こちらはスカラー3次元画像用フォーマットです。

 

入力するボリュームデータは16bit整数(short)のraw形式であるとします。dicomを使いたかったらライブラリを探してきましょう。

rawはこんな感じで始まるデータ



入力の擬似コード

バイナリリーダーを用意し、ファイルを読み込みます。

BinaryReader reader = new BinaryReader(new FileStream(filePath, FileMode.Open));

3DTextureを作る前に、ボリュームデータを配列(1次元)に入れます。

float[] data = new float[dimX * dimY * dimZ];

ここで、ボリュームデータのエンディアン(バイトオーダー)に注意!データとOSの組み合わせによっては変換が必要です。

ushort dataval = reader.ReadUInt16();

トルエンディアンからビッグエンディアンへの変換はこんな感じ

dataval = (ushort)(((dataval >> 8) | (dataval << 8)) & 0xFFFF);

データ値は、0~1の範囲の値に変換するとテクスチャアクセスの都合上扱いやすいです。シェーダーでは、テクスチャに0~1でアクセスします。

dataval = (dataval - min) / range;

 

forでもなんでもいいので回して配列に入れられたら実際にTexture3Dを作ります。

Texture3D texture3d = new Texture3D(dimX, dimY, dimZ, TextureFormat.RHalf, false);
texture3d.filterMode = FilterMode.Bilinear;
texture3d.wrapModeW = TextureWrapMode.Clamp;

作ったTexture3Dにデータを流し込みます。

texture3d.SetPixelData<float>(data, 0);
texture3d.Apply();

 

最後にキューブのシェーダーにデータを渡します。(シェーダー作った後で)
meshRenderer.sharedMaterial.SetTexture("_DataTex", 3dtexture);

 

これでボリュームデータの入力は完成です。

 

※発展

色設定はTexture2D(RGBA)で表現します。テクスチャの左端が0~1に変換したデータの0、右端がデータの1に対応する色に塗ります。

 

カスタムシェーダー編

Unityで3Dレンダリングするためには、基本的にMeshRendererが必要です。

また、ボリュームは3Dオブジェクトで、データの形も立方体(直方体は一旦無視)なので、プリミティブのキューブ立方体を用意しましょう。

Unityでボリュームレンダリングする際は、キューブのメッシュレンダリングを、ボリュームレンダリングに書き換えます。

 

カスタムシェーダー

では、カスタムシェーダーを作成します。かなり難しいのでUnityドキュメントをよく読む必要があります。

 

Unlit Shaderを改変していきます。

まず、Properties。マテリアルの項目設定にあたります。この場合3Dテクスチャを1つと、0~1の値を2つ設定できるマテリアルになります。

Properties
{
    //マテリアル
    _DataTex("Data Texture", 3D) = "" {}
    _MinVal("Min val", Range(0.0, 1.0)) = 0.0
    _MaxVal("Max val", Range(0.0, 1.0)) = 1.0
}

次はシェーダー設定。透過シェーダーであることを設定しています。途中のブロックコメントはコードのフォーマッタで壊されないためにあります。

//レンダリングステート
Tags { "Queue" = "Transparent" "RenderType" = "Transparent"}
//カリングモード
Cull /* */ Front
//デプスバッファテストモード
ZTest /* */ LEqual
//デプスバッファ書き込みモード
ZWrite /* */ On
//ブレンドモード アルファ値乗算 フレームバッファの(1 - Source Alpha)を乗算 # 昔ながらの透明
Blend /* */ SrcAlpha /* */OneMinusSrcAlpha
そしてCGPROGRAMの中です。
Unity組み込み関数を入れておきます。
    //頂点シェーダとしてvertをコンパイル
    #pragma vertex vert
    //フラグメントシェーダとしてfragをコンパイル
    #pragma fragment frag
 
#include "UnityCG.cginc"
入出力の内容を設定しています。
struct vert_in
{
    UNITY_VERTEX_INPUT_INSTANCE_ID //GPUインスタンスID
        //命名 : セマンティクス
    float4 vertex : POSITION; // 頂点位置
    float4 normal : NORMAL; //法線ベクトル
    float2 uv : TEXCOORD0; // テクスチャ座標
};

struct frag_in
{
    UNITY_VERTEX_OUTPUT_STEREO //GPUインスタンスID
        //SV... システム値セマンティクス
    float4 vertex : SV_POSITION;
    float2 uv : TEXCOORD0; //0
    float3 vertexLocal : TEXCOORD1; //1
    float3 normal : NORMAL;
};

struct frag_out
{
    float4 colour : SV_TARGET; //レンダー ターゲットに格納される出力値
    float depth : SV_DEPTH; //深度バッファー データ。 ピクセル シェーダーで記述できます。
};
マテリアルにアクセスするための宣言。同じ名前なら勝手にリンクされます。
sampler3D _DataTex;

float _MinVal;
float _MaxVal;
ボリュームレンダリング用のレイ構造体です。
//レイ単体
struct RayInfo
{
    float3 startPos;
    float3 endPos;
    float3 direction;
        //交差距離 近・遠
    float2 aabbInters;
};
//レイマーチング(レイを含む)
struct RaymarchInfo
{
    RayInfo ray;
    int numSteps;
    float numStepsRecip;
    float stepSize;
};
頂点シェーダー。座標系変換が行われています。この先出てくる見慣れない関数は大体Unity組み込みかHLSLの組み込みです。
        // 頂点シェーダエントリ
frag_in vert(vert_in v)
{
  // フラグメントシェーダ入力
    frag_in o;
    UNITY_SETUP_INSTANCE_ID(v);
    UNITY_INITIALIZE_VERTEX_OUTPUT_STEREO(o);
    o.vertex = UnityObjectToClipPos(v.vertex);
    o.uv = v.uv;
    o.vertexLocal = v.vertex;
    o.normal = UnityObjectToWorldNormal(v.normal);
    return o;
}
レイ生成関数1。レイの方向がわかります。
float3 getViewRayDir(float3 vertexLocal)
{
    if (unity_OrthoParams.w == 0)
    {
            // Perspective
            // カメラへのオブジェクト空間方向 (正規化なし)
        return normalize(ObjSpaceViewDir(float4(vertexLocal, 0.0f)));
    }
    else
    {
            // Orthographic
        float3 camfwd = mul( (float3x3) unity_CameraToWorld, float3(0, 0, -1) );
        float4 camfwdobjspace = mul(unity_WorldToObject, camfwd);
        return normalize(camfwdobjspace);
    }
}
レイ生成関数2。ボックスとの交差地点がわかります。
// 軸平行境界ボックスより交差地点を見つける
// 近・遠の距離を返す
float2 intersectAABB(float3 rayOrigin, float3 rayDir, float3 boxMin, float3 boxMax)
{
    float3 tMin = (boxMin - rayOrigin) / rayDir;
    float3 tMax = (boxMax - rayOrigin) / rayDir;
    float3 t1 = min(tMin, tMax);
    float3 t2 = max(tMin, tMax);
    float tNear = max(max(t1.x, t1.y), t1.z);
    float tFar = min(min(t2.x, t2.y), t2.z);
    return float2(tNear, tFar);
}
レイ生成関数3 レイを生成します。プリミティブのキューブを使ったのは、キューブサイズを1に固定するためでした。これで、テクスチャ座標系の0~1でのレイができました。
RayInfo getRayFront2Back(float3 vertexLocal)
{
    RayInfo ray;
    ray.direction = getViewRayDir(vertexLocal);
    ray.endPos = vertexLocal + float3(0.5f, 0.5f, 0.5f);
        // Find intersections with axis aligned boundinng box (the volume)
        // (1, 1, 1)サイズボリューム
    ray.aabbInters = intersectAABB(ray.endPos, ray.direction, float3(0.0, 0.0, 0.0), float3(1.0f, 1.0f, 1.0));

        // Check if camera is inside AABB
        // ボリューム内カメラ調整
    const float3 farPos = ray.endPos + ray.direction * ray.aabbInters.y - float3(0.5f, 0.5f, 0.5f);
        // オブジェクト空間からカメラのクリップ空間へ点を変換
    float4 clipPos = UnityObjectToClipPos(float4(farPos, 1.0f));
    ray.aabbInters += min(clipPos.w, 0.0);
   
    ray.startPos = ray.endPos + ray.direction * ray.aabbInters.y;
    ray.direction = -ray.direction;
    return ray;
}
レイマーチングの構造体です。
        // レイマーチングを作成
RaymarchInfo initRaymarch(RayInfo ray, int maxNumSteps)
{
    RaymarchInfo raymarchInfo;
    raymarchInfo.stepSize = 1.0 / maxNumSteps;
    raymarchInfo.numSteps = (int) clamp(length(ray.startPos - ray.endPos) / raymarchInfo.stepSize, 1, maxNumSteps);
    raymarchInfo.numStepsRecip = 1.0 / raymarchInfo.numSteps;
    return raymarchInfo;
}
ボリュームデータ値へのアクセスです。0~1でアクセスします。
        // ボリュームスカラー値取得
float getDensity(float3 pos)
{
    return tex3Dlod(_DataTex, float4(pos.x, pos.y, pos.z, 0.0f));
}
深度取得の関数です。深度は他オブジェクトとの前後を決定するために使われます。
    // Converts local position to depth value
    // デプス取得
float localToDepth(float3 localPos)
{
    float4 clipPos = UnityObjectToClipPos(float4(localPos, 1.0f));

#if defined(SHADER_API_GLCORE) || defined(SHADER_API_OPENGL) || defined(SHADER_API_GLES) || defined(SHADER_API_GLES3)
        return (clipPos.z / clipPos.w) * 0.5 + 0.5;
#else
    return clipPos.z / clipPos.w;
#endif
}
ではコアの部分です。レイを作成し、ボリューム内をループで行進、カラーを積算し、最後に書き出しています。
frag_out frag(frag_in i)
{

    RayInfo ray = getRayFront2Back(i.vertexLocal);
   
    RaymarchInfo raymarchInfo = initRaymarch(ray, 128);

    float4 col = float4(0.0f, 0.0f, 0.0f, 0.0f);
    float tDepth = raymarchInfo.numStepsRecip * (raymarchInfo.numSteps - 1);
        // レイマーチング
        [loop]
    for (float iStep = 0; iStep < raymarchInfo.numSteps; iStep += 1)
    {
        const float t = iStep * raymarchInfo.numStepsRecip;
        const float3 currPos = lerp(ray.startPos, ray.endPos, t);
 
       const float density = getDensity(currPos);

            // Apply visibility window
            // 可視スカラー範囲
        if (density < _MinVal || density > _MaxVal) continue;

       float4 src = float4(1.0f, 1.0f, 1.0f, density);
// 前からブレンド
    src.rgb *= src.a; // サンプルカラー
        col = (1.0f - col.a) * src + col; // 積算値とブレンド
    }
   
    //透明時は0の深度
    float fragDepth = step(0.01f, col.a) * localToDepth(ray.startPos);
       
    //フラグメント書き出し
    frag_out output;
    output.colour = col;
    output.depth = fragDepth;
    return output;
}
以上でシェーダーの完成です。お好みでsrcカラーの値をdensityに応じて変えましょう。
 
作ったシェーダーは、まずマテリアルを作ってシェーダーを割り当て、そのマテリアルをキューブのMeshRendererに割り当てると使用可能です。
 
今回は陰影やライティングは全く考慮していない簡単なバージョンを紹介しました。
陰影やライティングを考慮する場合は、法線を生成し、ライト位置を取ってきてシェーディングしましょう。また、深度も正確なものにする必要がありますね。
 

今回つくったもの

シェーダー(マテリアル)
 

参考

(これの簡易化版だよ!)

github.com

 

コンピュートシェーダー編

ボリュームはデータ量が大きいので、全体処理するときはコンピュートシェーダーを使うと高速化できます。法線計算だったり0~1化だったり。ただしメインメモリとGPUメモリの複雑な操作が必要になります。