3.2.1節 複数のディレクトリに分散したデータの集約

1) 正答した試行の反応時間

以下のように真偽値のリストを作成すればよい。

tf_list = (df['target']==target) & (
    df['items']==items) & (
    df['correct']==1)

2) 下限、上限を指定した抽出

通常、xがxmin以上xmax未満という条件は xmin<=x<xmax と書けるが、抽出の際には (xmin<=x) & (x<xmax) と書かねばならない。 従って、以下のように真偽値のリストを作成すればよい。

tf_list = (df['target']==target) & (
    df['items']==items) & (
    0.2<=df['rt']) & (
    df['rt']<3.0)

3.2.2節 ndarrayオブジェクトからの要素の並べ替え

1) 最大値とそのインデックスの抽出

C条件でターゲット無しの平均反応時間は np_data[:,0,3:6] である。インデックス0の次元に沿ってmax()を実行するとアイテム数毎に平均反応時間の最大値を得ることが出来る。従ってアイテム数毎の最大の平均反応時間は np_data[:,0,3:6].max(axis=0) である。

平均反応時間が最大であった参加者のインデックスを求めるには、argsort()を用いればよい。 np_data[:,0,3:6].argsort(axis=0) とすると、アイテム数毎に平均反応時間を昇順に並べる時のインデックスが得られる。このインデックスの最後の要素が平均反応時間が最大となる参加者のインデックスになるので、 np_data[:,0,3:6].argsort(axis=0)[-1] が答えである。

2) 文字列の積を利用したフォーマットの作成

本文中に示した通りだが、改行文字を結合してからformat()メソッドを呼ぶ必要があるため、文字列全体を()で囲む必要がある。

fp.write((('{},'*12)[:-1]+'\n').format(*data))

3) 小数の出力桁数の指定

略

3.2.3節 ndarrayオブジェクトを用いた計

1) 単位の一括変換

データフレームに値を読み込んだ直後に、以下のように反応時間の値を1000倍すればよい。

df = pandas.read_csv(os.path.join(root,file))
df['rt'] *= 1000

2) 対数変換

やはりデータフレームに値を読み込んだ直後に変換を行えばよい。

df = pandas.read_csv(os.path.join(root,file))
df['rt'] = np.log10(1000*df['rt'])

3.3.1節 データ加工の基礎

1) wavファイルへの保存

問題文にあるように、fileOpenDlg()と同様の引数を指定できるfileSaveDlg()を使うとよい。

import psychopy.gui

output_filename = psychopy.gui.fileSaveDlg(
    prompt='wav ファイルを保存します',
    allowed='wav files (*.wav)|*.wav')

保存にはscipy.io.wavfile.write()を用いる。

scipy.io.wavfile.write(output_filename, rate, filtered_data)

2) マーカー付き折れ線グラフの描画

以下に例を示す。

t = np.arange(data.shape[0])/rate
s = int(data.shape[0]/2)
plt.plot(t[s:s+100], data[s:s+100,0], 'bx:')
plt.plot(t[s:s+100], delayed_data[s:s+100,0], 'ro-')
plt.show()

3) 散布図の描画

以下に例を示す。

t = np.arange(data.shape[0])/rate
s = int(data.shape[0]/2)
plt.plot(t[s:s+100], data[s:s+100,0], 'g.')
plt.plot(t[s:s+100], delayed_data[s:s+100,0], 'kx')
plt.show()

3.3.2節 線形フィルタ処理

1) フィルタ効果の確認

略

2) 波形のプロットとファイルへの出力

以下に例を示す。fileSaveDlg()はキャンセルされたときにNoneを返すので、戻り値がNoneでない場合にファイルを保存している。

t = np.arange(data.shape[0])/rate
plt.plot(t, data[:,0], 'b-')
plt.plot(t, filtered_data[:,0], 'r-')
plt.show()

output_filename = psychopy.gui.fileSaveDlg(
    prompt='wavファイルを保存します',
    allowed='wav files (*.wav)|*.wav')
if output_filename is not None:
    scipy.io.wavefile.write(output_filename, rate, filtered_data)

