旅行好きなソフトエンジニアの備忘録

プログラミングや技術関連のメモを始めました

【アイトラッキング】 I-DTアルゴリズムの実装

Salvucci & Goldberg著 "Identifying Fixations and Saccades in Eye-Tracking Protocols"に記述されているI-DTアルゴリズムの実装メモになります。

まずは視線データを保存するためのクラスを定義します。

public class GazePoint
{
    public double X; // 視線X座標
    public double Y; // 視線Y座標
    public long TimeStamp; // タイムスタンプ
    public int? FixationId; // 固視ID
    public uint FixationDuration_msec; // 固視時間

    public static double operator- (GazePoint p1, GazePoint p2)
    {
        return Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));
    }

    public GazePoint Clone()
    {
        return new GazePoint() { X = X, Y = Y, TimeStamp = TimeStamp, FixationId = FixationId, FixationDuration_msec = FixationDuration_msec };
    }
}


次にI-DTアルゴリズムをIDTAlgorithmクラスとして実装します。AnalyzeメソッドがI-DTアルゴリズムを適用している箇所であり、論文に書かれている疑似コードと対比できるよう、論文から疑似コードのコメントを引用しています。

public class IDTAlgorithm
{
    private readonly double _dispersionThreshold;
    private readonly uint _durationThreshold_msec;

    public IDTAlgorithm(double dispersionThreshold, uint durationThreshold_msec)
    {
        _dispersionThreshold = dispersionThreshold;
        _durationThreshold_msec = durationThreshold_msec;
    }

    /// <summary>
    /// IDTアルゴリズムで視線データを分析し、固視/それ以外を判別するメソッド
    /// </summary>
    /// <param name="points">視線データ</param>
    /// <returns>固視フラグが付加された視線データ</returns>
    public GazePoint[] Analyze(IList<GazePoint> points)
    {
        GazePoint[] filteredPoints = Enumerable.Repeat<GazePoint>(null, points.Count).ToArray();

        int fixationId = 0;
        for (int i = 0; i < points.Count; ++i)
        {
            // Initialize window over first points to cover the duration threshold
            int endIndex;
            List<GazePoint> window = InitializeWindow(points, i, _durationThreshold_msec, out endIndex);

            // If dispersion of window points <= threshold
            double dispersion = ComputeDispersion(window);
            if (dispersion <= _dispersionThreshold)
            {
                // Add additional points to the window until dispersion > threshold
                AddAdditionalPoints(points, _dispersionThreshold, ref window, ref endIndex);

                // Note a fixation at the centroid of the window points
                GazePoint[] filteredFixations = ComputeAverageFixation(window);
                ++fixationId;
                uint fixationDuration_msec = (uint)(window.Last().TimeStamp - window.First().TimeStamp);
                for (int j = 0; j < filteredFixations.Length; ++j)
                {
                    filteredPoints[i + j] = filteredFixations[j];
                    filteredPoints[i + j].FixationId = fixationId;
                    filteredPoints[i + j].FixationDuration_msec = fixationDuration_msec;
                }

                // Remove window points from points
                i = endIndex;
            }
        }

        // 固視以外の視線データをセットする
        for (int i = 0; i < points.Count; ++i)
        {
            if (filteredPoints[i] == null)
            {
                filteredPoints[i] = points[i];
                filteredPoints[i].FixationId = -1;
            }
        }

        return filteredPoints;
    }

    private List<GazePoint> InitializeWindow(IList<GazePoint> points, int startIndex, uint durationThreshold_msec, out int endIndex)
    {
        var window = new List<GazePoint>();
        long start = points[startIndex].TimeStamp;
        long end = start + durationThreshold_msec;
        int index = 0;
        while ((startIndex + index < points.Count) && (points[startIndex + index].TimeStamp <= end))
        {
            window.Add(points[startIndex + index]);
            ++index;
        }
        endIndex = startIndex + index - 1;
        return window;
    }

    private void AddAdditionalPoints(IList<GazePoint> points, double dispersionThreshold, ref List<GazePoint> window, ref int endIndex)
    {
        double dispersion;
        int i = 0;
        bool outOfIndex = false;
        do
        {
            ++i;
            if (endIndex + i >= points.Count)
            {
                outOfIndex = true;
                break;
            }
            window.Add(points[endIndex + i]);
            dispersion = ComputeDispersion(window);
        } while (dispersion <= dispersionThreshold);
        if (!outOfIndex) window.RemoveAt(window.Count - 1);
        endIndex += i - 1;
    }

