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

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

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

PONTUS OLSSON著, “Real-time and Offline Filters for Eye Tracking"に書かれている疑似コードをプログラムに落としました。この論文では視線データから固視データを抽出するための方法が書かれています。
まずは視線データをGazePointクラスとして実装します。

public class GazePoint
{
    public double X;
    public double Y;
    public double D { get { return Math.Sqrt(X * X + Y * Y); } }
    public long TimeStamp;
    public int? FixationId;
    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 };
    }
}


次に論文で提案されているアルゴリズムをOfflineGazeFilterクラスとして実装します。

/// <summary>
/// PONTUS OLSSON, "Real-time and Offline Filters for Eye Tracking"
/// に書かれている疑似コードを実行するクラスです。
/// 視線データから固視位置、固視時間を抽出します。
/// </summary>
public class OfflineGazeFilter
{
    private readonly int _samplingFrequency_Hz;
    private readonly int _windowSize;
    private readonly double _peakThreshold;
    private readonly double _radiusThreshold;

    public OfflineGazeFilter(int samplingFrequency_Hz, double peakThreshold, double radiusThreshold, int? windowSize = null)
    {
        _samplingFrequency_Hz = samplingFrequency_Hz;
        _peakThreshold = peakThreshold;
        _radiusThreshold = radiusThreshold;
        if (windowSize == null)
        {
            _windowSize = DecideWindowSize(samplingFrequency_Hz);
        }
        else
        {
            _windowSize = (int)(windowSize);
        }
    }

    private int DecideWindowSize(int samplingFrequency_Hz)
    {
        // この数字はpp.16で言及されています
        const double impossiblyShortFixation_msec = 80.0;
        // 論文にはありませんが、サンプリング周期の長いアイトラッカーを利用した時に
        // ウインドウサイズが小さくなり過ぎるのを防ぐためのロジックです
        const int minimumWindowSize = 3;
        return Math.Max(minimumWindowSize, (int)Math.Ceiling(samplingFrequency_Hz * impossiblyShortFixation_msec / 1000.0));
    }

    /// <summary>
    /// 視線データから固視位置/固視時間を抽出するメソッド
    /// </summary>
    /// <param name="points">視線データ</param>
    /// <returns>固視位置/固視時間</returns>
    public GazePoint[] Analyze(IList<GazePoint> points)
    {
        var ts = points.Select(i => i.TimeStamp).ToArray();
        // Step1
        List<GazePoint> interpolatedPoints = InterpolateMissingData(points, _samplingFrequency_Hz);
        // Step2
        List<double> differenceVector = CalculateDifferenceVector(interpolatedPoints, _windowSize);
        // Step3
        List<double> peaks = FindPeaks(differenceVector);
        // Step4
        RemovePeaksIfTooClose(_windowSize, ref peaks);
        // Step5
        List<int> peakIndices = CreatePeakIndices(peaks, _peakThreshold);
        // Step6
        return EstimateFixationSpatialPositions(interpolatedPoints, _windowSize, _radiusThreshold, ref peakIndices).ToArray();
    }

    private List<GazePoint> InterpolateMissingData(IList<GazePoint> points, int samplingFrequency_Hz)
    {
        int samplingPeriod_msec = (int)Math.Round(1.0 / samplingFrequency_Hz * 1000, 0, MidpointRounding.AwayFromZero);
        int threshold = (int)Math.Round(samplingPeriod_msec * 1.3, 0, MidpointRounding.AwayFromZero);
        var interpolatedPoints = new List<GazePoint>();
        foreach (var point in points)
        {
            interpolatedPoints.Add(point.Clone());
        }
        for (int i = points.Count - 1; i >= 1; --i)
        {
            long lastTimeStamp = points[i - 1].TimeStamp;
            int count = 0;
            while (points[i].TimeStamp - lastTimeStamp >= threshold)
            {
                GazePoint point = points[i - 1].Clone();
                lastTimeStamp += samplingPeriod_msec;
                point.TimeStamp = lastTimeStamp;
                interpolatedPoints.Insert(i + count, point);
                ++count;
            }
        }

        return interpolatedPoints;
    }

    private List<double> CalculateDifferenceVector(List<GazePoint> points, int windowSize)
    {
        var differences = new List<double>();
        for (int n = windowSize; n < points.Count - windowSize; ++n)
        {
            IEnumerable<GazePoint> pointsBefore = points.Skip(n - windowSize).Take(windowSize);
            double mBeforeX = pointsBefore.Average(i => i.X);
            double mBeforeY = pointsBefore.Average(i => i.Y);
            IEnumerable<GazePoint> pointsAfter = points.Skip(n + 1).Take(windowSize);
            double mAfterX = pointsAfter.Average(i => i.X);
            double mAfterY = pointsAfter.Average(i => i.Y);

            double difference = Math.Sqrt(Math.Pow(mAfterX - mBeforeX, 2) + Math.Pow(mAfterY - mBeforeY, 2));
            differences.Add(difference);
        }

        return differences;
    }

    private List<double> FindPeaks(List<double> diffrenceVector)
    {
        var peaks = new double[diffrenceVector.Count];
        for (int n = 1; n < diffrenceVector.Count - 1; ++n)
        {
            if (diffrenceVector[n] > diffrenceVector[n - 1] &&
                diffrenceVector[n] > diffrenceVector[n + 1])
            {
                peaks[n] = diffrenceVector[n];
            }
        }

        return peaks.ToList();
    }

