Re - ImageJで学ぶ！

2022年8月12日金曜日

第72回　ImageJでRadiomicsを体験しよう！

今回はレディオミクスに挑戦する内容です。

2022/4に、ImageJでRadiomic特徴が計算できる「RadiomicsJ」というJavaのライブラリを公開しました。 RadiomicsJはImageJやFijiのプラグインインターフェースがインプリメントされており、ImageJ/Fijiプラグインとして動作します。今回はこのRadiomicsJを使ったRadiomicsをご紹介いたします。

Radiomicsとは

レディオミクスは、「高信号/低信号」「不明瞭」「辺縁に造影効果」「均一/不均一」などの医用画像から得られる特徴が、分子生物学的特徴のパターンを捉えることができるという仮説に基づき、腫瘍遺伝子型や予後などの予測やイメージングバイオマーカの発見などを目的として行われる研究手法です。

類似した研究領域として、従来から、コンピュータ支援診断（computer-aided diagnosis, CAD）がありますが、レディオミクスはCADの中の一つというよりも、独立した研究として説明されるように思います。筆者の個人的な、CADとレディオミクスとの違いについての認識は次の通りです。CADは、良悪性鑑別予測、異常検知、予後予測などを主な目的として、医師の意思決定支援のために研究されてきました。これに対してレディオミクスは、目的変数をこれらに限定せず、分子生物学的特徴という概念に置き換えた上で、マクロな画像の世界と細胞やＤＮＡなどのミクロな世界の関係の探索を目的として研究されています。

レディオミクスはだれでも試行できます。例えば、TCIA（The Cancer Imaging data Archive）などからオープンデータセットを探して、脳MRI画像を使って画像特徴を計算し、この計算結果から脳腫瘍のWHOグレードを予測したり、腫瘍遺伝子型の異常の有無を予測したりなどです。 PubMedなどで検索をすれば、他にもさまざまな応用例がヒットします。

しかし、いざやってみようとなると、なにから行えばよいかわからなくなることがあります。

本ブログでは、「胸部レントゲン画像の肺野・縦隔から60歳以上かどうかを予測する」ことに挑戦してみましたので、手順を公開いたします。

環境

Windows（Mac、Linuxでも動作するはずです）

必要なもの

ImageJ