    private double ComputeDispersion(List<GazePoint> window)
    {
        IEnumerable<double> xs = window.Select(p => p.X);
        IEnumerable<double> ys = window.Select(p => p.Y);
        return xs.Max() - xs.Min() + ys.Max() - ys.Min();
    }

    private GazePoint[] ComputeAverageFixation(IEnumerable<GazePoint> fixations)
    {
        double x = fixations.Average(i => i.X);
        double y = fixations.Average(i => i.Y);
        var filteredFixations = new GazePoint[fixations.Count()];
        for (int i = 0; i < filteredFixations.Length; ++i)
        {
            filteredFixations[i] = new GazePoint
            {
                X = x,
                Y = y,
                TimeStamp = fixations.ElementAt(i).TimeStamp
            };
        }
        return filteredFixations;
    }
}


最後にIDTAlgorithmクラスの使い方を記述します。

var points = new List<GazePoint>();
points.Add(new GazePoint() { X = 630.0, Y = 0.0, TimeStamp = 0 });
points.Add(new GazePoint() { X = 635.0, Y = 0.0, TimeStamp = 33 });
points.Add(new GazePoint() { X = 645.0, Y = 0.0, TimeStamp = 66 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 100 });
points.Add(new GazePoint() { X = 635.0, Y = 0.0, TimeStamp = 133 });
points.Add(new GazePoint() { X = 640.0, Y = 0.0, TimeStamp = 166 });
points.Add(new GazePoint() { X = 635.0, Y = 0.0, TimeStamp = 200 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 233 });
points.Add(new GazePoint() { X = 640.0, Y = 0.0, TimeStamp = 266 });
points.Add(new GazePoint() { X = 638.0, Y = 0.0, TimeStamp = 300 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 333 });
points.Add(new GazePoint() { X = 640.0, Y = 0.0, TimeStamp = 366 });
points.Add(new GazePoint() { X = 660.0, Y = 0.0, TimeStamp = 400 });
points.Add(new GazePoint() { X = 645.0, Y = 0.0, TimeStamp = 433 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 466 });
points.Add(new GazePoint() { X = 645.0, Y = 0.0, TimeStamp = 500 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 533 });
points.Add(new GazePoint() { X = 640.0, Y = 0.0, TimeStamp = 566 });
points.Add(new GazePoint() { X = 640.0, Y = 0.0, TimeStamp = 600 });
points.Add(new GazePoint() { X = 750.0, Y = 0.0, TimeStamp = 633 });
points.Add(new GazePoint() { X = 760.0, Y = 0.0, TimeStamp = 666 });
points.Add(new GazePoint() { X = 740.0, Y = 0.0, TimeStamp = 700 });
points.Add(new GazePoint() { X = 745.0, Y = 0.0, TimeStamp = 733 });
points.Add(new GazePoint() { X = 745.0, Y = 0.0, TimeStamp = 766 });
points.Add(new GazePoint() { X = 750.0, Y = 0.0, TimeStamp = 800 });
points.Add(new GazePoint() { X = 745.0, Y = 0.0, TimeStamp = 833 });
points.Add(new GazePoint() { X = 750.0, Y = 0.0, TimeStamp = 866 });
points.Add(new GazePoint() { X = 745.0, Y = 0.0, TimeStamp = 900 });
points.Add(new GazePoint() { X = 710.0, Y = 0.0, TimeStamp = 933 });
points.Add(new GazePoint() { X = 720.0, Y = 0.0, TimeStamp = 966 });
points.Add(new GazePoint() { X = 715.0, Y = 0.0, TimeStamp = 1000 });
points.Add(new GazePoint() { X = 600.0, Y = 0.0, TimeStamp = 1033 });
points.Add(new GazePoint() { X = 730.0, Y = 0.0, TimeStamp = 1066 });
points.Add(new GazePoint() { X = 740.0, Y = 0.0, TimeStamp = 1100 });
points.Add(new GazePoint() { X = 730.0, Y = 0.0, TimeStamp = 1133 });
points.Add(new GazePoint() { X = 735.0, Y = 0.0, TimeStamp = 1166 });
points.Add(new GazePoint() { X = 730.0, Y = 0.0, TimeStamp = 1200 });
points.Add(new GazePoint() { X = 735.0, Y = 0.0, TimeStamp = 1233 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 1266 });
points.Add(new GazePoint() { X = 645.0, Y = 0.0, TimeStamp = 1300 });
points.Add(new GazePoint() { X = 650.0, Y = 0.0, TimeStamp = 1333 });
points.Add(new GazePoint() { X = 640.0, Y = 0.0, TimeStamp = 1366 });
points.Add(new GazePoint() { X = 635.0, Y = 0.0, TimeStamp = 1400 });

