20180709

脳ドッグに行ってきた。→MRIの画像データをもらった。(DICOMファイル)→k-spaceは入ってなかった。

先に断っておく。

  • 医療情報っぽいことを書いたが、自分は医療従事者でもなんでもない。
  • 嘘を書いたつもりはない。しかし、間違っていても、責任は持てない。

脳ドッグに行ってきた。→MRIの画像データをもらった。(DICOMファイル)→残念ながら、k-spaceは入ってなかった。
しょうがないので、画像データを二次元フーリエ変換してみた。

MRIとは

MRI(Magnetic Resonance Imaging: 磁気共鳴画像)のこと。
装置は、google検索でトップだった、このサイトの画像のまさにこれ。これの中に入ってきた。

MRIの結果とフーリエ変換

MRIによって、脳とか、心臓の断面画像が手に入る。

この画像の元データがk-spaceとかいうデータで、
これを二次元フーリエ変換すると、断面画像ができるということで、これが欲しかった。
(逆二次元フーリエではないらしい。空間周波数のデータを、フーリエ変換するって何かいつもと違うのがいやだけど。)
脳ドッグでは、MRIを使うそうな。
だから、無保険(16万自腹)でちょっと高めの脳ドッグを受けてみた。

しかし、遊び目的ではない。

たぶんTIA(Transit Ischemic Attack:一過性脳虚血発作)ってやつにかかったかもしれないと思う出来事があった。
だから、脳の精密検査をしたかった。こっちが本当はメイン。

先日、ジムに行っていたときのこと。
週1で、トレッドミル(ランニングマシン)で10km走るのが習慣になっていた。10.5 km/h
換気が悪いなぁと思っていた。

そしたら汗をいつも以上にかいてしまったらしく、脱水症っぽい感じになった。
さらにひどいことに、走り終わった後、手がしびれ、ろれつが回らない状態になった。
こんなん初めて。

ポカリをたくさん飲んだら治った。
先月は、これについて書こうとしていたが、pythonの話とは毛色が違い、話がまとまらなかったため、投稿を断念した。上記のようなことを書こうとしていた。

脳が悪いのか!?と思って怖くなった。

#7119(神戸で救急車が必要か判断できない人のための電話番号)に電話したら、看護師さんにつないでくれて、事情を、上記以上に話して、何個か簡易テスト的なのを指示されて行ったら。「緊急性はないですね。けど、早いうちに医療機関に受診されたほうがいいですよ」とのこと。

もともと、k-space画像が欲しかったので、いわれなくても受けるつもりではあった。
しかし、もともと脳ドッグを予約していた日時(2か月くらい先)を言うと、
看護師さんがちょっと訂正するような口調で「っ、もうちょっと、早いほうが、いいですね」と言われてしまった。

で、直近の土曜日に、脳ドッグを受けた。

その結果が、3週間後(今)に届いた。
結果は、大丈夫だった。
どっか、血管が詰まっているとかはなかった。

ただ、そのほか健康面(コレステロールとか血圧とか)で、所見があったので、結果として行ってよかった。

MRIした感想

  • 狭い
  • 装置の音がすごかった。
  • ヘッドフォン装着→音楽は選べた。
  • 体は拘束された状態だった。

  • 「息を吸ってー、、、息を止めて。」と何回か指示が来た。

  • 「血管を拡げるげるお薬」なるスプレーを舌にかけられた。→寝すぎた時の頭痛と同じ状態に。
    スプレー1吹きで、体にあんな変化があるなんて、「医療すごい」ってなった。

DICOMファイル

医療情報の画像データの規格らしい。
画像のみならず、音声データ等も保存できるそうな。
その他、情報がたくさん詰まっている。

「これに、k-spaceが!?」と思い、
メモ帳で開くと、ところどころ文字化け。バイナリファイルなんでしょう。きっと。

DICOMファイルは、アップロードしない。
なぜならば、

  • 自分の個人情報
  • 受診した病院名 所在地
  • オペレーターの氏名

などが書いてあったから。
自分のは情報はいいとして、自分以外の情報を公開するのは気が引ける。
たとえやっても法律上問題ないとしても気が引ける。

なので、結果画像だけ読み込んでみた。
そしたら、すでに画像化されていた。
k-spaceのデータは含まれていないようです。

まぁそりゃそうだよなぁ。
医療従事者が二次元フーリエ変換なんで遊ぶ必要ないもん。
患者さんを助ける情報が必要なだけだもん。

というか、オペレーターの方が、「k-spaceはあげられないんですよぉー」(にこっ)
っておっしゃっていた気がする。いまさら思い出す。
MRIの技師の方は、二次元フーリエ変換知ってた。やはり使うんだろうか。

dicomファイルの中に、機械の製造者の情報(?)が、Philipsとなっていた。
Philipsといえばシェーバーの会社じゃないですか。
髭剃りだけじゃなかったんだな。
https://ja.wikipedia.org/w/index.php?title=%E3%83%95%E3%82%A3%E3%83%AA%E3%83%83%E3%83%97%E3%82%B9&oldid=69059425

pydicom

pythonでdicomファイルを開くには、このpydicomがいる。

  1. dicomファイルを読む(画像のnumpy arrayを含むデータの読み込み)
  2. matplotlibで、numpy arrayをimshow
    という、いつもの流れだった。

PythonでDICOM画像をなんとかするここが詳しい。

cmdで、全部小文字で、

