srecord noeq FastRegions_BWImage(BWImage image) is Runnable, IImageRegions<BWImage> {
  int w, h, runner;
  gettable int size;
  
  new IntBuffer stack; // locations as y*w+x
  int[] regionMatrix; // for each pixel: region index (starting at 1)
  new IntBuffer regionPixels; // collect all pixels for regions
  
  // initialize these to use them
  
  new IntBuffer regionFirstPixel; // for each region: index of first pixel found in regionPixels
  new IntBuffer regionSize; // for each region: number of pixels
  new IntBuffer regionBounds; // for each region: bounds (x1, y1, x2, y2)
  
  // index = dual log of region size, value = region index
  new L<IntBuffer> regionsBySize;
  
  int regionCounter;
  bool verbose;
  
  double regionStep = .1; // for rendering in regionsImage
  
  int x(int pos) { ret pos % w; }
  int y(int pos) { ret pos / w; }
  int pos(int x, int y) { ret y*w+x; }
  Pt pt(int pos) { ret Pt(x(pos), y(pos)); }
  bool validPos(int x, int y) { ret x >= 0 && y >= 0 && x < w && y < h; }
  
  int getColor(int pos) { ret image.getInt(x(pos), y(pos)); }
  
  *(BufferedImage img) { this(toBWImage(img)); }

  run {
    w = image.getWidth(); h = image.getHeight();
    size = w*h;
    regionMatrix = new int[size];
    
    // 0 entries are unused
    regionFirstPixel?.add(0); 
    regionSize?.add(0);
    
    regionPixels?.setSize(size);
    
    while (runner < size) {
      if (regionMatrix[runner] == 0) {
        // make a new region, get color
        int region = ++regionCounter;
        regionFirstPixel?.add(regionPixels != null ? l(regionPixels) : runner);
        stack.add(runner);
        int color = getColor(runner);
        int rsize = 0, x1 = w, y1 = h, x2 = 0, y2 = 0;
        ifdef FastRegions_debug
          printVars FastRegions(+region, +runner);
        endifdef

        // flood-fill region
        while (nempty(stack)) {
          int pos = stack.popLast();
          if (regionMatrix[pos] != 0) continue; // touching myself (or someone else)
          if (getColor(pos) != color) continue; // wrong color
        
          // new pixel found, mark as ours
          regionMatrix[pos] = region;
          ++rsize;
          regionPixels?.add(pos);
          int x = x(pos), y = y(pos);
          
          ifdef FastRegions_debug
            printVars FastRegions(+x, +y);
          endifdef
          
          if (x < x1) x1 = x;
          if (x > x2) x2 = x;
          if (y < y1) y1 = y;
          if (y > y2) y2 = y;

          // explore neighborhood
          if (x > 0)   stack.add(pos-1);
          if (x < w-1) stack.add(pos+1);
          if (y > 0)   stack.add(pos-w);
          if (y < h-1) stack.add(pos+w);
        }
        
        regionSize?.add(rsize);
        regionBounds?.addAll(x1, y1, x2+1, y2+1);
        if (regionsBySize != null) {
          int iBucket = dualLog(rsize);
          var buffer = listGetOrCreate(regionsBySize, iBucket, -> new IntBuffer);
          buffer.add(region);
        }
      }
      
      ++runner;
    }
  }
  
  IBWImage regionsImage() {
    ret iBWImageFromFunction((x, y) -> {
      var region = regionMatrix[pos(x, y)];
      ret ((region-1)*regionStep) % (1.0+regionStep-0.0001);
    }, w, h);
  }
  
  int regionCount aka nRegions() { ret regionCounter; }
  
  abstract class RegionIterator {
    int pos;
    
    abstract bool next();
    
    int pos() { ret pos; }
    int x() { ret BWImage_FastRegions.this.x(pos); }
    int y() { ret BWImage_FastRegions.this.y(pos); }
  }
  
  // returns points in no particular order
  class FloodRegionIterator > RegionIterator {
    int region;
    new IntBuffer stack; // locations as y*w+x
    BitSet seen = new(size);
    
    *(int *region) {
      int pos = regionFirstPixel.get(region);
      printVars(+region, +pos);
      seen.set(pos);
      stack.add(pos);
    }
    
    // flood-fill region
    bool next() {
      if (empty(stack)) false;
      
      pos = stack.popLast();

      // explore neighborhood
      int x = x(), y = y();
      if (x > 0)   tryPosition(pos-1);
      if (x < w-1) tryPosition(pos+1);
      if (y > 0)   tryPosition(pos-w);
      if (y < h-1) tryPosition(pos+w);

      true;
    }
    
    private void tryPosition(int p) {
      if (!seen.get(p) && regionMatrix[p] == region) {
        seen.set(p);
        stack.add(p);
      }
    }
  }
  
  class CachedRegionIterator > RegionIterator {
    int i, to;

    *(int region) {
      i = regionFirstPixel.get(region);
      to = region+1 < l(regionFirstPixel) ? regionFirstPixel.get(region+1) : l(regionPixels);
      
      ifdef FastRegions_debug
        printVars CachedRegionIterator(+region, +i, +to);
      endifdef
    }
    
    bool next() {
      if (i >= to) false;
      pos = regionPixels.get(i++);
      true;
    }
  }
  
  int regionSize(int iRegion) { ret regionSize.get(iRegion); }
  
  Pt samplePixel aka firstPixel(int iRegion) {
    ret pt(firstPixelPos(iRegion));
  }
  
  int firstPixelPos(int iRegion) {
    int i = regionFirstPixel.get(iRegion);
    ret regionPixels != null ? regionPixels.get(i) : i;
  }
  
  bool inRegion(int iRegion, int x, int y) {
    ret validPos(x, y) && regionMatrix[pos(x, y)] == iRegion;
  }
  Rect regionBounds(int iRegion) { ret rectFromPoints(
    regionBounds.get((iRegion-1)*4),
    regionBounds.get((iRegion-1)*4+1),
    regionBounds.get((iRegion-1)*4+2),
    regionBounds.get((iRegion-1)*4+3),
  ); }
  
  int regionAt(Pt p) { ret regionAt(p.x, p.y); }
  
  int regionAt(int x, int y) {
    ret !validPos(x, y) ? 0 : regionMatrix[pos(x, y)];
  }
  
  int regionAt(int pos) {
    ret regionMatrix[pos];
  }
  
  RegionIterator regionIterator(int iRegion) {
    ret regionPixels != null
      ? new CachedRegionIterator(iRegion)
      : new FloodRegionIterator(iRegion);
  }
  
  L<Pt> regionPixels(int iRegion) {
    var it = regionIterator(iRegion);
    new L<Pt> l;
    while (it.next())
      l.add(Pt(it.x(), it.y()));
    ret l;
  }
  
  // select extra features before regions are made
  // (not necessary anymore)
  
  void collectFirstPixels { /*regionFirstPixel = new IntBuffer;*/ }
  void collectBounds { /*regionBounds = new IntBuffer;*/ }
  
  BitMatrix regionBitMatrix(int iRegion) {
    ret new AbstractBitMatrix(w, h) {
      public Bool get(int x, int y) {
        ret inRegion(iRegion, x, y);
      }
    };
  }
  
  void markRegionInPixelArray(int[] pixels, int iRegion, int rgba) {
    if (iRegion <= 0) ret;
    for i over pixels:
      if (regionAt(i) == iRegion)
        pixels[i] = rgba;
  }
  
  L<Int> regionIndices() { ret virtualCountList(1, nRegions()+1); }
  
  ItIt<Int> regionsRoughlyByDecreasingSize() {
    ret nestedIterator(countIterator_inclusive_backwards(regionsBySize.size()-1, 0),
      iBucket -> iterator(regionsBySize.get(iBucket)));
  }
  
  IImageRegion<BWImage> getRegion aka get(int iRegion) {
    ret new ImageRegion(iRegion);
  }
  
  public L<IImageRegion<BWImage>> regions() {
    ret listFromFunction(i -> getRegion(i+1), nRegions());
  }
  
  record ImageRegion(int iRegion) is IImageRegion<BWImage> {
    public BWImage image() { ret image; }
    public O creator() { ret FastRegions_BWImage.this; }
    public int indexInCreator() { ret iRegion; }
  
    public Rect bounds() { ret regionBounds(iRegion); }
    
    public int numberOfPixels() { ret regionSize(iRegion); }
    
    public Pt firstPixel() { ret pt(firstPixelPos()); }
    public int firstPixelPos() { ret FastRegions_BWImage.this.firstPixelPos(iRegion); }
    
    public Iterator<Pt> pixelIterator() {
      var it = regionIterator(iRegion);
      ret iteratorFromFunction(-> it.next() ? pt(it.pos()) : null);
    }
    
    public bool contains(int x, int y) { ret inRegion(iRegion, x, y); }
    
    public RGB color() {
      ret rgbFromGrayscale(brightness());
    }
    
    public int brightness() {
      ret getColor(firstPixelPos());
    }
    
    toString {
      ret renderRecordVars("Region", +brightness(), +color(), pixels := numberOfPixels(), +bounds());
    }
  }
}

Began life as a copy of #1033761