医用機能画像による診断
- Mathematicaの適用可能性を探る

中嶋憲一, Kenichi Nakajima, MD, PhD
金沢大学先進予防医学研究科・機能画像人工知能学

Dec 6/2019

Join the Conversation @WolframJapan

- 核医学の画像の考え方とその例
- 統計結果のグラフの表示も簡単：多変量データ解析の特徴
- 画像の取り込み，色を付ける，3次元的な断面
- 画像の分解能を計算する（FWHM）
- 関心領域内の取り込み値を知る
- 時間放射能曲線を作る
- 機械学習の適用の例

(発表時の画像はDICOMですが，個人情報保護のため画像化したデータを本ファイル内に置いて修正しています．
画像の転載利用については，nakajima@med.kanazawa-u.ac.jpまでご連絡ください．)

心筋の血流

AnatomyPlot3D[]...

心筋の血流

心筋血流と冠動脈の融合画像

骨転移を見つける

甲状腺

統計データからグラフの表示も簡単：多変量解析の特徴

統計解析（ロジスティックモデル）から求められる関数：心臓の収縮から重症心事故を予測する

(*心機能と重症心事故の関係...入力を心臓の収縮、出力を心事故とする*)

例：統計解析の結果
-------------------------------------------------------------------------------
		Estimate			SE		Chi2		p
Intercept	0.192			0.340	0.32		0.5712
EF		0.0493			0.00627	61.70	<.0001
-------------------------------------------------------------------------------
データを{variable, outcome}のフォーマットでImportとして，以下のようにモデルフィットで計算することもできる．
NonlinearModelFit[data, 1/(1+Exp[-(b +a*x)]),{a,b},x]

estimateIntercept=-0.193;
estimateEf=-0.0493;
event[ef_]:=100/(1+Exp[-(estimateIntercept+estimateEf*ef)])

plot1=Plot[event[x],{x,20,80},AxesOrigin->{20,0},
PlotRange->All,AxesLabel->{"EF(%)","Events (%/3y)"},LabelStyle->Directive[Black,18],PlotTheme->"Classic"]

統計解析（ロジスティック解析）から求められる４変数の係数を視覚化する... 年齢、心臓の収縮、心筋血流欠損、合併症（糖尿病）

(*4変数を使った予測式の例　Nakajima K, et al Circ J 2012;76:168, J Nucl Cardiol 2018 (Epub ahead of print) *)
{b0,b1,b2,b3,b4}={-4.8,0.89,0.056,0.19,-0.048};
risk[dm_,age_,s_,ef_]:=100/(1+Exp[-(b0+b1*dm+b2*age+b3*s+b4*ef)])

Manipulate[risk[dm,age,s,ef],{age,40,80,10},{s,{0,1,2,3}},{ef,20,60,10},{dm,{0,1}},FrameLabel->Style["Events (%/3y)",Large],LabelStyle->Directive[Black,14]]

Plot3D[{risk[1,age,0,ef],risk[1,age,3,ef]},{ef,30,60},{age,40,80},PlotRange->{0,50},PlotTheme->"Scientific",BoxRatios->{1, 1, 1},AxesLabel->{"Cardiac function (%)","Age","Risk (%/3y)"},LabelStyle->Directive[Black,24]]

画像の処理に必要な処理方法

画像の表示方法
・DICOM, JPEG, TIFF, PNG, etc. 
・白黒（グレイスケール）、カラー表示
・３次元表示

画像の解析方法
・関心領域をとる：集積量を知る
・時間的な経過を見る：機能的な変化を知る

得られたデータの解析
・診断用の指標の計算

画像の表示

DICOMデータの読み込み
(脳の線条体のドパミンﾄランスポータの画像)

img1=;

画像の情報を見る

ImageDimensions[img1];
ImageMeasurements[img1,{"Max","Min","Mean","StandardDeviation"},"Dataset"]

imgA=ImageAdjust[img1]