const double dispersionThreshold = 30.0;
// 固視時間閾値(通常100~200msecあたりを指定)
const uint durationThreshold_msec = 150;
var algorithm = new IDTAlgorithm(dispersionThreshold, durationThreshold_msec);
GazePoint[] analyzedPoints = algorithm.Analyze(points);


Analyzeメソッドの出力をグラフ化して元のデータと比較すると以下のようになります。GazePointクラスのFixationIdが1から4まで割り当てられているため、固視回数は4回と分かります(グラフからも明らか)。 f:id:ni4muraano:20170628211212p:plain

【料理】 ゴーヤチャンプルー

ここを元に作成。 www.sirogohan.com

材料(二人分)

  • ゴーヤ   1/2本

  • 木綿豆腐  300g

  • バラ肉  100g

  • 卵     1個

  • 塩コショウ 適量

  • サラダ油/ごま油

  • 醤油

作り方

1. 木綿豆腐をキッチンペーパーでくるんで水切りする
2. 卵1個をといておく
3. ゴーヤを縦半分に切る
4. スプーンを使って種を取り除く
5. 4~5mm幅に切る
6. 塩小さじ1/4をまぶし、軽く混ぜ合わせる
7. 豚バラ肉を3~4cm幅に切り、塩をうっすら、粗挽き胡椒を多めにふる
8. フライパンにサラダ油をひき、温まったら豆腐を入れる
9. 豆腐を取り出し、ごま油小さじ1を追加し、ゴーヤを炒める
10. 豚バラ肉を入れて炒める
11. 塩ひとつまみ強、胡椒少々で味付けし、溶き卵を流し入れる
12. 10秒待ってフライパンを振る。最後に醤油を小さじ2回しかける

【数学】 フィッシャーの線形判別を理解するのに有効なサイトのメモ

フィッシャーの線形判別について調べていたのですが、分かりやすく書いてくれているサイトを見つけたのでメモします。

理論についてはこちらが分かりやすいです。

s0sem0y.hatenablog.com

コードが必要であればこちらにコード付きで記事を書かれています。

aidiary.hatenablog.com

【Python】 LSTMによる時系列データの予測

前回SimpleRNNによる時系列データの予測を行いましたが、今回はLSTMを用いて時系列データの予測を行ってみます。

ni4muraano.hatenablog.com

LSTMはSimpleRNNと比較すると長期依存性の高いデータに有効とのことなので、50回に一回パルスが発生する信号に対する予測をSimpleRNNとLSTMで行ってみました。

import numpy as np
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.layers.recurrent import SimpleRNN, LSTM
from keras.optimizers import Adam
from keras.callbacks import EarlyStopping

np.random.seed(0)

def generate_periodic_data():
    pulse = []
    for i in range(1000):
        if i%25 == 0:
            pulse.append(1.0)
        else:
            pulse.append(0.0)

    return np.array(pulse)


def generate_data(pulse, length_per_unit, dimension):
    sequences = []
    target = []
    for i in range(0, pulse.size - length_per_unit):
        sequences.append(pulse[i:i + length_per_unit])
        target.append(pulse[i + length_per_unit])

    X = np.array(sequences).reshape(len(sequences), length_per_unit, dimension)
    Y = np.array(target).reshape(len(sequences), dimension)

    N_train = int(len(sequences)*0.8)
    X_train = X[:N_train]
    X_validation = X[N_train:]
    Y_train = Y[:N_train]
    Y_validation = Y[N_train:]

    return (X_train, X_validation, Y_train, Y_validation)


