void TMKFeatureVectors::findPairOffsetsModuloPeriods()

in tmk/cpp/algo/tmkfv.cpp [570:676]

• 76 lines of code
• 12 McCabe index (conditional complexity)

void TMKFeatureVectors::findPairOffsetsModuloPeriods(
    const TMKFeatureVectors& fva,
    const TMKFeatureVectors& fvb,
    BestOffsets& bestOffsets,
    ValuesAtBestOffsets& valuesAtBestOffsets,
    bool printDetails) {
  int numPeriods = fva._periods.size();
  int numFourierCoefficients = fva._fourierCoefficients.size();

  // We assume areCompatible has already been called on fva and fvb.
  bestOffsets = BestOffsets(numPeriods);
  valuesAtBestOffsets = ValuesAtBestOffsets(numPeriods);

  if (numPeriods == 0 || numFourierCoefficients == 0) {
    return;
  }

  //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // We assume all cos/sin features are already L2-normalized.

  //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  std::vector<std::vector<float>> dotCosCos =
      facebook::tmk::libvec::allocateRank2(numPeriods, numFourierCoefficients);
  std::vector<std::vector<float>> dotSinSin =
      facebook::tmk::libvec::allocateRank2(numPeriods, numFourierCoefficients);
  std::vector<std::vector<float>> dotSinCos =
      facebook::tmk::libvec::allocateRank2(numPeriods, numFourierCoefficients);
  std::vector<std::vector<float>> dotCosSin =
      facebook::tmk::libvec::allocateRank2(numPeriods, numFourierCoefficients);

  for (int i = 0; i < numPeriods; i++) {
    for (int j = 0; j < numFourierCoefficients; j++) {
      const std::vector<float>& fvaCos = fva._cosFeatures[i][j];
      const std::vector<float>& fvaSin = fva._sinFeatures[i][j];
      const std::vector<float>& fvbCos = fvb._cosFeatures[i][j];
      const std::vector<float>& fvbSin = fvb._sinFeatures[i][j];
      dotCosCos[i][j] = facebook::tmk::libvec::computeDot(fvaCos, fvbCos);
      dotSinSin[i][j] = facebook::tmk::libvec::computeDot(fvaSin, fvbSin);
      dotSinCos[i][j] = facebook::tmk::libvec::computeDot(fvaSin, fvbCos);
      dotCosSin[i][j] = facebook::tmk::libvec::computeDot(fvaCos, fvbSin);
    }
  }

  //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  for (int i = 0; i < numPeriods; i++) {
    int period = fva._periods[i];
    std::vector<float> K_deltas(period);
    for (int offset = 0; offset < period; offset++) {
      float delta = 2.0 * M_PI * (float)offset / (float)period;
      float K_delta = dotCosCos[i][0];
#if 0
      // Direct evaluation
      for (int j = 1; j < numFourierCoefficients; j++) {
        K_delta += cos(j * delta) * (dotCosCos[i][j] + dotSinSin[i][j]);
        K_delta += sin(j * delta) * (dotSinCos[i][j] - dotCosSin[i][j]);
      }
#else
      /**
       * We noticed that sin and cos were the dominant cost to this function.
       * We can try using an incremental algorithm for equally spaced exp(1j)
       * as recommended in Numerical Recipes. A good demonstration of precision:
       * http://steve.hollasch.net/cgindex/math/inccos.html
       * this is based on the addition formulae.  We replace cos(d) with
       * sin(d/2) to avoid truncation error using the double angle formula.
       *
       * Initial tests show 1/3 improvement in performance.
       **/

      float cos_jd = 1;
      float sin_jd = 0;

      float beta = sin(delta);
      float alpha = sin(delta / 2);
      alpha = 2 * alpha * alpha;

      for (int j = 1; j < numFourierCoefficients; j++) {
        // a bit magical because we're starting at j = 1
        float n_cos_jd = (alpha * cos_jd) + (beta * sin_jd);
        float n_sin_jd = (alpha * sin_jd) - (beta * cos_jd);

        cos_jd -= n_cos_jd;
        sin_jd -= n_sin_jd;

        K_delta += cos_jd * (dotCosCos[i][j] + dotSinSin[i][j]);
        K_delta += sin_jd * (dotSinCos[i][j] - dotCosSin[i][j]);
      }
#endif
      if (printDetails) {
        printf("TODK %d %d %.6f %.6f\n", period, offset, delta, K_delta);
      }
      K_deltas[offset] = K_delta;
    }

    // Now find the offset with the largest K_delta for this period.
    // This is the best alignment of the two videos modulo this period.
    float max_K_delta = K_deltas[0];
    float maxIndex = 0;
    for (int offset = 1; offset < period; offset++) {
      if (K_deltas[offset] > max_K_delta) {
        max_K_delta = K_deltas[offset];
        maxIndex = offset;
      }
    }
    bestOffsets[i] = maxIndex;
    valuesAtBestOffsets[i] = max_K_delta;
  }
}