// Approximate pixel-by-pixel similarity of 2 images very quickly

// (probabilistically)

// Note: If my math is correct - and it always is -, floating-point

// rounding errors should not occur even for the largest images.

// Note 2: The data structures can get pretty big if you go to a

// very high detail level.

srecord noeq ImageSimilarityAtScale(

  IIntegralImage img1,

  IIntegralImage img2,

  int featureSize // smallest block size to be looked at (w/h in pixels)

) extends Probabilistic implements IF0<Double> {

  int channels;

  LookAtClip root;

  int blocksCalculated;

  bool verbose;

  bool updateInnerNodes;

  double factor;

  double rootDiffs;

  // key = area in pixels

  // Hmm. stil working this out. How to get the proper diff at

  // different feature sizes?

  //Map<Int, RatioAccumulator> diffAtScale = autoMap(() -> new RatioAccumulator);

  /*void updateDiffAtScale(int area, Double change) {

    addToDoubleValueMap(diffAtScale, area, change, area);

  }*/

  record LookAtClip(LookAtClip parent, Rect r) implements Runnable {

    L<LookAtClip> children;

    // diffs are in range [0;sqrt(channels)*pixelsCovered]

    double initialDiff, diffsFromChildren, areaCoveredByChildren;

    double latestDiff;

    double latestDiff() { ret latestDiff; }

    double calculateLatestDiff() {

      double area = area();

      double areaLeft = 1-areaCoveredByChildren/area;

      latestDiff = initialDiff*areaLeft+diffsFromChildren;

      if (verbose) printVars(+r, +latestDiff, +initialDiff, +areaLeft, +diffsFromChildren);

      ret latestDiff;

    }

    run {

      // make initial guess

      double sum = 0;

      for channel to channels: {

        double sum1 = img1.rectSum(r, channel);

        double sum2 = img2.rectSum(r, channel);

        sum += sqr(sum1-sum2);

        //if (verbose) printVars(+r, +sum1, +sum2, +sum);

      }

      setInitialDiff(sqrt(sum)*factor);

      ++blocksCalculated;

      if (r.w > r.h) {

        if (r.w > featureSize)

          splitHorizontally();

      } else {

        if (r.h > featureSize)

          splitVertically();

      }

    }

    void splitHorizontally {

      split(asList(splitRectInHorizontalHalves(r)));

    }

    void splitVertically {

      split(asList(splitRectInVerticalHalves(r)));

    }

    void split(L<Rect> parts) {

      double p = 0.99;

      children = map(parts, r -> new LookAtClip(this, r));

      scheduleAllRelative(p, children);

    }

    void setInitialDiff(double initialDiff) {    

      this.initialDiff = initialDiff;

      //updateDiffAtScale(area(), initialDiff);

      calculateLatestDiff();

      if (verbose) printVars setInitialDiff(+r, +initialDiff, +latestDiff);

      if (this == root) rootDiffs = latestDiff;

      if (parent != null) {

        parent.areaCoveredByChildren += area();

        parent.diffsFromChildren += latestDiff;

        parent.update();

      }

    }

    void update() {

      double oldValue = latestDiff;

      calculateLatestDiff();

      if (verbose) printVars update(+r, +latestDiff);

      rootDiffs += latestDiff-oldValue;

      if (updateInnerNodes && parent != null) {

        parent.diffsFromChildren += latestDiff-oldValue;

        if (verbose)

          printVars update2(+r, +oldValue, +latestDiff,

            dfc := parent.diffsFromChildren);

        parent.update();

      }

    }

    simplyCached int area() { ret rectArea(r); }

  }

  run {

    assertEquals("Sizes of compared images have to be equal", img1.getSize(), img2.getSize());

    assertEquals("Channel counts of compared images have to be equal", channels = img1.nChannels(), img2.nChannels());

    factor = 1/(255.0*sqrt(channels));

    root = new LookAtClip(null, imageRect(img1));

    root.run();

  }

  // similarity between 0 and 1

  // it's a best guess until the calculation is complete

  public Double similarity aka get() {

    ret 1-diff();

  }

  // difference (1-similarity) between 0 and 1

  // (this can be more precise in floating point then similarity)

  public double diff() {

    ret /*root.latestDiff()*/rootDiffs/root.area();

  }

  int blocksCalculated() { ret blocksCalculated; }

}