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

    }

  }

}