def build_model(input_shape, hidden_layer_count):
    model = Sequential()
    model.add(LSTM(hidden_layer_count, input_shape=input_shape))
    model.add(Dense(input_shape[1]))
    model.add(Activation('linear'))
    model.compile(loss='mse', optimizer=Adam())
    return model


# 一つの時系列データの長さ
LENGTH_PER_UNIT = 201
# 一次元データを扱う
DIMENSION = 1
# パルスの生成
pulse = generate_periodic_data()
# トレーニング、バリデーション用データの生成
X_train, X_validation, Y_train, Y_validation = generate_data(pulse, LENGTH_PER_UNIT, DIMENSION)

# LSTM隠れ層の数
HIDDEN_LAYER_COUNT = 25
# 入力の形状
input_shape=(LENGTH_PER_UNIT, DIMENSION)
# モデルの生成
model = build_model(input_shape, HIDDEN_LAYER_COUNT)

# モデルのトレーニング
epochs = 500
batch_size = 10
early_stopping = EarlyStopping(monitor='val_loss', patience=10, verbose=1)
model.fit(X_train, Y_train,
          batch_size=batch_size,
          epochs=epochs,
          validation_data=(X_validation, Y_validation),
          callbacks=[early_stopping])

# 予測を行う
part_of_sequence = np.array([pulse[i] for i in range(LENGTH_PER_UNIT)])
predicted = [None for i in range(LENGTH_PER_UNIT)]
Z = X_train[:1, :, :]
for i in range(pulse.size - LENGTH_PER_UNIT + 1):
    y_ = model.predict(Z)
    # 予測結果を入力として利用するため、第0項を削除し予測結果をひっつける
    Z = np.concatenate(
            (Z.reshape(LENGTH_PER_UNIT, DIMENSION)[1:], y_),
            axis=0).reshape(1, LENGTH_PER_UNIT, DIMENSION)
    predicted.append(y_.reshape(-1))
predicted = np.array(predicted)

# 予測結果の描画
plt.rc('font', family='serif')
plt.figure()
plt.ylim([0.0, 1.1])
plt.plot(pulse, linestyle='dotted', color='#aaaaaa')
plt.plot(part_of_sequence, linestyle='dashed', color='black')
plt.plot(predicted, color='black')
plt.show()

SimpleRNNの予測結果がこちらでf:id:ni4muraano:20170617162126p:plain

LSTMの予測結果がこうなりました。f:id:ni4muraano:20170617162147p:plain

LSTMはもう少し綺麗に予測できることを期待していたのですが、パラメータチューニングを行っていないので上図のような結果になったのかなと思います。とはいえ閾値決めてあげることで後処理でどうにでもなりそうなのでSimpleRNNよりは良い結果なのではと思います。

詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

【Python】 SimpleRNNで月平均気温を予測する

画像だけでなく時系列データにも手を出してみたい、ということで書籍「詳解ディープラーニング」を購入しました。書籍第5章から時系列データを扱っているのですが、そこで紹介されているSimpleRNNの例を写経します。書籍ではノイズの入ったサイン波の予測を行っているのですが、ここでは気象庁から取得した年別の月平均気温を使って予測を行います。

www.data.jma.go.jp


import numpy as np
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.layers.recurrent import SimpleRNN
from keras.optimizers import Adam
from keras.callbacks import EarlyStopping

np.random.seed(0)