pip install pydicom

だった気がする。間違ってたらごめん。
当然ですがpipが入ってないと、効かない。

環境

python2.7
pydicom, numpy, matplotlib, tkFileDialog

コード

  • このプログラムは、DICOMファイルを開いて、画像化するプログラムです。
  • 二次元フーリエ変換した画像も併せて載せていますが、これは(多分)k-spaceではありません。
  • DICOM以外のファイルを開いた時の挙動は意図していません。
import pydicom as dc
import numpy as np
import matplotlib.pyplot as plt


import
tkFileDialog as tk
name = tk.askopenfilename(filetypes=[("dicom",("all","*.*"))])
d = dc.read_file(name)


mx = np.max(np.real(np.fft.fft2(d.pixel_array)))/100000.

f,((ax1,ax2),(ax3,ax4))=plt.subplots(2,2)

ax1.imshow(d.pixel_array, cmap="gray")
ax2.imshow(np.abs(np.fft.fft2(d.pixel_array)), vmin=0, vmax=mx*1000, cmap="gray")
ax3.imshow(np.real(np.fft.fft2(d.pixel_array)), vmin=-mx, vmax=mx, cmap="gray")
ax4.imshow(np.imag(np.fft.fft2(d.pixel_array)), vmin=-mx, vmax=mx, cmap="gray")

ax1.set_title("pydicom pixel_data")
ax2.set_title("2DFFT absolute value")
ax3.set_title("2DFFT real part")
ax4.set_title("2DFFT imaginary part")

plt.show()

プログラムの実行方法

上記コードを、テキスト帳にコピーし、「a.py」というファイル名で保存。

そのファイルが表示されているブラウザ上で、普段アドレスが表示されているところに、

cmd

と入力し、Enterキーを押す。
(もしくは、ブラウザ上の何もないところを、shift+右クリックして、「PowerShellWindowをここに開く」をする。)

そうすると、黒い画面(青い画面)が出てくる。

その画面上で、

python a.py

と入力すると、実行できる。ファイルダイアログが出るはず。
(pythonのフォルダが、環境変数で適切に設定してあればの話。)dicomファイルを選ぶと、画像が表示される。


figure_1.png

  • 画像データはnumpy arrayで返ってくる
  • 画像は、mm感覚で撮影されているらしい
    なので、3次元でプロットすれば、立体的に見れるはず。
    大変そうだからまだやらない。
posted by yuchan at 07:00 | Comment(0) | python

20180719

Femapで円筒座標系を使って荷重をかける時の注意点。

有料ソフトなので、画像は載せなかった。
絵がないと、分かりにくい内容だけど。

なぜFemapを使っているのか。

計算力学技術者試験の受験資格を手に入れたい。
計算力学技術者の初級(2級の受験資格になる)は、講習で手に入る(結構高い)。

その講習で、Femapを使った。
Femapは300ノード制限バージョンを使用した。

Femapをいじっていて、使ったことのない円筒座標系を使ってみようと思った。
円筒座標系の荷重でつまづいたので、書き残す。

先に結論を書いておくと、
円筒座標系の中心は、原点。任意の場所を円筒の中心にできない。

円筒座標系で荷重をかけたかった。

こうすりゃいいんじゃない?と思って、

  • 座標系を「全体直交座標系」→「全体円筒座標系」に変更。
  • 方向を「ベクトル」にしdR, dT, dZ(中心からの距離Radius、角度Theta, 高さ座標Z)を(1., 0., 0.)に変更。
  • 荷重を入力。

これだけでできると思っていた。
(通常は、これだけでできる。)

[あらぬ方向に荷重方向の矢印が伸びている絵をここに乗せたかった。]

うまくいっていません。
明らかに、円筒の中心から伸びてきた矢印ではない。
この矢印の根元は、どこをさしているのかというと、(0, 0, 0)

問題なのは、モデルの原点。

最初に書かなかったけど、すでにモデルはできていた。
そのとき、円筒座標系のことを何も考えずに作ってしまっていた。
円筒の中心は、原点から離れた場所に作っていた。

これが原因だった。

円筒座標系とベクトルを使って荷重をかけると、原点から線が伸びる。
だから、荷重方向がおかしなことになった。

これは、絵がないとわからないでしょう。
なのに絵は載せられない(怖くて)。

解決策

  1. 円筒の中心点を、原点に移動する。
  2. 座標系を「全体直交座標系」→「全体円筒座標系」に変更。
  3. 方向を「ベクトル」にしdR, dT, dZ(中心からの距離Radius、角度Theta, 高さ座標Z)を(1., 0., 0.)に変更。
  4. 荷重を入力。

2以降は、上と同じ。

局所円筒座標系はないのか

円筒座標系の中心は、原点。任意の場所を円筒の中心にできない。
「全体円筒座標系」の中心は、絶対原点。
なら、「局所円筒座標系」とか作って、任意の場所を中心に荷重をかけられるようにしてくれたらいいのに。
あっても多分使わないけど。

  • 2個以上の筒状のタンクか何かが、連結されていて
  • 筒の内部で、圧力ではない荷重が放射状にかかっている
    こういうのを模擬するときに使えるかもしれないと思った。

まぁ、そんなん無くても、三角関数(円周上の座標(cosΘ, sinΘ))を使って値を出して、ノード一点一点にその値を当てはめればできますが。

posted by yuchan at 07:00 | Comment(0) | FEMAP