3.3.3節 Fourier変換

1) フィルタの効果の確認

略

2) 4000Hzに最も近い成分の抽出

最初から4000Hz成分のみを抽出すれば簡単だが、問題文に「変数powerから取り出して」とあるのでそれに従うことにする。 変数powerの内容を生成する際に、変数resultsで周波数成分を逆順に並べているので、この点を考慮する必要がある。 以下の例ではreversed_freq_listという変数に周波数成分の一覧を逆順にしたものを用意して、ここから最も4000に近い値のインデックスを探している。 このインデックスさえ求めることが出来れば、あとは難しくないはずである。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import scipy.io.wavfile
 5import numpy as np
 6import matplotlib.pyplot as plt
 7
 8rate, data = scipy.io.wavfile.read('test.wav')
 9N = 1024     # 切り出す要素数 2の累乗数が望ましい
10step = 256   # 切り出し位置の移動ステップ
11duration = 3 # 処理の終了時間(秒)
12s = 0        # 切り出し開始位置を保持する変数
13hamming_window = np.hamming(N) # 窓関数(ハミング窓)
14freq_list = np.fft.fftfreq(N, d=1.0/rate) # 周波数リストを作成
15max_freq_index = int(6000*N/rate) # 6000Hzに対応するインデックス
16results = [] # 結果を保持するリスト
17
18while s+N < rate*duration: # durationを超えて切り出していないか？
19    wd = hamming_window * data[s:s+N,0] # 窓関数との積をとる
20    W = np.fft.fft(wd)     # Fourier変換を実行
21
22    # 6000Hzまでの振幅の対数を計算してresultsに追加
23    results.append(np.log(np.abs(W[max_freq_index:1:-1])))
24    s+=step  # 切り出し開始位置を更新
25
26power = np.array(results).transpose() # 縦軸が周波数になるよう転置
27
28#resultsで周波数成分を逆順にしているので逆順のリストを用意する
29reversed_freq_list = freq_list[max_freq_index:1:-1]
30#4000との差の絶対値が最小となるインデックスを求める
31index4kHz = np.argsort(np.abs(reversed_freq_list-4000))[0]
32#時刻のリストを作る
33t = np.arange(power.shape[1])*(step/rate)
34#プロットする
35plt.plot(t, power[index4kHz,:],'b-')
36plt.show()

変数powerの生成過程自体を変更してよいのであれば、以下のように書くことも出来る。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import scipy.io.wavfile
 5import numpy as np
 6import matplotlib.pyplot as plt
 7
 8rate, data = scipy.io.wavfile.read('test.wav')
 9N = 1024     # 切り出す要素数 2の累乗数が望ましい
10step = 256   # 切り出し位置の移動ステップ
11duration = 3 # 処理の終了時間(秒)
12s = 0        # 切り出し開始位置を保持する変数
13hamming_window = np.hamming(N) # 窓関数(ハミング窓)
14freq_list = np.fft.fftfreq(N, d=1.0/rate) # 周波数リストを作成
15max_freq_index = int(6000*N/rate) # 6000Hzに対応するインデックス
16results = [] # 結果を保持するリスト
17
18while s+N < rate*duration: # durationを超えて切り出していないか？
19    wd = hamming_window * data[s:s+N,0] # 窓関数との積をとる
20    W = np.fft.fft(wd)     # Fourier変換を実行
21
22    # 6000Hzまでの振幅の対数を計算してresultsに追加
23    results.append(np.log(np.abs(W[1:max_freq_index])))
24    s+=step  # 切り出し開始位置を更新
25
26power = np.array(results)
27
28#4000との差の絶対値が最小となるインデックスを求める
29index4kHz = np.argsort(np.abs(freq_list[1:max_freq_index]-4000))[0]
30#時刻のリストを作る
31t = np.arange(power.shape[0])*(step/rate)
32#プロットする
33plt.plot(t, power[:,index4kHz],'b-')
34plt.show()

3.4.2節 PIL を用いた画像処理