    private void RemovePeaksIfTooClose(int windowSize, ref List<double> peaks)
    {
        for (int n = windowSize; n < peaks.Count - windowSize; ++n)
        {
            if (peaks[n] != 0.0)
            {
                for (int i = n - windowSize; i < n; ++i)
                {
                    if (peaks[i] < peaks[n]) peaks[i] = 0.0;
                }
                for (int i = n + 1; i < n + windowSize + 1; ++i)
                {
                    if (peaks[i] < peaks[n]) peaks[i] = 0.0;
                }
            }
        }
    }

    private List<int> CreatePeakIndices(List<double> peaks, double peakThreshold)
    {
        var peakIndices = new List<int>();
        peakIndices.Add(0);
        for (int n = 0; n < peaks.Count; ++n)
        {
            if (peaks[n] >= peakThreshold)
            {
                peakIndices.Add(n);
            }
        }
        peakIndices.Add(peaks.Count - 1);

        return peakIndices;
    }

    private List<GazePoint> EstimateFixationSpatialPositions(List<GazePoint> points, int windowSize, double radiusThrehsold, ref List<int> peakIndices)
    {
        double shortestDistance = 0.0;
        var fixations = new List<GazePoint>();
        while (shortestDistance < radiusThrehsold)
        {
            /* Clear Fixation list and calculate new estimates for each iteration. */
            fixations.Clear();

            int id = 1;
            for (int n = 1; n < peakIndices.Count; ++n)
            {
                /* Add median estimate to the Fixation list. Use all rawdata samples
                    * between PeakIndices(n - 1) and PeakIndices(n) when calculating the
                    * estimate. */
                var data = new List<GazePoint>();
                for (int i = peakIndices[n - 1]; i <= peakIndices[n]; ++i)
                {
                    data.Add(points[i + windowSize]);
                }
                GazePoint median = Median(data);
                median.FixationId = id;
                ++id;
                median.FixationDuration_msec = (uint)(data.Last().TimeStamp - data.First().TimeStamp - 1);
                fixations.Add(median);
            }

            /* Set shortestdistance to Inf. If this is not done, it will remember the shortest
                * distance from the previous iteration and no new shortest distance between
                * fixation estimates will be found. */
            shortestDistance = double.PositiveInfinity;

            int index = -1;
            for (int n = 1; n < fixations.Count; ++n)
            {
                /* Calculate the distance as teh Euclidean norm between the current and
                    * previous fixation estimate. If it is shorter than the current shortest
                    * distance, mark it as the shortest and remember its index. */
                double distance = fixations[n] - fixations[n - 1];

                if (distance < shortestDistance)
                {
                    shortestDistance = distance;
                    index = n;
                }
            }

            /* If the shortest intermediate distance between two fixation estimates is shorter
                * than the threshold radius, remove the peak separating them. During the next
                * iteration, they will be joined together. */
            if (shortestDistance < radiusThrehsold)
            {
                peakIndices.RemoveAt(index);
            }
        }

        return fixations;
    }

    private GazePoint Median(List<GazePoint> points)
    {
        IEnumerable<GazePoint> orderedPointsX = points.OrderBy(i => i.X);
        IEnumerable<GazePoint> orderedPointsY = points.OrderBy(i => i.Y);
        GazePoint median;
        if (points.Count % 2 == 0)
        {
            GazePoint p1X = orderedPointsX.ElementAt(points.Count / 2 - 1);
            GazePoint p1Y = orderedPointsY.ElementAt(points.Count / 2 - 1);
            GazePoint p2X = orderedPoints.ElementAt(points.Count / 2);
            GazePoint p2Y = orderedPoints.ElementAt(points.Count / 2);
            median = new GazePoint()
            {
                X = (p1X.X + p2X.X) / 2.0,
                Y = (p1Y.Y + p2Y.Y) / 2.0
            };
        }
        else
        {
            medianX = orderedPointsX.ElementAt(points.Count / 2);
            medianY = orderedPointsY.ElementAt(points.Count / 2);
            median = new GazePoint()
            {
                X = medianX.X,
                Y = medianY.Y
            };
        }
        median.TimeStamp = points.First().TimeStamp;
        return median;
    }
}


最後にOfflineGazeFilterクラスの利用例です。

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 int samplingFrequency = 30;
const double peakThreshold = 30.0;
const double radiusThreshold = 20.0;
const int windowSize = 4;
var algorithm = new OfflineGazeFilter(samplingFrequency, peakThreshold, radiusThreshold, windowSize);
GazePoint[] analyzedPoints = algorithm.Analyze(points);

下記グラフは入力データになりますが、これに対してAnalyzeメソッドの結果は下記のように出力されました。当然ハイパーパラメータにより結果は変わるため、実際には自然な結果に近くなるようなパラメータを探すことが必要です(具体的には上記コードのpeakThreshold/radiusThreshold/windowSizeがハイパーパラメータに該当します)。

  • 固視1 座標:642.5
  • 固視2 座標:737.5

f:id:ni4muraano:20170709002457p:plain