凹みTips

C++、JavaScript、Unity、ガジェット等の Tips について雑多に書いています。

Unity で DICOM 画像をレンダリングしてみる

はじめに

先日、Twitter で手持ちの DICOM 画像を以前書いたボリュームレンダリングの記事のように描画したい、というご相談をいただき試してみましたので内容を共有します。

tips.hecomi.com

DICOM(ダイコム)とは「Digital Imaging and COmmunications in Medicine」の略で、CT や MRI などのデータを保存する画像や通信プロトコルといった医用の標準規格です。今回の場合は、MRI のスライス画像が DICOM という形式のもと保存されているといった状況で、それをどうレンダリングすればよいのか、というお話でした。

ja.wikipedia.org

デモ

f:id:hecomi:20210927001258g:plain

f:id:hecomi:20210927001311g:plain

ダウンロード

github.com

情報収集

Unity で DICOM を扱うことに関してはインターネット上にいくつか情報がありました。

Unite Tokyo 2018

カンファレンスイベントである Unite Tokyo 2018 では、DICOM に関連した講演を東京大学医学部脳神経外科の金太一助教がされていました。

cgworld.jp

この中で、Kompath 社と協力して作られた DICOM を読み込める Simple DICOM Loader や、マーチングキューブ法でメッシュ化する High Speed CPU-based Marching cubes が紹介されていました。

assetstore.unity.com

assetstore.unity.com

fo-dicom の利用

また、Q&A サイトの teratail では正に似たような質問がされていました。

teratail.com

こちらでは回答者さんが fo-dicom という無料の DICOM ローダーを使って Texture3D 化をされていました。

方針

上記 fo-dicom を使って同様に Texture3D 化を行う方針で進めることにしました。

私の記事では PVM という形式で配布されているデータを RAW 形式に変換、更にそれをパースするみたいな部分に記事の多くの分量が割かれてしまっていましたが、本質はいくつかの工程を経てデータを Texture3D 化したあとの、それをシェーダでどのようにレンダリングするかというところにあるので、大きな方針としては変わりません。

つまり、DICOM データを fo-dicom を利用して Texture3D に変換するスクリプトを書き、その Texture3D を入力として以前の記事と同じコードでレンダリングを行います。

データのダウンロード

ご相談頂いた方のデータは諸事情でブログでは使えないため、本記事では配布されているデータで試してみることにします。上記質問回答の文中にもあるページから「Visible Female CT Datasets > Regional Tar Files (download) > Head」をダウンロード・展開しておきます。

mri.radiology.uiowa.edu

fo-dicom のインポート

fo-dicom 本体は GitHub にてオープンソースとして配布されています。

github.com

Unity 向けのパッケージは Asset Store にて配布されています。

assetstore.unity.com

Asset Store からパッケージをインストールし、プロジェクトに展開しましょう。これで fo-dicom をスクリプトから使えるようになります。

インポートスクリプトの記述とデータのインポート

方針

次にデータのインポートスクリプトを書きます。先程の展開したデータの中身を見ると大量の .dcm ファイルが有ることが分かります。それぞれのファイルがスライスした画になるので、1 枚 1 枚を Texture2D となるようインポートし、これを後で結合して Texture3D になるようにします。

teratail で回答者さんが便利な Editor スクリプトを書いていらっしゃいましたが、DICOM 以外のユースケースでもわかりやすいように最小限のコードで変換をしてみたいと思います。なので、.dcm を Texture2D としてインポートするスクリプトおよび、複数の Texture2D をまとめて Texture3D にするスクリプトの 2 つを別々に最小限のコードで書いてみます。

.dcm ファイルのインポート

まずは、.dcm を Texture2D としてインポートしてみます。

using UnityEngine;
using System;

[UnityEditor.AssetImporters.ScriptedImporter(1, "dcm")]
public class DicomImporter: UnityEditor.AssetImporters.ScriptedImporter
{
    public override void OnImportAsset(UnityEditor.AssetImporters.AssetImportContext ctx)
    {
        try
        {
            var image = new Dicom.Imaging.DicomImage(ctx.assetPath);
            var tex2d = image.RenderImage().As<Texture2D>();
            ctx.AddObjectToAsset("Volume", tex2d);
        }
        catch (Exception e)
        {
            Debug.LogException(e);
        }
    }
}

fo-dicom のドキュメントを見ると DicomImage を作成し、これを RenderImage().As<Texture2D>() するだけで Unity の Texture2D に変換できるようになっています。これを AssetImportContext に登録することで .dcm ファイルが Texture2D として使えるようになります。

Texture3D へ結合

ウィンドウを作っても良いのですが、シンプルに MonoBehaviourスクリプトとそのエディタ拡張の組み合わせで片付けます。まず、次のようなコンポーネントを作ります。

using UnityEngine;
using System.Collections.Generic;

public class Texture2DArrayToTexture3DConverter : MonoBehaviour
{
    public List<Texture2D> texture2DArray;
}

そして、このエディタ拡張でインスペクタから登録された Texture2D の配列を結合して Texture3D にします。

using UnityEngine;
using UnityEditor;

[CustomEditor(typeof(Texture2DArrayToTexture3DConverter))]
public class Texture2DArrayToTexture3DConverterEditor : Editor
{
    string error;

    public override void OnInspectorGUI()
    {
        base.OnInspectorGUI();

        if (GUILayout.Button("Convert"))
        {
            error = "";
            Convert();
        }

        if (!string.IsNullOrEmpty(error))
        {
            EditorGUILayout.HelpBox(error, MessageType.Error);
        }
    }

