velocity_tracker.cc revision a1401311d1ab56c4ed0a474bd38c108f75cb0cd9
1bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant// Copyright 2014 The Chromium Authors. All rights reserved.
2bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant// Use of this source code is governed by a BSD-style license that can be
3f5256e16dfc425c1d466f6308d4026d529ce9e0bHoward Hinnant// found in the LICENSE file.
4bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
5b64f8b07c104c6cc986570ac8ee0ed16a9f23976Howard Hinnant#include "ui/events/gesture_detection/velocity_tracker.h"
6b64f8b07c104c6cc986570ac8ee0ed16a9f23976Howard Hinnant
7bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant#include <cmath>
8bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
9bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant#include "base/logging.h"
10bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant#include "ui/events/gesture_detection/motion_event.h"
11bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
12bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnantusing base::TimeDelta;
13bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnantusing base::TimeTicks;
14bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
15bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnantnamespace ui {
16bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
17bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant// Implements a particular velocity tracker algorithm.
18061d0cc4db18d17bf01ed14c5db0be098205bd47Marshall Clowclass VelocityTrackerStrategy {
1981381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant public:
20bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  virtual ~VelocityTrackerStrategy() {}
21bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
22bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  virtual void Clear() = 0;
23bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  virtual void ClearPointers(BitSet32 id_bits) = 0;
24bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  virtual void AddMovement(const base::TimeTicks& event_time,
25bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant                           BitSet32 id_bits,
26bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant                           const Position* positions) = 0;
27bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  virtual bool GetEstimator(uint32_t id, Estimator* out_estimator) const = 0;
28bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
29bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant protected:
30bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  VelocityTrackerStrategy() {}
31bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant};
32bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
33bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnantnamespace {
34bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
35bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard HinnantCOMPILE_ASSERT(MotionEvent::MAX_POINTER_ID < 32, max_pointer_id_too_large);
36bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
37bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnantstruct Position {
38bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  float x, y;
39bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant};
40bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
41bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnantstruct Estimator {
42bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  enum { MAX_DEGREE = 4 };
43bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
44bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  // Estimator time base.
45bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  TimeTicks time;
46bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
47bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  // Polynomial coefficients describing motion in X and Y.
48bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  float xcoeff[MAX_DEGREE + 1], ycoeff[MAX_DEGREE + 1];
49bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
50bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  // Polynomial degree (number of coefficients), or zero if no information is
51bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  // available.
52bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  uint32_t degree;
53bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
54bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  // Confidence (coefficient of determination), between 0 (no fit)
55bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  // and 1 (perfect fit).
56bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  float confidence;
57bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
58bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  inline void Clear() {
59bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant    time = TimeTicks();
60bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant    degree = 0;
61bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant    confidence = 0;
62bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant    for (size_t i = 0; i <= MAX_DEGREE; i++) {
63bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant      xcoeff[i] = 0;
64bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant      ycoeff[i] = 0;
65bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant    }
66bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  }
67bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant};
68bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant
6981381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant// Threshold for determining that a pointer has stopped moving.
7081381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant// Some input devices do not send ACTION_MOVE events in the case where a pointer
7181381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant// hasstopped.  We need to detect this case so that we can accurately predict
7281381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant// the velocity after the pointer starts moving again.
7381381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnantconst TimeDelta ASSUME_POINTER_STOPPED_TIME = TimeDelta::FromMilliseconds(40);
7481381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
7581381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnantstatic float VectorDot(const float* a, const float* b, uint32_t m) {
7681381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  float r = 0;
7781381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  while (m--) {
7881381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    r += *(a++) * *(b++);
7981381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  }
8081381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  return r;
8181381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant}
8281381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
8381381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnantstatic float VectorNorm(const float* a, uint32_t m) {
8481381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  float r = 0;
8581381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  while (m--) {
8681381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    float t = *(a++);
8781381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    r += t * t;
8881381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  }
8981381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  return sqrtf(r);
9081381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant}
9181381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
9281381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant// Velocity tracker algorithm based on least-squares linear regression.
9381381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnantclass LeastSquaresVelocityTrackerStrategy : public VelocityTrackerStrategy {
9481381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant public:
9581381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  enum Weighting {
9681381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    // No weights applied.  All data points are equally reliable.
9781381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    WEIGHTING_NONE,
9881381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
9981381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    // Weight by time delta.  Data points clustered together are weighted less.
10081381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    WEIGHTING_DELTA,
10181381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
10281381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    // Weight such that points within a certain horizon are weighed more than
10381381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    // those outside of that horizon.
10481381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    WEIGHTING_CENTRAL,
10581381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
10681381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    // Weight such that points older than a certain amount are weighed less.
10781381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant    WEIGHTING_RECENT,
10881381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  };
10981381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
11081381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  // Number of samples to keep.
11181381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  enum { HISTORY_SIZE = 20 };
11281381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
11381381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  // Degree must be no greater than Estimator::MAX_DEGREE.
11481381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  LeastSquaresVelocityTrackerStrategy(uint32_t degree,
11581381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant                                      Weighting weighting = WEIGHTING_NONE);
11681381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant  virtual ~LeastSquaresVelocityTrackerStrategy();
11781381a932fbebb384adfe5c0116d45b37479efdeHoward Hinnant
118bc8d3f97eb5c958007f2713238472e0c1c8fe02Howard Hinnant  virtual void Clear() OVERRIDE;
119  virtual void ClearPointers(BitSet32 id_bits) OVERRIDE;
120  virtual void AddMovement(const TimeTicks& event_time,
121                           BitSet32 id_bits,
122                           const Position* positions) OVERRIDE;
123  virtual bool GetEstimator(uint32_t id,
124                            Estimator* out_estimator) const OVERRIDE;
125
126 private:
127  // Sample horizon.
128  // We don't use too much history by default since we want to react to quick
129  // changes in direction.
130  static const TimeDelta HORIZON;
131
132  struct Movement {
133    TimeTicks event_time;
134    BitSet32 id_bits;
135    Position positions[VelocityTracker::MAX_POINTERS];
136
137    inline const Position& GetPosition(uint32_t id) const {
138      return positions[id_bits.get_index_of_bit(id)];
139    }
140  };
141
142  float ChooseWeight(uint32_t index) const;
143
144  const uint32_t degree_;
145  const Weighting weighting_;
146  uint32_t index_;
147  Movement movements_[HISTORY_SIZE];
148};
149
150// Velocity tracker algorithm that uses an IIR filter.
151class IntegratingVelocityTrackerStrategy : public VelocityTrackerStrategy {
152 public:
153  // Degree must be 1 or 2.
154  explicit IntegratingVelocityTrackerStrategy(uint32_t degree);
155  virtual ~IntegratingVelocityTrackerStrategy();
156
157  virtual void Clear() OVERRIDE;
158  virtual void ClearPointers(BitSet32 id_bits) OVERRIDE;
159  virtual void AddMovement(const TimeTicks& event_time,
160                           BitSet32 id_bits,
161                           const Position* positions) OVERRIDE;
162  virtual bool GetEstimator(uint32_t id,
163                            Estimator* out_estimator) const OVERRIDE;
164
165 private:
166  // Current state estimate for a particular pointer.
167  struct State {
168    TimeTicks update_time;
169    uint32_t degree;
170
171    float xpos, xvel, xaccel;
172    float ypos, yvel, yaccel;
173  };
174
175  const uint32_t degree_;
176  BitSet32 pointer_id_bits_;
177  State mPointerState[MotionEvent::MAX_POINTER_ID + 1];
178
179  void InitState(State& state,
180                 const TimeTicks& event_time,
181                 float xpos,
182                 float ypos) const;
183  void UpdateState(State& state,
184                   const TimeTicks& event_time,
185                   float xpos,
186                   float ypos) const;
187  void PopulateEstimator(const State& state, Estimator* out_estimator) const;
188};
189
190VelocityTrackerStrategy* CreateStrategy(VelocityTracker::Strategy strategy) {
191  switch (strategy) {
192    case VelocityTracker::LSQ1:
193      return new LeastSquaresVelocityTrackerStrategy(1);
194    case VelocityTracker::LSQ2:
195      return new LeastSquaresVelocityTrackerStrategy(2);
196    case VelocityTracker::LSQ3:
197      return new LeastSquaresVelocityTrackerStrategy(3);
198    case VelocityTracker::WLSQ2_DELTA:
199      return new LeastSquaresVelocityTrackerStrategy(
200          2, LeastSquaresVelocityTrackerStrategy::WEIGHTING_DELTA);
201    case VelocityTracker::WLSQ2_CENTRAL:
202      return new LeastSquaresVelocityTrackerStrategy(
203          2, LeastSquaresVelocityTrackerStrategy::WEIGHTING_CENTRAL);
204    case VelocityTracker::WLSQ2_RECENT:
205      return new LeastSquaresVelocityTrackerStrategy(
206          2, LeastSquaresVelocityTrackerStrategy::WEIGHTING_RECENT);
207    case VelocityTracker::INT1:
208      return new IntegratingVelocityTrackerStrategy(1);
209    case VelocityTracker::INT2:
210      return new IntegratingVelocityTrackerStrategy(2);
211  }
212  NOTREACHED() << "Unrecognized velocity tracker strategy: " << strategy;
213  return CreateStrategy(VelocityTracker::STRATEGY_DEFAULT);
214}
215
216}  // namespace
217
218// --- VelocityTracker ---
219
220VelocityTracker::VelocityTracker()
221    : current_pointer_id_bits_(0),
222      active_pointer_id_(-1),
223      strategy_(CreateStrategy(STRATEGY_DEFAULT)) {}
224
225VelocityTracker::VelocityTracker(Strategy strategy)
226    : current_pointer_id_bits_(0),
227      active_pointer_id_(-1),
228      strategy_(CreateStrategy(strategy)) {}
229
230VelocityTracker::~VelocityTracker() {}
231
232void VelocityTracker::Clear() {
233  current_pointer_id_bits_.clear();
234  active_pointer_id_ = -1;
235  strategy_->Clear();
236}
237
238void VelocityTracker::ClearPointers(BitSet32 id_bits) {
239  BitSet32 remaining_id_bits(current_pointer_id_bits_.value & ~id_bits.value);
240  current_pointer_id_bits_ = remaining_id_bits;
241
242  if (active_pointer_id_ >= 0 && id_bits.has_bit(active_pointer_id_)) {
243    active_pointer_id_ = !remaining_id_bits.is_empty()
244                             ? remaining_id_bits.first_marked_bit()
245                             : -1;
246  }
247
248  strategy_->ClearPointers(id_bits);
249}
250
251void VelocityTracker::AddMovement(const TimeTicks& event_time,
252                                  BitSet32 id_bits,
253                                  const Position* positions) {
254  while (id_bits.count() > MAX_POINTERS)
255    id_bits.clear_last_marked_bit();
256
257  if ((current_pointer_id_bits_.value & id_bits.value) &&
258      event_time >= (last_event_time_ + ASSUME_POINTER_STOPPED_TIME)) {
259    // We have not received any movements for too long.  Assume that all
260    // pointers
261    // have stopped.
262    strategy_->Clear();
263  }
264  last_event_time_ = event_time;
265
266  current_pointer_id_bits_ = id_bits;
267  if (active_pointer_id_ < 0 || !id_bits.has_bit(active_pointer_id_))
268    active_pointer_id_ = id_bits.is_empty() ? -1 : id_bits.first_marked_bit();
269
270  strategy_->AddMovement(event_time, id_bits, positions);
271}
272
273void VelocityTracker::AddMovement(const MotionEvent& event) {
274  int32_t actionMasked = event.GetAction();
275
276  switch (actionMasked) {
277    case MotionEvent::ACTION_DOWN:
278      // case MotionEvent::HOVER_ENTER:
279      // Clear all pointers on down before adding the new movement.
280      Clear();
281      break;
282    case MotionEvent::ACTION_POINTER_DOWN: {
283      // Start a new movement trace for a pointer that just went down.
284      // We do this on down instead of on up because the client may want to
285      // query the final velocity for a pointer that just went up.
286      BitSet32 downIdBits;
287      downIdBits.mark_bit(event.GetPointerId(event.GetActionIndex()));
288      ClearPointers(downIdBits);
289      break;
290    }
291    case MotionEvent::ACTION_MOVE:
292      // case MotionEvent::ACTION_HOVER_MOVE:
293      break;
294    default:
295      // Ignore all other actions because they do not convey any new information
296      // about pointer movement.  We also want to preserve the last known
297      // velocity of the pointers.
298      // Note that ACTION_UP and ACTION_POINTER_UP always report the last known
299      // position of the pointers that went up.  ACTION_POINTER_UP does include
300      // the new position of pointers that remained down but we will also
301      // receive an ACTION_MOVE with this information if any of them actually
302      // moved.  Since we don't know how many pointers will be going up at once
303      // it makes sense to just wait for the following ACTION_MOVE before adding
304      // the movement.
305      return;
306  }
307
308  size_t pointer_count = event.GetPointerCount();
309  if (pointer_count > MAX_POINTERS) {
310    pointer_count = MAX_POINTERS;
311  }
312
313  BitSet32 id_bits;
314  for (size_t i = 0; i < pointer_count; i++) {
315    id_bits.mark_bit(event.GetPointerId(i));
316  }
317
318  uint32_t pointer_index[MAX_POINTERS];
319  for (size_t i = 0; i < pointer_count; i++) {
320    pointer_index[i] = id_bits.get_index_of_bit(event.GetPointerId(i));
321  }
322
323  Position positions[MAX_POINTERS];
324  size_t historySize = event.GetHistorySize();
325  for (size_t h = 0; h < historySize; h++) {
326    for (size_t i = 0; i < pointer_count; i++) {
327      uint32_t index = pointer_index[i];
328      positions[index].x = event.GetHistoricalX(i, h);
329      positions[index].y = event.GetHistoricalY(i, h);
330    }
331    AddMovement(event.GetHistoricalEventTime(h), id_bits, positions);
332  }
333
334  for (size_t i = 0; i < pointer_count; i++) {
335    uint32_t index = pointer_index[i];
336    positions[index].x = event.GetX(i);
337    positions[index].y = event.GetY(i);
338  }
339  AddMovement(event.GetEventTime(), id_bits, positions);
340}
341
342bool VelocityTracker::GetVelocity(uint32_t id,
343                                  float* out_vx,
344                                  float* out_vy) const {
345  Estimator estimator;
346  if (GetEstimator(id, &estimator) && estimator.degree >= 1) {
347    *out_vx = estimator.xcoeff[1];
348    *out_vy = estimator.ycoeff[1];
349    return true;
350  }
351  *out_vx = 0;
352  *out_vy = 0;
353  return false;
354}
355
356void LeastSquaresVelocityTrackerStrategy::AddMovement(
357    const TimeTicks& event_time,
358    BitSet32 id_bits,
359    const Position* positions) {
360  if (++index_ == HISTORY_SIZE) {
361    index_ = 0;
362  }
363
364  Movement& movement = movements_[index_];
365  movement.event_time = event_time;
366  movement.id_bits = id_bits;
367  uint32_t count = id_bits.count();
368  for (uint32_t i = 0; i < count; i++) {
369    movement.positions[i] = positions[i];
370  }
371}
372
373bool VelocityTracker::GetEstimator(uint32_t id,
374                                   Estimator* out_estimator) const {
375  return strategy_->GetEstimator(id, out_estimator);
376}
377
378// --- LeastSquaresVelocityTrackerStrategy ---
379
380const TimeDelta LeastSquaresVelocityTrackerStrategy::HORIZON =
381    TimeDelta::FromMilliseconds(100);
382
383LeastSquaresVelocityTrackerStrategy::LeastSquaresVelocityTrackerStrategy(
384    uint32_t degree,
385    Weighting weighting)
386    : degree_(degree), weighting_(weighting) {
387  DCHECK_LT(degree_, static_cast<uint32_t>(Estimator::MAX_DEGREE));
388  Clear();
389}
390
391LeastSquaresVelocityTrackerStrategy::~LeastSquaresVelocityTrackerStrategy() {}
392
393void LeastSquaresVelocityTrackerStrategy::Clear() {
394  index_ = 0;
395  movements_[0].id_bits.clear();
396}
397
398/**
399 * Solves a linear least squares problem to obtain a N degree polynomial that
400 * fits the specified input data as nearly as possible.
401 *
402 * Returns true if a solution is found, false otherwise.
403 *
404 * The input consists of two vectors of data points X and Y with indices 0..m-1
405 * along with a weight vector W of the same size.
406 *
407 * The output is a vector B with indices 0..n that describes a polynomial
408 * that fits the data, such the sum of W[i] * W[i] * abs(Y[i] - (B[0] + B[1]
409 * X[i] * + B[2] X[i]^2 ... B[n] X[i]^n)) for all i between 0 and m-1 is
410 * minimized.
411 *
412 * Accordingly, the weight vector W should be initialized by the caller with the
413 * reciprocal square root of the variance of the error in each input data point.
414 * In other words, an ideal choice for W would be W[i] = 1 / var(Y[i]) = 1 /
415 * stddev(Y[i]).
416 * The weights express the relative importance of each data point.  If the
417 * weights are* all 1, then the data points are considered to be of equal
418 * importance when fitting the polynomial.  It is a good idea to choose weights
419 * that diminish the importance of data points that may have higher than usual
420 * error margins.
421 *
422 * Errors among data points are assumed to be independent.  W is represented
423 * here as a vector although in the literature it is typically taken to be a
424 * diagonal matrix.
425 *
426 * That is to say, the function that generated the input data can be
427 * approximated by y(x) ~= B[0] + B[1] x + B[2] x^2 + ... + B[n] x^n.
428 *
429 * The coefficient of determination (R^2) is also returned to describe the
430 * goodness of fit of the model for the given data.  It is a value between 0
431 * and 1, where 1 indicates perfect correspondence.
432 *
433 * This function first expands the X vector to a m by n matrix A such that
434 * A[i][0] = 1, A[i][1] = X[i], A[i][2] = X[i]^2, ..., A[i][n] = X[i]^n, then
435 * multiplies it by w[i]./
436 *
437 * Then it calculates the QR decomposition of A yielding an m by m orthonormal
438 * matrix Q and an m by n upper triangular matrix R.  Because R is upper
439 * triangular (lower part is all zeroes), we can simplify the decomposition into
440 * an m by n matrix Q1 and a n by n matrix R1 such that A = Q1 R1.
441 *
442 * Finally we solve the system of linear equations given by
443 * R1 B = (Qtranspose W Y) to find B.
444 *
445 * For efficiency, we lay out A and Q column-wise in memory because we
446 * frequently operate on the column vectors.  Conversely, we lay out R row-wise.
447 *
448 * http://en.wikipedia.org/wiki/Numerical_methods_for_linear_least_squares
449 * http://en.wikipedia.org/wiki/Gram-Schmidt
450 */
451static bool SolveLeastSquares(const float* x,
452                              const float* y,
453                              const float* w,
454                              uint32_t m,
455                              uint32_t n,
456                              float* out_b,
457                              float* out_det) {
458  // MSVC does not support variable-length arrays (used by the original Android
459  // implementation of this function).
460#if defined(COMPILER_MSVC)
461  enum {
462    M_ARRAY_LENGTH = LeastSquaresVelocityTrackerStrategy::HISTORY_SIZE,
463    N_ARRAY_LENGTH = Estimator::MAX_DEGREE
464  };
465  DCHECK_LE(m, static_cast<uint32_t>(M_ARRAY_LENGTH));
466  DCHECK_LE(n, static_cast<uint32_t>(N_ARRAY_LENGTH));
467#else
468  const uint32_t M_ARRAY_LENGTH = m;
469  const uint32_t N_ARRAY_LENGTH = n;
470#endif
471
472  // Expand the X vector to a matrix A, pre-multiplied by the weights.
473  float a[N_ARRAY_LENGTH][M_ARRAY_LENGTH];  // column-major order
474  for (uint32_t h = 0; h < m; h++) {
475    a[0][h] = w[h];
476    for (uint32_t i = 1; i < n; i++) {
477      a[i][h] = a[i - 1][h] * x[h];
478    }
479  }
480
481  // Apply the Gram-Schmidt process to A to obtain its QR decomposition.
482
483  // Orthonormal basis, column-major order.
484  float q[N_ARRAY_LENGTH][M_ARRAY_LENGTH];
485  // Upper triangular matrix, row-major order.
486  float r[N_ARRAY_LENGTH][N_ARRAY_LENGTH];
487  for (uint32_t j = 0; j < n; j++) {
488    for (uint32_t h = 0; h < m; h++) {
489      q[j][h] = a[j][h];
490    }
491    for (uint32_t i = 0; i < j; i++) {
492      float dot = VectorDot(&q[j][0], &q[i][0], m);
493      for (uint32_t h = 0; h < m; h++) {
494        q[j][h] -= dot * q[i][h];
495      }
496    }
497
498    float norm = VectorNorm(&q[j][0], m);
499    if (norm < 0.000001f) {
500      // vectors are linearly dependent or zero so no solution
501      return false;
502    }
503
504    float invNorm = 1.0f / norm;
505    for (uint32_t h = 0; h < m; h++) {
506      q[j][h] *= invNorm;
507    }
508    for (uint32_t i = 0; i < n; i++) {
509      r[j][i] = i < j ? 0 : VectorDot(&q[j][0], &a[i][0], m);
510    }
511  }
512
513  // Solve R B = Qt W Y to find B.  This is easy because R is upper triangular.
514  // We just work from bottom-right to top-left calculating B's coefficients.
515  float wy[M_ARRAY_LENGTH];
516  for (uint32_t h = 0; h < m; h++) {
517    wy[h] = y[h] * w[h];
518  }
519  for (uint32_t i = n; i-- != 0;) {
520    out_b[i] = VectorDot(&q[i][0], wy, m);
521    for (uint32_t j = n - 1; j > i; j--) {
522      out_b[i] -= r[i][j] * out_b[j];
523    }
524    out_b[i] /= r[i][i];
525  }
526
527  // Calculate the coefficient of determination as 1 - (SSerr / SStot) where
528  // SSerr is the residual sum of squares (variance of the error),
529  // and SStot is the total sum of squares (variance of the data) where each
530  // has been weighted.
531  float ymean = 0;
532  for (uint32_t h = 0; h < m; h++) {
533    ymean += y[h];
534  }
535  ymean /= m;
536
537  float sserr = 0;
538  float sstot = 0;
539  for (uint32_t h = 0; h < m; h++) {
540    float err = y[h] - out_b[0];
541    float term = 1;
542    for (uint32_t i = 1; i < n; i++) {
543      term *= x[h];
544      err -= term * out_b[i];
545    }
546    sserr += w[h] * w[h] * err * err;
547    float var = y[h] - ymean;
548    sstot += w[h] * w[h] * var * var;
549  }
550  *out_det = sstot > 0.000001f ? 1.0f - (sserr / sstot) : 1;
551  return true;
552}
553
554void LeastSquaresVelocityTrackerStrategy::ClearPointers(BitSet32 id_bits) {
555  BitSet32 remaining_id_bits(movements_[index_].id_bits.value & ~id_bits.value);
556  movements_[index_].id_bits = remaining_id_bits;
557}
558
559bool LeastSquaresVelocityTrackerStrategy::GetEstimator(
560    uint32_t id,
561    Estimator* out_estimator) const {
562  out_estimator->Clear();
563
564  // Iterate over movement samples in reverse time order and collect samples.
565  float x[HISTORY_SIZE];
566  float y[HISTORY_SIZE];
567  float w[HISTORY_SIZE];
568  float time[HISTORY_SIZE];
569  uint32_t m = 0;
570  uint32_t index = index_;
571  const Movement& newest_movement = movements_[index_];
572  do {
573    const Movement& movement = movements_[index];
574    if (!movement.id_bits.has_bit(id))
575      break;
576
577    TimeDelta age = newest_movement.event_time - movement.event_time;
578    if (age > HORIZON)
579      break;
580
581    const Position& position = movement.GetPosition(id);
582    x[m] = position.x;
583    y[m] = position.y;
584    w[m] = ChooseWeight(index);
585    time[m] = -age.InSecondsF();
586    index = (index == 0 ? HISTORY_SIZE : index) - 1;
587  } while (++m < HISTORY_SIZE);
588
589  if (m == 0)
590    return false;  // no data
591
592  // Calculate a least squares polynomial fit.
593  uint32_t degree = degree_;
594  if (degree > m - 1)
595    degree = m - 1;
596
597  if (degree >= 1) {
598    float xdet, ydet;
599    uint32_t n = degree + 1;
600    if (SolveLeastSquares(time, x, w, m, n, out_estimator->xcoeff, &xdet) &&
601        SolveLeastSquares(time, y, w, m, n, out_estimator->ycoeff, &ydet)) {
602      out_estimator->time = newest_movement.event_time;
603      out_estimator->degree = degree;
604      out_estimator->confidence = xdet * ydet;
605      return true;
606    }
607  }
608
609  // No velocity data available for this pointer, but we do have its current
610  // position.
611  out_estimator->xcoeff[0] = x[0];
612  out_estimator->ycoeff[0] = y[0];
613  out_estimator->time = newest_movement.event_time;
614  out_estimator->degree = 0;
615  out_estimator->confidence = 1;
616  return true;
617}
618
619float LeastSquaresVelocityTrackerStrategy::ChooseWeight(uint32_t index) const {
620  switch (weighting_) {
621    case WEIGHTING_DELTA: {
622      // Weight points based on how much time elapsed between them and the next
623      // point so that points that "cover" a shorter time span are weighed less.
624      //   delta  0ms: 0.5
625      //   delta 10ms: 1.0
626      if (index == index_) {
627        return 1.0f;
628      }
629      uint32_t next_index = (index + 1) % HISTORY_SIZE;
630      float delta_millis =
631          static_cast<float>((movements_[next_index].event_time -
632                              movements_[index].event_time).InMillisecondsF());
633      if (delta_millis < 0)
634        return 0.5f;
635      if (delta_millis < 10)
636        return 0.5f + delta_millis * 0.05;
637
638      return 1.0f;
639    }
640
641    case WEIGHTING_CENTRAL: {
642      // Weight points based on their age, weighing very recent and very old
643      // points less.
644      //   age  0ms: 0.5
645      //   age 10ms: 1.0
646      //   age 50ms: 1.0
647      //   age 60ms: 0.5
648      float age_millis =
649          static_cast<float>((movements_[index_].event_time -
650                              movements_[index].event_time).InMillisecondsF());
651      if (age_millis < 0)
652        return 0.5f;
653      if (age_millis < 10)
654        return 0.5f + age_millis * 0.05;
655      if (age_millis < 50)
656        return 1.0f;
657      if (age_millis < 60)
658        return 0.5f + (60 - age_millis) * 0.05;
659
660      return 0.5f;
661    }
662
663    case WEIGHTING_RECENT: {
664      // Weight points based on their age, weighing older points less.
665      //   age   0ms: 1.0
666      //   age  50ms: 1.0
667      //   age 100ms: 0.5
668      float age_millis =
669          static_cast<float>((movements_[index_].event_time -
670                              movements_[index].event_time).InMillisecondsF());
671      if (age_millis < 50) {
672        return 1.0f;
673      }
674      if (age_millis < 100) {
675        return 0.5f + (100 - age_millis) * 0.01f;
676      }
677      return 0.5f;
678    }
679
680    case WEIGHTING_NONE:
681    default:
682      return 1.0f;
683  }
684}
685
686// --- IntegratingVelocityTrackerStrategy ---
687
688IntegratingVelocityTrackerStrategy::IntegratingVelocityTrackerStrategy(
689    uint32_t degree)
690    : degree_(degree) {
691  DCHECK_LT(degree_, static_cast<uint32_t>(Estimator::MAX_DEGREE));
692}
693
694IntegratingVelocityTrackerStrategy::~IntegratingVelocityTrackerStrategy() {}
695
696void IntegratingVelocityTrackerStrategy::Clear() { pointer_id_bits_.clear(); }
697
698void IntegratingVelocityTrackerStrategy::ClearPointers(BitSet32 id_bits) {
699  pointer_id_bits_.value &= ~id_bits.value;
700}
701
702void IntegratingVelocityTrackerStrategy::AddMovement(
703    const TimeTicks& event_time,
704    BitSet32 id_bits,
705    const Position* positions) {
706  uint32_t index = 0;
707  for (BitSet32 iter_id_bits(id_bits); !iter_id_bits.is_empty();) {
708    uint32_t id = iter_id_bits.clear_first_marked_bit();
709    State& state = mPointerState[id];
710    const Position& position = positions[index++];
711    if (pointer_id_bits_.has_bit(id))
712      UpdateState(state, event_time, position.x, position.y);
713    else
714      InitState(state, event_time, position.x, position.y);
715  }
716
717  pointer_id_bits_ = id_bits;
718}
719
720bool IntegratingVelocityTrackerStrategy::GetEstimator(
721    uint32_t id,
722    Estimator* out_estimator) const {
723  out_estimator->Clear();
724
725  if (pointer_id_bits_.has_bit(id)) {
726    const State& state = mPointerState[id];
727    PopulateEstimator(state, out_estimator);
728    return true;
729  }
730
731  return false;
732}
733
734void IntegratingVelocityTrackerStrategy::InitState(State& state,
735                                                   const TimeTicks& event_time,
736                                                   float xpos,
737                                                   float ypos) const {
738  state.update_time = event_time;
739  state.degree = 0;
740  state.xpos = xpos;
741  state.xvel = 0;
742  state.xaccel = 0;
743  state.ypos = ypos;
744  state.yvel = 0;
745  state.yaccel = 0;
746}
747
748void IntegratingVelocityTrackerStrategy::UpdateState(
749    State& state,
750    const TimeTicks& event_time,
751    float xpos,
752    float ypos) const {
753  const base::TimeDelta MIN_TIME_DELTA = TimeDelta::FromMicroseconds(2);
754  const float FILTER_TIME_CONSTANT = 0.010f;  // 10 milliseconds
755
756  if (event_time <= state.update_time + MIN_TIME_DELTA)
757    return;
758
759  float dt = static_cast<float>((event_time - state.update_time).InSecondsF());
760  state.update_time = event_time;
761
762  float xvel = (xpos - state.xpos) / dt;
763  float yvel = (ypos - state.ypos) / dt;
764  if (state.degree == 0) {
765    state.xvel = xvel;
766    state.yvel = yvel;
767    state.degree = 1;
768  } else {
769    float alpha = dt / (FILTER_TIME_CONSTANT + dt);
770    if (degree_ == 1) {
771      state.xvel += (xvel - state.xvel) * alpha;
772      state.yvel += (yvel - state.yvel) * alpha;
773    } else {
774      float xaccel = (xvel - state.xvel) / dt;
775      float yaccel = (yvel - state.yvel) / dt;
776      if (state.degree == 1) {
777        state.xaccel = xaccel;
778        state.yaccel = yaccel;
779        state.degree = 2;
780      } else {
781        state.xaccel += (xaccel - state.xaccel) * alpha;
782        state.yaccel += (yaccel - state.yaccel) * alpha;
783      }
784      state.xvel += (state.xaccel * dt) * alpha;
785      state.yvel += (state.yaccel * dt) * alpha;
786    }
787  }
788  state.xpos = xpos;
789  state.ypos = ypos;
790}
791
792void IntegratingVelocityTrackerStrategy::PopulateEstimator(
793    const State& state,
794    Estimator* out_estimator) const {
795  out_estimator->time = state.update_time;
796  out_estimator->confidence = 1.0f;
797  out_estimator->degree = state.degree;
798  out_estimator->xcoeff[0] = state.xpos;
799  out_estimator->xcoeff[1] = state.xvel;
800  out_estimator->xcoeff[2] = state.xaccel / 2;
801  out_estimator->ycoeff[0] = state.ypos;
802  out_estimator->ycoeff[1] = state.yvel;
803  out_estimator->ycoeff[2] = state.yaccel / 2;
804}
805
806}  // namespace ui
807