以下に例を示す。切り出し範囲の座標計算は少々冗長だが分かり易い方法を用いた。

1) 画像の切り出し

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4from PIL import Image
 5from PIL import ImageFilter
 6import os
 7import sys
 8source_dir = 'image'
 9output_dir = 'output_image'
10output_size = 640
11
12if not os.path.exists(output_dir): # output_imageがない
13    os.mkdir(output_dir)           # output_imageを作る
14elif not os.path.isdir(output_dir): # output_imageがディレクトリでない
15    sys.exit()                      # スクリプトを終了
16
17for file in os.listdir(source_dir):
18    fileroot, ext = os.path.splitext(file)  # 拡張子を分離
19    if ext.lower() in ('.jpg', '.jpeg'): # 拡張子が.jpgか.jpeg
20        image = Image.open(os.path.join(source_dir, file)) # 画像を開く
21
22        if image.size[0] == image.size[1]: # 正方形ならそのまま保存
23            image.save(os.path.join(output_dir, 'c_'+file))
24
25        else: # 長方形の場合、切り出し範囲を計算してcrop
26            if image.size[0] > image.size[1]:
27                diff = image.size[0]-image.size[1]
28                if diff % 2 == 1:  # 2の剰余が1であれば奇数である
29                    left = int((diff-1)/2)
30                    right = image.size[0]-(left+1)
31                else:
32                    left = int(diff/2)
33                    right = image.size[0]-left
34                top = 0
35                bottom = image.size[1]
36            else:
37                diff = image.size[1]-image.size[0]
38                if diff % 2 == 1:
39                    top = int((diff-1)/2)
40                    bottom = image.size[1]-(top+1)
41                else:
42                    top = int(diff/2)
43                    bottom = image.size[1]-top
44                left = 0
45                right = image.size[0]
46
47            image.crop((left, top, right, bottom)).save(
48                os.path.join(output_dir, 'c_'+file))

2) 画像の縦横比を考慮した処理

以下に例を示す。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4from PIL import Image
 5import numpy as np
 6# 画像を開いて透明度チャネルを追加
 7image = Image.open('test.jpg').convert('RGBA')
 8w,h = image.size
 9# Imageオブジェクトからndarrayオブジェクトを作成
10image_np = np.array(image)
11
12# 画像サイズと同じ要素数のndarrayを作成
13if w>h:
14    r = w/h
15    x = np.arange(-2.0*r, 2.0*r, 4.0*r/w)
16    y = np.arange(-2.0, 2.0, 4.0/h)
17else:
18    r = h/w
19    x = np.arange(-2.0, 2.0, 4.0/w)
20    y = np.arange(-2.0*r, 2.0*r, 4.0*r/h)
21
22# 2次元ガウス関数を計算するためのメッシュを作成
23x_mesh, y_mesh = np.meshgrid(x, y)
24
25# 透明度チャネルに値を書き込む
26image_np[:,:,3] = np.exp(-x_mesh**2 -y_mesh**2)*255
27
28# Imageオブジェクトを作成しpngで保存
29Image.fromarray(image_np).save('ouput.png')

3.4.3節 OpenCVを用いた画像処理

1) フォントの変更

OpenCV 2.4.13のドキュメントに記載されている、putText()で使用できるfontFaceは以下のとおりである。コード3.10のFONT_HERSHEY_SIMPLEXを他のものに書き換えればよい。

  • FONT_HERSHEY_SIMPLEX

  • FONT_HERSHEY_PLAIN

  • FONT_HERSHEY_DUPLEX

  • FONT_HERSHEY_COMPLEX

  • FONT_HERSHEY_TRIPLEX

  • FONT_HERSHEY_COMPLEX_SMALL

  • FONT_HERSHEY_SCRIPT_SIMPLEX

  • FONT_HERSHEY_SCRIPT_COMPLEX

2) 検出する顔の大きさの範囲を指定

OpenCV 2.4.13のドキュメントによると、detectMultiScaleの引数は以下のとおりである。

cv2.CascadeClassifier.detectMultiScale(
    image[, scaleFactor[, minNeighbors[, flags[, minSize[, maxSize]]]]]) → objects