    bool Convert()
    {
        var converter = target as Texture2DArrayToTexture3DConverter;
        var tex2dArray = converter.texture2DArray;
        if (tex2dArray.Count == 0)
        {
            error = "no image";
            return false;
        }
        tex2dArray.Reverse();

        var w = tex2dArray[0].width;
        var h = tex2dArray[0].height;
        var d = tex2dArray.Count;
        var format = tex2dArray[0].format;
        var colors = new Color32[w * h * d];

        for (int i = 0; i < d; ++i)
        {
            var tex2d = tex2dArray[i];
            if (tex2d.width != w || tex2d.height != h)
            {
                error = "texture size error";
                return false;
            }
            if (tex2d.format != format)
            {
                error = "texture format error";
                return false;
            }
            tex2d.GetPixels32().CopyTo(colors, w * h * i);
        }

        var tex3d = new Texture3D(w, h, d, format, false);
        tex3d.SetPixels32(colors);
        tex3d.Apply();

        var path = EditorUtility.SaveFilePanelInProject(
            "Texture3D",
            $"New Texture3D Asset",
            "asset",
            "Save Texture3D");
        AssetDatabase.CreateAsset(tex3d, path);
        AssetDatabase.SaveAssets();

        return true;
    }
}

こうして次のように適当なシーンを作成、DICOM ファイルを先ほど作成したコンポーネントにドラッグして登録し、エディタ拡張で作成した Convert ボタンを押下します。

f:id:hecomi:20210926230117p:plain

これで、Texture3D が生成されます。

レンダリング

レンダリングは以下の記事のものを使います。

tips.hecomi.com

少しだけ改変しました。おさらいのためコードを貼っておきます。

Shader "VolumeRendering/Basic"
{

Properties
{
    [Header(Rendering)]
    _Volume("Volume", 3D) = "" {}
    _Transfer("Transfer", 2D) = "" {}
    _Iteration("Iteration", Int) = 10
    _Intensity("Intensity", Range(0.0, 1.0)) = 0.1
    [Enum(UnityEngine.Rendering.BlendMode)] _BlendSrc ("Blend Src", Float) = 5
    [Enum(UnityEngine.Rendering.BlendMode)] _BlendDst ("Blend Dst", Float) = 10

    [Header(Ranges)]
    _MinX("MinX", Range(0, 1)) = 0.0
    _MaxX("MaxX", Range(0, 1)) = 1.0
    _MinY("MinY", Range(0, 1)) = 0.0
    _MaxY("MaxY", Range(0, 1)) = 1.0
    _MinZ("MinZ", Range(0, 1)) = 0.0
    _MaxZ("MaxZ", Range(0, 1)) = 1.0
}

CGINCLUDE

#include "UnityCG.cginc"

struct appdata
{
    float4 vertex : POSITION;
};

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

sampler3D _Volume;
sampler2D _Transfer;
int _Iteration;
float _Intensity;
float _MinX, _MaxX, _MinY, _MaxY, _MinZ, _MaxZ;

struct Ray
{
    float3 from;
    float3 dir;
    float tmax;
};

void intersection(inout Ray ray)
{
    float3 invDir = 1.0 / ray.dir;
    float3 t1 = (-0.5 - ray.from) * invDir;
    float3 t2 = (+0.5 - ray.from) * invDir;
    float3 tmax3 = max(t1, t2);
    float2 tmax2 = min(tmax3.xx, tmax3.yz);
    ray.tmax = min(tmax2.x, tmax2.y);
}

inline float sampleVolume(float3 pos)
{
    float x = step(pos.x, _MaxX) * step(_MinX, pos.x);
    float y = step(pos.y, _MaxY) * step(_MinY, pos.y);
    float z = step(pos.z, _MaxZ) * step(_MinZ, pos.z);
    return tex3D(_Volume, pos).r * (x * y * z);
}

inline float4 transferFunction(float t)
{
    return tex2D(_Transfer, float2(t, 0));
}

v2f vert(appdata v)
{
    v2f o;
    o.vertex = UnityObjectToClipPos(v.vertex);
    o.localPos = v.vertex;
    o.worldPos = mul(unity_ObjectToWorld, v.vertex);
    return o;
}

float4 frag(v2f i) : SV_Target
{
    float3 worldDir = i.worldPos - _WorldSpaceCameraPos;
    float3 localDir = normalize(mul(unity_WorldToObject, worldDir));

    Ray ray;
    ray.from = i.localPos;
    ray.dir = localDir;
    intersection(ray);

    int n = _Iteration * ray.tmax / sqrt(3);
    float3 localStep = localDir * ray.tmax / n;
    float3 localPos = i.localPos;
    float4 output = 0;

    [loop]
    for (int i = 0; i < n; ++i)
    {
        float volume = sampleVolume(localPos + 0.5);
        float4 color = transferFunction(volume) * volume * _Intensity;
        output += (1.0 - output.a) * color;
        localPos += localStep;
    }

    return output;
}

ENDCG

SubShader
{

Tags 
{ 
    "Queue" = "Transparent"
    "RenderType" = "Transparent" 
}

Pass
{
    Cull Back
    ZWrite Off
    ZTest LEqual
    Blend [_BlendSrc] [_BlendDst]
    Lighting Off

    CGPROGRAM
    #pragma vertex vert
    #pragma fragment frag
    ENDCG
}

}

}

これで Transfer Function で色を設定すると狙った密度の場所に指定した色が付く形になります。詳しくはデモプロジェクトをダウンロードしてご参照ください。

おわりに

密度を表す Texture3D が作れさえすれば形式はどんなものでも同じ用に描画出来ます。そのうちマーチングキューブによるメッシュ化の方も試してみたいです。