Jeff Brown | 5912f95 | 2013-07-01 19:10:31 -0700 | [diff] [blame] | 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 | |
| 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 | |
| 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 |
| 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 | |
| 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 | |
| 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 | |
| 133 | VelocityTracker::~VelocityTracker() { |
| 134 | delete mStrategy; |
| 135 | } |
| 136 | |
| 137 | bool VelocityTracker::configureStrategy(const char* strategy) { |
| 138 | mStrategy = createStrategy(strategy); |
| 139 | return mStrategy != NULL; |
| 140 | } |
| 141 | |
| 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 | |
| 202 | void VelocityTracker::clear() { |
| 203 | mCurrentPointerIdBits.clear(); |
| 204 | mActivePointerId = -1; |
| 205 | |
| 206 | mStrategy->clear(); |
| 207 | } |
| 208 | |
| 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 | |
| 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 | |
| 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 | |
| 335 | 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 | |
| 347 | 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 | |
| 357 | LeastSquaresVelocityTrackerStrategy::LeastSquaresVelocityTrackerStrategy( |
| 358 | uint32_t degree, Weighting weighting) : |
| 359 | mDegree(degree), mWeighting(weighting) { |
| 360 | clear(); |
| 361 | } |
| 362 | |
| 363 | LeastSquaresVelocityTrackerStrategy::~LeastSquaresVelocityTrackerStrategy() { |
| 364 | } |
| 365 | |
| 366 | void LeastSquaresVelocityTrackerStrategy::clear() { |
| 367 | mIndex = 0; |
| 368 | mMovements[0].idBits.clear(); |
| 369 | } |
| 370 | |
| 371 | void LeastSquaresVelocityTrackerStrategy::clearPointers(BitSet32 idBits) { |
| 372 | BitSet32 remainingIdBits(mMovements[mIndex].idBits.value & ~idBits.value); |
| 373 | mMovements[mIndex].idBits = remainingIdBits; |
| 374 | } |
| 375 | |
| 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 | */ |
| 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 | |
| 557 | 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 | |
| 625 | 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 | |
| 695 | IntegratingVelocityTrackerStrategy::IntegratingVelocityTrackerStrategy(uint32_t degree) : |
| 696 | mDegree(degree) { |
| 697 | } |
| 698 | |
| 699 | IntegratingVelocityTrackerStrategy::~IntegratingVelocityTrackerStrategy() { |
| 700 | } |
| 701 | |
| 702 | void IntegratingVelocityTrackerStrategy::clear() { |
| 703 | mPointerIdBits.clear(); |
| 704 | } |
| 705 | |
| 706 | void IntegratingVelocityTrackerStrategy::clearPointers(BitSet32 idBits) { |
| 707 | mPointerIdBits.value &= ~idBits.value; |
| 708 | } |
| 709 | |
| 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 | |
| 727 | 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 | |
| 740 | 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 | |
| 753 | 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 | |
| 795 | 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 | |
| 815 | LegacyVelocityTrackerStrategy::LegacyVelocityTrackerStrategy() { |
| 816 | clear(); |
| 817 | } |
| 818 | |
| 819 | LegacyVelocityTrackerStrategy::~LegacyVelocityTrackerStrategy() { |
| 820 | } |
| 821 | |
| 822 | void LegacyVelocityTrackerStrategy::clear() { |
| 823 | mIndex = 0; |
| 824 | mMovements[0].idBits.clear(); |
| 825 | } |
| 826 | |
| 827 | void LegacyVelocityTrackerStrategy::clearPointers(BitSet32 idBits) { |
| 828 | BitSet32 remainingIdBits(mMovements[mIndex].idBits.value & ~idBits.value); |
| 829 | mMovements[mIndex].idBits = remainingIdBits; |
| 830 | } |
| 831 | |
| 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 | |
| 847 | 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 |