def generate_temperatures():
    # 月平均気温
    temperatures = [[7.6, 6.0, 9.4, 14.5, 19.8, 22.5, 27.7, 28.3, 25.6, 18.8, 13.3, 8.8], #2000年
                    [4.9, 6.6, 9.8, 15.7, 19.5, 23.1, 28.5, 26.4, 23.2, 18.7, 13.1, 8.4], #2001年
                    [7.4, 7.9, 12.2, 16.1, 18.4, 21.6, 28.0, 28.0, 23.1, 19.0, 11.6, 7.2],#2002年
                    [5.5, 6.4, 8.7, 15.1, 18.8, 23.2, 22.8, 26.0, 24.2, 17.8, 14.4, 9.2], #2003年
                    [6.3, 8.5, 9.8, 16.4, 19.6, 23.7, 28.5, 27.2, 25.1, 17.5, 15.6, 9.9], #2004年
                    [6.1, 6.2, 9.0, 15.1, 17.7, 23.2, 25.6, 28.1, 24.7, 19.2, 13.3, 6.4],
                    [5.1, 6.7, 9.8, 13.6, 19.0, 22.5, 25.6, 27.5, 23.5, 19.5, 14.4, 9.5],
                    [7.6, 8.6, 10.8, 13.7, 19.8, 23.2, 24.4, 29.0, 25.2, 19.0, 13.3, 9.0],
                    [5.9, 5.5, 10.7, 14.7, 18.5, 21.3, 27.0, 26.8, 24.4, 19.4, 13.1, 9.8],
                    [6.8, 7.8, 10.0, 15.7, 20.1, 22.5, 26.3, 26.6, 23.0, 19.0, 13.5, 9.0],
                    [7.0, 6.5, 9.1, 12.4, 19.0, 23.6, 28.0, 29.6, 25.1, 18.9, 13.5, 9.9],
                    [5.1, 7.0, 8.1, 14.5, 18.5, 22.8, 27.3, 27.5, 25.1, 19.5, 14.9, 7.5],
                    [4.8, 5.4, 8.8, 14.5, 19.6, 21.4, 26.4, 29.1, 26.2, 19.4, 12.7, 7.3],
                    [5.5, 6.2, 12.1, 15.2, 19.8, 22.9, 27.3, 29.2, 25.2, 19.8, 13.5, 8.3],
                    [6.3, 5.9, 10.4, 15.0, 20.3, 23.4, 26.8, 27.7, 23.2, 19.1, 14.2, 6.7],
                    [5.8, 5.7, 10.3, 14.5, 21.1, 22.1, 26.2, 26.7, 22.6, 18.4, 13.9, 9.3],
                    [6.1, 7.2, 10.1, 15.4, 20.2, 22.4, 25.4, 27.1, 24.4, 18.7, 11.4, 8.9]]#2016年
    temperatures = np.array(temperatures)
    temperatures = np.reshape(temperatures, (temperatures.size))
    return temperatures


def generate_data(temperatures, length_per_unit, dimension):
    sequences = []
    target = []
    for i in range(0, temperatures.size - length_per_unit):
        sequences.append(temperatures[i:i + length_per_unit])
        target.append(temperatures[i + length_per_unit])

    X = np.array(sequences).reshape(len(sequences), length_per_unit, dimension)
    Y = np.array(target).reshape(len(sequences), dimension)

    N_train = int(len(sequences) * 0.9)
    X_train = X[:N_train]
    X_validation = X[N_train:]
    Y_train = Y[:N_train]
    Y_validation = Y[N_train:]

    return (X_train, X_validation, Y_train, Y_validation)


def build_model(input_shape, hidden_layer_count):
    model = Sequential()
    model.add(SimpleRNN(hidden_layer_count, input_shape=input_shape))
    model.add(Dense(input_shape[1]))
    model.add(Activation('linear'))
    model.compile(loss='mse', optimizer=Adam())
    return model


# 一つの時系列データの長さ
LENGTH_PER_UNIT = 24
# 一次元データを扱う
DIMENSION = 1
# 年別月平均気温の生成
temperatures = generate_temperatures()
# トレーニング、バリデーション用データの生成
X_train, X_validation, Y_train, Y_validation = generate_data(temperatures, LENGTH_PER_UNIT, DIMENSION)

# SimpleRNN隠れ層の数
HIDDEN_LAYER_COUNT = 25
# 入力の形状
input_shape=(LENGTH_PER_UNIT, DIMENSION)
# モデルの生成
model = build_model(input_shape, HIDDEN_LAYER_COUNT)

# モデルのトレーニング
epochs = 500
batch_size = 10
early_stopping = EarlyStopping(monitor='val_loss', patience=10, verbose=1)
model.fit(X_train, Y_train,
          batch_size=batch_size,
          epochs=epochs,
          validation_data=(X_validation, Y_validation),
          callbacks=[early_stopping])

