// see https://nanonets.com/blog/optical-flow/

srecord noeq GoodThingsToTrack(BWImage image) {
  settable int windowSize = 3;
  settable double scoreFactor = 1/256.0;
  gettable BWImage outputImage;
  
  gettable double maxScore;
  
  run {
    int w = image.getWidth(), h = image.getHeight();
    int wn = windowSize*windowSize;
    int ofs = (windowSize+1)/2;
    outputImage = new BWImage(w, h);
    for (int y = 1; y < h-windowSize-1; y++) {
      for (int x = 1; x < w-windowSize-1; x++) {
        double[] ix = new[wn];
        double[] iy = new[wn];
        int iw = 0;
        for yy to windowSize:
          for xx to windowSize: {
            //double value = image.getInt(x+xx, y+yy);
            double left = image.getInt(x+xx-1, y+yy);
            double right = image.getInt(x+xx+1, y+yy);
            double above = image.getInt(x+xx, y+yy-1);
            double below = image.getInt(x+xx, y+yy+1);
            //ix[iw] = right-value;
            //iy[iw] = below-value;
            ix[iw] = right-left;
            iy[iw] = below-above;
            iw++;
          }
        
        double sideSum = 0;
        for i to wn: sideSum += ix[i]*iy[i];
        double xSum = 0;
        for i to wn: xSum += sqr(ix[i]);
        double ySum = 0;
        for i to wn: ySum += sqr(iy[i]);
        double[][] matrix = new double[][]
          { { xSum, sideSum }, { sideSum, ySum } };
        Complex[] lambdas = eigenvaluesOfSymmetricMatrix(matrix);
        if (l(lambdas) == 2) {
          double score = min(abs(lambdas[0]), abs(lambdas[1]));
          if (score > maxScore) maxScore = score;
          outputImage.setInt(x+ofs, y+ofs, iround(score*scoreFactor));
        }
      }
    }
  }
}