.. title:: 「心理学実験プログラミング: Python/PsychoPyによる実験作成・データ処理」 第3章練習問題の解答


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

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

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

.. code-block:: python

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


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

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

.. code-block:: python

    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()メソッドを呼ぶ必要があるため、文字列全体を()で囲む必要がある。

.. code-block:: python

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


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

略


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

1) 単位の一括変換
------------------------------

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

.. code-block:: python

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

2) 対数変換
------------------------------

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

.. code-block:: python

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


3.3.1節 データ加工の基礎
===========================

1) wavファイルへの保存
----------------------------

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

.. code-block:: python

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

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

.. code-block:: python

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


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

以下に例を示す。

.. code-block:: python

    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) 散布図の描画
----------------------------

以下に例を示す。

.. code-block:: python

    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でない場合にファイルを保存している。

.. code-block:: python

    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に近い値のインデックスを探している。
このインデックスさえ求めることが出来れば、あとは難しくないはずである。

.. code-block:: python
    :linenos:
    :emphasize-lines: 28-35

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import scipy.io.wavfile
    import numpy as np
    import matplotlib.pyplot as plt

    rate, data = scipy.io.wavfile.read('test.wav')
    N = 1024     # 切り出す要素数 2の累乗数が望ましい
    step = 256   # 切り出し位置の移動ステップ
    duration = 3 # 処理の終了時間(秒)
    s = 0        # 切り出し開始位置を保持する変数
    hamming_window = np.hamming(N) # 窓関数(ハミング窓)
    freq_list = np.fft.fftfreq(N, d=1.0/rate) # 周波数リストを作成
    max_freq_index = int(6000*N/rate) # 6000Hzに対応するインデックス
    results = [] # 結果を保持するリスト

    while s+N < rate*duration: # durationを超えて切り出していないか？
        wd = hamming_window * data[s:s+N,0] # 窓関数との積をとる
        W = np.fft.fft(wd)     # Fourier変換を実行

        # 6000Hzまでの振幅の対数を計算してresultsに追加
        results.append(np.log(np.abs(W[max_freq_index:1:-1])))
        s+=step  # 切り出し開始位置を更新

    power = np.array(results).transpose() # 縦軸が周波数になるよう転置

    #resultsで周波数成分を逆順にしているので逆順のリストを用意する
    reversed_freq_list = freq_list[max_freq_index:1:-1]
    #4000との差の絶対値が最小となるインデックスを求める
    index4kHz = np.argsort(np.abs(reversed_freq_list-4000))[0]
    #時刻のリストを作る
    t = np.arange(power.shape[1])*(step/rate)
    #プロットする
    plt.plot(t, power[index4kHz,:],'b-')
    plt.show()

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

.. code-block:: python
    :linenos:
    :emphasize-lines: 28-33

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import scipy.io.wavfile
    import numpy as np
    import matplotlib.pyplot as plt

    rate, data = scipy.io.wavfile.read('test.wav')
    N = 1024     # 切り出す要素数 2の累乗数が望ましい
    step = 256   # 切り出し位置の移動ステップ
    duration = 3 # 処理の終了時間(秒)
    s = 0        # 切り出し開始位置を保持する変数
    hamming_window = np.hamming(N) # 窓関数(ハミング窓)
    freq_list = np.fft.fftfreq(N, d=1.0/rate) # 周波数リストを作成
    max_freq_index = int(6000*N/rate) # 6000Hzに対応するインデックス
    results = [] # 結果を保持するリスト

    while s+N < rate*duration: # durationを超えて切り出していないか？
        wd = hamming_window * data[s:s+N,0] # 窓関数との積をとる
        W = np.fft.fft(wd)     # Fourier変換を実行

        # 6000Hzまでの振幅の対数を計算してresultsに追加
        results.append(np.log(np.abs(W[1:max_freq_index])))
        s+=step  # 切り出し開始位置を更新

    power = np.array(results)

    #4000との差の絶対値が最小となるインデックスを求める
    index4kHz = np.argsort(np.abs(freq_list[1:max_freq_index]-4000))[0]
    #時刻のリストを作る
    t = np.arange(power.shape[0])*(step/rate)
    #プロットする
    plt.plot(t, power[:,index4kHz],'b-')
    plt.show()



3.4.2節 PIL を用いた画像処理
=================================

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

1) 画像の切り出し
-----------------------

