エンコーダ用の階層型データを整える(vim-2)

ここでは2種類のデータを扱うことになる。視覚刺激としての自然動画像と、fMRIで測定したBOLD信号(血中酸素濃度に依存する信号)の2種類である。前者はいわば「入力」で、後者は脳活動の指標で、「応答」である。

Stimuli Image

fMRI Image

主要の目的として、入力画像から脳活動を正確に予測することを掲げる。イメージとして、画像を応答信号へと「符号化」していることから、このシステムを「エンコーダ」と呼ぶことが多い。

データセットの概要

今回使うデータの名称は「Gallant Lab Natural Movie 4T fMRI Data set」で、通称はvim-2である。その中身を眺めてみると:

視覚刺激

Stimuli.tar.gz 4089110463 (3.8 GB)

BOLD信号

VoxelResponses_subject1.tar.gz 3178624411 (2.9 GB)
VoxelResponses_subject2.tar.gz 3121761551 (2.9 GB)
VoxelResponses_subject3.tar.gz 3216874972 (2.9 GB)

視覚刺激はすべてStimuli.matというファイルに格納されている(形式:Matlab v.7.3)。その中身は訓練・検証データから構成される。

  • st: training stimuli. 128x128x3x108000 matrix (108000 128x128 rgb frames).
  • sv: validation stimuli. 128x128x3x8100 matrix (8100 128x128 rgb frames).

訓練データに関しては、視覚刺激は15fpsで120分間提示されたため、7200の時点で合計108000枚のフレームから成る。検証データについては、同様に15fpsで9分間提示されたため、540の時点で合計8100枚のフレームから成る。検証用の視覚刺激は10回被験者に提示されたが、今回使う応答信号は、その10回の試行から算出した平均値である。平均を取る前の「生」データは公開されているが、ここでは使わない。

あと、データを並べ替える必要は特にない。作者の説明:

"The order of the stimuli in the "st" and "sv" variables matches the order of the stimuli in the "rt" and "rv" variables in the response files."

前へ進むには、これらのファイルを解凍しておく必要がある。

$ tar -xzf Stimuli.tar.gz
$ tar -xzf VoxelResponses_subject1.tar.gz
$ tar -xzf VoxelResponses_subject2.tar.gz
$ tar -xzf VoxelResponses_subject3.tar.gz

