はじめに
先日、Twitter で手持ちの DICOM 画像を以前書いたボリュームレンダリングの記事のように描画したい、というご相談をいただき試してみましたので内容を共有します。
DICOM(ダイコム)とは「Digital Imaging and COmmunications in Medicine」の略で、CT や MRI などのデータを保存する画像や通信プロトコルといった医用の標準規格です。今回の場合は、MRI のスライス画像が DICOM という形式のもと保存されているといった状況で、それをどうレンダリングすればよいのか、というお話でした。
デモ
ダウンロード
情報収集
Unity で DICOM を扱うことに関してはインターネット上にいくつか情報がありました。
Unite Tokyo 2018
カンファレンスイベントである Unite Tokyo 2018 では、DICOM に関連した講演を東京大学医学部脳神経外科の金太一助教がされていました。
この中で、Kompath 社と協力して作られた DICOM を読み込める Simple DICOM Loader や、マーチングキューブ法でメッシュ化する High Speed CPU-based Marching cubes が紹介されていました。
fo-dicom の利用
また、Q&A サイトの teratail では正に似たような質問がされていました。
こちらでは回答者さんが fo-dicom という無料の DICOM ローダーを使って Texture3D
化をされていました。
方針
上記 fo-dicom を使って同様に Texture3D
化を行う方針で進めることにしました。
私の記事では PVM という形式で配布されているデータを RAW 形式に変換、更にそれをパースするみたいな部分に記事の多くの分量が割かれてしまっていましたが、本質はいくつかの工程を経てデータを Texture3D
化したあとの、それをシェーダでどのようにレンダリングするかというところにあるので、大きな方針としては変わりません。
つまり、DICOM データを fo-dicom を利用して Texture3D に変換するスクリプトを書き、その Texture3D
を入力として以前の記事と同じコードでレンダリングを行います。
データのダウンロード
ご相談頂いた方のデータは諸事情でブログでは使えないため、本記事では配布されているデータで試してみることにします。上記質問回答の文中にもあるページから「Visible Female CT Datasets > Regional Tar Files (download) > Head」をダウンロード・展開しておきます。
fo-dicom のインポート
fo-dicom 本体は GitHub にてオープンソースとして配布されています。
Unity 向けのパッケージは Asset Store にて配布されています。
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 ボタンを押下します。
これで、Texture3D
が生成されます。
レンダリング
レンダリングは以下の記事のものを使います。
少しだけ改変しました。おさらいのためコードを貼っておきます。
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
が作れさえすれば形式はどんなものでも同じ用に描画出来ます。そのうちマーチングキューブによるメッシュ化の方も試してみたいです。