Parametersの欄に

  • minSize – Minimum possible object size. Objects smaller than that are ignored.

  • maxSize – Maximum possible object size. Objects larger than that are ignored.

とあるので、これらの引数を指定すればよい。問題はオブジェクトサイズの指定方法である。OpenCVではサイズはSizeというクラスを用いて指定し、幅と高さを指定する必要がある。このことを知っていれば要素数2のシーケンスで指定すればよいことは想像がつく。だが、筆者が確認する限りシーケンスではなんでも良いわけではなくタプルでならないようだ。本文でも述べたようにOpenCVのドキュメントの記述は非常に簡潔なので、こういった情報を読み取ることは難しい。検索を活用して実際にPythonで使用している例を探すと参考になるはずである。

以上より、この問題の解答は以下のようにdetectMultiScale()の引数を変更すればよい。サイズは各自が使用する画像に応じて変更すること。

faces = classifier.detectMultiScale(image,
    maxSize=(250,250), minSize=(100,100))

3.5.1節 USBカメラの活用

1) 左右反転映像の録画

以下に例を示す。reversedという変数を用いずに、write()とsetImage()の引数を直接imageから計算してもよい。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import numpy as np
 5import cv2
 6import psychopy.visual
 7import psychopy.event
 8
 9writer = cv2.VideoWriter('output.m4v',
10    cv2.cv.CV_FOURCC(b'M', b'P', b'4', b'V'), 30, (640,480))
11
12win = psychopy.visual.Window(units='pix', fullscr=False)
13stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
14capture = cv2.VideoCapture(0)
15#if not capture.isOpened():
16#    # 開けなかったときの処理はこのようにif文とisOpened()を使う
17
18while not 'escape' in psychopy.event.getKeys():
19    retval, image = capture.read() # grabとretrieveを一気に
20    if retval:
21        reversed = image[:,::-1,:]    # 動画用は左右反転のみ
22        writer.write(reversed)        # 動画として書き出し
23        stim.setImage(reversed[::-1,:,::-1]/128-1) # PsychoPy用に変更
24        stim.draw()
25        win.flip()
26
27capture.release()
28win.close()

2) 映像のディレイ処理

以下に例を示す。「framesにappend()してからlen(frames)>=11であれば切り詰める」のではなく、「len(frames)==10であれば最も古いものを破棄してからappend()する」方が良いかも知れない。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import numpy as np
 5import cv2
 6import psychopy.visual
 7import psychopy.event
 8
 9win = psychopy.visual.Window(units='pix', fullscr=False)
10stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
11capture = cv2.VideoCapture(0)
12
13frames = []
14
15while not 'escape' in psychopy.event.getKeys():
16    retval, image = capture.read() # grabとretrieveを一気に
17    if retval:
18        frames.append(image)
19        if len(frames) >= 11:
20            frames = frames[-10:]
21        if len(frames) == 10:
22            stim.setImage(frames[0][::-1,:,::-1]/128-1)
23            stim.draw()
24            win.flip()
25
26capture.release()
27win.close()

3) カメラ映像から顔検出

以下に例を示す。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import numpy as np
 5import cv2
 6import psychopy.visual
 7import psychopy.event
 8
 9win = psychopy.visual.Window(units='pix', fullscr=False)
10stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
11capture = cv2.VideoCapture(0)
12
13classifier = cv2.CascadeClassifier('haarcascade_frontalface_default.xml')
14
15while not 'escape' in psychopy.event.getKeys():
16    retval, image = capture.read()
17    if retval:
18        faces = classifier.detectMultiScale(image)
19
20        for rect in faces:
21            cv2.rectangle(image,
22            tuple(rect[0:2]), tuple(rect[0:2]+rect[2:4]), (0,0,255), 2)
23
24        cv2.putText(image, '{} face(s) detected'.format(len(faces)),
25            (10, image.shape[0]-10), cv2.FONT_HERSHEY_SIMPLEX, 1.0,
26            (0, 0, 255))
27
28        stim.setImage(image[::-1,:,::-1]/128-1)
29        stim.draw()
30        win.flip()
31
32capture.release()
33win.close()