ImageJサイト(https://imagej.nih.gov/ij/download.html)の説明に従ってダウンロードします。

RadiomicsJ

こちらからダウンロードします。
https://drive.google.com/file/d/148ToZ1knJOpi4pzIwSf8g23kn-d006tR/view?usp=sharing
ダウンロードしてzip解凍後します。

図　RadiomicsJを解凍して内容を展開した様子

「RadiomicsLib」フォルダと「My_RadiomicsJPlugin.class」ファイルを、インストール済みのImageJのpluginsフォルダにコピーします。

図　ImageJのpluginsフォルダへコピー（これでインストール完了）

JDK 1.8

OpenJDK1.8の最新バージョンをインストールして、環境変数にJAVA_HOMEパスを追加してください。 (JDK 1.8以外ではエラーが起こる可能性があります。)
参考：JDKのインストール方法（Windows）https://www.javadrive.jp/start/install/index4.html

データセット

こちらからダウンロードします。
https://drive.google.com/file/d/1Tet5rSuD1Zw1tbyIicGGU3dyYIO5OIht/view?usp=sharing
画像データ、Ageラベルデータ、解析処理のsettingsファイルが含まれています。
データセットの元は、JSRTから教育研究用に公開されているSegmentation02データセットおよびAgeデータセットです。ダウンロード用のデータセットは、元のデータから欠損ラベル除外、ファイル名統一を行ったデータになっています。

Weka

こちら(https://waikato.github.io/weka-wiki/downloading_weka/)からOSに合わせてダウンロードします。

レディオミクスのプロセスの全体像

レディオミクスのプロセスは、大きく二つのステップで構成されます。

画像特徴を計算する
機械学習で予測モデル作成する

①は、解析したい画像上の関心領域を決めて、画像特徴を計算します。

図　胸部レントゲン画像（case001）

今回の胸部レントゲンの例では、解析したい画像上の関心領域はマスク画像（ラベル画像とも呼びます）で指定されます。今回の例では、肺野領域:255、心臓:85、肺野外:170、体外:0の画素値になっています。肺野は、マスク画像上で255の値の白い領域です。この肺野から画像の特徴を計算したいときは、この領域を計算するように設定すればよいということになります。心臓の縦隔領域を解析したい場合は、85の領域（濃いグレーの部分）を指定すればいいということです。

②の機械学習は、一般的な機械学習の手順で、説明変数と目的変数を設定して、予測モデルを作成・テストします。

レディオミクスで一番大変なのは、このようなデータセットの準備ですが、ImageJを使って比較的簡便に作成できます。

ここで、マスク画像の作成方法を紹介します。（この操作は以降の機械学習モデル作成には不要です）

作成時のポイントは、ROIマネージャツールを使うことです。 ROIマネージャはAnalyze>Tools>ROI Managerから起動します。ROIマネージャにROIを追加するには、好きなROIツールで画像にROIを描いた後、「Add」します。

図　ROIを描いた様子

（描かれたROIはROIマネージャにAddされている）

ROIを更新したとき（形や位置を変えたとき）は、「Update」ボタンを必ず押すようにします。

ROIが設定出来たら、ROI以外の領域を０、ROIの内側を任意の値にします。 RadiomicsJは8-bitグレースケールのマスク画像（Tifがデフォルト）に対応していますので、ROIの内側の値は、0～255である必要があります。例えば、肺の領域を255、そのほかを0にしたい場合は、ImageJメニューのEdit>Clear OutsideでROI以外の領域を0にできます。

図　Edit>Clear Outside

Clear Outsideをしてから、Process＞Math＞Set...で、値を255にして実行します。

図　Process＞Math＞Set...（値を255にした例）

この操作を繰り返してマスク画像ができます。

図　出来上がったマスク画像

完成したら、解析画像とはフォルダを分けて、Tif形式で保存しておきます。保存するときの画像ファイル名は、解析画像と同じ画像名などにしておくとわかりやすいです。RadiomicsJはTif形式が標準画像フォーマットです（ImageJで読み込めるファイルならばOKです）。保存の際、元の解析画像を上書きしないように注意します。 ROIデータの結合や分離もROI Manager上で行うことができます。ROIデータは保存しておくことができるので、やり直しも可能です。

Radiomic特徴を計算する

RadiomicsJによるRadiomic特徴の計算には、３つのデータが必要です。

画像データ
マスクデータ
セッティングファイル

RadiomicsJはグラフィカルユーザーインターフェースからの操作や、コマンドプロンプトから一括計算処理することもできます。セッティングファイルはコマンドプロンプトによる操作の時に必要になります。順に解説します。

ImageJのプラグインとして利用する方法（一例ずつ処理）

画像を一枚一枚確かめながら解析を行いたい場合は、解析画像とマスク画像をIJで表示した状態で、ImageJのプラグインを起動し、画像を設定して、計算を実行します。

図　Plugins>My_RadiomicsJPlugin (解析画像とマスク画像はプラグインを起動する前にIJに表示しておく)

図　特徴量の計算結果がResultTableに表示される（少し時間かかります）

追加で解析をしていく場合も同様の操作を繰り返します。このとき、ResultTableを閉じずにそのままにしておくと、行が追加されていきます。

図　連続して解析していく例（画像は表示しなおしてから再度プラグインを起動、設定（図中矢印）を適用する）

図　ResultTableに新しい解析結果が追加される

結果を保存したいときは、このResult Tableの機能（File>Save）を使って保存できます。

しかし、一般に、大量のデータを自動的に繰り返し処理したいことの方が多いと思います。そのような場合は、次に説明するコマンドプロンプトからの操作が可能です。

コマンドプロンプトから利用する方法（自動で全処理）

コマンドプロンプトから利用する場合、データセットをフォルダに綺麗に分けておく必要があります。綺麗にというのは、画像とマスクがセットになるように、症例ごとに分けておくという意味です。この際、画像とマスクのファイル枚数は同じになる必要があります。

今回の画像データセットはあらかじめorgとlabelにフォルダが分かれており、対応する画像ファイルがそれぞれに準備されています。ファイル数も同じになっています。

（マトリクスサイズが異なる場合は先に整えておくか、症例ごとにフォルダを分けておく必要があります。また、CTやMRIのようにシリーズが複数の画像で構成される場合もフォルダで分けておく必要があります。）

あとは、計算に必要なセッティングをします。

計算の実行には、必要に応じて、セッティング情報（.propertiesファイル）を指定できます。 RadiomicsJは3Dがデフォルトです。そのため、今回の胸部レントゲンのような2D画像を解析するときは2Dで解析するように設定します。

CTやMRIのシリーズ画像を3Dで処理する場合は、セッティングファイルは省略可能です。

今回のデータセットはすべての画像の縦横のマトリクスサイズが同じですので、一連として一括処理します。

セッティングファイルの内容

セッティングファイルには、RadiomicsJで計算可能な特徴ファミリーの計算設定をセットします。よく使うのは、マスクのラベル値の指定（INT_label）、２D解析とするかどうか（BOOL_force2D）、2Dの画像特徴を計算するか（BOOL_enableShape2D）、デフォルトでない特徴を有効にするか（BOOL_activate_no_default_features）です。

有効無効の設定は、1（True）、0（False）で指定します。

今回は、2Dでの処理に合わせるために、BOOL_force2D=1, BOOL_enableShape2D=1としています。

セッティングファイル内で、「!」「#」が行の先頭についている行はコメントアウト（その行を無視する）を意味しています。

他の項目は、特徴計算の専門家が設定するような項目なので、最初は気にしなくてもよいです（例えば、32-bitのPET画像を解析するときはBinWidthで実行したほうがよいなど、モダリティによっても考慮すべきことがあります。詳しくはIBSIリファレンスを参照ください）。

図　settings.propertiesの例

特徴計算を実行していきます。

コマンドプロンプトを開き、カレントディレクトリをRadiomicsJLibフォルダにします。

下のようにコマンドを入力し、Enterで実行します。

>cd "RadiomicsJLibフォルダまでのパス"

図　カレントディレクトリの移動（ダウンロードフォルダにあるRadiomicsJLibフォルダへ移動している例）

次に、処理を実行するためのコマンドを入力し、Enterで実行します。

>java -jar RadiomicsJ.jar -i "orgフォルダまでの絶対パス" -m "labelフォルダまでの絶対パス" -s "settings.propertiesまでの絶対パス" -d

図　デスクトップにあるデータセットフォルダ内の各ファイル（org, label, settings.properties）を指定した例

オプションの説明

-i：解析画像フォルダパス
-m：マスク画像フォルダパス
-s：セッティングファイルパス
-o：解析結果CSV保存先（.csvまで必要）
-d：処理経過の表示

パスというのは、PC内のフォルダやファイルの場所のことです。

解析結果CSVはデフォルトでは自動で保存されません。自動的に保存したい場合は、-oオプションでcsvファイルの保存場所を指定しておきます。

処理が始まると、計算中のプロセスが表示されます。

245件の処理が終わるまで継続されるので、小一時間かかります。

図　処理中画面

処理が終わると、Rediomics Features結果テーブルが表示されます。この結果テーブルのFileメニューから結果をCSV保存します。

図　Radimomics Features 結果テーブル例（Fileメニューから結果をCSVしておく）

保存したCSVをエクセルで開きます。

図　処理結果例（今回はDICOMデータではないので、基本情報はNaNになっている）

これで、画像特徴の計算は完了です。

機械学習モデルを作成する

先に、学習データを完成させます。

今回は年齢を推定するのですが、データも少ないので、分類モデルとして扱う例で示します。60歳以上かどうかを予測するためのモデルを作成してみます。

Radiomic特徴の計算結果CSVをエクセルなどで開き、データセットに含まれているage.csvファイル内にある目的変数「elder_class」を、計算結果CSVの一番左のカラムに追加します。 Ageの情報はファイル名順に並んでおり、解析結果もファイル名順に並んでいますので、そのまま列を挿入して追加します。

次に、不要なカラムを削除しておきます。RadiomicsJはデータの基本的な情報なども取得してテーブルにまとめてくれますが、これらは学習には不要なデータですので、先に削除しておきます。

ここでは、OperationalInfo_ではじまるカラム、Diagnostics_ではじまるカラムはすべて削除します。

この他にも、①相関の強い、または完全相関する変数は削除、②すべてのデータで同じ数になっている（例えばすべて０であるなど）変数は削除、などしてデータを調整しますが、ここでは、この手順を省略します。

図　elder_class（目的変数）を学習データに追加し、OperationalInfo_ではじまるカラム、Diagnostics_ではじまるカラムはすべて削除

ここまでできたら、CSVを別名で保存しておきます（今回はRadiomics Features lung ML-2.csvとして保存しました）。

これで、練習用の学習データの準備が整いました。 Wekaを起動しましょう。

図　WEKA

「Explorer」を起動します。Preprocessタブ内の「Open file...」から学習データCSVを読み込みます。 CSVを読み込むために、「Files of Type:」で読み込むファイルの形式をCSVに変更します。

図　CSVのロード

図　データを開いた結果

データのサマリーは「Visualize」タブから確認できます。

図　データサマリー

次に、「Preprocess」タブで、前処理を実行していきます。今回は、定数のみを含むの説明変数の除外のために「RemoveUseless」とデータの標準化「Standardization」を適用してみます。（外れ値の除外も可能ですが、説明が煩雑になるため省略します。方法はこちら：https://youtu.be/WrjpO7CmUoQ）

「RemoveUseless」は、PreprocessタブのFilterボタンから、unsupervised>attributes>RemoveUselessを選択し、右側のApplyボタンを押して実行します。

この操作で、すべて「0」などの値となっていた変数がドロップされます。

図　RemoveUseless

次に、同様の操作で、標準化を実施します。unsupervised>attributes>Standardizeで選択できます。同様にApplyして実行します。文字列のカテゴリ変数でない数値の変数はすべて標準化されます。

図　標準化（Standardize）

最後に、学習方法を設定して学習を実行します。利用する分類器は今回はランダムフォレストにしました。

学習の設定

Classifier: RandomForest
Cross-validation: 10 folds-CV
Target class: elder_class (間違えやすいので注意)
Startで実行

図　学習の設定

分類精度は学習終了時に表示されます。

AUCは0.679でした。

Startボタン下にある「Result list」に解析結果が保存されます。このリストを右クリックして、「Visualize thresholdCurve」(TRUEの方)を実行すると、ROC曲線を確認できます。

図　ROC-AUC (TRUEクラス)

今回は、分類精度が低くなってしまったのですが、これがAUC 0.9などの予測モデルであれば、結構いい感じの予測モデルになりそうということがわかります。

訓練した予測モデルは、modelファイルとして保存(Result listを右クリックしてSave model)しておき、Weka APIを通じてアプリケーションから利用できます。

Knowledge Flowを使う

Wekaには、複雑な機械学習プロセスを可視化しながら実行できる機能があります。

「Knowledge Flow」です。

この機能は、Explorerでできることをアイコンを使ってパイプラインを描くことができます。Explorerでは、途中でどこまでやったかわからなくなってしまったりすることもあるので、直感的にわかりやすいインターフェースは重宝されます。

今回の解析に特徴選択というプロセスを追加した例を示します。

（Wekaはいろいろな特徴選択が可能です：https://youtu.be/Pf9yKjSiVnw）

図　Knowledge Flowの構築例（訓練データとテストデータを分割し、訓練データで特徴選択を実行。訓練データのみの交差検証と、特徴選択済みのモデルをテストデータで評価した結果を同時に得る）

利用したい処理をツリーから選び、キャンバスに追加していきます。
右クリックで処理の流れに使うデータセットを選び、繋げていきます。
処理の結果や、マップは、Visualize機能を使います。
各アイコンには詳細設定を追加できます。
フローが完成したら、実行ボタンで実行できます。

図　訓練データによる交差検証結果（訓練データによって特徴選択を行っているため少し精度が高くなる, AUC 0.73）

図　テストデータによる精度（AUC 0.65）

縦隔から60歳以上かどうかを予測する

これまでの例では、肺野から取得されたラジオミクス特徴で予測モデルをテストしていました。同じ手順で縦隔（主に心臓）のラジオミクス特徴から60歳以上かどうかを予測することもできます。ラジオミクス計算時に、settings.propertiesの「INT_label」を85にして計算するだけです。

また、AgeラベルCSVには10歳刻みのカテゴリを予測するためのラベルも用意しています。マルチクラス分類は、筆者が試した限りでは難しそうでしたが、ご興味のある方は、ラベルを入れ替えて試してみることができます。

図　ラベルを年代に変更したレディオミクスCSV例

他の秀逸なレディオミクスのツール

RadiomicsJ以外にもGUIから操作できるツールがあります。

LifeX
PyRadiomics + 3D Slicer
etc...

レディオミクスの詳しい話について

IBSIリファレンス(https://theibsi.github.io/documentation/)に詳しい手順が解説されていますので、専門家になりたい方は一読をおすすめします。

おわりに

自分の研究なんて無効な研究だと思っても、無効の組み合わせが有効なものへと形を変えることがあります。医用画像を使った研究でなんかやろうかなと思っている方は、レディオミクスなんていかがでしょうか。

医師

治療や診断の高品質化に直結するような研究のために。研究の成果を積み上げるために。

コメディカル

専門の領域で、治療や診断とは区別したうえで、生物学的な探索や診断や治療を支援するための研究に。（医師でなくても、オープンデータを使ってさまざまな検討が可能な時代になりました。）
最近は、予測モデルをクラウドから実行できるようなサービスもありますので、自分で作った予測モデルを市場に出すことも夢ではなくなりました。

情報系エンジニア

医療の知識は医療系の人たちから吸収すれば、あとは得意な情報処理。医療応用のための基礎技術特許などを通じた社会貢献のために。

バイオ系研究者

人と自然との関連の探求に。

RadiomicsJの改良（特徴の追加など）、バグ報告など、ご意見・コメントお待ちしております。

Visionary Imaging Services, Inc. Tatsuaki Kobayashi

References

元データセット：http://imgcom.jsrt.or.jp/minijsrtdb/
IBSIリファレンス：https://theibsi.github.io/documentation/
RadiomicsJ（論文引用をしていただけますと幸いです）：https://pubmed.ncbi.nlm.nih.gov/35792994/

2021年12月17日金曜日

第71回 ImageJで深層学習モデルを動かす！－NVIDIA pre-trained modelの例－

深層学習はアカデミックからの視点で見ればその歴史はその分野内では長いほうで、以前から脳の機能を模倣したアルゴリズムとして研究されてきました。しかし、計算量が膨大になるため、ハードウェア的に計算コストがかかりすぎるなどのリミテーションがあり（あったので）、大衆化までの道のりが長くなってしまいました。しかし、TensorFlowやChainer（現在のPyTorchにも影響を与えた秀逸な日本製・日本発機械学習向けの計算機ライブラリ、現在は開発は終了している）などの深層学習に親和性の高いオープンソースの計算処理ライブラリや、オープンデータ（例えば、MNISTやImageNetなど）を世界で共有しようというデータサイエンス界の文化もこの分野が急速に注目される追い風になりました。

医療応用も進んでいます。今日では、「医用AI」「医療AI」「医療機器プログラム（あるいはすでに申請済みの医療機器プログラムへの機能追加をするなど）」などという枠として承認申請ができるようになり、診療報酬制度下での運用方法も厚労省を中心に議論されています。例えば、CTスキャナーでは、深層学習を使った超解像画像を作ったり、内視鏡では、内視鏡画像からリアルタイムに病変を検出し、悪性度などを予測する、など、臨床でのクリニカル・クエッションに応えようとするモデル（いわゆる、医療AI）が販売されています。

このような中で、今回はImageJで深層学習モデルを動かす方法について、NVIDIA社がNGCという取り組みの中で研究者とともに公開している深層学習モデル（肺のセグメンテーション）を例に、その概要と簡単な使い方を説明します。

NVIDIA NGC

NGCは「NVIDIA GPU CLOUD」の略です。このNGCの取り組みでは、GPUを普及させるために、GPUを利用するとメリットの大きい深層学習モデル開発などを対象として、必要なツールを提供しています。訓練済みの深層学習モデルデータや、モデルを動かすために必要なコンピュータプログラム群などです。

図 NVIDIA NGC Catalogのページ例

これらのツールはカタログにまとめられており、ユーザーが必要とするサービスを選べるようになっています。クラウド上で深層学習を動かせるように、環境が整えられています。

とはいうものの、クラウドサービスとなると、個人の研究やちょっとした勉強レベルでは手を出しにくいですし、試してみようと思いづらいところもあるのではないでしょうか（難しそうだから）。また慣れが必要と感じてしまうのは私だけでしょうか。普段からクラウドが当然で、「あ、PCの電源つけよう！クラウド上の！」とすいすいできる人はクラウドサービスエンジニアくらいでしょう。私のような凡人には遠い感覚です。ふだんの使い慣れたパソコンで、なんならImageJで動かしたい（できるなら無料で）と思ってしまいます。

今回は、このような私が無料・深層学習・ImageJにトライしてみました。

以下、手順になります。

pre-trained modelのダウンロード

訓練済みモデル（pre-trained model）は、NVIDIA NGC > Catalog > modelsのページからダウンロードできます。キーワード検索できます。例として、「lung」で検索すると、このようなモデルのリストが表示されました。

図　「lung」検索結果

今回はこの「clara_pt_covid19_ct_lung_segmentation」を使ってみたいと思います。このモデルは、「A pre-trained model for volumetric (3D) segmentation of the lung from CT images.」ということで、3Dのボリュームデータから、肺の領域を予測してくれるようです。

では、このモデルをダウンロードし、解凍します。

ダウンロード時は、「files.zip」というファイル名になっていました。

解凍すると、中身はこのようになっています。

図　files.zip

モデルデータは、解凍したフォルダー内のmodelsフォルダーに入っています。

この例では、「model.pt」「model.ts」です（拡張子が見えないときは、右クリックからプロパティを開いて確認できます）。

図　modelファイル（files>modelsフォルダの中身）

さて、このptファイルは何のファイルでしょうか。私はTensorflow派なので最初は何のファイルだか分かりませんでした。このptファイルは、PyTorchで作成された深層学習モデル用のファイルです。tsファイルはトーチスクリプトファイルと呼ばれるファイルで、ptファイルと同様に、モデルデータを保存しているデータです。ptファイルとの違いは、モデルを推論実行させる際に高度な処理の設定を追加できることです。この「.pt」拡張子のファイルは、モデルの学習状態とモデルの構造の両方を保持できます。ただし、モデルの学習状態のみを保持することもできます。モデルの学習状態のみを保持している場合は、モデルの構造は保持しません。これはクリエイター事情です。なので、ptとtsのファイルが合わせてある場合は、モデルとして利用する際に両方のファイルが必要であることが多いです。今回のケースは2つとも必要です。

CPUかGPUか

深層学習モデルは、（特にPyTorchの場合）GPU用あるいはCPU用に出力することができます。NVIDIA NGCで公開されているモデルはすべてGPU用になっています。そのため、なにもしなければ、GPUを持っていない人はここから先の操作ができません。しかし、それでは「大衆化」とは呼べません（個人的に）。ここでは、CPUでも動かせるようにモデルに前処理を追加します。GPU前提の方は、ここの処理は不要です。CPUで動かしたい方は、この処理が必要です。

GoogleColabを使って、GPU用に保存されたモデル（先ほどのpt, tsファイル）を、CPU用にリトレースします。新しくColabノートブックを作成し、ランタイムをGPUに切り替えます。方法は以下の通りです。

      import torch
      # 自分のモデルが保存されている場所を指定します。
      # Colabノートブックに一時的に直接アップロードした場合は、"model.pt"や"model.ts"としてください
      p2state_dict = "/content/drive/MyDrive/NVIDIA_PreTrained_Models/clara_pt_covid19_ct_lung_segmentation/models/model.pt"
      p2m = "/content/drive/MyDrive/NVIDIA_PreTrained_Models/clara_pt_covid19_ct_lung_segmentation/models/model.ts"
      state_dict = torch.load(p2state_dict) # 学習重み
      model = torch.load(p2m) # モデルの構造など 
      model.load_state_dict(state_dict) # All keys matched successfully, モデルに学習状態をロード
      model = model.to("cpu")# gpuにトレースするときはmodel = model.to("cuda")
      traced_model = torch.jit.trace(model, (None,224,224,32))# traced_model = torch.jit.script(model) # こちらでもよさそう（試していません）
      torch.jit.save(traced_model, "model_cpu.pt")
      # model_cpu.ptというファイルがColab内のフォルダに作成されるのでダウンロードして使います。

CPU用にリトレースしたモデルデータは、こちらに公開しています。

ImageJ/Fijiで動かしてみる（プラグインを開発する）

ここから、Javaの強みである「Write once, run anywhere」を享受させてもらいます。

今回は、プラグインとしてファイルを作るのではなく、ImageJをライブラリとして使って、簡単に動作を確かめる方法で解説させていただきます。

最終的に動かすところまではできるので、プラグインとしてJar化して使いたいなどは、ImageJ公式のプラグイン開発方法をご参照ください。

深層学習モデルを推論実行させるために必要なライブラリとして、今回はAmazonが開発している「Deep Java Library」を使います。

開発環境には、Eclipseを使っていきます。

Javaは、私の環境では、AdoptOpenJDK11を使っています。

では、Mavenプロジェクトを新しく作成し、pom.xmlを次のように編集しています。

groupid、artifactid、versionなどは自分の設定に合わせてください。

一部抜粋していますが、<project></project>で挟まれる領域です。

<groupid>com.vis.machinelearning</groupid>

<version>0.0.1-SNAPSHOT</version>

</properties>

<groupid>net.imagej</groupid>

</dependency>

</dependency>

<groupid>ai.djl.pytorch</groupid>

<artifactid>pytorch-engine</artifactid>

</dependency>

<groupid>ai.djl.pytorch</groupid>

<artifactid>pytorch-native-auto</artifactid>

<scope>runtime</scope>

</dependency>

<groupid>ai.djl.tensorflow</groupid>

<artifactid>tensorflow-engine</artifactid>

</dependency>

<groupid>ai.djl.tensorflow</groupid>

<artifactid>tensorflow-native-auto</artifactid>

<scope>runtime</scope>

</dependency>

<artifactid>model-zoo</artifactid>

</dependency>

<groupid>ai.djl.pytorch</groupid>

<artifactid>pytorch-model-zoo</artifactid>

</dependency>

</dependencies>

</project>

今回は、Windowsで動かしていますが、MacでもLinuxでも動くFatなライブラリをインストールしています（そのはずです）。必要なものは、ImageJ、DJL、DJL用のPyTorchのラッパーです。Tesorflowのラッパーも加えればTensorflowモデルを動かすことができます。上記のpom.xmlにはこれらが含まれています。

うまくライブラリを参照できると、エラーのない状態になります。

図　うまくライブラリをMavenで参照できた状態（少し時間がかかります）

サンプルデータ

ここでは、サンプルCTシリーズ画像として、TCIAで公開されているCOVID-19データで試していきます。

モデルのdocs>READMEを見てみると、下記のような記載があります。

## Input

Input: 1 channel CT image with intensity in HU and arbitary spacing

1. Resampling spacing to (0.8, 0.8, 5) mm

2. Clipping intensity to [-1500, 500] HU

3. Converting to channel first

4. Randomly cropping the volume to a fixed size (224, 224, 32)

5. Randomly applying spatial flipping

6. Randomly applying spatial rotation

7. Randomly shifting intensity of the volume

## Output

Output: 2 channels

- Label 0: everything else

- Label 1: lung

入力は0.8*0.8*5.0ボクセルのCTシリーズ画像を前提としているようです。

入力画像をこのボクセルサイズになるようにリサンプリングしておきます。

ただし、後から肺の領域をクロップし、224*224*32、つまり、マトリクスサイズが224*224の32スライスのボリュームに整形する（4の手順）ので、ここはそこまで厳密でなくともよいかと思います(結局、クロップするときのサイズで拡大率が変わります。ただ、何もしないよりはばらつきが少なくなり、ある程度の整合性が担保されます。ここは開発者に詳細を伺わなければ解釈が難しいところです)。

リサンプリング済みのデータはこちらです。

オリジナルはこちらです（非圧縮済み）。

オリジナルからインプットデータを作成したい方のために、Fijiを使って作成する手順を簡単に示します。

例えば、ここで利用するオリジナルのサンプルCTシリーズは、Voxel size: 0.7031x0.7031x1.25 mm^3です（ImageJで開く>Image>Show Info...から確認）。
XY方向のリサンプリングは、縦横方向のマトリクスサイズを、0.7031/0.8倍して補正します。オリジナルの512*512マトリクスを450*450にします。
Z方向のリサンプリングは、簡易的ですが、Image>Stacks>ResliceZというFijiのプラグイン機能を利用します。ImageJではデフォルトで入っていないので、ここはFijiを使ったほうが楽です。Z軸方向に5.0mmにしたいので、5.0としてリスライスすればOKです。
次に、大まかに肺領域を32スライス選びます(Image>Duplicateからスライス選択)。
最後に、矩形のROIを肺を囲うように1つ設定し、ROIが設定された画像をTIFとして1つのファイルに保存します。ほかのフォーマットはお勧めしません。
これで入力画像の準備ができました。

図　サンプルデータ

（ボクセルサイズ0.8*0.8*5.0, 32スライスに限定）

モデルのロードから推論実行まで

ここではメソッド全体を記載しています。

各コードの解説はコード内のコメントに記載しています。

  
  public static void main(String[] args) throws IOException, ModelException, TranslateException {
	   
		/*
         * load model
         * load image
         * create input
         * predict
         * save results
		 */
        //リソースにモデルファイルを置き、モデルファイルまでのパスで設定している例
		String model_loc = "clara_pt_covid19_ct_lung_segmentation/models/model_cpu.pt";
		URI modelUri = null;
		try {
			modelUri = HelloPyTorch.class.getClassLoader().getResource(model_loc).toURI();
		} catch (URISyntaxException e) {
			e.printStackTrace();
		}
		Path modelPath = Paths.get(modelUri);
		
		//224*224*32 in HU unit.
		ImagePlus in = new ImagePlus("C:\\Users\\ユーザー\\Desktop\\Resampled.tif");//入力画像ファイルへのパス
		ij.gui.Roi rect = in.getRoi();//画像からROIを取得
		System.out.println("roi size:"+rect.getFloatWidth()+" "+rect.getFloatHeight());
		float[] vol = new float[224*224*32];
		NDList input = new NDList();
		NDManager manager = NDManager.newBaseManager();
		int vox_pos = 0;
		for(int i=0;i<in.getNSlices();i++) {
			in.setPosition(i+1);
			in.setRoi(rect);
			Calibration cal = in.getCalibration().copy();
			ImagePlus sliceCrop = in.crop();
			sliceCrop = sliceCrop.resize(224,224,"bilinear");
			ImageProcessor crop = sliceCrop.getProcessor();
			crop.setCalibrationTable(cal.getCTable());
			crop.setMinAndMax(-1500d+32768, 500d+32768);//ct 16-bit on IJ.
			for(int r =0;r<224;r++) {
				for(int h =0;h<224;h++) {
					float val = crop.getPixelValue(h, r);
					vol[vox_pos++] = val;
					if(i==0 && r < 10 && h < 10) {
						// System.out.println(val);//test print pixel val in HU.
					}
				}
			}
		}
		Shape s = new Shape(new long[] {1l,1l,32l,224l,224l});//batch,channel,d,h,w
		NDArray ndar = manager.create(vol, s);//create input to predict.
//		System.out.println(ndar.getShape());
		input.add(ndar);
		
        //入力から出力までのクライテリアを定義する入力も出力もNDListになる。
		Criteria<NDList, NDList> criteria =
                Criteria.builder()
                        .setTypes(NDList.class, NDList.class)
                        .optModelPath(modelPath)
                        .optTranslator(new NoopTranslator())
                        .build();
		//ZooModelクラスを利用して、訓練済みモデルをモデル化し、予測実行する
        try (ZooModel<NDList, NDList> model = ModelZoo.loadModel(criteria);
            Predictor<NDList, NDList> predictor = model.newPredictor()) {
                NDList result =  predictor.predict(input);
                NDArray res = result.get(0);
                NDArray bg = res.get(new NDIndex("0,0,:,:,:"));
                NDArray lung = res.get(new NDIndex("0,1,:,:,:"));
                System.out.println(res.getShape());//(1, 2, 32, 224, 224)
                System.out.println("min max : "+res.min()+", "+res.max());
                System.out.println(bg.getShape());//(32, 224, 224)
                System.out.println(lung.getShape());//(32, 224, 224)
                /*
                 * - Label 0: everything else
                 * - Label 1: lung
                 */
                //予測結果を画像へ
                float[] pred = res.toFloatArray();//1*2*32*224*224 = 3211264
                
                for(int c = 0;c<2;c++) {
                	ImageStack stack = new ImageStack(224, 224, 32);
                    for(int d=0;d<32;d++) {
                    	float[] pred_lbl = new float[224*224];//per slice
                    	int pred_lbl_pos = 0;
                    	for(int h=0;h<224;h++) {
                    		for(int w=0;w<224;w++) {
                    			if(c == 0) {
                					pred_lbl[pred_lbl_pos] = pred[(d*224*224)+pred_lbl_pos];
                				} else {
                					pred_lbl[pred_lbl_pos] = pred[c*(224*224*32)-1+((d)*224*224)+pred_lbl_pos];
                				}
                    			pred_lbl_pos++;
                            }
                        }
                    	FloatProcessor fp = new FloatProcessor(224, 224, pred_lbl);
                    	stack.setProcessor(fp, d+1);
                    }
                    ImagePlus imp = new ImagePlus("", stack);
                    imp.show();
                    IJ.saveAs(imp, "tif", System.getProperty("user.home") + "/Desktop/test_"+(c+1));
                }
        }
		
	}

実行すると、推論結果が得られます。結果の予測画像はデスクトップに保存されます。

図　推論結果（左：予測画像肺以外、右：肺）

図　予測画像の二値化結果

あとは、これを元のサイズに戻し、切り出した位置にもっていけば、どの領域が肺か、肺以外かを予想した結果を確認できます（この部分は省略します。ImageROIにしてオリジナル画像に重ねて表示するなどで実装できます）。

まとめ

今回は、ImageJを用いた深層学習モデルの利用例を紹介しました。Tesorflowモデルも動くので、いろいろなことができると思います。また、深層学習が必要ない場合は、WEKAという機械学習ライブラリを利用することもできます。
自分だけの・自分の施設だけの深層学習モデルをあなたの手で動かしてみることもできるのではないでしょうか。

References

NVIDIA NGC https://docs.nvidia.com/ngc/ngc-overview/index.html
https://xtech.nikkei.com/atcl/nxt/column/18/00001/03341/

Visionary Imaging Services, Inc.

2015年12月24日木曜日

第70回　ImageJを用いたBoneJの紹介で学ぶ！

“BoneJ”は、骨画像解析のためのImageJプラグインです。

骨の生成構造の解析の歴史は長く、昔から二次元スペクトル解析やフラクタル解析などの医用画像解析で骨粗鬆症の進行度の評価などに積極的に用いられてきました。
今回はImageJ のオープンソースコミュニティによって紹介されている、骨梁の形状と全体の骨の形状解析のためのプラグインBoneJについて、その概要と簡単な使い方を説明します。

BoneJとは

骨の三次元形状は、一般的にはCTとX線computed microtomography（μCT）の画像を用いて計測が行われます。これらCT、μCTなどの画像データを処理する画像処理や解析ソフトウェアは、装置が高額であることからソフトウェアも高額であることが多いです。
しかし、ImageJのBoneJを用いれば、高度な骨の構造解析が無料のプラグインで解析できます。

BoneJのインストール方法と主な機能

BoneJをインストールするには、あらかじめImageJの最新バージョンをダウンロードして、“BoneJ_.jar”をImageJフォルダの中のpluginsディレクトリにコピーします。
BoneJ_.jarは、（http://bonej.org/）からダウンロードできます。
また、BoneJはImageJ3DViewerが必要です。こちらがインストールされていない場合は合わせて取得しておきましょう。（ImageJ_3D_Viewer.jar）

BoneJには、骨梁の解析のために細線化解析（枝と分岐点の分類、カウント、計測）、非等方性解析、結合性解析（画像のオイラー特徴、1m㎥あたりの骨梁の解釈として濃淡結合性評価）フラクタル解析、平坦度解析（構造モデルインデックス）、等方性サーフフェス変換、purification（浄化：結合性解析のための前処理）、細線化3D、構造モデルインデックス、最適2値化処理、結合性度を最小化する自動閾値処理、厚み解析（Local Thickness pluginを用いた骨梁厚Tb.Thと分離度Tb.Sp の計算）、ボリュームフラクション（体積比：BV/TV を計算）が用意されています。

（BoneJの解析機能リスト）

また、骨の全体解析のために、楕円フィット（ROIマネージャにあるポイントROIのセットから、最適フィットする楕円体を検出）、球体フィット（上記楕円体フィットの球体バージョン）、モーメント3D（骨全体、構造体内部のモーメントの計算）、ネットシャフト角（大腿骨の曲率とネックシャフト角の計算）、スライスジオメトリ（面積とセクションモジュールの二次モーメントのようなクロスセクションパラメータの解析）が用意されています。
さらに粒子解析として、三次元2値化画像スタックの中から粒子を検出し、計測することが可能です。

BoneJの実践

BoneJ を用いたいくつかのImageJによる実践画像処理について解説します。（ここでご紹介するのは、ほんの一部です。）

ImageJのサンプル画像（bat cochlea volume.tif）のスタック画像を読み込んで、BoneJのthicknessプラグインを実施した例を示します。

このカラーは、3D Local Thickness機能によって厚みのカラーマップを表しており、その数値としては平均厚み、平均スペース距離が計算されています。

次は、同じ画像でボリュームフラクションを計算した例で、この三次元画像の密度が計算されています。

次は、2値化されたセルコロニーの画像を使って最適閾値処理を施してフラクタル解析
結果を表示した例です。

その他、Fly Brain画像を最適な2値化処理を行って異方性を計算し、三次元可視化表示することなども可能です。

ひとを対象した解析には、実際にはμCTで撮影した骨梁のDICOM画像や切除乳房の病巣部のμCTでの構造解析などいろいろな画像で定量評価が試みられています。

今回は、ImageJを用いたBoneJによる解析方法について説明を行いました。
骨梁の解析やそのほか何らかの緻密な内線構造物の評価などで画像処理手順が決まれば、BoneJで使用する機能をマクロでピックアップして自動処理すれば、画像処理・解析作業は効率が良くなると思われます。

参考記事：「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(27・4) 2012, p78-79」

第69回　ImageJを用いたDNAマイクロアレイの画像解析で学ぶ！

ImageJでは、特定のDNAの検出のために行われるDNAマイクロアレイの画像解析が可能です。いくつかのWebサイトでは、さまざまなパターンのDNAチップの解析プラグインなどが紹介されています。今回はDNAマイクロアレイの画像解析の基本的な説明と、ImageJによる実践例をご紹介します。

DNAマイクロアレイとは

DNAは生物の遺伝情報を担う高分子生体物質です。DNAはデオキシリボ核酸の略で、デオキシリボース（五炭糖）と燐酸、塩基から構成される核酸です。核酸は、塩基と酸、リン酸からなりますが、スクレオチド結合で連なった生体高分子です。糖の違い（2位が水素基（DNA）か水酸基（RNA）であるか）によって、2-デオキシリボースを持つデオキシリボ核酸（DNA）とリボースを持つリボ核酸（RNA）とがあります。
塩基には、アデニン、グアニン、チミン、シトシンの4種類があり、それぞれA、G、T、Cと略します。

このDNAについて特定のDNAを抽出し解析する方法のひとつがDNAマイクロアレイ解析と呼ばれています。DNAマイクロアレイ（もしくはDNAチップ）とは、細胞内の遺伝子発現量を測定するために、多数のDNA断片をスライドガラスなどの基板上に高密度に配置したものです。

あらかじめ、塩基配列が既知である1本鎖のDNAを多種、基板上に配置しておき、これに検体を反応させれば、検体のDNA配列と相補的な塩基配列の部分にのみ、検体のDNA鎖が結合します。これを解析することによって、目的の細胞・組織でどの遺伝子が働いていたのかを調べることができます。また、アレイ上での結合位置を蛍光や電流によって検出し、最初の配置から検体に含まれるDNA配列を知ることができます。検体の塩基配列が予測できる場合には、効率的にその配列が特定できます。

次の図は、Gilles Carpentier氏が公開しているウェブサイトのDNA（プロテイン）マイクロアレイ表示例です。

The Protein Array Analyzer

Computer Image Analysis & Biochemistry Practice Works

（http://image.bio.methods.free.fr/ImageJ/?Protein-Array-Analyzer-for-ImageJ.html&lang=en）

ImageJによるDNAマイクロアレイの画像解析

ImageJでは、画像化されたDNAマイクロアレイの解析プラグインがいくつかのWebサイトで紹介されています。
シンプルな解析プラグインとしては、Bob Dougherty and Wayne Rasband氏作成によるMicroArray Profileプラグインがあります。

http://www.optinav.com/MicroArray_Profile.htm

次に、ImageJのサンプル画像"Dot_Blot.jpg"画像をロードし、MicroArray_Profileを実行した例を示します。

デフォルトでは、6行8列の円形ROIが画像のサイズに合わせて表示されるため、これをプラグインメニューの<Reset Grid>ボタンを押して、4行7列に変更してROIを表示しています。

それぞれのROIは、ImageJで操作できますので、位置を微調整可能です。微調整した場合は、<Save Grid>で記録しておけば、次回も同じ位置にROI設定できます（<Read ROI>を利用）。

ROIの計測データは、Analyze>Mesureで結果を参照できます。Mesureの項目は、Analyze>Set Mesurements...で設定可能です。

次の図は、ROIを再配置して、プラグインメニューの<Histgrams>ボタンを押して、各ROIの配列座標ごとに統計値を表示した例です。

マクロを用いたDNAマイクロアレイの画像解析

ImageJのマクロを用いたDNAマイクロアレイの画像解析でユニークなマクロ（こちら）が、上記のGilles Carpentier Research Web Siteで紹介されています。

ここからダウンロードできるマクロ（"Dot Blot Analyzer.txt"）を、ImageJフォルダ下のmacrosフォルダに置いてImageJを起動し、<Plugin>Macros>Install...から選択すれば、インストールできます。
インストールすると、ImageJメニューバーにアイコンが現れます。

（マクロインストール後）

その中の<Dot Blot>の絵が描かれたアイコンをクリックして、Analysis a Dot Blotを選択すると、自動的にROI設定のユーザーインターフェースがついたアプリケーションが起動します。サンプル画像もロードできます。

サンプル画像

新しく出現するインターフェース

その後は、マニュアルを参照しながら、ROIを設置（Add）し、計測が実行できます。
先述のプラグインのように、設定ROIを保存することや、そのロード、計測結果の出力などが可能です。

今回はImageJを用いたDNAマイクロアレイの解析方法について説明を行いました。これは、DNAの解析だけではなく、複数の画像像上に複数のROIを設置して自動で一挙に統計解析するプラグインやマクロとしても使いやすいので機会があれば、いろいろな画像で試してみてはいかがでしょうか。

参考記事：「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(26・11) 2011, p104-105」

登録: 投稿 (Atom)

2022年8月12日金曜日

第72回 ImageJでRadiomicsを体験しよう！

Radiomicsとは

環境

必要なもの

レディオミクスのプロセスの全体像

Radiomic特徴を計算する

ImageJのプラグインとして利用する方法（一例ずつ処理）

コマンドプロンプトから利用する方法（自動で全処理）

機械学習モデルを作成する

Knowledge Flowを使う

縦隔から60歳以上かどうかを予測する

他の秀逸なレディオミクスのツール

レディオミクスの詳しい話について

おわりに

2021年12月17日金曜日

第71回 ImageJで深層学習モデルを動かす！－NVIDIA pre-trained modelの例－

NVIDIA NGC

pre-trained modelのダウンロード

CPUかGPUか

ImageJ/Fijiで動かしてみる（プラグインを開発する）

サンプルデータ

モデルのロードから推論実行まで

まとめ

References

2015年12月24日木曜日

第70回 ImageJを用いたBoneJの紹介で学ぶ！

BoneJとは

BoneJのインストール方法と主な機能

BoneJの実践

第69回 ImageJを用いたDNAマイクロアレイの画像解析で学ぶ！

DNAマイクロアレイとは

ImageJによるDNAマイクロアレイの画像解析

マクロを用いたDNAマイクロアレイの画像解析

第72回　ImageJでRadiomicsを体験しよう！

第70回　ImageJを用いたBoneJの紹介で学ぶ！

第69回　ImageJを用いたDNAマイクロアレイの画像解析で学ぶ！