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