.. code-block:: python
    :linenos:

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    from PIL import Image
    from PIL import ImageFilter
    import os
    import sys
    source_dir = 'image'
    output_dir = 'output_image'
    output_size = 640

    if not os.path.exists(output_dir): # output_imageがない
        os.mkdir(output_dir)           # output_imageを作る
    elif not os.path.isdir(output_dir): # output_imageがディレクトリでない
        sys.exit()                      # スクリプトを終了

    for file in os.listdir(source_dir):
        fileroot, ext = os.path.splitext(file)  # 拡張子を分離
        if ext.lower() in ('.jpg', '.jpeg'): # 拡張子が.jpgか.jpeg
            image = Image.open(os.path.join(source_dir, file)) # 画像を開く

            if image.size[0] == image.size[1]: # 正方形ならそのまま保存
                image.save(os.path.join(output_dir, 'c_'+file))

            else: # 長方形の場合、切り出し範囲を計算してcrop
                if image.size[0] > image.size[1]:
                    diff = image.size[0]-image.size[1]
                    if diff % 2 == 1:  # 2の剰余が1であれば奇数である
                        left = int((diff-1)/2)
                        right = image.size[0]-(left+1)
                    else:
                        left = int(diff/2)
                        right = image.size[0]-left
                    top = 0
                    bottom = image.size[1]
                else:
                    diff = image.size[1]-image.size[0]
                    if diff % 2 == 1:
                        top = int((diff-1)/2)
                        bottom = image.size[1]-(top+1)
                    else:
                        top = int(diff/2)
                        bottom = image.size[1]-top
                    left = 0
                    right = image.size[0]
                
                image.crop((left, top, right, bottom)).save(
                    os.path.join(output_dir, 'c_'+file))


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

以下に例を示す。

.. code-block:: python
    :linenos:
    :emphasize-lines: 12-20

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    from PIL import Image
    import numpy as np
    # 画像を開いて透明度チャネルを追加
    image = Image.open('test.jpg').convert('RGBA')
    w,h = image.size
    # Imageオブジェクトからndarrayオブジェクトを作成
    image_np = np.array(image)

    # 画像サイズと同じ要素数のndarrayを作成
    if w>h:
        r = w/h
        x = np.arange(-2.0*r, 2.0*r, 4.0*r/w)
        y = np.arange(-2.0, 2.0, 4.0/h)
    else:
        r = h/w
        x = np.arange(-2.0, 2.0, 4.0/w)
        y = np.arange(-2.0*r, 2.0*r, 4.0*r/h)
        
    # 2次元ガウス関数を計算するためのメッシュを作成
    x_mesh, y_mesh = np.meshgrid(x, y)

    # 透明度チャネルに値を書き込む
    image_np[:,:,3] = np.exp(-x_mesh**2 -y_mesh**2)*255

    # Imageオブジェクトを作成しpngで保存
    Image.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の引数は以下のとおりである。

.. code-block:: none

    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()の引数を変更すればよい。サイズは各自が使用する画像に応じて変更すること。

.. code-block:: python

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


3.5.1節 USBカメラの活用
===========================

1) 左右反転映像の録画
---------------------------------

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

.. code-block:: python
    :linenos:
    :emphasize-lines: 9-10,21-23

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import numpy as np
    import cv2
    import psychopy.visual
    import psychopy.event

    writer = cv2.VideoWriter('output.m4v',
        cv2.cv.CV_FOURCC(b'M', b'P', b'4', b'V'), 30, (640,480))

    win = psychopy.visual.Window(units='pix', fullscr=False)
    stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
    capture = cv2.VideoCapture(0)
    #if not capture.isOpened():
    #    # 開けなかったときの処理はこのようにif文とisOpened()を使う

    while not 'escape' in psychopy.event.getKeys():
        retval, image = capture.read() # grabとretrieveを一気に
        if retval:
            reversed = image[:,::-1,:]    # 動画用は左右反転のみ
            writer.write(reversed)        # 動画として書き出し
            stim.setImage(reversed[::-1,:,::-1]/128-1) # PsychoPy用に変更
            stim.draw()
            win.flip()

    capture.release()
    win.close()


2) 映像のディレイ処理
-------------------------

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

.. code-block:: python
    :linenos:
    :emphasize-lines: 13,18-24

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import numpy as np
    import cv2
    import psychopy.visual
    import psychopy.event

    win = psychopy.visual.Window(units='pix', fullscr=False)
    stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
    capture = cv2.VideoCapture(0)

    frames = []

    while not 'escape' in psychopy.event.getKeys():
        retval, image = capture.read() # grabとretrieveを一気に
        if retval:
            frames.append(image)
            if len(frames) >= 11:
                frames = frames[-10:]
            if len(frames) == 10:
                stim.setImage(frames[0][::-1,:,::-1]/128-1)
                stim.draw()
                win.flip()

    capture.release()
    win.close()



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

