専攻実験をPythonでやってみた（その２）

前回の続きです。
smon.hatenablog.com

今回も，実験テキストでサンプルコードとして載っている部分をPythonで書いてみます。
自分でコーディングしてもらう部分や，考察に関わる部分などは一切触れていません。（授業・実験の公平性を期すためです）

正方形を描く

import numpy as np
import matplotlib.pyplot as plt

nx = 256
ny = 256

ar = np.zeros((nx, ny))

for j in range(ny):
    for i in range(nx):
        if 64 <= i <= 192 and 64 <= j <= 192:
            ar[j][i] = 100
        else:
            ar[j][i] = 0

plt.imshow(ar)

実行すると正方形が表示されます。

最初にnumpyとmatplotlibというのをインポートしていますが，numpyは数値計算ライブラリ，matplotlibはグラフ描画などに使うライブラリです。
とくにnumpyは関数や配列などを扱う際には必須のとても重要かつ便利なライブラリ。これのおかげで学術分野でPythonが使いやすくなっています。

また，

a < b < c  # a < b and b < c と同義

のように2つ以上の不等式を一発で書けるのも地味に便利。

サイノグラムの作成

逆投影を行う前のサイノグラムを数値データで作るという作業があるのですが，そのサンプルコードをPythonで書き直すとこんな感じ。

import numpy as np
import matplotlib.pyplot as plt

nx = 256
ny = 256
ntheta = 256

x1 = 32.0
y1 = 32.0
r  = 24.0

sino = np.zeros((nx, ny))

for k in range(ntheta):
    theta = 2.0 * np.pi * k / ntheta
    x0 = x1 * np.cos(theta) + y1 * np.sin(theta)
    
    for i in range(nx):
        x = i - nx/2
        
        if -r <= x-x0 <= r:
            sino[k][i] = np.sqrt(r**2 - (x0-x)**2)
        else:
            sino[k][i] = 0.0

plt.imshow(sino)

実行するとサイノグラムが表示されます。

さらにこの配列（sino）をバイナリ形式で保存したければ，

sino.tofile('sinogram.flt')

と最後に1行書くだけでfltファイルとして保存できてしまいます。すばらしい。