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