以下に例を示す。

.. code-block:: python
    :linenos:

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import numpy as np
    import cv2
    import psychopy.visual
    import psychopy.event

    win = psychopy.visual.Window(units='pix', fullscr=False)
    stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
    capture = cv2.VideoCapture(0)

    classifier = cv2.CascadeClassifier('haarcascade_frontalface_default.xml')

    while not 'escape' in psychopy.event.getKeys():
        retval, image = capture.read()
        if retval:
            faces = classifier.detectMultiScale(image)

            for rect in faces:
                cv2.rectangle(image,
                tuple(rect[0:2]), tuple(rect[0:2]+rect[2:4]), (0,0,255), 2)

            cv2.putText(image, '{} face(s) detected'.format(len(faces)),
                (10, image.shape[0]-10), cv2.FONT_HERSHEY_SIMPLEX, 1.0,
                (0, 0, 255))

            stim.setImage(image[::-1,:,::-1]/128-1)
            stim.draw()
            win.flip()

    capture.release()
    win.close()

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

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

.. code-block:: python
    :linenos:
    :emphasize-lines: 8-9,18-20,23,25,29,31,39-50

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import numpy as np
    import cv2
    import psychopy.visual
    import psychopy.event
    import psychopy.core
    import matplotlib.pyplot as plt

    writer = cv2.VideoWriter('output.m4v',
        cv2.cv.CV_FOURCC(b'M', b'P', b'4', b'V'), 30, (640,480))

    win = psychopy.visual.Window(units='pix', fullscr=False)
    stim = psychopy.visual.ImageStim(win, size=(640,480), contrast=0.3)
    capture = cv2.VideoCapture(0)

    clock = psychopy.core.Clock()
    read_time = []
    write_time = []

    for frames in range(1000):
        clock.reset()
        retval, image = capture.read()
        read_time.append(1000*clock.getTime())

        if retval:
            reversed = image[:,::-1,:]
            clock.reset()
            writer.write(reversed)
            write_time.append(1000*clock.getTime())
            stim.setImage(reversed[::-1,:,::-1]/128-1)
            stim.draw()
            win.flip()

    capture.release()
    win.close()

    bins = np.arange(0,32,2)
    plt.subplot(1,2,1)
    plt.hist(read_time, bins=bins, color='blue')
    plt.title('read time')
    plt.ylabel('frequency')
    plt.xlabel('ms')
    plt.subplot(1,2,2)
    plt.hist(write_time, bins=bins, color='red')
    plt.title('write time')
    plt.ylabel('frequency')
    plt.xlabel('ms')
    plt.show()


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

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

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

.. code-block:: python
    :linenos:

    #coding:utf-8
    from __future__ import division
    from __future__ import unicode_literals
    import numpy as np
    import matplotlib.pyplot as plt
    import numpy.random as random
    import scipy.optimize as optimize  # 最適化の数値計算を行うモジュール

    rt = 1000/random.normal(5.0,1.5,size=100) # ダミーデータの作成

    n, bins = np.histogram(rt, bins=np.arange(100,650,50)) # ヒストグラムを得る
    bin_center = (bins[:-1]+bins[1:])/2 # 各階級の中央の値を計算

    def pdf_LATER(t, k, mu, sigma):  # 当てはめに用いる関数を定義する
        return k/(t**2*np.sqrt(2*np.pi)*sigma)*np.exp(
            -(1-mu*t)**2/(2*sigma**2*t**2))

    residuals_list = []
    params_list = []
    for i in range(1000):
        k = 100 * random.random()
        mu = 1 + random.random()
        sigma = 1 + random.random()
        params, pcov = optimize.curve_fit(
            pdf_LATER, bin_center, n, p0=(k, mu, sigma))
        params_list.append(params)
        residuals_list.append(np.sum((n - pdf_LATER(bin_center, *params))**2))

    i_min = np.argmin(residuals_list)
    i_max = np.argmax(residuals_list)

    plt.bar(bins[:-1], n, width=50, color='0.75')
    t = np.arange(1,650)
    plt.plot(t, pdf_LATER(t, *params_list[i_min]), 'k-', linewidth=3)
    plt.plot(t, pdf_LATER(t, *params_list[i_max]), 'k:', linewidth=3)
    plt.text(400, 20, 'MIN k:{}\nmu:{}\nsigma:{}'.format(*params_list[i_min]) )
    plt.text(400, 30, 'MAX k:{}\nmu:{}\nsigma:{}'.format(*params_list[i_max]) )
    plt.show()