srecord noeq ImageSimplifier(BufferedImage inputImage) {
  replace Channel with int.
  
  // input image (mandatory)

  CompactIntegralImage mainImage;
  
  // Be verbose?
  
  bool verbose, verboseImageSize, debug;
    
  // the big internal cache

  new Map<Rect, IIntegralImage> clipCache;
  
  // channel definitions (r, g, b and grayscale are channels)
  
  static final int grayscale = 3; // channel number for grayscale
  static final int channels = 4;
  
  // user-settable maximums
  
  long maxSteps = 1000;

  // very important value. size of smallest features in relation to min(image width, image height)
  
  double featureSize = 0.1;

  long steps;
  double lowestExecutedProbability;
  
  // OUTPUT of the algorithm (found points)
  
  new ProbabilisticList<IIntegralImage> liveliestPoints;

  // INTERNAL (probabilistic scheduling, memory)
  
  ProbabilisticList<IIntegralImage> scheduler;
  Set<IIntegralImage> lookedAt;

  abstract class IIntegralImage implements main IIntegralImage {
    // width and height of image
    int w, h;
    
    public int getWidth() { ret w; }
    public int getHeight() { ret h; }

    int liveliness_cachedChannel = -1;
    double liveliness_cache;
    
    long discoveredInStep;
    
    public abstract double getIntegralValue(int x, int y, Channel channel);
    
    BufferedImage render() {
      ret imageFromFunction(w, h, (x, y) -> rgbPixel(x, y, x+1, y+1) | fullAlphaMask());
    }

    double getPixel(Rect r, int channel) {
      ret getPixel(r.x, r.y, r.x2(), r.y2(), channel);
    }

    double getPixel(int channel) { ret getPixel(0, 0, w, h, channel); }
    
    // return value ranges from 0 to 1 (usually)
    double getPixel(int x1, int y1, int x2, int y2, int channel) {
      ret doubleRatio(rectSum(x1, y1, x2, y2, channel), (x2-x1)*(y2-y1)*255.0);
    }
    
    public double rectSum(Rect r, int channel) {
      ret rectSum(r.x, r.y, r.x2(), r.y2(), channel);
    }
    
    public double rectSum(int x1, int y1, int x2, int y2, int channel) {
      double bottomRight = getIntegralValue(x2-1, y2-1, channel);
      double topRight    = getIntegralValue(x2-1, y1-1, channel);
      double bottomLeft  = getIntegralValue(x1-1, y2-1, channel);
      double topLeft     = getIntegralValue(x1-1, y1-1, channel);
      ret bottomRight-topRight-bottomLeft+topLeft;
    }

    int rgbPixel(int x1, int y1, int x2, int y2) {
      int r = iround(clampZeroToOne(getPixel(x1, y1, x2, y2, 0))*255);
      int g = iround(clampZeroToOne(getPixel(x1, y1, x2, y2, 1))*255);
      int b = iround(clampZeroToOne(getPixel(x1, y1, x2, y2, 2))*255);
      ret rgbInt(r, g, b);
    }
    
    double liveliness(int channel) {
      if (liveliness_cachedChannel != channel) {
        // optimization (but no change in semantics):
        // if (w <= 1 && h <= 1) ret 0; // liveliness of single pixel is 0
        liveliness_cache = standardDeviation(map(q -> q.getPixel(channel), quadrants()));
        liveliness_cachedChannel = channel;
      }
      ret liveliness_cache;
    }

    // no duplicates, without full image
    L<IIntegralImage> descentShapes_cleaned() {
      ret uniquify(listMinus(descentShapes(), this));
    }

    L<IIntegralImage> descentShapes() {
      ret centerPlusQuadrants();
    }
    
    L<IIntegralImage> centerPlusQuadrants() {
      int midX = w/2, midY = h/2;
      Rect r = rectAround(iround(midX), iround(midY), max(midX, 1), max(midY, 1));
      ret itemPlusList(clip(r), quadrants());
    }
    
    L<IIntegralImage> quadrants() {
      if (w <= 1 && h <= 1) null; // let's really not have quadrants of a single pixel
      int midX = w/2, midY = h/2;
      ret mapLL clip(
        rect(0, 0, max(midX, 1), max(midY, 1)),
        rect(midX, 0, w-midX, max(midY, 1)),
        rect(0, midY, max(midX, 1), h-midY),
        rect(midX, midY, w-midX, h-midY)
      );
    }

    IIntegralImage liveliestSubshape(int channel) {
      ret highestBy(q -> q.liveliness(channel), quadrants());
    }
    
    ProbabilisticList<IIntegralImage> liveliestSubshape_probabilistic(int channel) {
      ret new ProbabilisticList<IIntegralImage>(map(descentShapes(), shape ->
        withProbability(shape.liveliness(channel), shape)));
    }

    public IIntegralImage clip(Rect r) {
      Rect me = rect(0, 0, w, h);
      r = intersectRects(r, 0, 0, w, h);
      if (r.x == 0 && r.y == 0 && r.w == w && r.h == h) this;
      ret actuallyClip(r);
    }

    IIntegralImage actuallyClip(Rect r) {
      ret newClip(this, r);
    }    
    
    public IIntegralImage clip(int x1, int y1, int w, int h) { ret clip(rect(x1, y1, w, h)); }
    
    Rect positionInImage(IIntegralImage mainImage) {
      ret this == mainImage ? positionInImage() : null;
    }

    Rect positionInImage() {
      ret rect(0, 0, w, h);
    }
    
    Pt center() {
      ret main center(positionInImage());
    }
    
    double area() { ret w*h; }
    double relativeArea() { ret area()/mainImage.area(); }

    bool singlePixel() { ret w <= 1 && h <= 1; }

    toString { ret w + "*" + h; }
  }

  // virtual clip of an integral image
  class Clip extends IIntegralImage {
    IIntegralImage fullImage;
    int x1, y1;

    *(IIntegralImage *fullImage, Rect r) {
      x1 = r.x; y1 = r.y; w = r.w; h = r.h;
    }
     
    *(IIntegralImage *fullImage, int *x1, int *y1, int *w, int *h) {}
    
    public double getIntegralValue(int x, int y, int channel) {
      ret fullImage.getIntegralValue(x+x1, y+y1, channel);
    }

    // don't clip a clip - be smarter than that!
    IIntegralImage actuallyClip(Rect r) {
      ret newClip(fullImage, translateRect(r, x1, y1));
    }

    Rect positionInImage() {
      ret rect(x1, y1, w, h);
    }

    Rect positionInImage(IIntegralImage mainImage) {
      try object Rect r = super.positionInImage(mainImage);
      if (fullImage == mainImage) ret rect(x1, y1, w, h);
      null;
    }

    toString { ret positionInImage() + " in " + fullImage; }
  }

  class CompactIntegralImage extends IIntegralImage {
    short[] data;
    
    // for each 8*8 block, 4 shorts
    // We store the bits 15 to 30 so we have one overlapping
    // bit with the lower words (only way to get the addition right)
    short[] highWords;
    static final int blockShift = 3;
    static final int blockSize = 1 << 3;
    int gridW, gridH;

    *(CompactIntegralImage img) {
      w = img.w;
      h = img.h;
      data = img.data;
      highWords = img.highWords;
      gridW = img.gridW;
      gridH = img.gridH;
    }
    
    *(BufferedImage img) {
      w = img.getWidth(); h = img.getHeight();
      int safety = 64;
      if (longMul(w, h)*channels*256 > Int.MAX_VALUE-safety) fail("Image too big: " + w + "*" + h);
      gridW = iceil_divideByPowerOfTwo(blockShift, w);
      gridH = iceil_divideByPowerOfTwo(blockShift, h);
      highWords = new short[gridW*gridH*channels];
      
      int[] pixels = pixelsOfBufferedImage(img);
      data = new short[w*h*channels];
      int i = 0, j = 0;
      int[] sum = new[channels];
      int[] lastRow = new[w*channels];
      int iHighWords = 0;
      for y to h: {
        for c to channels: sum[c] = 0;
        int iLastRow = 0;
        for x to w: {
          int rgb = pixels[j++] & 0xFFFFFF;
          for c to channels: {
            int current;
            if (c == grayscale)
              current = iround((sum[0]+sum[1]+sum[2])/3);
            else {
              current = (sum[c] += rgb >> 16);
              rgb = (rgb << 8) & 0xFFFFFF;
            }
            int value = current + lastRow[iLastRow];
            short low = (short) value;
            data[i] = low;
            lastRow[iLastRow++] = value;
            i++;
            
            if (debug)
              printVars(+x, +y, +c, rgb := intToHex(rgb), +current, +value, low := ushortToInt(low), +iLastRow, +iHighWords);

            // if top-left corner of a block, set high word
            if ((x & (blockSize-1)) == 0 && (y & (blockSize-1)) == 0) {
              short high = (short) (value >> 15);
              highWords[iHighWords++] = high;
              if (debug)
                printVars(+high);
              ifdef CompactIntegralImage_safetyTests
                int blockX = x >> blockShift, blockY = y >> blockShift;
                assertEquals(iHighWords-1, (blockY*gridW+blockX)*channels+c);
                assertEquals(ushortToInt(high), getHighWord(x, y, c));
                if (blockX == 141) {
                  printVars(+x, +y, +c, rgb := intToHex(rgb), +current, +value, low := ushortToInt(low), +iLastRow, +iHighWords);
                  printVars(+blockX, +blockY, +high);
                }
              endifdef
            }              
          }
        }
      }
    }

    int getHighWord(int x, int y, int channel) {    
      int blockX = x >> blockShift, blockY = y >> blockShift;
      ret ushortToInt(highWords[(blockY*gridW+blockX)*channels+channel]);
    }
    
    public double getIntegralValue(int x, int y, Channel channel) {
      if (x < 0 || y < 0) ret 0;
      y = min(y, h-1);
      x = min(x, w-1);
      int blockX = x >> blockShift, blockY = y >> blockShift;
      int low = data[(y*w+x)*channels+channel];
      int high = ushortToInt(highWords[(blockY*gridW+blockX)*channels+channel]);
      
      int value = ((high & 1) == 0)
        //   - need unsigned expansion of low
        // bit 15 in block is 0 in top left corner
        ? (high << 15) + (low & 0xFFFF)
        
        // bit 15 in block is 1 in top left corner
        //   - need signed expansion of low
        : ((high+1) << 15) + low;
        
      if (debug)
        printVars getIntegralValue(+x, +y, +channel, +blockX, +blockY, +high, low := low + "/" + ushortToInt((short) low), +value);
        
      ret value;
    }
  }

  IIntegralImage newClip(IIntegralImage fullImage, Rect r) {
    assertSame(fullImage, mainImage);
    ret getOrCreate(clipCache, r, () -> new Clip(fullImage, r));
  }

  IIntegralImage liveliestPointIn(IIntegralImage image) {
    ret applyUntilEqual_goOneBackOnNull(c -> c.liveliestSubshape(grayscale), image);
  }

  // level++ <=> a fourth the area
  double level(IIntegralImage image) {
    ret -log(image.relativeArea(), 4);
  }
  
  // featureSize = relative to smaller image dimension
  double actualFeatureSize() {
    ret featureSize*min(mainImage.w, mainImage.h);
  }

  void clearCaches {
    clipCache.clear();
  }

  void prepareImage {
    if (mainImage == null)
      // make integral image
      mainImage = new CompactIntegralImage(inputImage);
    else
      // not sure why we are doing this one or whether we should do it
      mainImage = new CompactIntegralImage(mainImage);
    inputImage = null;
    
    if (verbose || verboseImageSize) print("Full image size: " + mainImage.w + "*" + mainImage.h);
  }
  
  run {
    prepareImage();

    time "Simplification" {
    }
  }
  
  void setInputImage aka setImage(BufferedImage image) {
    inputImage = image;
    mainImage = null;
  }
}