// 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));

        }

      }

    }

  }

}