Pythonプログラミング(ステップ7・配列・CT)
このページでは、典型的な逆問題である、CT (Computed Tomography)の動作原理について考える(ここは作成中)。
CTを使った画像診断
X線のように指向性と透過性の強い電磁波を物体に照射すると、内部の組成や密度に応じて吸収され、「反対側」で検出される信号の強度が変化する。
線源を並行移動しながらスキャンすることで、ある方向から見た、吸収の程度の分布を得ることができる。
センサーの方向角を
写真は筆者の脳画像
センサーを回転させて、あるいは物体を回転させさがら、様々な
これに似た画像診断装置にMRIがあるが、MRIの場合は、空間的に磁場の勾配を設けることで、特定の面内からの核磁気共鳴の信号を検出し、 水の密度分布を推定している。 上図との対応では、紫色で描いた「線」が、奥行きのある「面」に変わる。
図のようなセンサーを基準にした座標上で、計測対象となる物体の位置
もしも物体内部の吸収率分布
ところが、計測可能なのは
CTのシミュレーション
モノクロ画像を読み込んで、その画素値を
サンプル画像データ: lena.png
ここで用いる画像は
[ [角度1, 位置1, 吸収量1], [角度2, 位置2, 吸収量2], [角度3, 位置3, 吸収量3], ... ]
として構成し、NumPyの配列データとしてファイルに保存する(ファイル名:scanned_data.npy
)。
具体的には、画素値が二次元の正方格子上の点に割当られていると考え、
「X線」を想定した直線とそれら格子点との距離を求め、格子点からの距離
# coding: utf-8 import numpy as np from PIL import Image import math def calc_signal(img, angle, step): NX = img.shape[0] NY = img.shape[1] x_c = NX/2 y_c = NY/2 nx = math.cos(angle*math.pi/180) ny = math.sin(angle*math.pi/180) x0 = x_c + nx * step y0 = y_c + ny * step c = -nx*x0 - ny*y0 sig=0 for x in range(NX): for y in range(NY): dist = abs(nx*x + ny*y + c) if dist<5: wt = math.exp(-dist*2/3.0) sig = sig + wt*img[x,y] return sig def scan_image(img): data=[] L = int(np.max(img.shape)/2*math.sqrt(2)+1) for angle in np.arange(0,180,5): # step angle of 5 degrees print("angle=",angle) for step in range(-L,L,1): sig = calc_signal(img, angle, step) if sig>0: data.append([angle,step,sig]) return np.array(data) # img = np.array(Image.open('lena.png')) data = scan_image(img) np.save('scanned_data.npy',data)
画像の再構成
では、前節で生成されたデータを読み込んで、画像を再構成してみよう。
もともとの画像のサイズは既知であるとして、未知である画素値を1次元配列(あるいはベクトル)として表現しておく(配列 X
)。
「X線」を表す直線が与えられれば、
格子点と添字の対応関係は任意であるので、ここでは、格子の左から右に、および上から下に順に番号を振ることにする。
以下のコードでは、
格子の幅をNX
とすれば、格子点上の座標x,y
の点の添字を k = x + NX*y
としている。
モデルに従って係数を決め、格子点(画素)の値を未知数とする連立方程式を解く問題に帰着させる
行列
未知数の数
さらに、ここで示したPythonのコードでは、計算の効率化のため、直線(ビーム)までの距離が
サンプリングするデータ点が問題の規模に対して少ない、すなわち
ここでは、SciPyライブラリの中の疎行列用のソルバーである sparse.linalg.lsqr を使って、この逆問題を解かせてみよう。
以下に具体的なコードの例を示す。
前節のコードが生成したデータ(scanned_data.npy
)を読み込み、
計算結果から再構成した画像をファイル(reconstructed_image.png
)に出力する。
なお、大きなサイズの行列計算を伴うので、一般的なパソコンでは処理が完了するまで数分程度以上を要するので注意すること。
# coding: utf-8 import numpy as np from scipy.sparse import lil_matrix,csc_matrix from scipy.sparse.linalg import lsqr from PIL import Image import math def gen_col_index(NX,NY,angle,step): L = int(max(NX,NY)*1.42) nx = math.cos(angle*math.pi/180) ny = math.sin(angle*math.pi/180) x0 = NX/2 + nx * step y0 = NY/2 + ny * step c = -nx*x0 - ny*y0 cols=[ ] if abs(nx)>abs(ny): for y in range(NY): x1 = int(-(ny*y+c)/nx+0.5) for x in range(x1-5,x1+6): if x<0 or x>=NX: continue dist = abs(nx*x+ny*y+c) if dist<5: wt = math.exp(-dist**2/3.0) k = x + y*NX cols.append([k,wt]) else: for x in range(NX): y1 = int(-(nx*x+c)/ny + 0.5) for y in range(y1-5,y1+6): if y<0 or y>=NY: continue dist = abs(nx*x+ny*y+c) if dist<5: wt = math.exp(-dist**2/3.0) k = x + y*NX cols.append([k,wt]) return cols def gen_matrix(NX,NY,data): N = data.shape[0] M = NX*NY W = lil_matrix((N,M)) F = np.zeros((N,)) for j,x in enumerate(data): angle = x[0] step = x[1] val = x[2] cols = gen_col_index(NX,NY,angle,step) for k,wt in cols: W[j,k]=wt F[j]=val return W,F ### data = np.load('scanned_data.npy') NX=128 NY=128 W,F = gen_matrix(NX,NY,data) W = csc_matrix(W) res = lsqr(W,F,show=True,iter_lim=100) X = res[0] image = np.zeros((NX,NY)) for x in range(NX): for y in range(NY): val = X[x + y*NX] image[x,y] = val pil_img = Image.fromarray(image.astype(np.uint8)) pil_img.show() pil_img.save('reconstructed_image.png')
上記のコードで再構成された画像を以下に示す。
練習:画像を使った実験
グレースケールの小さめの画像ファイルを用意し、それを使って上記のシミュレーションを追試してみなさい。
- 上の例では、方向角を5度毎にサンプルを取得しているが、それを倍程度、あるいは半分程度にすると結果はどのように変わるか。
- それぞれの格子点の影響度を与えるモデル(式(2))が結果にどのように影響するか。