1 /*
2  * Copyright (C) 2012 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *      http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #define LOG_TAG "VelocityTracker"
18 //#define LOG_NDEBUG 0
19 
20 // Log debug messages about velocity tracking.
21 #define DEBUG_VELOCITY 0
22 
23 // Log debug messages about the progress of the algorithm itself.
24 #define DEBUG_STRATEGY 0
25 
26 #include <math.h>
27 #include <limits.h>
28 
29 #include <cutils/properties.h>
30 #include <input/VelocityTracker.h>
31 #include <utils/BitSet.h>
32 #include <utils/String8.h>
33 #include <utils/Timers.h>
34 
35 namespace android {
36 
37 // Nanoseconds per milliseconds.
38 static const nsecs_t NANOS_PER_MS = 1000000;
39 
40 // Threshold for determining that a pointer has stopped moving.
41 // Some input devices do not send ACTION_MOVE events in the case where a pointer has
42 // stopped.  We need to detect this case so that we can accurately predict the
43 // velocity after the pointer starts moving again.
44 static const nsecs_t ASSUME_POINTER_STOPPED_TIME = 40 * NANOS_PER_MS;
45 
46 
vectorDot(const float * a,const float * b,uint32_t m)47 static float vectorDot(const float* a, const float* b, uint32_t m) {
48     float r = 0;
49     while (m--) {
50         r += *(a++) * *(b++);
51     }
52     return r;
53 }
54 
vectorNorm(const float * a,uint32_t m)55 static float vectorNorm(const float* a, uint32_t m) {
56     float r = 0;
57     while (m--) {
58         float t = *(a++);
59         r += t * t;
60     }
61     return sqrtf(r);
62 }
63 
64 #if DEBUG_STRATEGY || DEBUG_VELOCITY
vectorToString(const float * a,uint32_t m)65 static String8 vectorToString(const float* a, uint32_t m) {
66     String8 str;
67     str.append("[");
68     while (m--) {
69         str.appendFormat(" %f", *(a++));
70         if (m) {
71             str.append(",");
72         }
73     }
74     str.append(" ]");
75     return str;
76 }
77 
matrixToString(const float * a,uint32_t m,uint32_t n,bool rowMajor)78 static String8 matrixToString(const float* a, uint32_t m, uint32_t n, bool rowMajor) {
79     String8 str;
80     str.append("[");
81     for (size_t i = 0; i < m; i++) {
82         if (i) {
83             str.append(",");
84         }
85         str.append(" [");
86         for (size_t j = 0; j < n; j++) {
87             if (j) {
88                 str.append(",");
89             }
90             str.appendFormat(" %f", a[rowMajor ? i * n + j : j * m + i]);
91         }
92         str.append(" ]");
93     }
94     str.append(" ]");
95     return str;
96 }
97 #endif
98 
99 
100 // --- VelocityTracker ---
101 
102 // The default velocity tracker strategy.
103 // Although other strategies are available for testing and comparison purposes,
104 // this is the strategy that applications will actually use.  Be very careful
105 // when adjusting the default strategy because it can dramatically affect
106 // (often in a bad way) the user experience.
107 const char* VelocityTracker::DEFAULT_STRATEGY = "lsq2";
108 
VelocityTracker(const char * strategy)109 VelocityTracker::VelocityTracker(const char* strategy) :
110         mLastEventTime(0), mCurrentPointerIdBits(0), mActivePointerId(-1) {
111     char value[PROPERTY_VALUE_MAX];
112 
113     // Allow the default strategy to be overridden using a system property for debugging.
114     if (!strategy) {
115         int length = property_get("debug.velocitytracker.strategy", value, NULL);
116         if (length > 0) {
117             strategy = value;
118         } else {
119             strategy = DEFAULT_STRATEGY;
120         }
121     }
122 
123     // Configure the strategy.
124     if (!configureStrategy(strategy)) {
125         ALOGD("Unrecognized velocity tracker strategy name '%s'.", strategy);
126         if (!configureStrategy(DEFAULT_STRATEGY)) {
127             LOG_ALWAYS_FATAL("Could not create the default velocity tracker strategy '%s'!",
128                     strategy);
129         }
130     }
131 }
132 
~VelocityTracker()133 VelocityTracker::~VelocityTracker() {
134     delete mStrategy;
135 }
136 
configureStrategy(const char * strategy)137 bool VelocityTracker::configureStrategy(const char* strategy) {
138     mStrategy = createStrategy(strategy);
139     return mStrategy != NULL;
140 }
141 
createStrategy(const char * strategy)142 VelocityTrackerStrategy* VelocityTracker::createStrategy(const char* strategy) {
143     if (!strcmp("lsq1", strategy)) {
144         // 1st order least squares.  Quality: POOR.
145         // Frequently underfits the touch data especially when the finger accelerates
146         // or changes direction.  Often underestimates velocity.  The direction
147         // is overly influenced by historical touch points.
148         return new LeastSquaresVelocityTrackerStrategy(1);
149     }
150     if (!strcmp("lsq2", strategy)) {
151         // 2nd order least squares.  Quality: VERY GOOD.
152         // Pretty much ideal, but can be confused by certain kinds of touch data,
153         // particularly if the panel has a tendency to generate delayed,
154         // duplicate or jittery touch coordinates when the finger is released.
155         return new LeastSquaresVelocityTrackerStrategy(2);
156     }
157     if (!strcmp("lsq3", strategy)) {
158         // 3rd order least squares.  Quality: UNUSABLE.
159         // Frequently overfits the touch data yielding wildly divergent estimates
160         // of the velocity when the finger is released.
161         return new LeastSquaresVelocityTrackerStrategy(3);
162     }
163     if (!strcmp("wlsq2-delta", strategy)) {
164         // 2nd order weighted least squares, delta weighting.  Quality: EXPERIMENTAL
165         return new LeastSquaresVelocityTrackerStrategy(2,
166                 LeastSquaresVelocityTrackerStrategy::WEIGHTING_DELTA);
167     }
168     if (!strcmp("wlsq2-central", strategy)) {
169         // 2nd order weighted least squares, central weighting.  Quality: EXPERIMENTAL
170         return new LeastSquaresVelocityTrackerStrategy(2,
171                 LeastSquaresVelocityTrackerStrategy::WEIGHTING_CENTRAL);
172     }
173     if (!strcmp("wlsq2-recent", strategy)) {
174         // 2nd order weighted least squares, recent weighting.  Quality: EXPERIMENTAL
175         return new LeastSquaresVelocityTrackerStrategy(2,
176                 LeastSquaresVelocityTrackerStrategy::WEIGHTING_RECENT);
177     }
178     if (!strcmp("int1", strategy)) {
179         // 1st order integrating filter.  Quality: GOOD.
180         // Not as good as 'lsq2' because it cannot estimate acceleration but it is
181         // more tolerant of errors.  Like 'lsq1', this strategy tends to underestimate
182         // the velocity of a fling but this strategy tends to respond to changes in
183         // direction more quickly and accurately.
184         return new IntegratingVelocityTrackerStrategy(1);
185     }
186     if (!strcmp("int2", strategy)) {
187         // 2nd order integrating filter.  Quality: EXPERIMENTAL.
188         // For comparison purposes only.  Unlike 'int1' this strategy can compensate
189         // for acceleration but it typically overestimates the effect.
190         return new IntegratingVelocityTrackerStrategy(2);
191     }
192     if (!strcmp("legacy", strategy)) {
193         // Legacy velocity tracker algorithm.  Quality: POOR.
194         // For comparison purposes only.  This algorithm is strongly influenced by
195         // old data points, consistently underestimates velocity and takes a very long
196         // time to adjust to changes in direction.
197         return new LegacyVelocityTrackerStrategy();
198     }
199     return NULL;
200 }
201 
clear()202 void VelocityTracker::clear() {
203     mCurrentPointerIdBits.clear();
204     mActivePointerId = -1;
205 
206     mStrategy->clear();
207 }
208 
clearPointers(BitSet32 idBits)209 void VelocityTracker::clearPointers(BitSet32 idBits) {
210     BitSet32 remainingIdBits(mCurrentPointerIdBits.value & ~idBits.value);
211     mCurrentPointerIdBits = remainingIdBits;
212 
213     if (mActivePointerId >= 0 && idBits.hasBit(mActivePointerId)) {
214         mActivePointerId = !remainingIdBits.isEmpty() ? remainingIdBits.firstMarkedBit() : -1;
215     }
216 
217     mStrategy->clearPointers(idBits);
218 }
219 
addMovement(nsecs_t eventTime,BitSet32 idBits,const Position * positions)220 void VelocityTracker::addMovement(nsecs_t eventTime, BitSet32 idBits, const Position* positions) {
221     while (idBits.count() > MAX_POINTERS) {
222         idBits.clearLastMarkedBit();
223     }
224 
225     if ((mCurrentPointerIdBits.value & idBits.value)
226             && eventTime >= mLastEventTime + ASSUME_POINTER_STOPPED_TIME) {
227 #if DEBUG_VELOCITY
228         ALOGD("VelocityTracker: stopped for %0.3f ms, clearing state.",
229                 (eventTime - mLastEventTime) * 0.000001f);
230 #endif
231         // We have not received any movements for too long.  Assume that all pointers
232         // have stopped.
233         mStrategy->clear();
234     }
235     mLastEventTime = eventTime;
236 
237     mCurrentPointerIdBits = idBits;
238     if (mActivePointerId < 0 || !idBits.hasBit(mActivePointerId)) {
239         mActivePointerId = idBits.isEmpty() ? -1 : idBits.firstMarkedBit();
240     }
241 
242     mStrategy->addMovement(eventTime, idBits, positions);
243 
244 #if DEBUG_VELOCITY
245     ALOGD("VelocityTracker: addMovement eventTime=%lld, idBits=0x%08x, activePointerId=%d",
246             eventTime, idBits.value, mActivePointerId);
247     for (BitSet32 iterBits(idBits); !iterBits.isEmpty(); ) {
248         uint32_t id = iterBits.firstMarkedBit();
249         uint32_t index = idBits.getIndexOfBit(id);
250         iterBits.clearBit(id);
251         Estimator estimator;
252         getEstimator(id, &estimator);
253         ALOGD("  %d: position (%0.3f, %0.3f), "
254                 "estimator (degree=%d, xCoeff=%s, yCoeff=%s, confidence=%f)",
255                 id, positions[index].x, positions[index].y,
256                 int(estimator.degree),
257                 vectorToString(estimator.xCoeff, estimator.degree + 1).string(),
258                 vectorToString(estimator.yCoeff, estimator.degree + 1).string(),
259                 estimator.confidence);
260     }
261 #endif
262 }
263 
addMovement(const MotionEvent * event)264 void VelocityTracker::addMovement(const MotionEvent* event) {
265     int32_t actionMasked = event->getActionMasked();
266 
267     switch (actionMasked) {
268     case AMOTION_EVENT_ACTION_DOWN:
269     case AMOTION_EVENT_ACTION_HOVER_ENTER:
270         // Clear all pointers on down before adding the new movement.
271         clear();
272         break;
273     case AMOTION_EVENT_ACTION_POINTER_DOWN: {
274         // Start a new movement trace for a pointer that just went down.
275         // We do this on down instead of on up because the client may want to query the
276         // final velocity for a pointer that just went up.
277         BitSet32 downIdBits;
278         downIdBits.markBit(event->getPointerId(event->getActionIndex()));
279         clearPointers(downIdBits);
280         break;
281     }
282     case AMOTION_EVENT_ACTION_MOVE:
283     case AMOTION_EVENT_ACTION_HOVER_MOVE:
284         break;
285     default:
286         // Ignore all other actions because they do not convey any new information about
287         // pointer movement.  We also want to preserve the last known velocity of the pointers.
288         // Note that ACTION_UP and ACTION_POINTER_UP always report the last known position
289         // of the pointers that went up.  ACTION_POINTER_UP does include the new position of
290         // pointers that remained down but we will also receive an ACTION_MOVE with this
291         // information if any of them actually moved.  Since we don't know how many pointers
292         // will be going up at once it makes sense to just wait for the following ACTION_MOVE
293         // before adding the movement.
294         return;
295     }
296 
297     size_t pointerCount = event->getPointerCount();
298     if (pointerCount > MAX_POINTERS) {
299         pointerCount = MAX_POINTERS;
300     }
301 
302     BitSet32 idBits;
303     for (size_t i = 0; i < pointerCount; i++) {
304         idBits.markBit(event->getPointerId(i));
305     }
306 
307     uint32_t pointerIndex[MAX_POINTERS];
308     for (size_t i = 0; i < pointerCount; i++) {
309         pointerIndex[i] = idBits.getIndexOfBit(event->getPointerId(i));
310     }
311 
312     nsecs_t eventTime;
313     Position positions[pointerCount];
314 
315     size_t historySize = event->getHistorySize();
316     for (size_t h = 0; h < historySize; h++) {
317         eventTime = event->getHistoricalEventTime(h);
318         for (size_t i = 0; i < pointerCount; i++) {
319             uint32_t index = pointerIndex[i];
320             positions[index].x = event->getHistoricalX(i, h);
321             positions[index].y = event->getHistoricalY(i, h);
322         }
323         addMovement(eventTime, idBits, positions);
324     }
325 
326     eventTime = event->getEventTime();
327     for (size_t i = 0; i < pointerCount; i++) {
328         uint32_t index = pointerIndex[i];
329         positions[index].x = event->getX(i);
330         positions[index].y = event->getY(i);
331     }
332     addMovement(eventTime, idBits, positions);
333 }
334 
getVelocity(uint32_t id,float * outVx,float * outVy) const335 bool VelocityTracker::getVelocity(uint32_t id, float* outVx, float* outVy) const {
336     Estimator estimator;
337     if (getEstimator(id, &estimator) && estimator.degree >= 1) {
338         *outVx = estimator.xCoeff[1];
339         *outVy = estimator.yCoeff[1];
340         return true;
341     }
342     *outVx = 0;
343     *outVy = 0;
344     return false;
345 }
346 
getEstimator(uint32_t id,Estimator * outEstimator) const347 bool VelocityTracker::getEstimator(uint32_t id, Estimator* outEstimator) const {
348     return mStrategy->getEstimator(id, outEstimator);
349 }
350 
351 
352 // --- LeastSquaresVelocityTrackerStrategy ---
353 
354 const nsecs_t LeastSquaresVelocityTrackerStrategy::HORIZON;
355 const uint32_t LeastSquaresVelocityTrackerStrategy::HISTORY_SIZE;
356 
LeastSquaresVelocityTrackerStrategy(uint32_t degree,Weighting weighting)357 LeastSquaresVelocityTrackerStrategy::LeastSquaresVelocityTrackerStrategy(
358         uint32_t degree, Weighting weighting) :
359         mDegree(degree), mWeighting(weighting) {
360     clear();
361 }
362 
~LeastSquaresVelocityTrackerStrategy()363 LeastSquaresVelocityTrackerStrategy::~LeastSquaresVelocityTrackerStrategy() {
364 }
365 
clear()366 void LeastSquaresVelocityTrackerStrategy::clear() {
367     mIndex = 0;
368     mMovements[0].idBits.clear();
369 }
370 
clearPointers(BitSet32 idBits)371 void LeastSquaresVelocityTrackerStrategy::clearPointers(BitSet32 idBits) {
372     BitSet32 remainingIdBits(mMovements[mIndex].idBits.value & ~idBits.value);
373     mMovements[mIndex].idBits = remainingIdBits;
374 }
375 
addMovement(nsecs_t eventTime,BitSet32 idBits,const VelocityTracker::Position * positions)376 void LeastSquaresVelocityTrackerStrategy::addMovement(nsecs_t eventTime, BitSet32 idBits,
377         const VelocityTracker::Position* positions) {
378     if (++mIndex == HISTORY_SIZE) {
379         mIndex = 0;
380     }
381 
382     Movement& movement = mMovements[mIndex];
383     movement.eventTime = eventTime;
384     movement.idBits = idBits;
385     uint32_t count = idBits.count();
386     for (uint32_t i = 0; i < count; i++) {
387         movement.positions[i] = positions[i];
388     }
389 }
390 
391 /**
392  * Solves a linear least squares problem to obtain a N degree polynomial that fits
393  * the specified input data as nearly as possible.
394  *
395  * Returns true if a solution is found, false otherwise.
396  *
397  * The input consists of two vectors of data points X and Y with indices 0..m-1
398  * along with a weight vector W of the same size.
399  *
400  * The output is a vector B with indices 0..n that describes a polynomial
401  * that fits the data, such the sum of W[i] * W[i] * abs(Y[i] - (B[0] + B[1] X[i]
402  * + B[2] X[i]^2 ... B[n] X[i]^n)) for all i between 0 and m-1 is minimized.
403  *
404  * Accordingly, the weight vector W should be initialized by the caller with the
405  * reciprocal square root of the variance of the error in each input data point.
406  * In other words, an ideal choice for W would be W[i] = 1 / var(Y[i]) = 1 / stddev(Y[i]).
407  * The weights express the relative importance of each data point.  If the weights are
408  * all 1, then the data points are considered to be of equal importance when fitting
409  * the polynomial.  It is a good idea to choose weights that diminish the importance
410  * of data points that may have higher than usual error margins.
411  *
412  * Errors among data points are assumed to be independent.  W is represented here
413  * as a vector although in the literature it is typically taken to be a diagonal matrix.
414  *
415  * That is to say, the function that generated the input data can be approximated
416  * by y(x) ~= B[0] + B[1] x + B[2] x^2 + ... + B[n] x^n.
417  *
418  * The coefficient of determination (R^2) is also returned to describe the goodness
419  * of fit of the model for the given data.  It is a value between 0 and 1, where 1
420  * indicates perfect correspondence.
421  *
422  * This function first expands the X vector to a m by n matrix A such that
423  * A[i][0] = 1, A[i][1] = X[i], A[i][2] = X[i]^2, ..., A[i][n] = X[i]^n, then
424  * multiplies it by w[i]./
425  *
426  * Then it calculates the QR decomposition of A yielding an m by m orthonormal matrix Q
427  * and an m by n upper triangular matrix R.  Because R is upper triangular (lower
428  * part is all zeroes), we can simplify the decomposition into an m by n matrix
429  * Q1 and a n by n matrix R1 such that A = Q1 R1.
430  *
431  * Finally we solve the system of linear equations given by R1 B = (Qtranspose W Y)
432  * to find B.
433  *
434  * For efficiency, we lay out A and Q column-wise in memory because we frequently
435  * operate on the column vectors.  Conversely, we lay out R row-wise.
436  *
437  * http://en.wikipedia.org/wiki/Numerical_methods_for_linear_least_squares
438  * http://en.wikipedia.org/wiki/Gram-Schmidt
439  */
solveLeastSquares(const float * x,const float * y,const float * w,uint32_t m,uint32_t n,float * outB,float * outDet)440 static bool solveLeastSquares(const float* x, const float* y,
441         const float* w, uint32_t m, uint32_t n, float* outB, float* outDet) {
442 #if DEBUG_STRATEGY
443     ALOGD("solveLeastSquares: m=%d, n=%d, x=%s, y=%s, w=%s", int(m), int(n),
444             vectorToString(x, m).string(), vectorToString(y, m).string(),
445             vectorToString(w, m).string());
446 #endif
447 
448     // Expand the X vector to a matrix A, pre-multiplied by the weights.
449     float a[n][m]; // column-major order
450     for (uint32_t h = 0; h < m; h++) {
451         a[0][h] = w[h];
452         for (uint32_t i = 1; i < n; i++) {
453             a[i][h] = a[i - 1][h] * x[h];
454         }
455     }
456 #if DEBUG_STRATEGY
457     ALOGD("  - a=%s", matrixToString(&a[0][0], m, n, false /*rowMajor*/).string());
458 #endif
459 
460     // Apply the Gram-Schmidt process to A to obtain its QR decomposition.
461     float q[n][m]; // orthonormal basis, column-major order
462     float r[n][n]; // upper triangular matrix, row-major order
463     for (uint32_t j = 0; j < n; j++) {
464         for (uint32_t h = 0; h < m; h++) {
465             q[j][h] = a[j][h];
466         }
467         for (uint32_t i = 0; i < j; i++) {
468             float dot = vectorDot(&q[j][0], &q[i][0], m);
469             for (uint32_t h = 0; h < m; h++) {
470                 q[j][h] -= dot * q[i][h];
471             }
472         }
473 
474         float norm = vectorNorm(&q[j][0], m);
475         if (norm < 0.000001f) {
476             // vectors are linearly dependent or zero so no solution
477 #if DEBUG_STRATEGY
478             ALOGD("  - no solution, norm=%f", norm);
479 #endif
480             return false;
481         }
482 
483         float invNorm = 1.0f / norm;
484         for (uint32_t h = 0; h < m; h++) {
485             q[j][h] *= invNorm;
486         }
487         for (uint32_t i = 0; i < n; i++) {
488             r[j][i] = i < j ? 0 : vectorDot(&q[j][0], &a[i][0], m);
489         }
490     }
491 #if DEBUG_STRATEGY
492     ALOGD("  - q=%s", matrixToString(&q[0][0], m, n, false /*rowMajor*/).string());
493     ALOGD("  - r=%s", matrixToString(&r[0][0], n, n, true /*rowMajor*/).string());
494 
495     // calculate QR, if we factored A correctly then QR should equal A
496     float qr[n][m];
497     for (uint32_t h = 0; h < m; h++) {
498         for (uint32_t i = 0; i < n; i++) {
499             qr[i][h] = 0;
500             for (uint32_t j = 0; j < n; j++) {
501                 qr[i][h] += q[j][h] * r[j][i];
502             }
503         }
504     }
505     ALOGD("  - qr=%s", matrixToString(&qr[0][0], m, n, false /*rowMajor*/).string());
506 #endif
507 
508     // Solve R B = Qt W Y to find B.  This is easy because R is upper triangular.
509     // We just work from bottom-right to top-left calculating B's coefficients.
510     float wy[m];
511     for (uint32_t h = 0; h < m; h++) {
512         wy[h] = y[h] * w[h];
513     }
514     for (uint32_t i = n; i-- != 0; ) {
515         outB[i] = vectorDot(&q[i][0], wy, m);
516         for (uint32_t j = n - 1; j > i; j--) {
517             outB[i] -= r[i][j] * outB[j];
518         }
519         outB[i] /= r[i][i];
520     }
521 #if DEBUG_STRATEGY
522     ALOGD("  - b=%s", vectorToString(outB, n).string());
523 #endif
524 
525     // Calculate the coefficient of determination as 1 - (SSerr / SStot) where
526     // SSerr is the residual sum of squares (variance of the error),
527     // and SStot is the total sum of squares (variance of the data) where each
528     // has been weighted.
529     float ymean = 0;
530     for (uint32_t h = 0; h < m; h++) {
531         ymean += y[h];
532     }
533     ymean /= m;
534 
535     float sserr = 0;
536     float sstot = 0;
537     for (uint32_t h = 0; h < m; h++) {
538         float err = y[h] - outB[0];
539         float term = 1;
540         for (uint32_t i = 1; i < n; i++) {
541             term *= x[h];
542             err -= term * outB[i];
543         }
544         sserr += w[h] * w[h] * err * err;
545         float var = y[h] - ymean;
546         sstot += w[h] * w[h] * var * var;
547     }
548     *outDet = sstot > 0.000001f ? 1.0f - (sserr / sstot) : 1;
549 #if DEBUG_STRATEGY
550     ALOGD("  - sserr=%f", sserr);
551     ALOGD("  - sstot=%f", sstot);
552     ALOGD("  - det=%f", *outDet);
553 #endif
554     return true;
555 }
556 
getEstimator(uint32_t id,VelocityTracker::Estimator * outEstimator) const557 bool LeastSquaresVelocityTrackerStrategy::getEstimator(uint32_t id,
558         VelocityTracker::Estimator* outEstimator) const {
559     outEstimator->clear();
560 
561     // Iterate over movement samples in reverse time order and collect samples.
562     float x[HISTORY_SIZE];
563     float y[HISTORY_SIZE];
564     float w[HISTORY_SIZE];
565     float time[HISTORY_SIZE];
566     uint32_t m = 0;
567     uint32_t index = mIndex;
568     const Movement& newestMovement = mMovements[mIndex];
569     do {
570         const Movement& movement = mMovements[index];
571         if (!movement.idBits.hasBit(id)) {
572             break;
573         }
574 
575         nsecs_t age = newestMovement.eventTime - movement.eventTime;
576         if (age > HORIZON) {
577             break;
578         }
579 
580         const VelocityTracker::Position& position = movement.getPosition(id);
581         x[m] = position.x;
582         y[m] = position.y;
583         w[m] = chooseWeight(index);
584         time[m] = -age * 0.000000001f;
585         index = (index == 0 ? HISTORY_SIZE : index) - 1;
586     } while (++m < HISTORY_SIZE);
587 
588     if (m == 0) {
589         return false; // no data
590     }
591 
592     // Calculate a least squares polynomial fit.
593     uint32_t degree = mDegree;
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, outEstimator->xCoeff, &xdet)
601                 && solveLeastSquares(time, y, w, m, n, outEstimator->yCoeff, &ydet)) {
602             outEstimator->time = newestMovement.eventTime;
603             outEstimator->degree = degree;
604             outEstimator->confidence = xdet * ydet;
605 #if DEBUG_STRATEGY
606             ALOGD("estimate: degree=%d, xCoeff=%s, yCoeff=%s, confidence=%f",
607                     int(outEstimator->degree),
608                     vectorToString(outEstimator->xCoeff, n).string(),
609                     vectorToString(outEstimator->yCoeff, n).string(),
610                     outEstimator->confidence);
611 #endif
612             return true;
613         }
614     }
615 
616     // No velocity data available for this pointer, but we do have its current position.
617     outEstimator->xCoeff[0] = x[0];
618     outEstimator->yCoeff[0] = y[0];
619     outEstimator->time = newestMovement.eventTime;
620     outEstimator->degree = 0;
621     outEstimator->confidence = 1;
622     return true;
623 }
624 
chooseWeight(uint32_t index) const625 float LeastSquaresVelocityTrackerStrategy::chooseWeight(uint32_t index) const {
626     switch (mWeighting) {
627     case WEIGHTING_DELTA: {
628         // Weight points based on how much time elapsed between them and the next
629         // point so that points that "cover" a shorter time span are weighed less.
630         //   delta  0ms: 0.5
631         //   delta 10ms: 1.0
632         if (index == mIndex) {
633             return 1.0f;
634         }
635         uint32_t nextIndex = (index + 1) % HISTORY_SIZE;
636         float deltaMillis = (mMovements[nextIndex].eventTime- mMovements[index].eventTime)
637                 * 0.000001f;
638         if (deltaMillis < 0) {
639             return 0.5f;
640         }
641         if (deltaMillis < 10) {
642             return 0.5f + deltaMillis * 0.05;
643         }
644         return 1.0f;
645     }
646 
647     case WEIGHTING_CENTRAL: {
648         // Weight points based on their age, weighing very recent and very old points less.
649         //   age  0ms: 0.5
650         //   age 10ms: 1.0
651         //   age 50ms: 1.0
652         //   age 60ms: 0.5
653         float ageMillis = (mMovements[mIndex].eventTime - mMovements[index].eventTime)
654                 * 0.000001f;
655         if (ageMillis < 0) {
656             return 0.5f;
657         }
658         if (ageMillis < 10) {
659             return 0.5f + ageMillis * 0.05;
660         }
661         if (ageMillis < 50) {
662             return 1.0f;
663         }
664         if (ageMillis < 60) {
665             return 0.5f + (60 - ageMillis) * 0.05;
666         }
667         return 0.5f;
668     }
669 
670     case WEIGHTING_RECENT: {
671         // Weight points based on their age, weighing older points less.
672         //   age   0ms: 1.0
673         //   age  50ms: 1.0
674         //   age 100ms: 0.5
675         float ageMillis = (mMovements[mIndex].eventTime - mMovements[index].eventTime)
676                 * 0.000001f;
677         if (ageMillis < 50) {
678             return 1.0f;
679         }
680         if (ageMillis < 100) {
681             return 0.5f + (100 - ageMillis) * 0.01f;
682         }
683         return 0.5f;
684     }
685 
686     case WEIGHTING_NONE:
687     default:
688         return 1.0f;
689     }
690 }
691 
692 
693 // --- IntegratingVelocityTrackerStrategy ---
694 
IntegratingVelocityTrackerStrategy(uint32_t degree)695 IntegratingVelocityTrackerStrategy::IntegratingVelocityTrackerStrategy(uint32_t degree) :
696         mDegree(degree) {
697 }
698 
~IntegratingVelocityTrackerStrategy()699 IntegratingVelocityTrackerStrategy::~IntegratingVelocityTrackerStrategy() {
700 }
701 
clear()702 void IntegratingVelocityTrackerStrategy::clear() {
703     mPointerIdBits.clear();
704 }
705 
clearPointers(BitSet32 idBits)706 void IntegratingVelocityTrackerStrategy::clearPointers(BitSet32 idBits) {
707     mPointerIdBits.value &= ~idBits.value;
708 }
709 
addMovement(nsecs_t eventTime,BitSet32 idBits,const VelocityTracker::Position * positions)710 void IntegratingVelocityTrackerStrategy::addMovement(nsecs_t eventTime, BitSet32 idBits,
711         const VelocityTracker::Position* positions) {
712     uint32_t index = 0;
713     for (BitSet32 iterIdBits(idBits); !iterIdBits.isEmpty();) {
714         uint32_t id = iterIdBits.clearFirstMarkedBit();
715         State& state = mPointerState[id];
716         const VelocityTracker::Position& position = positions[index++];
717         if (mPointerIdBits.hasBit(id)) {
718             updateState(state, eventTime, position.x, position.y);
719         } else {
720             initState(state, eventTime, position.x, position.y);
721         }
722     }
723 
724     mPointerIdBits = idBits;
725 }
726 
getEstimator(uint32_t id,VelocityTracker::Estimator * outEstimator) const727 bool IntegratingVelocityTrackerStrategy::getEstimator(uint32_t id,
728         VelocityTracker::Estimator* outEstimator) const {
729     outEstimator->clear();
730 
731     if (mPointerIdBits.hasBit(id)) {
732         const State& state = mPointerState[id];
733         populateEstimator(state, outEstimator);
734         return true;
735     }
736 
737     return false;
738 }
739 
initState(State & state,nsecs_t eventTime,float xpos,float ypos) const740 void IntegratingVelocityTrackerStrategy::initState(State& state,
741         nsecs_t eventTime, float xpos, float ypos) const {
742     state.updateTime = eventTime;
743     state.degree = 0;
744 
745     state.xpos = xpos;
746     state.xvel = 0;
747     state.xaccel = 0;
748     state.ypos = ypos;
749     state.yvel = 0;
750     state.yaccel = 0;
751 }
752 
updateState(State & state,nsecs_t eventTime,float xpos,float ypos) const753 void IntegratingVelocityTrackerStrategy::updateState(State& state,
754         nsecs_t eventTime, float xpos, float ypos) const {
755     const nsecs_t MIN_TIME_DELTA = 2 * NANOS_PER_MS;
756     const float FILTER_TIME_CONSTANT = 0.010f; // 10 milliseconds
757 
758     if (eventTime <= state.updateTime + MIN_TIME_DELTA) {
759         return;
760     }
761 
762     float dt = (eventTime - state.updateTime) * 0.000000001f;
763     state.updateTime = eventTime;
764 
765     float xvel = (xpos - state.xpos) / dt;
766     float yvel = (ypos - state.ypos) / dt;
767     if (state.degree == 0) {
768         state.xvel = xvel;
769         state.yvel = yvel;
770         state.degree = 1;
771     } else {
772         float alpha = dt / (FILTER_TIME_CONSTANT + dt);
773         if (mDegree == 1) {
774             state.xvel += (xvel - state.xvel) * alpha;
775             state.yvel += (yvel - state.yvel) * alpha;
776         } else {
777             float xaccel = (xvel - state.xvel) / dt;
778             float yaccel = (yvel - state.yvel) / dt;
779             if (state.degree == 1) {
780                 state.xaccel = xaccel;
781                 state.yaccel = yaccel;
782                 state.degree = 2;
783             } else {
784                 state.xaccel += (xaccel - state.xaccel) * alpha;
785                 state.yaccel += (yaccel - state.yaccel) * alpha;
786             }
787             state.xvel += (state.xaccel * dt) * alpha;
788             state.yvel += (state.yaccel * dt) * alpha;
789         }
790     }
791     state.xpos = xpos;
792     state.ypos = ypos;
793 }
794 
populateEstimator(const State & state,VelocityTracker::Estimator * outEstimator) const795 void IntegratingVelocityTrackerStrategy::populateEstimator(const State& state,
796         VelocityTracker::Estimator* outEstimator) const {
797     outEstimator->time = state.updateTime;
798     outEstimator->confidence = 1.0f;
799     outEstimator->degree = state.degree;
800     outEstimator->xCoeff[0] = state.xpos;
801     outEstimator->xCoeff[1] = state.xvel;
802     outEstimator->xCoeff[2] = state.xaccel / 2;
803     outEstimator->yCoeff[0] = state.ypos;
804     outEstimator->yCoeff[1] = state.yvel;
805     outEstimator->yCoeff[2] = state.yaccel / 2;
806 }
807 
808 
809 // --- LegacyVelocityTrackerStrategy ---
810 
811 const nsecs_t LegacyVelocityTrackerStrategy::HORIZON;
812 const uint32_t LegacyVelocityTrackerStrategy::HISTORY_SIZE;
813 const nsecs_t LegacyVelocityTrackerStrategy::MIN_DURATION;
814 
LegacyVelocityTrackerStrategy()815 LegacyVelocityTrackerStrategy::LegacyVelocityTrackerStrategy() {
816     clear();
817 }
818 
~LegacyVelocityTrackerStrategy()819 LegacyVelocityTrackerStrategy::~LegacyVelocityTrackerStrategy() {
820 }
821 
clear()822 void LegacyVelocityTrackerStrategy::clear() {
823     mIndex = 0;
824     mMovements[0].idBits.clear();
825 }
826 
clearPointers(BitSet32 idBits)827 void LegacyVelocityTrackerStrategy::clearPointers(BitSet32 idBits) {
828     BitSet32 remainingIdBits(mMovements[mIndex].idBits.value & ~idBits.value);
829     mMovements[mIndex].idBits = remainingIdBits;
830 }
831 
addMovement(nsecs_t eventTime,BitSet32 idBits,const VelocityTracker::Position * positions)832 void LegacyVelocityTrackerStrategy::addMovement(nsecs_t eventTime, BitSet32 idBits,
833         const VelocityTracker::Position* positions) {
834     if (++mIndex == HISTORY_SIZE) {
835         mIndex = 0;
836     }
837 
838     Movement& movement = mMovements[mIndex];
839     movement.eventTime = eventTime;
840     movement.idBits = idBits;
841     uint32_t count = idBits.count();
842     for (uint32_t i = 0; i < count; i++) {
843         movement.positions[i] = positions[i];
844     }
845 }
846 
getEstimator(uint32_t id,VelocityTracker::Estimator * outEstimator) const847 bool LegacyVelocityTrackerStrategy::getEstimator(uint32_t id,
848         VelocityTracker::Estimator* outEstimator) const {
849     outEstimator->clear();
850 
851     const Movement& newestMovement = mMovements[mIndex];
852     if (!newestMovement.idBits.hasBit(id)) {
853         return false; // no data
854     }
855 
856     // Find the oldest sample that contains the pointer and that is not older than HORIZON.
857     nsecs_t minTime = newestMovement.eventTime - HORIZON;
858     uint32_t oldestIndex = mIndex;
859     uint32_t numTouches = 1;
860     do {
861         uint32_t nextOldestIndex = (oldestIndex == 0 ? HISTORY_SIZE : oldestIndex) - 1;
862         const Movement& nextOldestMovement = mMovements[nextOldestIndex];
863         if (!nextOldestMovement.idBits.hasBit(id)
864                 || nextOldestMovement.eventTime < minTime) {
865             break;
866         }
867         oldestIndex = nextOldestIndex;
868     } while (++numTouches < HISTORY_SIZE);
869 
870     // Calculate an exponentially weighted moving average of the velocity estimate
871     // at different points in time measured relative to the oldest sample.
872     // This is essentially an IIR filter.  Newer samples are weighted more heavily
873     // than older samples.  Samples at equal time points are weighted more or less
874     // equally.
875     //
876     // One tricky problem is that the sample data may be poorly conditioned.
877     // Sometimes samples arrive very close together in time which can cause us to
878     // overestimate the velocity at that time point.  Most samples might be measured
879     // 16ms apart but some consecutive samples could be only 0.5sm apart because
880     // the hardware or driver reports them irregularly or in bursts.
881     float accumVx = 0;
882     float accumVy = 0;
883     uint32_t index = oldestIndex;
884     uint32_t samplesUsed = 0;
885     const Movement& oldestMovement = mMovements[oldestIndex];
886     const VelocityTracker::Position& oldestPosition = oldestMovement.getPosition(id);
887     nsecs_t lastDuration = 0;
888 
889     while (numTouches-- > 1) {
890         if (++index == HISTORY_SIZE) {
891             index = 0;
892         }
893         const Movement& movement = mMovements[index];
894         nsecs_t duration = movement.eventTime - oldestMovement.eventTime;
895 
896         // If the duration between samples is small, we may significantly overestimate
897         // the velocity.  Consequently, we impose a minimum duration constraint on the
898         // samples that we include in the calculation.
899         if (duration >= MIN_DURATION) {
900             const VelocityTracker::Position& position = movement.getPosition(id);
901             float scale = 1000000000.0f / duration; // one over time delta in seconds
902             float vx = (position.x - oldestPosition.x) * scale;
903             float vy = (position.y - oldestPosition.y) * scale;
904             accumVx = (accumVx * lastDuration + vx * duration) / (duration + lastDuration);
905             accumVy = (accumVy * lastDuration + vy * duration) / (duration + lastDuration);
906             lastDuration = duration;
907             samplesUsed += 1;
908         }
909     }
910 
911     // Report velocity.
912     const VelocityTracker::Position& newestPosition = newestMovement.getPosition(id);
913     outEstimator->time = newestMovement.eventTime;
914     outEstimator->confidence = 1;
915     outEstimator->xCoeff[0] = newestPosition.x;
916     outEstimator->yCoeff[0] = newestPosition.y;
917     if (samplesUsed) {
918         outEstimator->xCoeff[1] = accumVx;
919         outEstimator->yCoeff[1] = accumVy;
920         outEstimator->degree = 1;
921     } else {
922         outEstimator->degree = 0;
923     }
924     return true;
925 }
926 
927 } // namespace android
928