4) read()とwrite()の速度

1)の解答を元にした例を示す。時刻を測るにはいろいろな方法があるが、ここではPsychoPyのClockオブジェクトを利用した。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import numpy as np
 5import cv2
 6import psychopy.visual
 7import psychopy.event
 8import psychopy.core
 9import matplotlib.pyplot as plt
10
11writer = cv2.VideoWriter('output.m4v',
12    cv2.cv.CV_FOURCC(b'M', b'P', b'4', b'V'), 30, (640,480))
13
14win = psychopy.visual.Window(units='pix', fullscr=False)
15stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
16capture = cv2.VideoCapture(0)
17
18clock = psychopy.core.Clock()
19read_time = []
20write_time = []
21
22for frames in range(1000):
23    clock.reset()
24    retval, image = capture.read()
25    read_time.append(1000*clock.getTime())
26
27    if retval:
28        reversed = image[:,::-1,:]
29        clock.reset()
30        writer.write(reversed)
31        write_time.append(1000*clock.getTime())
32        stim.setImage(reversed[::-1,:,::-1]/128-1)
33        stim.draw()
34        win.flip()
35
36capture.release()
37win.close()
38
39bins = np.arange(0,32,2)
40plt.subplot(1,2,1)
41plt.hist(read_time, bins=bins, color='blue')
42plt.title('read time')
43plt.ylabel('frequency')
44plt.xlabel('ms')
45plt.subplot(1,2,2)
46plt.hist(write_time, bins=bins, color='red')
47plt.title('write time')
48plt.ylabel('frequency')
49plt.xlabel('ms')
50plt.show()

3.5.2節 データへの非線形当てはめ

1) 残差平方和を用いた処理

以下に例を示す。この例では残差平方和が最も大きい結果と最も小さい結果を最後にプロットしている。 強いて挙げるならば、argmin()、argmax()を使って最大、最小の残差平方和のインデックスを得て、このインデックスを使って対応するパラメータを取り出しているのがポイントである。

 1#coding:utf-8
 2from __future__ import division
 3from __future__ import unicode_literals
 4import numpy as np
 5import matplotlib.pyplot as plt
 6import numpy.random as random
 7import scipy.optimize as optimize  # 最適化の数値計算を行うモジュール
 8
 9rt = 1000/random.normal(5.0,1.5,size=100) # ダミーデータの作成
10
11n, bins = np.histogram(rt, bins=np.arange(100,650,50)) # ヒストグラムを得る
12bin_center = (bins[:-1]+bins[1:])/2 # 各階級の中央の値を計算
13
14def pdf_LATER(t, k, mu, sigma):  # 当てはめに用いる関数を定義する
15    return k/(t**2*np.sqrt(2*np.pi)*sigma)*np.exp(
16        -(1-mu*t)**2/(2*sigma**2*t**2))
17
18residuals_list = []
19params_list = []
20for i in range(1000):
21    k = 100 * random.random()
22    mu = 1 + random.random()
23    sigma = 1 + random.random()
24    params, pcov = optimize.curve_fit(
25        pdf_LATER, bin_center, n, p0=(k, mu, sigma))
26    params_list.append(params)
27    residuals_list.append(np.sum((n - pdf_LATER(bin_center, *params))**2))
28
29i_min = np.argmin(residuals_list)
30i_max = np.argmax(residuals_list)
31
32plt.bar(bins[:-1], n, width=50, color='0.75')
33t = np.arange(1,650)
34plt.plot(t, pdf_LATER(t, *params_list[i_min]), 'k-', linewidth=3)
35plt.plot(t, pdf_LATER(t, *params_list[i_max]), 'k:', linewidth=3)
36plt.text(400, 20, 'MIN k:{}\nmu:{}\nsigma:{}'.format(*params_list[i_min]) )
37plt.text(400, 30, 'MAX k:{}\nmu:{}\nsigma:{}'.format(*params_list[i_max]) )
38plt.show()