# 予測を行う
part_of_sequence = np.array([temperatures[i] for i in range(LENGTH_PER_UNIT)])
predicted = [None for i in range(LENGTH_PER_UNIT)]
Z = X_train[:1, :, :]
for i in range(temperatures.size - LENGTH_PER_UNIT + 1):
    y_ = model.predict(Z)
    # 予測結果を入力として利用するため、第0項を削除し予測結果をひっつける
    Z = np.concatenate(
            (Z.reshape(LENGTH_PER_UNIT, DIMENSION)[1:], y_),
            axis=0).reshape(1, LENGTH_PER_UNIT, DIMENSION)
    predicted.append(y_.reshape(-1))
predicted = np.array(predicted)

# 予測結果の描画
plt.rc('font', family='serif')
plt.figure()
plt.ylim([0.0, 35.0])
plt.plot(temperatures, linestyle='dotted', color='#aaaaaa')
plt.plot(part_of_sequence, linestyle='dashed', color='black')
plt.plot(predicted, color='black')
plt.show()


割と無難な結果を得ることが出来ました。 f:id:ni4muraano:20170615213133p:plain

詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

【料理】 たらこパスタの実装

下記を参考に作りました。 cookpad.com

用意するもの(二名分)

  • たらこ   2腹
  • バター   30g
  • 牛乳    大さじ4
  • 醤油    大さじ1
  • きざみのり 適量
  • パスタ   180g

1 お湯を沸かす

2 ボールにたらこ、バター、牛乳、醤油を入れてソースを作る

3 パスタを茹でる

4 お湯を捨ててそこにソースを入れて良く混ぜる

5 容器に移してきざみのりをかける

【アイトラッキング】 オンラインスムージングフィルタの実装

アイトラッカーから取得したデータにオンラインでスムージングフィルタをかけようとしています。Manu Kumar著 “GUIDe SACCADE DETECTION AND SMOOTHING ALGORITHM"に書いてある疑似コードを参考に実装しました。プログラムの注意点としては論文の疑似コードが分かりにくいため間違っているかもしれないこと、スレッドセーフでないためサンプリングレートの高いアイトラッカーにコードを適用した場合意図した通りに動作しないことです。

まずはアイトラッカーからのデータを格納するGazePointクラスの実装です。-演算子オーバーロードし、二点間の距離を計算できるようにしています。

public class GazePoint
{
    public double X; // X座標
    public double Y; // Y座標
    public long TimeStamp; // タイムスタンプ

    public static double operator- (GazePoint p1, GazePoint p2)
    {
        return Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));
    }

    public GazePoint Clone()
    {
        return new GazePoint() { X = X, Y = Y, TimeStamp = TimeStamp };
    }
}


次にオンラインスムージングフィルタを行うクラスをOnlineGazeFilterクラスとして実装します(2017/06/14追記:OnlineGazeFilterクラスにバグがあったため修正)。

public class OnlineGazeFilter
{
    private readonly uint _windowSize;
    private readonly double _saccadeThreshold;
    private readonly long _timeout;
    private List<GazePoint> _points = new List<GazePoint>();
    private GazePoint _currentFixation = null;
    private GazePoint _potentialFixation = null;
    private GazePoint _output = null;

    public OnlineGazeFilter(uint windowSize, double saccadeThreshold, long timeout)
    {
        _windowSize = windowSize;
        _saccadeThreshold = saccadeThreshold;
        _timeout = timeout;
    }

    public void AddPoint(GazePoint point)
    {
        // 新しいデータが前回のデータから大分時間が経っている場合は古いデータを消去する
        Reset(point.TimeStamp);
        if (_currentFixation == null)
        {
            AddPoint(point, _windowSize, ref _points);
            _points.Add(point);
            _currentFixation = ComputeFixation(_points);
            return;
        }
        if (_potentialFixation != null)
        {
            double distanceFromCurrent = point - _currentFixation;
            double distanceFromPotential = point - _potentialFixation;
            if (distanceFromCurrent <= distanceFromPotential)
            {
                if (distanceFromCurrent < _saccadeThreshold)
                {
                    // 現在の固視位置を固視し続けている
                    AddPoint(point, _windowSize, ref _points);
                    _potentialFixation = null;
                    _currentFixation = ComputeFixation(_points);
                    _output = _currentFixation.Clone();
                }
                else
                {
                    // サッカード、もしくは固視位置の変更の可能性有り
                    _potentialFixation = point;
                    _output = _currentFixation.Clone();
                }
            }
            else
            {
                if (distanceFromPotential < _saccadeThreshold)
                {
                    // 固視位置が変わった
                    _potentialFixation = point;
                    _points.Clear();
                    AddPoint(point, _windowSize, ref _points);
                    _currentFixation = _potentialFixation;
                    _output = _currentFixation.Clone();
                }
                else
                {
                    // サッカード中
                    GazePoint savedPotentialFixation = _potentialFixation.Clone();
                    _potentialFixation = point;
                    _output = savedPotentialFixation;
                }
            }
        }
        else
        {
            double distance = point - _currentFixation;
            if (distance > _saccadeThreshold)
            {
                // サッカード、もしくは固視位置の変更の可能性有り
                _potentialFixation = point;
                _output = _currentFixation.Clone();
            }
            else
            {
                // 現在の固視位置を固視し続けている
                AddPoint(point, _windowSize, ref _points);
                _currentFixation = ComputeFixation(_points);
                _output = _currentFixation.Clone();
            }
        }
        if (_points.Count > _windowSize)
        {
            _points.RemoveAt(0);
        }
    }