するとStimuli.matおよびVoxelResponses_subject{1,2,3}.matが得られる。階層的な構造を持つデータなので、開閉、読み書き、編集等を楽にしてくれるPyTables( http://www.pytables.org/usersguide/index.html )というライブラリを使う。シェルからファイルの中身とドキュメンテーションを照らし合わせてみると、下記のような結果が出てくる。

$ ptdump Stimuli.mat
/ (RootGroup) ''
/st (EArray(108000, 3, 128, 128), zlib(3)) ''
/sv (EArray(8100, 3, 128, 128), zlib(3)) ''

かなり単純な「階層」ではあるが、RootGroupにはフォルダー(stsv)が2つある。それぞれの座標軸の意味を確認すると、1つ目は時点、2つ目は色チャネル(RGB)、3つ目と4つ目のペアは2次元配列における位置を示す。

次に応答信号のデータに注視すると、もう少し豊かな階層構造が窺える。

$ ptdump VoxelResponses_subject1.mat 
/ (RootGroup) ''
/rt (EArray(73728, 7200), zlib(3)) ''
/rv (EArray(73728, 540), zlib(3)) ''
/rva (EArray(73728, 10, 540), zlib(3)) ''
(...Warnings...)
/ei (Group) ''
/ei/TRsec (Array(1, 1)) ''
/ei/datasize (Array(3, 1)) ''
/ei/imhz (Array(1, 1)) ''
/ei/valrepnum (Array(1, 1)) ''
/roi (Group) ''
/roi/FFAlh (EArray(18, 64, 64), zlib(3)) ''
/roi/FFArh (EArray(18, 64, 64), zlib(3)) ''
/roi/IPlh (EArray(18, 64, 64), zlib(3)) ''
/roi/IPrh (EArray(18, 64, 64), zlib(3)) ''
/roi/MTlh (EArray(18, 64, 64), zlib(3)) ''
/roi/MTplh (EArray(18, 64, 64), zlib(3)) ''
/roi/MTprh (EArray(18, 64, 64), zlib(3)) ''
/roi/MTrh (EArray(18, 64, 64), zlib(3)) ''
/roi/OBJlh (EArray(18, 64, 64), zlib(3)) ''
/roi/OBJrh (EArray(18, 64, 64), zlib(3)) ''
/roi/PPAlh (EArray(18, 64, 64), zlib(3)) ''
/roi/PPArh (EArray(18, 64, 64), zlib(3)) ''
/roi/RSCrh (EArray(18, 64, 64), zlib(3)) ''
/roi/STSrh (EArray(18, 64, 64), zlib(3)) ''
/roi/VOlh (EArray(18, 64, 64), zlib(3)) ''
/roi/VOrh (EArray(18, 64, 64), zlib(3)) ''
/roi/latocclh (EArray(18, 64, 64), zlib(3)) ''
/roi/latoccrh (EArray(18, 64, 64), zlib(3)) ''
/roi/v1lh (EArray(18, 64, 64), zlib(3)) ''
/roi/v1rh (EArray(18, 64, 64), zlib(3)) ''
/roi/v2lh (EArray(18, 64, 64), zlib(3)) ''
/roi/v2rh (EArray(18, 64, 64), zlib(3)) ''
/roi/v3alh (EArray(18, 64, 64), zlib(3)) ''
/roi/v3arh (EArray(18, 64, 64), zlib(3)) ''
/roi/v3blh (EArray(18, 64, 64), zlib(3)) ''
/roi/v3brh (EArray(18, 64, 64), zlib(3)) ''
/roi/v3lh (EArray(18, 64, 64), zlib(3)) ''
/roi/v3rh (EArray(18, 64, 64), zlib(3)) ''
/roi/v4lh (EArray(18, 64, 64), zlib(3)) ''
/roi/v4rh (EArray(18, 64, 64), zlib(3)) ''

応答のデータでは、RootGroupのなかには、まずrtrvrvaという3つの配列がある。これらはBOLD信号の測定値を格納している。また、roieiと名付けられたgroupがある。前者は応答信号の時系列とボクセルを結びつけるためのインデックスである。後者は実験の条件等を示す数値が格納されている。ここでroiのほうに注目すると、計測された脳の領域全体を分割して(構成要素:ボクセル)、それを生理的・解剖学的に関心を持つべき「関心領域」(ROI)に振り分けていくのだが、このroiなるグループは、各ROIの名が付いた配列を含む。たとえば、v4rhとはV4という領域で、左半球(left hemisphere)に限定したROIのことである。明らかなように、ROIの数はボクセル数($18 \times 64 \times 64 = 73728$)よりも遥かに少ないので、各ROIが多数のボクセルを含むことがわかる。

入力パターン(視覚刺激):元のデータを調べる

それでは、蓋を開けて画像データを見てみよう。

In [1]:
import numpy as np
import tables
import matplotlib
from matplotlib import pyplot as plt
import pprint as pp

# Open file connection.
f = tables.open_file("data/vim-2/Stimuli.mat", mode="r")

# Get object and array.
stimulus_object = f.get_node(where="/", name="st")
print("stimulus_object:")
print(stimulus_object)
print(type(stimulus_object))
print("----")

stimulus_array = stimulus_object.read()
print("stimulus_array:")
#print(stimulus_array)
print(type(stimulus_array))
print("----")

# Close the connection.
f.close()

# Check that it is closed.
if not f.isopen:
    print("Successfully closed.")
else:
    print("File connection is still open.")
stimulus_object:
/st (EArray(108000, 3, 128, 128), zlib(3)) ''
<class 'tables.earray.EArray'>
----
stimulus_array:
<class 'numpy.ndarray'>
----
Successfully closed.

指定した配列の中身が予想通りのものを格納していること、正しく読み込めていることなどを確かめる必要がある。

In [2]:
# Open file connection.
f = tables.open_file("data/vim-2/Stimuli.mat", mode="r")
In [3]:
# Get object and array.
stimulus_object = f.get_node(where="/", name="st")
print("stimulus_object:")
print(stimulus_object)
print(type(stimulus_object))
print("----")
stimulus_object:
/st (EArray(108000, 3, 128, 128), zlib(3)) ''
<class 'tables.earray.EArray'>
----
In [4]:
# Print out basic attributes.
stimulus_array = stimulus_object.read()
#print("stimulus_array:")
#print(stimulus_array)
print("(type)")
print(type(stimulus_array))
print("(dtype)")
print(stimulus_array.dtype)
print("----")
(type)
<class 'numpy.ndarray'>
(dtype)
uint8
----
In [5]:
# Print out some frames.
num_frames = stimulus_array.shape[0]
num_channels = stimulus_array.shape[1]
frame_w = stimulus_array.shape[2]
frame_h = stimulus_array.shape[3]

frames_to_play = 5

oneframe = np.zeros(num_channels*frame_h*frame_w, dtype=np.uint8).reshape((frame_h, frame_w, num_channels))
im = plt.imshow(oneframe)

for t in range(frames_to_play):
    oneframe[:,:,0] = stimulus_array[t,0,:,:] # red
    oneframe[:,:,1] = stimulus_array[t,1,:,:] # green
    oneframe[:,:,2] = stimulus_array[t,2,:,:] # blue
    plt.imshow(oneframe)
    plt.show()
In [6]:
# Close file connection.
f.close()

十分なフレーム数を重ねてみると、読み込んだデータが確かに動画像であることが確認できた。但し、向きがおかしいことと、フレームレートがかなり高いこと、この2点は後ほど対応する。

練習問題 (A):

  1. 座標軸を適当に入れ替えて(help(np.swapaxes)を参照)、動画の向きを修正してみてください。

  2. 上では最初の数枚のフレームを取り出して確認したのだが、次は最後の数枚のフレームを確認しておくこと。

  3. 上記は検証データであったが、同様な操作(取り出して可視化すること)を訓練データstに対しても行うこと。

  4. 上記の操作で、最初の100枚のフレームを取り出してください。フレームレートが15fpsなので、7秒弱の提示時間に相当する。上記のコードのforループに着目し、rangeに代わって適当にnp.arangeを使うことで、時間的なダウンサンプリングが容易にできる。15枚につき1枚だけ取り出して、1fpsに変えてから、最初の100枚を見て、変更前との違いを確認してください。

入力パターン(視覚刺激):空間的ダウンサンプリング

画素数が多すぎると、後の計算が大変になるので、空間的なダウンサンプリング、つまり画像を縮小することが多い。色々なやり方はあるが、ここではscikit-imageというライブラリのtransformモジュールから、resizeという関数を使う。まずは動作確認。

In [2]:
from skimage import transform
from matplotlib import pyplot as plt
import imageio
import tables
import numpy as np

im = imageio.imread("img/bishop.png") # a 128px x 128px image
In [3]:
med_h = 96 # desired height in pixels
med_w = 96 # desired width in pixels
im_med = transform.resize(image=im, output_shape=(med_h,med_w), mode="reflect")

small_h = 32
small_w = 32
im_small = transform.resize(image=im, output_shape=(small_h,small_w), mode="reflect")

tiny_h = 16
tiny_w = 16
im_tiny = transform.resize(image=im, output_shape=(tiny_h,tiny_w), mode="reflect")
In [4]:
myfig = plt.figure(figsize=(18,4))

ax_im = myfig.add_subplot(1,4,1)
plt.imshow(im)
plt.title("Original image")
ax_med = myfig.add_subplot(1,4,2)
plt.imshow(im_med)
plt.title("Resized image (Medium)")
ax_small = myfig.add_subplot(1,4,3)
plt.imshow(im_small)
plt.title("Resized image (Small)")
ax_small = myfig.add_subplot(1,4,4)
plt.imshow(im_tiny)
plt.title("Resized image (Tiny)")

plt.show()

予想するように動いているようなので、視覚刺激の全フレームに対して同様な操作を行う。

In [5]:
# Open file connection.
f = tables.open_file("data/vim-2/Stimuli.mat", mode="r")
In [6]:
# Specify which subset to use.
touse_subset = "sv" 
In [7]:
# Get object, array, and properties.
stimulus_object = f.get_node(where="/", name=touse_subset)
stimulus_array = stimulus_object.read()
num_frames = stimulus_array.shape[0]
num_channels = stimulus_array.shape[1]
In [8]:
# Swap the axes.
print("Stimulus array (before):", stimulus_array.shape)
stimulus_array = np.swapaxes(a=stimulus_array, axis1=0, axis2=3)
stimulus_array = np.swapaxes(a=stimulus_array, axis1=1, axis2=2)
print("Stimulus array (after):", stimulus_array.shape)
Stimulus array (before): (8100, 3, 128, 128)
Stimulus array (after): (128, 128, 3, 8100)
In [9]:
# Prepare new array of downsampled stimulus.
ds_h = 64
ds_w = 64

new_array = np.zeros(num_frames*num_channels*ds_h*ds_w, dtype=np.float32)
new_array = new_array.reshape((ds_h, ds_w, num_channels, num_frames))
In [10]:
# Iterate over frames to be resized.
for t in range(num_frames):
    new_array[:,:,:,t] = transform.resize(image=stimulus_array[:,:,:,t],
                                          output_shape=(ds_h,ds_w),
                                          mode="reflect")
    if t % 1000 == 0:
        print("Update: t =", t)
        
f.close()
Update: t = 0
Update: t = 1000
Update: t = 2000
Update: t = 3000
Update: t = 4000
Update: t = 5000
Update: t = 6000
Update: t = 7000
Update: t = 8000

ここで一つ注意すべきことは、縮小したあとの画像は、$\{0,1,\ldots,255\}$ではなく、$[0,1]$の実数値を取ることである。そのため、dtypenp.float32に変えている。

In [12]:
print("(Pre-downsize) max:", np.max(stimulus_array),
      "min:", np.min(stimulus_array),
      "ave:", np.mean(stimulus_array))
print("(Post-downsize) max:", np.max(new_array),
      "min:", np.min(new_array),
      "ave:", np.mean(new_array))
(Pre-downsize) max: 255 min: 0 ave: 90.50942391854745
(Post-downsize) max: 1.0 min: 0.0 ave: 0.3549391

バイナリ形式にして、ディスクに保存しておくことにする(PyTablesと比べて、速いかどうかは各自で試すと良い)。

In [13]:
fname = "data/vim-2/" + str(touse_subset) + "_downsampled.dat"
with open(fname, mode="bw") as fbin:
    new_array.tofile(fbin)

fname = "data/vim-2/" + str(touse_subset) + "_downsampled.info"
with open(fname, mode="w") as f:
    f.write("dtype: "+str(new_array.dtype)+"\n")
    f.write("shape: "+str(new_array.shape)+"\n")

練習問題 (B):

  1. 上記の空間的ダウンサンプリングを訓練データ(st)に対しても行なうこと。

  2. さらに思い切ったダウンサンプリングを試してみること。32x32、16x16などと小さくすることで、どのような利点があると思われるか。デメリットはどのようなものがあるか。

期待通りに作動しているかどうか確認するために、先ほど保存したファイルを読み込み、その中身を確かめる。まず、保存したデータの寸法を読み込む:

In [14]:
! cat data/vim-2/st_downsampled.info
! cat data/vim-2/sv_downsampled.info
dtype: float32
shape: (64, 64, 3, 108000)
dtype: float32
shape: (64, 64, 3, 8100)

これらの値を下記のように使うと、確かに期待通りのデータが保存されていることがわかる。

In [23]:
print("Training data:")
fname = "data/vim-2/" + "st" + "_downsampled.dat"
with open(fname, mode="br") as fbin:
    
    # Read array.
    arr_tocheck = np.fromfile(file=fbin, dtype=np.float32).reshape((64, 64, 3, 108000))

    # Check a few frames.
    num_frames = arr_tocheck.shape[3]
    frames_to_play = 10
    for t in np.arange(0, frames_to_play*15, 15):
        plt.imshow(arr_tocheck[:,:,:,t])
        plt.show()
Training data:
In [25]:
print("Testing data:")
fname = "data/vim-2/" + "sv" + "_downsampled.dat"
with open(fname, mode="br") as fbin:
    
    # Read array.
    arr_tocheck = np.fromfile(file=fbin, dtype=np.float32).reshape((64, 64, 3, 8100))

    # Check a few frames.
    num_frames = arr_tocheck.shape[3]
    frames_to_play = 10
    for t in np.arange(0, frames_to_play*15, 15):
        plt.imshow(arr_tocheck[:,:,:,t])
        plt.show()
Testing data: