blob: 17cefbe836bf3add465603b442cc3ad67df2fa55 [file] [log] [blame]
Jeff Brown8a90e6e2012-05-11 12:24:35 -07001/*
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
Jeff Brown9eb7d862012-06-01 12:39:25 -070023// Log debug messages about the progress of the algorithm itself.
24#define DEBUG_STRATEGY 0
Jeff Brown8a90e6e2012-05-11 12:24:35 -070025
26#include <math.h>
27#include <limits.h>
28
29#include <androidfw/VelocityTracker.h>
30#include <utils/BitSet.h>
31#include <utils/String8.h>
32#include <utils/Timers.h>
33
Jeff Brown9eb7d862012-06-01 12:39:25 -070034#include <cutils/properties.h>
35
Jeff Brown8a90e6e2012-05-11 12:24:35 -070036namespace android {
37
Jeff Brown90729402012-05-14 18:46:18 -070038// Nanoseconds per milliseconds.
39static const nsecs_t NANOS_PER_MS = 1000000;
40
41// Threshold for determining that a pointer has stopped moving.
42// Some input devices do not send ACTION_MOVE events in the case where a pointer has
43// stopped. We need to detect this case so that we can accurately predict the
44// velocity after the pointer starts moving again.
45static const nsecs_t ASSUME_POINTER_STOPPED_TIME = 40 * NANOS_PER_MS;
46
47
Jeff Brown85bd0d62012-05-13 15:30:42 -070048static float vectorDot(const float* a, const float* b, uint32_t m) {
Jeff Brown8a90e6e2012-05-11 12:24:35 -070049 float r = 0;
50 while (m--) {
51 r += *(a++) * *(b++);
52 }
53 return r;
54}
55
Jeff Brown85bd0d62012-05-13 15:30:42 -070056static float vectorNorm(const float* a, uint32_t m) {
Jeff Brown8a90e6e2012-05-11 12:24:35 -070057 float r = 0;
58 while (m--) {
59 float t = *(a++);
60 r += t * t;
61 }
62 return sqrtf(r);
63}
64
Jeff Brown9eb7d862012-06-01 12:39:25 -070065#if DEBUG_STRATEGY || DEBUG_VELOCITY
Jeff Brown8a90e6e2012-05-11 12:24:35 -070066static String8 vectorToString(const float* a, uint32_t m) {
67 String8 str;
68 str.append("[");
69 while (m--) {
70 str.appendFormat(" %f", *(a++));
71 if (m) {
72 str.append(",");
73 }
74 }
75 str.append(" ]");
76 return str;
77}
78
79static String8 matrixToString(const float* a, uint32_t m, uint32_t n, bool rowMajor) {
80 String8 str;
81 str.append("[");
82 for (size_t i = 0; i < m; i++) {
83 if (i) {
84 str.append(",");
85 }
86 str.append(" [");
87 for (size_t j = 0; j < n; j++) {
88 if (j) {
89 str.append(",");
90 }
91 str.appendFormat(" %f", a[rowMajor ? i * n + j : j * m + i]);
92 }
93 str.append(" ]");
94 }
95 str.append(" ]");
96 return str;
97}
98#endif
99
Jeff Brown85bd0d62012-05-13 15:30:42 -0700100
101// --- VelocityTracker ---
102
Jeff Brown9eb7d862012-06-01 12:39:25 -0700103// The default velocity tracker strategy.
104// Although other strategies are available for testing and comparison purposes,
105// this is the strategy that applications will actually use. Be very careful
106// when adjusting the default strategy because it can dramatically affect
107// (often in a bad way) the user experience.
108const char* VelocityTracker::DEFAULT_STRATEGY = "lsq2";
109
110VelocityTracker::VelocityTracker(const char* strategy) :
111 mLastEventTime(0), mCurrentPointerIdBits(0), mActivePointerId(-1) {
112 char value[PROPERTY_VALUE_MAX];
113
114 // Allow the default strategy to be overridden using a system property for debugging.
115 if (!strategy) {
116 int length = property_get("debug.velocitytracker.strategy", value, NULL);
117 if (length > 0) {
118 strategy = value;
119 } else {
120 strategy = DEFAULT_STRATEGY;
121 }
122 }
123
124 // Configure the strategy.
125 if (!configureStrategy(strategy)) {
126 ALOGD("Unrecognized velocity tracker strategy name '%s'.", strategy);
127 if (!configureStrategy(DEFAULT_STRATEGY)) {
128 LOG_ALWAYS_FATAL("Could not create the default velocity tracker strategy '%s'!",
129 strategy);
130 }
131 }
Jeff Brown85bd0d62012-05-13 15:30:42 -0700132}
133
Jeff Brown85bd0d62012-05-13 15:30:42 -0700134VelocityTracker::~VelocityTracker() {
135 delete mStrategy;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700136}
137
Jeff Brown9eb7d862012-06-01 12:39:25 -0700138bool VelocityTracker::configureStrategy(const char* strategy) {
139 mStrategy = createStrategy(strategy);
140 return mStrategy != NULL;
141}
142
143VelocityTrackerStrategy* VelocityTracker::createStrategy(const char* strategy) {
144 if (!strcmp("lsq1", strategy)) {
145 // 1st order least squares. Quality: POOR.
146 // Frequently underfits the touch data especially when the finger accelerates
147 // or changes direction. Often underestimates velocity. The direction
148 // is overly influenced by historical touch points.
149 return new LeastSquaresVelocityTrackerStrategy(1);
150 }
151 if (!strcmp("lsq2", strategy)) {
152 // 2nd order least squares. Quality: VERY GOOD.
153 // Pretty much ideal, but can be confused by certain kinds of touch data,
154 // particularly if the panel has a tendency to generate delayed,
155 // duplicate or jittery touch coordinates when the finger is released.
156 return new LeastSquaresVelocityTrackerStrategy(2);
157 }
158 if (!strcmp("lsq3", strategy)) {
159 // 3rd order least squares. Quality: UNUSABLE.
160 // Frequently overfits the touch data yielding wildly divergent estimates
161 // of the velocity when the finger is released.
162 return new LeastSquaresVelocityTrackerStrategy(3);
163 }
Jeff Brown18f329e2012-06-03 14:18:26 -0700164 if (!strcmp("wlsq2-delta", strategy)) {
165 // 2nd order weighted least squares, delta weighting. Quality: EXPERIMENTAL
166 return new LeastSquaresVelocityTrackerStrategy(2,
167 LeastSquaresVelocityTrackerStrategy::WEIGHTING_DELTA);
168 }
169 if (!strcmp("wlsq2-central", strategy)) {
170 // 2nd order weighted least squares, central weighting. Quality: EXPERIMENTAL
171 return new LeastSquaresVelocityTrackerStrategy(2,
172 LeastSquaresVelocityTrackerStrategy::WEIGHTING_CENTRAL);
173 }
174 if (!strcmp("wlsq2-recent", strategy)) {
175 // 2nd order weighted least squares, recent weighting. Quality: EXPERIMENTAL
176 return new LeastSquaresVelocityTrackerStrategy(2,
177 LeastSquaresVelocityTrackerStrategy::WEIGHTING_RECENT);
178 }
Jeff Brown53dd12a2012-06-01 13:24:04 -0700179 if (!strcmp("int1", strategy)) {
180 // 1st order integrating filter. Quality: GOOD.
181 // Not as good as 'lsq2' because it cannot estimate acceleration but it is
182 // more tolerant of errors. Like 'lsq1', this strategy tends to underestimate
183 // the velocity of a fling but this strategy tends to respond to changes in
184 // direction more quickly and accurately.
185 return new IntegratingVelocityTrackerStrategy();
186 }
Jeff Brown9eb7d862012-06-01 12:39:25 -0700187 return NULL;
188}
189
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700190void VelocityTracker::clear() {
Jeff Brown85bd0d62012-05-13 15:30:42 -0700191 mCurrentPointerIdBits.clear();
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700192 mActivePointerId = -1;
Jeff Brown85bd0d62012-05-13 15:30:42 -0700193
194 mStrategy->clear();
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700195}
196
197void VelocityTracker::clearPointers(BitSet32 idBits) {
Jeff Brown85bd0d62012-05-13 15:30:42 -0700198 BitSet32 remainingIdBits(mCurrentPointerIdBits.value & ~idBits.value);
199 mCurrentPointerIdBits = remainingIdBits;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700200
201 if (mActivePointerId >= 0 && idBits.hasBit(mActivePointerId)) {
202 mActivePointerId = !remainingIdBits.isEmpty() ? remainingIdBits.firstMarkedBit() : -1;
203 }
Jeff Brown85bd0d62012-05-13 15:30:42 -0700204
205 mStrategy->clearPointers(idBits);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700206}
207
208void VelocityTracker::addMovement(nsecs_t eventTime, BitSet32 idBits, const Position* positions) {
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700209 while (idBits.count() > MAX_POINTERS) {
210 idBits.clearLastMarkedBit();
211 }
212
Jeff Brown90729402012-05-14 18:46:18 -0700213 if ((mCurrentPointerIdBits.value & idBits.value)
214 && eventTime >= mLastEventTime + ASSUME_POINTER_STOPPED_TIME) {
215#if DEBUG_VELOCITY
216 ALOGD("VelocityTracker: stopped for %0.3f ms, clearing state.",
217 (eventTime - mLastEventTime) * 0.000001f);
218#endif
219 // We have not received any movements for too long. Assume that all pointers
220 // have stopped.
221 mStrategy->clear();
222 }
223 mLastEventTime = eventTime;
224
Jeff Brown85bd0d62012-05-13 15:30:42 -0700225 mCurrentPointerIdBits = idBits;
226 if (mActivePointerId < 0 || !idBits.hasBit(mActivePointerId)) {
227 mActivePointerId = idBits.isEmpty() ? -1 : idBits.firstMarkedBit();
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700228 }
229
Jeff Brown85bd0d62012-05-13 15:30:42 -0700230 mStrategy->addMovement(eventTime, idBits, positions);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700231
232#if DEBUG_VELOCITY
233 ALOGD("VelocityTracker: addMovement eventTime=%lld, idBits=0x%08x, activePointerId=%d",
234 eventTime, idBits.value, mActivePointerId);
235 for (BitSet32 iterBits(idBits); !iterBits.isEmpty(); ) {
236 uint32_t id = iterBits.firstMarkedBit();
237 uint32_t index = idBits.getIndexOfBit(id);
238 iterBits.clearBit(id);
239 Estimator estimator;
Jeff Brown85bd0d62012-05-13 15:30:42 -0700240 getEstimator(id, &estimator);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700241 ALOGD(" %d: position (%0.3f, %0.3f), "
242 "estimator (degree=%d, xCoeff=%s, yCoeff=%s, confidence=%f)",
243 id, positions[index].x, positions[index].y,
244 int(estimator.degree),
Jeff Browndcab1902012-05-14 18:44:17 -0700245 vectorToString(estimator.xCoeff, estimator.degree + 1).string(),
246 vectorToString(estimator.yCoeff, estimator.degree + 1).string(),
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700247 estimator.confidence);
248 }
249#endif
250}
251
252void VelocityTracker::addMovement(const MotionEvent* event) {
253 int32_t actionMasked = event->getActionMasked();
254
255 switch (actionMasked) {
256 case AMOTION_EVENT_ACTION_DOWN:
257 case AMOTION_EVENT_ACTION_HOVER_ENTER:
258 // Clear all pointers on down before adding the new movement.
259 clear();
260 break;
261 case AMOTION_EVENT_ACTION_POINTER_DOWN: {
262 // Start a new movement trace for a pointer that just went down.
263 // We do this on down instead of on up because the client may want to query the
264 // final velocity for a pointer that just went up.
265 BitSet32 downIdBits;
266 downIdBits.markBit(event->getPointerId(event->getActionIndex()));
267 clearPointers(downIdBits);
268 break;
269 }
270 case AMOTION_EVENT_ACTION_MOVE:
271 case AMOTION_EVENT_ACTION_HOVER_MOVE:
272 break;
273 default:
274 // Ignore all other actions because they do not convey any new information about
275 // pointer movement. We also want to preserve the last known velocity of the pointers.
276 // Note that ACTION_UP and ACTION_POINTER_UP always report the last known position
277 // of the pointers that went up. ACTION_POINTER_UP does include the new position of
278 // pointers that remained down but we will also receive an ACTION_MOVE with this
279 // information if any of them actually moved. Since we don't know how many pointers
280 // will be going up at once it makes sense to just wait for the following ACTION_MOVE
281 // before adding the movement.
282 return;
283 }
284
285 size_t pointerCount = event->getPointerCount();
286 if (pointerCount > MAX_POINTERS) {
287 pointerCount = MAX_POINTERS;
288 }
289
290 BitSet32 idBits;
291 for (size_t i = 0; i < pointerCount; i++) {
292 idBits.markBit(event->getPointerId(i));
293 }
294
Jeff Browndcab1902012-05-14 18:44:17 -0700295 uint32_t pointerIndex[MAX_POINTERS];
296 for (size_t i = 0; i < pointerCount; i++) {
297 pointerIndex[i] = idBits.getIndexOfBit(event->getPointerId(i));
298 }
299
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700300 nsecs_t eventTime;
301 Position positions[pointerCount];
302
303 size_t historySize = event->getHistorySize();
304 for (size_t h = 0; h < historySize; h++) {
305 eventTime = event->getHistoricalEventTime(h);
306 for (size_t i = 0; i < pointerCount; i++) {
Jeff Browndcab1902012-05-14 18:44:17 -0700307 uint32_t index = pointerIndex[i];
308 positions[index].x = event->getHistoricalX(i, h);
309 positions[index].y = event->getHistoricalY(i, h);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700310 }
311 addMovement(eventTime, idBits, positions);
312 }
313
314 eventTime = event->getEventTime();
315 for (size_t i = 0; i < pointerCount; i++) {
Jeff Browndcab1902012-05-14 18:44:17 -0700316 uint32_t index = pointerIndex[i];
317 positions[index].x = event->getX(i);
318 positions[index].y = event->getY(i);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700319 }
320 addMovement(eventTime, idBits, positions);
321}
322
Jeff Brown85bd0d62012-05-13 15:30:42 -0700323bool VelocityTracker::getVelocity(uint32_t id, float* outVx, float* outVy) const {
324 Estimator estimator;
325 if (getEstimator(id, &estimator) && estimator.degree >= 1) {
326 *outVx = estimator.xCoeff[1];
327 *outVy = estimator.yCoeff[1];
328 return true;
329 }
330 *outVx = 0;
331 *outVy = 0;
332 return false;
333}
334
335bool VelocityTracker::getEstimator(uint32_t id, Estimator* outEstimator) const {
336 return mStrategy->getEstimator(id, outEstimator);
337}
338
339
340// --- LeastSquaresVelocityTrackerStrategy ---
341
Jeff Brown85bd0d62012-05-13 15:30:42 -0700342const nsecs_t LeastSquaresVelocityTrackerStrategy::HORIZON;
343const uint32_t LeastSquaresVelocityTrackerStrategy::HISTORY_SIZE;
344
Jeff Brown18f329e2012-06-03 14:18:26 -0700345LeastSquaresVelocityTrackerStrategy::LeastSquaresVelocityTrackerStrategy(
346 uint32_t degree, Weighting weighting) :
347 mDegree(degree), mWeighting(weighting) {
Jeff Brown85bd0d62012-05-13 15:30:42 -0700348 clear();
349}
350
351LeastSquaresVelocityTrackerStrategy::~LeastSquaresVelocityTrackerStrategy() {
352}
353
354void LeastSquaresVelocityTrackerStrategy::clear() {
355 mIndex = 0;
356 mMovements[0].idBits.clear();
357}
358
359void LeastSquaresVelocityTrackerStrategy::clearPointers(BitSet32 idBits) {
360 BitSet32 remainingIdBits(mMovements[mIndex].idBits.value & ~idBits.value);
361 mMovements[mIndex].idBits = remainingIdBits;
362}
363
364void LeastSquaresVelocityTrackerStrategy::addMovement(nsecs_t eventTime, BitSet32 idBits,
365 const VelocityTracker::Position* positions) {
366 if (++mIndex == HISTORY_SIZE) {
367 mIndex = 0;
368 }
369
370 Movement& movement = mMovements[mIndex];
371 movement.eventTime = eventTime;
372 movement.idBits = idBits;
373 uint32_t count = idBits.count();
374 for (uint32_t i = 0; i < count; i++) {
375 movement.positions[i] = positions[i];
376 }
377}
378
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700379/**
380 * Solves a linear least squares problem to obtain a N degree polynomial that fits
381 * the specified input data as nearly as possible.
382 *
383 * Returns true if a solution is found, false otherwise.
384 *
Jeff Brown18f329e2012-06-03 14:18:26 -0700385 * The input consists of two vectors of data points X and Y with indices 0..m-1
386 * along with a weight vector W of the same size.
387 *
Jeff Brown9eb7d862012-06-01 12:39:25 -0700388 * The output is a vector B with indices 0..n that describes a polynomial
Jeff Brown18f329e2012-06-03 14:18:26 -0700389 * that fits the data, such the sum of W[i] * W[i] * abs(Y[i] - (B[0] + B[1] X[i]
390 * + B[2] X[i]^2 ... B[n] X[i]^n)) for all i between 0 and m-1 is minimized.
391 *
392 * Accordingly, the weight vector W should be initialized by the caller with the
393 * reciprocal square root of the variance of the error in each input data point.
394 * In other words, an ideal choice for W would be W[i] = 1 / var(Y[i]) = 1 / stddev(Y[i]).
395 * The weights express the relative importance of each data point. If the weights are
396 * all 1, then the data points are considered to be of equal importance when fitting
397 * the polynomial. It is a good idea to choose weights that diminish the importance
398 * of data points that may have higher than usual error margins.
399 *
400 * Errors among data points are assumed to be independent. W is represented here
401 * as a vector although in the literature it is typically taken to be a diagonal matrix.
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700402 *
403 * That is to say, the function that generated the input data can be approximated
404 * by y(x) ~= B[0] + B[1] x + B[2] x^2 + ... + B[n] x^n.
405 *
406 * The coefficient of determination (R^2) is also returned to describe the goodness
407 * of fit of the model for the given data. It is a value between 0 and 1, where 1
408 * indicates perfect correspondence.
409 *
410 * This function first expands the X vector to a m by n matrix A such that
Jeff Brown18f329e2012-06-03 14:18:26 -0700411 * A[i][0] = 1, A[i][1] = X[i], A[i][2] = X[i]^2, ..., A[i][n] = X[i]^n, then
412 * multiplies it by w[i]./
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700413 *
414 * Then it calculates the QR decomposition of A yielding an m by m orthonormal matrix Q
415 * and an m by n upper triangular matrix R. Because R is upper triangular (lower
416 * part is all zeroes), we can simplify the decomposition into an m by n matrix
417 * Q1 and a n by n matrix R1 such that A = Q1 R1.
418 *
Jeff Brown18f329e2012-06-03 14:18:26 -0700419 * Finally we solve the system of linear equations given by R1 B = (Qtranspose W Y)
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700420 * to find B.
421 *
422 * For efficiency, we lay out A and Q column-wise in memory because we frequently
423 * operate on the column vectors. Conversely, we lay out R row-wise.
424 *
425 * http://en.wikipedia.org/wiki/Numerical_methods_for_linear_least_squares
426 * http://en.wikipedia.org/wiki/Gram-Schmidt
427 */
Jeff Brown18f329e2012-06-03 14:18:26 -0700428static bool solveLeastSquares(const float* x, const float* y,
429 const float* w, uint32_t m, uint32_t n, float* outB, float* outDet) {
Jeff Brown9eb7d862012-06-01 12:39:25 -0700430#if DEBUG_STRATEGY
Jeff Brown18f329e2012-06-03 14:18:26 -0700431 ALOGD("solveLeastSquares: m=%d, n=%d, x=%s, y=%s, w=%s", int(m), int(n),
432 vectorToString(x, m).string(), vectorToString(y, m).string(),
433 vectorToString(w, m).string());
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700434#endif
435
Jeff Brown18f329e2012-06-03 14:18:26 -0700436 // Expand the X vector to a matrix A, pre-multiplied by the weights.
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700437 float a[n][m]; // column-major order
438 for (uint32_t h = 0; h < m; h++) {
Jeff Brown18f329e2012-06-03 14:18:26 -0700439 a[0][h] = w[h];
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700440 for (uint32_t i = 1; i < n; i++) {
441 a[i][h] = a[i - 1][h] * x[h];
442 }
443 }
Jeff Brown9eb7d862012-06-01 12:39:25 -0700444#if DEBUG_STRATEGY
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700445 ALOGD(" - a=%s", matrixToString(&a[0][0], m, n, false /*rowMajor*/).string());
446#endif
447
448 // Apply the Gram-Schmidt process to A to obtain its QR decomposition.
449 float q[n][m]; // orthonormal basis, column-major order
450 float r[n][n]; // upper triangular matrix, row-major order
451 for (uint32_t j = 0; j < n; j++) {
452 for (uint32_t h = 0; h < m; h++) {
453 q[j][h] = a[j][h];
454 }
455 for (uint32_t i = 0; i < j; i++) {
456 float dot = vectorDot(&q[j][0], &q[i][0], m);
457 for (uint32_t h = 0; h < m; h++) {
458 q[j][h] -= dot * q[i][h];
459 }
460 }
461
462 float norm = vectorNorm(&q[j][0], m);
463 if (norm < 0.000001f) {
464 // vectors are linearly dependent or zero so no solution
Jeff Brown9eb7d862012-06-01 12:39:25 -0700465#if DEBUG_STRATEGY
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700466 ALOGD(" - no solution, norm=%f", norm);
467#endif
468 return false;
469 }
470
471 float invNorm = 1.0f / norm;
472 for (uint32_t h = 0; h < m; h++) {
473 q[j][h] *= invNorm;
474 }
475 for (uint32_t i = 0; i < n; i++) {
476 r[j][i] = i < j ? 0 : vectorDot(&q[j][0], &a[i][0], m);
477 }
478 }
Jeff Brown9eb7d862012-06-01 12:39:25 -0700479#if DEBUG_STRATEGY
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700480 ALOGD(" - q=%s", matrixToString(&q[0][0], m, n, false /*rowMajor*/).string());
481 ALOGD(" - r=%s", matrixToString(&r[0][0], n, n, true /*rowMajor*/).string());
482
483 // calculate QR, if we factored A correctly then QR should equal A
484 float qr[n][m];
485 for (uint32_t h = 0; h < m; h++) {
486 for (uint32_t i = 0; i < n; i++) {
487 qr[i][h] = 0;
488 for (uint32_t j = 0; j < n; j++) {
489 qr[i][h] += q[j][h] * r[j][i];
490 }
491 }
492 }
493 ALOGD(" - qr=%s", matrixToString(&qr[0][0], m, n, false /*rowMajor*/).string());
494#endif
495
Jeff Brown18f329e2012-06-03 14:18:26 -0700496 // Solve R B = Qt W Y to find B. This is easy because R is upper triangular.
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700497 // We just work from bottom-right to top-left calculating B's coefficients.
Jeff Brown18f329e2012-06-03 14:18:26 -0700498 float wy[m];
499 for (uint32_t h = 0; h < m; h++) {
500 wy[h] = y[h] * w[h];
501 }
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700502 for (uint32_t i = n; i-- != 0; ) {
Jeff Brown18f329e2012-06-03 14:18:26 -0700503 outB[i] = vectorDot(&q[i][0], wy, m);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700504 for (uint32_t j = n - 1; j > i; j--) {
505 outB[i] -= r[i][j] * outB[j];
506 }
507 outB[i] /= r[i][i];
508 }
Jeff Brown9eb7d862012-06-01 12:39:25 -0700509#if DEBUG_STRATEGY
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700510 ALOGD(" - b=%s", vectorToString(outB, n).string());
511#endif
512
513 // Calculate the coefficient of determination as 1 - (SSerr / SStot) where
Jeff Brown18f329e2012-06-03 14:18:26 -0700514 // SSerr is the residual sum of squares (variance of the error),
515 // and SStot is the total sum of squares (variance of the data) where each
516 // has been weighted.
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700517 float ymean = 0;
518 for (uint32_t h = 0; h < m; h++) {
519 ymean += y[h];
520 }
521 ymean /= m;
522
523 float sserr = 0;
524 float sstot = 0;
525 for (uint32_t h = 0; h < m; h++) {
526 float err = y[h] - outB[0];
527 float term = 1;
528 for (uint32_t i = 1; i < n; i++) {
529 term *= x[h];
530 err -= term * outB[i];
531 }
Jeff Brown18f329e2012-06-03 14:18:26 -0700532 sserr += w[h] * w[h] * err * err;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700533 float var = y[h] - ymean;
Jeff Brown18f329e2012-06-03 14:18:26 -0700534 sstot += w[h] * w[h] * var * var;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700535 }
536 *outDet = sstot > 0.000001f ? 1.0f - (sserr / sstot) : 1;
Jeff Brown9eb7d862012-06-01 12:39:25 -0700537#if DEBUG_STRATEGY
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700538 ALOGD(" - sserr=%f", sserr);
539 ALOGD(" - sstot=%f", sstot);
540 ALOGD(" - det=%f", *outDet);
541#endif
542 return true;
543}
544
Jeff Brown85bd0d62012-05-13 15:30:42 -0700545bool LeastSquaresVelocityTrackerStrategy::getEstimator(uint32_t id,
546 VelocityTracker::Estimator* outEstimator) const {
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700547 outEstimator->clear();
548
549 // Iterate over movement samples in reverse time order and collect samples.
550 float x[HISTORY_SIZE];
551 float y[HISTORY_SIZE];
Jeff Brown18f329e2012-06-03 14:18:26 -0700552 float w[HISTORY_SIZE];
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700553 float time[HISTORY_SIZE];
554 uint32_t m = 0;
555 uint32_t index = mIndex;
556 const Movement& newestMovement = mMovements[mIndex];
557 do {
558 const Movement& movement = mMovements[index];
559 if (!movement.idBits.hasBit(id)) {
560 break;
561 }
562
563 nsecs_t age = newestMovement.eventTime - movement.eventTime;
Jeff Brown85bd0d62012-05-13 15:30:42 -0700564 if (age > HORIZON) {
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700565 break;
566 }
567
Jeff Brown85bd0d62012-05-13 15:30:42 -0700568 const VelocityTracker::Position& position = movement.getPosition(id);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700569 x[m] = position.x;
570 y[m] = position.y;
Jeff Brown18f329e2012-06-03 14:18:26 -0700571 w[m] = chooseWeight(index);
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700572 time[m] = -age * 0.000000001f;
573 index = (index == 0 ? HISTORY_SIZE : index) - 1;
574 } while (++m < HISTORY_SIZE);
575
576 if (m == 0) {
577 return false; // no data
578 }
579
580 // Calculate a least squares polynomial fit.
Jeff Brown9eb7d862012-06-01 12:39:25 -0700581 uint32_t degree = mDegree;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700582 if (degree > m - 1) {
583 degree = m - 1;
584 }
585 if (degree >= 1) {
586 float xdet, ydet;
587 uint32_t n = degree + 1;
Jeff Brown18f329e2012-06-03 14:18:26 -0700588 if (solveLeastSquares(time, x, w, m, n, outEstimator->xCoeff, &xdet)
589 && solveLeastSquares(time, y, w, m, n, outEstimator->yCoeff, &ydet)) {
Jeff Brown90729402012-05-14 18:46:18 -0700590 outEstimator->time = newestMovement.eventTime;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700591 outEstimator->degree = degree;
592 outEstimator->confidence = xdet * ydet;
Jeff Brown9eb7d862012-06-01 12:39:25 -0700593#if DEBUG_STRATEGY
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700594 ALOGD("estimate: degree=%d, xCoeff=%s, yCoeff=%s, confidence=%f",
595 int(outEstimator->degree),
596 vectorToString(outEstimator->xCoeff, n).string(),
597 vectorToString(outEstimator->yCoeff, n).string(),
598 outEstimator->confidence);
599#endif
600 return true;
601 }
602 }
603
604 // No velocity data available for this pointer, but we do have its current position.
605 outEstimator->xCoeff[0] = x[0];
606 outEstimator->yCoeff[0] = y[0];
Jeff Brown90729402012-05-14 18:46:18 -0700607 outEstimator->time = newestMovement.eventTime;
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700608 outEstimator->degree = 0;
609 outEstimator->confidence = 1;
610 return true;
611}
612
Jeff Brown18f329e2012-06-03 14:18:26 -0700613float LeastSquaresVelocityTrackerStrategy::chooseWeight(uint32_t index) const {
614 switch (mWeighting) {
615 case WEIGHTING_DELTA: {
616 // Weight points based on how much time elapsed between them and the next
617 // point so that points that "cover" a shorter time span are weighed less.
618 // delta 0ms: 0.5
619 // delta 10ms: 1.0
620 if (index == mIndex) {
621 return 1.0f;
622 }
623 uint32_t nextIndex = (index + 1) % HISTORY_SIZE;
624 float deltaMillis = (mMovements[nextIndex].eventTime- mMovements[index].eventTime)
625 * 0.000001f;
626 if (deltaMillis < 0) {
627 return 0.5f;
628 }
629 if (deltaMillis < 10) {
630 return 0.5f + deltaMillis * 0.05;
631 }
632 return 1.0f;
633 }
634
635 case WEIGHTING_CENTRAL: {
636 // Weight points based on their age, weighing very recent and very old points less.
637 // age 0ms: 0.5
638 // age 10ms: 1.0
639 // age 50ms: 1.0
640 // age 60ms: 0.5
641 float ageMillis = (mMovements[mIndex].eventTime - mMovements[index].eventTime)
642 * 0.000001f;
643 if (ageMillis < 0) {
644 return 0.5f;
645 }
646 if (ageMillis < 10) {
647 return 0.5f + ageMillis * 0.05;
648 }
649 if (ageMillis < 50) {
650 return 1.0f;
651 }
652 if (ageMillis < 60) {
653 return 0.5f + (60 - ageMillis) * 0.05;
654 }
655 return 0.5f;
656 }
657
658 case WEIGHTING_RECENT: {
659 // Weight points based on their age, weighing older points less.
660 // age 0ms: 1.0
661 // age 50ms: 1.0
662 // age 100ms: 0.5
663 float ageMillis = (mMovements[mIndex].eventTime - mMovements[index].eventTime)
664 * 0.000001f;
665 if (ageMillis < 50) {
666 return 1.0f;
667 }
668 if (ageMillis < 100) {
669 return 0.5f + (100 - ageMillis) * 0.01f;
670 }
671 return 0.5f;
672 }
673
674 case WEIGHTING_NONE:
675 default:
676 return 1.0f;
677 }
678}
679
Jeff Brown53dd12a2012-06-01 13:24:04 -0700680
681// --- IntegratingVelocityTrackerStrategy ---
682
683IntegratingVelocityTrackerStrategy::IntegratingVelocityTrackerStrategy() {
684}
685
686IntegratingVelocityTrackerStrategy::~IntegratingVelocityTrackerStrategy() {
687}
688
689void IntegratingVelocityTrackerStrategy::clear() {
690 mPointerIdBits.clear();
691}
692
693void IntegratingVelocityTrackerStrategy::clearPointers(BitSet32 idBits) {
694 mPointerIdBits.value &= ~idBits.value;
695}
696
697void IntegratingVelocityTrackerStrategy::addMovement(nsecs_t eventTime, BitSet32 idBits,
698 const VelocityTracker::Position* positions) {
699 uint32_t index = 0;
700 for (BitSet32 iterIdBits(idBits); !iterIdBits.isEmpty();) {
701 uint32_t id = iterIdBits.clearFirstMarkedBit();
702 State& state = mPointerState[id];
703 const VelocityTracker::Position& position = positions[index++];
704 if (mPointerIdBits.hasBit(id)) {
705 updateState(state, eventTime, position.x, position.y);
706 } else {
707 initState(state, eventTime, position.x, position.y);
708 }
709 }
710
711 mPointerIdBits = idBits;
712}
713
714bool IntegratingVelocityTrackerStrategy::getEstimator(uint32_t id,
715 VelocityTracker::Estimator* outEstimator) const {
716 outEstimator->clear();
717
718 if (mPointerIdBits.hasBit(id)) {
719 const State& state = mPointerState[id];
720 populateEstimator(state, outEstimator);
721 return true;
722 }
723
724 return false;
725}
726
727void IntegratingVelocityTrackerStrategy::initState(State& state,
728 nsecs_t eventTime, float xpos, float ypos) {
729 state.updateTime = eventTime;
730 state.first = true;
731
732 state.xpos = xpos;
733 state.xvel = 0;
734 state.ypos = ypos;
735 state.yvel = 0;
736}
737
738void IntegratingVelocityTrackerStrategy::updateState(State& state,
739 nsecs_t eventTime, float xpos, float ypos) {
740 const nsecs_t MIN_TIME_DELTA = 2 * NANOS_PER_MS;
741 const float FILTER_TIME_CONSTANT = 0.010f; // 10 milliseconds
742
743 if (eventTime <= state.updateTime + MIN_TIME_DELTA) {
744 return;
745 }
746
747 float dt = (eventTime - state.updateTime) * 0.000000001f;
748 state.updateTime = eventTime;
749
750 float xvel = (xpos - state.xpos) / dt;
751 float yvel = (ypos - state.ypos) / dt;
752 if (state.first) {
753 state.xvel = xvel;
754 state.yvel = yvel;
755 state.first = false;
756 } else {
757 float alpha = dt / (FILTER_TIME_CONSTANT + dt);
758 state.xvel += (xvel - state.xvel) * alpha;
759 state.yvel += (yvel - state.yvel) * alpha;
760 }
761 state.xpos = xpos;
762 state.ypos = ypos;
763}
764
765void IntegratingVelocityTrackerStrategy::populateEstimator(const State& state,
766 VelocityTracker::Estimator* outEstimator) {
767 outEstimator->time = state.updateTime;
768 outEstimator->degree = 1;
769 outEstimator->confidence = 1.0f;
770 outEstimator->xCoeff[0] = state.xpos;
771 outEstimator->xCoeff[1] = state.xvel;
772 outEstimator->yCoeff[0] = state.ypos;
773 outEstimator->yCoeff[1] = state.yvel;
774}
775
Jeff Brown8a90e6e2012-05-11 12:24:35 -0700776} // namespace android