    private void Reset(long newTimeStamp)
    {
        if (_currentFixation != null)
        {
            long lastTimeStamp;
            if (_potentialFixation != null)
            {
                lastTimeStamp = Math.Max(_currentFixation.TimeStamp, _potentialFixation.TimeStamp);
            }
            else
            {
                lastTimeStamp = _currentFixation.TimeStamp;
            }
            if (newTimeStamp - lastTimeStamp > _timeout)
            {
                _currentFixation = null;
                _potentialFixation = null;
                _points.Clear();
            }
        }
    }

    private void AddPoint(GazePoint point, uint windowSize, ref List<GazePoint> points)
    {
        points.Add(point);
        if (points.Count > windowSize)
        {
            points.RemoveAt(0);
        }
    }

    /// <summary>
    /// スムージングフィルタをかけるメソッド
    /// </summary>
    /// <param name="fixations">視線データ群</param>
    /// <returns>フィルタ後の視線位置</returns>
    private GazePoint ComputeFixation(List<GazePoint> fixations)
    {
        double sumX = 0.0;
        double sumY = 0.0;
        double denominator = 0.0;
        for (int i = 0; i < fixations.Count; ++i)
        {
            denominator += (i + 1);
            sumX += (i + 1)*fixations[i].X;
            sumY += (i + 1)*fixations[i].Y;
        }
        return new GazePoint { X = sumX / denominator, Y = sumY / denominator, TimeStamp = fixations.Last().TimeStamp };
    }

    /// <summary>
    /// フィルタ後の視線位置を取得する
    /// </summary>
    /// <returns>フィルタ後の視線位置</returns>
    public GazePoint GetCurrentPoint()
    {
        return _output;
    }
}


最後にOnlineGazeFilterクラスの使い方を示します。

class Program
{
    static void Main(string[] args)
    {
        uint windowSize = 6;
        double saccadeThreshold = 50.0;
        long timeout = 1000;
        var filter = new OnlineGazeFilter(windowSize, saccadeThreshold, timeout);
        var data = new List<GazePoint>();
        data.Add(new GazePoint() { X = 630.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 635.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 645.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 635.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 640.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 635.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 640.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 638.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 640.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 660.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 645.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 645.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 640.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 640.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 750.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 760.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 740.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 745.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 745.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 750.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 745.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 750.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 745.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 710.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 720.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 715.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 600.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 730.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 740.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 730.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 735.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 730.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 735.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 645.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 650.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 640.0, Y = 0.0 });
        data.Add(new GazePoint() { X = 635.0, Y = 0.0 });
        var filtered = new List<GazePoint>();
        foreach (var point in data)
        {
            filter.AddPoint(point);
            var filteredPoint = filter.GetCurrentPoint();
            if (filteredPoint != null)
            {
                filtered.Add(filteredPoint);
            }
        }
        var xs = data.Select(i => i.X).ToArray();
        var filteredXs = filtered.Select(i => i.X).ToArray();
    }
}


下記図が上記コードの実行結果になります(青が生データ、オレンジがフィルタ後のデータ)。アルゴリズムの性質上どうしても遅延が避けられませんが、アイトラッカーにありがちなスパイクノイズを防いでいることが読み取れます。 f:id:ni4muraano:20170611214551p:plain