widget/gonk/libui/VelocityTracker.cpp

Wed, 31 Dec 2014 07:22:50 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Wed, 31 Dec 2014 07:22:50 +0100
branch
TOR_BUG_3246
changeset 4
fc2d59ddac77
permissions
-rw-r--r--

Correct previous dual key logic pending first delivery installment.

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

mercurial