Re: filling area by color atack safety - worst memory size

On Fri, 3 May 2024 18:33:05 +0300
Michael S <already5chosen@yahoo.com> wrote:

On Thu, 25 Apr 2024 17:56:06 +0300
Michael S <already5chosen@yahoo.com> wrote:
 
 
A solution (sort of) is in line with the famous quite of David Wheeler
- to turn todo lists from bit maps into arrays of
abscesses-or-ordinates of contact points.
 
The cost is a memory footprint - 4x bigger than the previous version,
32 times bigger than above-mentioned "packed" variant of the previous
version. But in BigO sense it's the same.
 
In my tests it reduced the worst case time from O(max(A,B)**3) to
O(A*B*log(max(A,B)). Which is non-ideal, but probably acceptable,
because the bad cases should be very rare in practice.
 
The real trouble is different - I don't know if my "worst case" is
really the worst.
 
The code below is for presentation of algorithm in both clear and
compact manner, with emphasis on symmetry between x and y directions.
It is not optimal in any sense and can be made no-trivially faster
both by algorithm enhancements an by specialization of critical loops.
 
>

Following code improves on ideas from the previous post.
Unlike the previous one, it is purely iterative, with no recursion.
The algorithm is simpler and access storage in more compact manner, i.e.
all accessed memory area starts from beginning and grows according to
need. Previous attempt did not have this property.
It's still longer and less simple than I would like.

// try+split algorithm with flat storage
// - horizontal intervals
// - two stacks: main stack for intervals,
//   auxiliary stack of areas of interest (AoI)
// - both stacks implemented as arrays
#include <stddef.h>
#include <string.h>
#include <stdlib.h>
#include <stdint.h>
#define NDEBUG
#include <assert.h>
#include <stdio.h>

typedef unsigned char Color;

typedef struct {
  size_t n_intervals;
  size_t n_splits;
} stack_sizes_t;

static
stack_sizes_t floodfill4_calc_stack_size(int width, int height)
{
  stack_sizes_t sz = { .n_intervals = 1, .n_splits = 1 };
  for (;;) {
    ptrdiff_t len;
    if (width > height) { // split vertically
      len = height;
      width = (width + 1)/2;
    } else { // split horizontally
      len = width;
      height = (height + 1)/2;
    }
    if (len <= 1)
      break;
    sz.n_intervals += len*2 + 4;
    sz.n_splits += 1;
  }
  return sz;
}

int floodfill4(
  Color* image,
  int width, int height,
  int x,     int y,
  Color old_color, Color new_color)
{
  if (width <= 0 || height <= 0)
    return 0;

  if (x < 0 || x >= width || y < 0 || y >= height)
    return 0;

  size_t w = width;
  Color* row = &image[w*y];
  if (row[x] != old_color)
    return 0;

  enum coordinate_axes {
    x_i = 0, y_i, // index of pos[] MS bit of index of limits[][]
  };
  #define X2Y(axis) ((axis) ^ 1)
  enum beg_or_end {
    beg_i = 0, end_i // LS bit of index of limits[],
                     // I use 0 and 1 more commonly
  };
  enum limits_idx { // index of limits[]
    x0_i = x_i*2+beg_i,
    x1_i = x_i*2+end_i,
    y0_i = y_i*2+beg_i,
    y1_i = y_i*2+end_i,
  };

  typedef struct {
    int x0, x1, y;
    int from; // 0 => from y-1, 1 => from y+1
  } interval_t;
  typedef struct {
    interval_t* parent_todo;
    int         saved_limit_val;
    uint8_t     saved_limit_idx; // axis*2+beg_or_end
    int         frame_capacity_deficit;
  } parent_info_t;

  stack_sizes_t stacks_len = floodfill4_calc_stack_size(width, height);
  const size_t parent_info_sz = stacks_len.n_splits *
  sizeof(parent_info_t); const size_t todo_sz = stacks_len.n_intervals
  * sizeof(interval_t); void* stacks = malloc(parent_info_sz + todo_sz);
  if (!stacks)
    return -1;

  parent_info_t* parents_stack = stacks;
  parent_info_t* parents_stack_end =
  &parents_stack[stacks_len.n_splits]; interval_t* todo_stack =
  (interval_t*)parents_stack_end; interval_t* todo = todo_stack;
  #ifndef NDEBUG
  interval_t* todo_stack_end = &todo[stacks_len.n_intervals];
  #endif

  int limits[2*2] = { 0, width-1, 0, height-1}; // {x0, x1, y0, y1};

  // recolor initial horizontal interval
  row[x] = new_color;
  // look backward
  int x00;
  for (x00 = x-1; x00 >= 0 && row[x00]==old_color; --x00)
    row[x00] = new_color;
  x00 += 1;
  // look forward
  int x01;
  for (x01 = x+1; x01 < width && row[x01]==old_color; ++x01)
    row[x01] = new_color;
  x01 -= 1;
  // push neighbors of initial interval on todo stack
  for (enum beg_or_end from = beg_i; from <= end_i; ++from) {
    unsigned next_y = y+1-from*2;
    if (next_y < (unsigned)height) {
      todo->x0 = x00;
      todo->x1 = x01;
      todo->y  = next_y;
      todo->from = from;
      ++todo;
    }
  }

  parent_info_t* parent_aoi = parents_stack;
  interval_t* parent_todo = todo_stack;
  ptrdiff_t frame_capacity = width < height ? width : height;
  for (;;) {
    while (todo != parent_todo) {
      assert(todo_stack_end != todo);
      assert(parent_todo >= todo_stack && parent_todo <=
  todo_stack_end); // Get interval from top of todo stack
      --todo; // pop interval from todo stack
      int xBeg = todo->x0;
      int xEnd = todo->x1;
      int y = todo->y;
      int from = todo->from;
      // check range
      if ((unsigned)(y-limits[y0_i]) >
  (unsigned)(limits[y1_i]-limits[y0_i]) || xEnd < limits[x0_i] || xBeg
  > limits[x1_i]) { // Whole interval belongs to parent
        // Bring value from the bottom of todo stack to the top
        // freeing stack slot for parent stack
        assert(todo_stack_end != todo);
        *todo = *parent_todo; ++todo;
        // Store interval on top of parent stack
        parent_todo->x0   = xBeg;
        parent_todo->x1   = xEnd;
        parent_todo->y    = y;
        parent_todo->from = from;
        ++parent_todo;
        continue;
      }
      // At least a part of the interval is in current rectangle
      if (xBeg < limits[x0_i]) {
        // left part of interval belongs to parent
        // Store left part of interval on todo stack
        // for later demotion to parent's stack
        assert(todo_stack_end != todo);
        todo->x0 = xBeg;
        todo->x1 = limits[x0_i]-1;
        todo->y = y;
        todo->from = from;
        ++todo;
        xBeg = limits[x0_i]; // adjust xBeg
      }
      if (xEnd > limits[x1_i] ) {
        // right part of interval belongs to parent
        // Store right part of interval on todo stack
        // for later demotion to parent's stack
        assert(todo_stack_end != todo);
        todo->x0 = limits[x1_i]+1;
        todo->x1 = xEnd;
        todo->y = y;
        todo->from = from;
        ++todo;
        xEnd = limits[x1_i]; // adjust xEnd
      }
      // remaining part of interval is within limits

      // look for target points
      Color* row = &image[y*w];
      int x = xBeg;
      do {
        if (row[x] == old_color) { // target found
          if (todo-parent_todo > frame_capacity) {
            // can't complete floodfill of current rectangle
            // due to space constraints.
            // Split
            const int dLim[] = {
              limits[x1_i]-limits[x0_i],
              limits[y1_i]-limits[y0_i]};
            const enum coordinate_axes axis =
              dLim[x_i] > dLim[y_i] ?
                x_i :  // split vertically
                y_i ;  // split horizontally
            // select half
            const int hpos0 = (limits[axis*2+0] +
                          limits[axis*2+1])/2; // lower split point
            const int pos[2] = { [x_i] = x, [y_i] = y};
            enum beg_or_end src_i = pos[axis] > hpos0;

            // preserve state of current rectangle on parents stack
            assert(parent_aoi != parents_stack_end);
            parent_aoi->parent_todo = parent_todo;
            enum beg_or_end save_i = 1 - src_i;
            enum limits_idx saved_limit_idx = axis*2+save_i;
            parent_aoi->saved_limit_idx = saved_limit_idx;
            parent_aoi->saved_limit_val = limits[saved_limit_idx];
            parent_aoi->frame_capacity_deficit =
                                todo-parent_todo - frame_capacity;
            ++parent_aoi;

            // switch processing to selected half of rectangle
            frame_capacity = dLim[X2Y(axis)] + 1;
            limits[saved_limit_idx] = hpos0+src_i;
            parent_todo = todo;
            // push interval on fresh todo stack
            assert(todo_stack_end != todo);
            todo->x0 = x;
            todo->x1 = xEnd;
            todo->y  = y;
            todo->from = from;
            ++todo;
            break;
          }

          row[x] = new_color;
          // look forward
          int x1;
          for (x1 = x+1; x1 < width && row[x1]==old_color; ++x1)
            row[x1] = new_color;
          x1 -= 1;

          int x0 = x;
          if (x == xBeg) {
            // look backward
            for (x0 = x-1; x0 >= 0 && row[x0]==old_color; --x0)
              row[x0] = new_color;
            x0 += 1;

            // bounce
            if (x0 <= xBeg-2) {
              assert(todo_stack_end != todo);
              todo->x0 = x0;
              todo->x1 = xBeg-2;
              todo->y  = y+from*2-1;
              todo->from = 1-from;
              ++todo;
            }
          }
          // bounce
          if (x1 >= xEnd+2) {
            assert(todo_stack_end != todo);
            todo->x1 = x1;
            todo->x0 = xEnd+2;
            todo->y  = y+from*2-1;
            todo->from = 1-from;
            ++todo;
          }
          unsigned next_y = y+1-from*2;
          if (next_y < (unsigned)height) {
            // continuation
            #if 1
            // The following if is not necessary for correction
            // It is here to speed up few test cases
            if (y != limits[y0_i+1-from] &&
                x0 >= limits[x0_i]       &&
                x1 <= limits[x1_i]       &&
                x1+2 > xEnd              &&
                todo-parent_todo <= frame_capacity)
            { // Bypass stack
              // Advance vertically in the same direction
              xBeg = x = x0;
              xEnd = x1;
              y = next_y;
              row = &image[y*w];
              continue;
            }
            #endif
            // put new interval on current stack
            assert(todo_stack_end != todo);
            todo->x0 = x0;
            todo->x1 = x1;
            todo->y = next_y;
            todo->from = from;
            ++todo;
          }
          x = x1+1;
        }
        ++x;
      } while (x <= xEnd);
    }

    if (parent_aoi == parents_stack)
      break; // top AOI finished

    // back to parent rectangle
    --parent_aoi;
    parent_todo = parent_aoi->parent_todo;
    limits[parent_aoi->saved_limit_idx]= parent_aoi->saved_limit_val;
    frame_capacity = todo - parent_todo
                     - parent_aoi->frame_capacity_deficit;
  }
  assert((void*)todo == (void*)parents_stack_end);

  free(stacks);
  return 1;
}