画像に色をつける

imgA//Colorize

(*画像のスムージング*)
imgB=Blur[,3]
imgC=Colorize[imgB,ColorFunction->"Rainbow"];

cm1=ComponentMeasurements[imgB,{"Image","Count","Max","Mean","StandardDeviation"},All,"Dataset"]
Dimensions[cm1];
imgP=cm1[[1,1]];

(*画像の演算も容易*)

imgPsub=(imgP-0.5)*2;
imgCrp=ImageCrop[imgPsub,{90,90}]
ComponentMeasurements[imgCrp,{"Image","Count","Max","Mean","StandardDeviation"},All,"Dataset"]

3次元的分布をみる

SetDirectory[NotebookDirectory[]];
ClearAll;
(*元データはDICOMですが今回はサンプルとして.mov型式から，dat=Import["datsAsym.mov","ImageList"]で読み込んだデータとして以下に含めます*)

dat={,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,};

ndata=Dimensions[dat][[1]];
maxTable=Table[Max[dat[[i]]],{i,1,ndata}];
coeff=1/Max[maxTable];
dat10=dat*coeff;
(*表示のため一部を選択*)
dat10s=dat10[[10;;54]];
dat10st=ImageTrim[dat10s[[#]],{{30,30},{100,100}}]&/@Range[45];
dat10str=Reverse[dat10st];
number=Dimensions[dat10st][[1]];

Manipulate[Image3D[dat10str,ClipRange->{All,All,{z,number}},ColorFunction->"Rainbow"],{{z,number},1,number,1}]

画像分解能を調べる

ClearAll["Global`*"];
(*線状の画像の断面図*)
h1[x_]:=If[And[x>-.5,x<.5],1,0];
k1=Integrate[h1[x],{x,-100,100}]//N;
(*平均＝0，標準偏差＝1の曲線*)
sigma=1;mu=0;
h2[x_]:=1/Sqrt[2Pi*sigma^2]*Exp[-(x-mu)^2/(2*sigma^2)]
Integrate[h2[x],{x,-100,100}]//N;

(*断面のプロフィール*)
pl=Plot[{h1[x],h2[x]},{x,-4,4},PlotRange->All,Filling->Bottom]

(*Full-width at half maximum, FWHM，高さの半分の幅を計算して分解能の指標とする*)
max=1/Sqrt[2Pi]*Exp[-0^2/2]//N;
NSolve[1/Sqrt[2Pi]*Exp[-x^2/2]==0.5*max,x]

l=Graphics[{Arrowheads[{-0.05,.05}],Style[Arrow[{{-1.177,0.5/Sqrt[2 Pi]},{1.177,0.5/Sqrt[2 Pi]}}],Brown]}];
Show[{pl,l},PlotLabel->Style["Full-width at half maximum",Brown,FontSize->24]]

(*実際に撮った線の画像を読み込む*)
imgN=;
ImageDimensions[imgN];
(*mm/pixelに換算: 500mmで512 pixel*)
pix=500/512//N;

(*画像のトリミング*)
imgT=ImageTrim[imgN,{{136,221},{234,319}}]

ImageDimensions[imgT];
(*k行目のプロフィールを作る*)
profile[k_]:=ImageData[imgT][[k,1;;100]]
ListPlot[profile[40],PlotRange->All]

pr=profile[40];n=Range[100];pr30=pr[[40;;69]];n30=Range[0,29];
combined={n30,pr30};
prof=Transpose[combined];

model:=k/s*Exp[-(x-m)^2/(2*s^2)];
fit=FindFit[prof,model,{k,m,s},x];
k1=fit[[1,2]];
m1=fit[[2,2]];
s1=fit[[3,2]];
plot1=ListPlot[prof];
g[x_]:=k1/s1*Exp[-(x-m1)^2/(2*s1^2)];
plot2=Plot[g[x],{x,0,30},PlotRange->All];

Show[{plot1,plot2}]

g[x_]:=k1/s1*Exp[-(x-m1)^2/(2*s1^2)]
findmx=FindMaximum[g[x],{x,0,30}];
max=findmx[[1]];
half=max/2;

nsr1=FindRoot[k1/s1*Exp[-(x-m1)^2/(2*s1^2)]==half,{x,14}];
nsr2=FindRoot[k1/s1*Exp[-(x-m1)^2/(2*s1^2)]==half,{x,20}];
left=nsr1[[1,2]];
right=nsr2[[1,2]];
fwhm=pix*(right-left)

Region of interest: 甲状腺の関心領域

ClearAll["Global`*"];
SetDirectory[NotebookDirectory[]];
thyroid0=;
thyroidT=ImageTrim[thyroid0,{{90,90},{180,180}}]

mask=;

im1=ImageMeasurements[thyroidT,{"Mean","Min","Max","Total"},Masking->mask]
im2=ImageMeasurements[mask,{"Mean","Min","Max","Total"}]
totalCount=im1[[4]]

(*もし投与量が1500なら、摂取率は...*)
totalCount/1500

Regions of Interest using Dynamic and Locator

この画像は123I-MIBGによる交感神経の画像です．心臓と縦隔に関心領域をとって，指標として用います．

img=;
coeff=1/Max[img];
igs=ImageAdjust[img];
igsdata=ImageData[img];

ImageDimensions[img];
dk=ColorNegate[Graphics@Disk[]];
p10={155,120};p11={160,138};
p20={123,160}; p21={133,180};Dynamic[trm1rct=ImageTrim[igs,{p10-Norm[p10-p11],p10+Norm[p10-p11]}];
(* 矩形の心臓領域画像から円形領域を抽出 trm1*)
trm1=ImageMultiply[trm1rct,ImageResize[dk,ImageDimensions[trm1rct][[1]]]]; 
(* 矩形の縦隔領域を抽出 trm2 *)
trm2=ImageTrim[igs,{p20,p21}];
(* 画素値のうち0より大きいもののみを抽出 *)
mat1=Select[Flatten@ImageData[trm1],#>0&];mat2=Select[Flatten@ImageData[trm2],#>0&];
mean1=Mean[mat1];
mean2=Mean[mat2];
std1=StandardDeviation[mat1];std2=StandardDeviation[mat2];
Column[{Show[HighlightImage[Image[igs,ImageSize->600],{EdgeForm[{Green,Thickness[Medium]}],Graphics[Rectangle[p20,p21]]}],Graphics[{Line[{p10,p11}],{White,Thickness[Small],Circle[p10,Norm[p10-p11]]},Locator[Dynamic[p10],Appearance->Small,Background->LightPink],Locator[Dynamic[p11],Appearance->Small,Background->LightBlue],Locator[Dynamic[p20],Appearance->Small,Background->LightPink],Locator[Dynamic[p21],Appearance->Small,Background->LightBlue]},PlotRange->{0,152},ImageSize->600]],
Grid[{{"","Mean","Std dev"},{"Heart",mean1,std1},{"Mediastinum",mean2,std2},{"Heart/Mediastinum ratio",mean1/mean2}},Frame->All]}]]

Time activity curve ...臓器内の時間的経過を見る

腎臓の機能をみる. 元画像はDICOMだが，ここではkidney.movから読み込むように修正．データは最初が1枚2秒，後が1枚20秒なので，グラフ作成時に調整する（40 frames x 2 sec = 80 sec, 56 frames x 20 sec = 1120 sec）．

(*kidney=Import["kidney.mov","ImageList"]*)

kidney={,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,};

tabmax=Table[Max[kidney[[k]]],{k,1,96}];
coeff=1/Max[tabmax];
numberData=Dimensions[kid];
kid1=kidney*coeff;

(*初期80秒と以後の20分間と画像の時間が異なるので加算して調整*)
p1=kidney[[1;;40]];p2=kidney[[41;;96]];
p1a=ImageAdd[p1[[#]],p1[[#+1]],p1[[#+2]],p1[[#+3]],p1[[#+4]],p1[[#+5]],p1[[#+6]],p1[[#+7]],p1[[#+8]],p1[[#+9]]]&/@ Range[1,39,10];

pAdd={p1a,p2}//Flatten;
ndata=Dimensions[pAdd][[1]];

Motion images

f=pAdd;
Manipulate[f[[i]],{i,1,ndata,1}];

fc=Colorize[f[[#]],ColorFunction->"Rainbow"]&/@Range[60];

Manipulate[fc[[i]],{i,1,ndata,1}]

(*時間変化を解析する／ROIをとるために適宜加算してとりやすくしてから*)
fadd=ImageAdd[f[[10]],f[[30]]]

maskB=;
maskL=;
maskR=;

ImageData[maskB]//Flatten;
pxB=Total[%];
ImageData[maskL]//Flatten;
pxL=Total[%];
ImageData[maskR]//Flatten;
pxR=Total[%];

imB=ImageMultiply[#,maskB]&/@f;
imMeasB=ImageMeasurements[imB,"Total"];
datB=imMeasB;
maxB=Max[datB];
imR=ImageMultiply[#,maskR]&/@f;
imL=ImageMultiply[#,maskL]&/@f;
imMeasR=ImageMeasurements[imR,"Total"];
imMeasL=ImageMeasurements[imL,"Total"];
datR=imMeasR-(pxR/pxB)*imMeasB;
maxR=Max[datR];
datL=imMeasL-(pxL/pxB)*imMeasB;
maxL=Max[datL];
ListPlot[{datL,datR},Filling->Bottom]

(*時間軸をつけてグラフ化*)

rg={Range[0,1180,20]}//Flatten;
dataR=Transpose[{rg,datR}];
dataL=Transpose[{rg,datL}];

ListPlot[{dataR,dataL},Filling->Bottom,Joined->True,PlotLegends->{"Right","Left"},PlotTheme->"Detailed"]

(*左右の腎機能の比率や、腎臓の糸球体濾過量も推定可能（略）*)
addR=Total[datR[[8;;12]]];
addL=Total[datL[[8;;12]]];
addRL=addR+addL;
Grid[{{"","R","L","Total"},{"Count",addR,addL,addRL},{"%",addR 100/addRL,addL 100/addRL,addRL 100/addRL}},Dividers->{{False,True},{False,True}}]

新しい関数：pandamon （kumamonとpandaを識別するclassifier）

kumamon   Panda

trainingDataKp={,,,,,,,,,,,,,,,,,,,,,,,};
{k,p}={"kumamon","panda"};
answer={Table[k,12],Table[p,12]}//Flatten;

pandamon=Classify[trainingDataKp->answer,Method->"NaiveBayes"]
Information[pandamon]

trialR1={,,}
pandamon[trialR1,{"Decision","Probabilities"}]//TableForm

消化管の中の食物の通過を見てみよう

{->大腸，->胃, ー＞小腸}のようなセットを作って，トレーニングも可能．
トレーニング用のデータセットを作る(trainingSet)．
テスト用のデータセットを作る(testSet)．
cs=classify[traningSet]
Information[cs]
cs[testSet,{“Decision”,”TopProbabilities”}]
cm=ClassifierMeasurements[cs,comp,”ConfusionMatrixPlot”]

Gogh and Monet painting

まとめ

・画像に由来する様々なデータの視覚化
・放射性医薬品の集積量の定量が可能
・放射性医薬品の経時的変化の解析が可能
・３次元画像への適用
・画像専用ソフトでのみ解析可能であったデータの新規解析の可能性: 新たな関数とツールの作成に期待
・機械学習による画像診断

・小道具として使う〜本格的な画像処理として使う