After trying my luck at implementing it in Fortran 90, I turned to C. I'm not sure whether it's my queue implementation or the actual flood fill algorithm but it doesn't work as expected and tends to fill the whole array, irrespective of any border I'm putting in it's way. My goal is to create a Python extension.
Below is the code, which is based on scikit-image's implementation of this algorithm.
The initial intent was to implement the queue as a dynamic array and borrow the code from rosettacode but I couldn't make it work either so I've tried to make it more like scikit-image's QueueWithHistory implementation, although with minor adjustments.
queue.h
#ifndef QUEUE_H
#define QUEUE_H
typedef size_t DATA; /* type of data to store in queue */
typedef struct {
DATA *buf;
size_t head;
size_t tail;
size_t alloc;
} queue_t, *queue;
extern queue init_queue();
extern void free_queue(queue q);
extern void push_queue(queue q, DATA *n);
extern unsigned char pop_queue(queue q, DATA *n);
#endif /* QUEUE_H */
queue.c
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "queue.h"
#define QUEUE_INIT_SIZE 64
#define QUEUE_THRESHOLD 512
static void _handle_memory_allocation_failure(const char *message, queue q) {
fprintf(stderr, "Error: Failed to (re)allocate memory for %s\n", message);
free_queue(q); // Ensure proper cleanup
exit(EXIT_FAILURE);
}
queue init_queue() {
queue q = malloc(sizeof(queue_t));
if (!q) {
_handle_memory_allocation_failure("queue structure", NULL);
}
q->buf = malloc(sizeof(DATA) * (q->alloc = QUEUE_INIT_SIZE));
if (!q->buf) {
_handle_memory_allocation_failure("queue buffer", q);
}
q->head = q->tail = 0;
return q;
}
void free_queue(queue q) {
if (!q) return;
free(q->buf);
free(q);
}
static void _grow_queue(queue q) {
q->alloc *= 2;
DATA *new_buf = realloc(q->buf, sizeof(DATA) * q->alloc);
if (!new_buf) {
_handle_memory_allocation_failure("growing queue buffer", q);
}
q->buf = new_buf;
}
void push_queue(queue q, DATA *n) {
if (q->alloc <= q->head) {
_grow_queue(q);
}
q->buf[q->head++] = *n;
}
unsigned char pop_queue(queue q, DATA *n) {
if (q->tail < q->head) {
*n = q->buf[q->tail++];
return 1;
}
return 0;
}
floodfill.h
#ifndef FLOODFILL_H
#define FLOODFILL_H
// Function declarations
int flood_fill(int* image, int* offsets, size_t n_offsets,
size_t start_index, int new_value, size_t image_size);
#endif /* FLOODFILL_H */
floodfill.c
#include <stdlib.h>
#include <stdio.h>
#include "floodfill.h"
#include "queue.h"
// C implementation of floodfill algorithm
int flood_fill(int* image, int* offsets, size_t n_offsets,
size_t start_index, int new_value, size_t image_size) {
size_t current_index, neighbor_index;
// Get the seed value (image value at start location)
int seed_value = image[start_index];
// Short circuit to avoid infinite loop
if (seed_value == new_value) {
return 0;
}
image[start_index] = new_value;
// Initialize the queue and push start position
queue q = init_queue();
push_queue(q, &start_index);
// Break loop if all queued positions were evaluated
while (pop_queue(q, ¤t_index)) {
// Look at all neighbors
for (int i = 0; i < n_offsets; ++i) {
neighbor_index = current_index + offsets[i];
if (neighbor_index < image_size) {
if (image[neighbor_index] == seed_value) {
// Process neighbor and push to queue
image[neighbor_index] = new_value;
push_queue(q, &neighbor_index);
}
}
}
}
// Ensure memory is released
free_queue(q);
return 0;
}
The above C code is wrapped in a Python module, below is the main function that calls the C code. It's called flood_fill in the module API.
#define RAVEL_INDEX(tuple, h) ( \
(int)PyLong_AsLong(PyTuple_GetItem(tuple, 0)) + \
h * (int)PyLong_AsLong(PyTuple_GetItem(tuple, 1)) \
)
#define N_OFFSETS 8
static PyObject *algorithms_flood_fill(PyObject *self, PyObject *args) {
PyObject *image_obj;
PyObject *start_obj;
int seed_value;
if (!PyArg_ParseTuple(args, "O!O!i", &PyArray_Type, &image_obj,
&PyTuple_Type, &start_obj,
&seed_value)) {
PyErr_SetString(PyExc_TypeError, "Error parsing input");
return NULL;
}
PyArrayObject *image_array;
image_array = (PyArrayObject*)PyArray_FROM_OTF(image_obj,
NPY_INT,
NPY_ARRAY_IN_FARRAY);
if (image_array == NULL) {
// Py_XDECREF(image_array);
PyErr_SetString(PyExc_TypeError, "Couldn't parse the input array");
return NULL;
}
int h = (int)PyArray_DIM(image_array, 0);
int w = (int)PyArray_DIM(image_array, 1);
intptr_t start_index = (intptr_t)RAVEL_INDEX(start_obj, h);
// it should be safe to use PyArray_DATA as we ensured aligned contiguous Fortran array
int *image = (int *)PyArray_DATA(image_array);
size_t image_size = h * w;
int offsets[N_OFFSETS] = {-(h+1), -h, -(h-1), -1, +1, h-1, h, h+1};
flood_fill(image, offsets, N_OFFSETS, start_index, seed_value, image_size);
Py_DECREF(image_array);
Py_RETURN_NONE;
}
If I give it a 9x9 array of integers initialized to 0 with a row of 1 in the middle (row index 4), int offsets[] = {-10, -9, -8, -1, 1, 8, 9, 10}, size_t n_offsets = 8 size_t start_index = 0, int new_value = 3, size_t image_size = 81, all the 0 in the array get replaced by 3, which isn't the intended behavior.
import numpy as np
from mymodule import flood_fill
a = np.zeros((9, 9), order='F', dtype=np.int32)
a[4, :] = 1
flood_fill(a, (0, 0), 3)
# output :
# array([[3, 3, 3, 3, 3, 3, 3, 3, 3],
# [3, 3, 3, 3, 3, 3, 3, 3, 3],
# [3, 3, 3, 3, 3, 3, 3, 3, 3],
# [3, 3, 3, 3, 3, 3, 3, 3, 3],
# [1, 1, 1, 1, 1, 1, 1, 1, 1],
# [3, 3, 3, 3, 3, 3, 3, 3, 3],
# [3, 3, 3, 3, 3, 3, 3, 3, 3],
# [3, 3, 3, 3, 3, 3, 3, 3, 3],
# [3, 3, 3, 3, 3, 3, 3, 3, 3]])
Is there something wrong with either the queue implementation or the flood fill algorithm logic itself or maybe the way it is called by the Python wrapper that could cause this problem ?
When working with an image stored as a 1D array, you need more complex logic to check whether the neighbors are valid.
Try something like this:
// C implementation of floodfill algorithm
//offset indexes are [upleft,up,upright,left,right,downleft,down,downright]
const int UP_LEFT=0,UP=1,UP_RIGHT=2,LEFT=3,RIGHT=4,DOWN_LEFT=5,DOWN=6,DOWN_RIGHT=7;
int valid_neighbour(int current_index, int dim, int offset_index){
// first column
if(current_index%dim==0
&& (offset_index==UP_LEFT
||offset_index==LEFT
||offset_index==DOWN_LEFT)){
return 0;
}
//last column
if(current_index%dim==dim-1
&& (offset_index==UP_RIGHT
||offset_index==RIGHT
||offset_index==DOWN_RIGHT)){
return 0;
}
//first row
if(current_index/dim==0
&& (offset_index==UP_LEFT
||offset_index==UP
||offset_index==UP_RIGHT)){
return 0;
}
//last row
if(current_index/dim==dim-1
&& (offset_index==DOWN_LEFT
||offset_index==DOWN
||offset_index==DOWN_RIGHT)){
return 0;
}
return 1;
}
int flood_fill(int* image, int* offsets, int n_offsets,
int start_index, int new_value, int image_size,int image_dim) {
int current_index, neighbor_index;
// Get the seed value (image value at start location)
int seed_value = image[start_index];
// Short circuit to avoid infinite loop
if (seed_value == new_value) {
return 0;
}
//image[start_index] = new_value;
// Initialize the queue and push start position
queue q = init_queue();
push_queue(q, &start_index);
// Break loop if all queued positions were evaluated
while (pop_queue(q, ¤t_index)) {
//Don't process already processed pixels
if(image[current_index]==new_value){
continue;
}
image[current_index]=new_value;
printf("curr index: %d\n",current_index);
// Look at all neighbors
for (int i = 0; i < n_offsets; ++i) {
neighbor_index = current_index + offsets[i];
if (neighbor_index < image_size && neighbor_index>=0) {
if (image[neighbor_index] == seed_value && valid_neighbour(current_index,image_dim,i)) {
// Process neighbor and push to queue
//image[neighbor_index] = new_value;
push_queue(q, &neighbor_index);
}
}
}
}
// Ensure memory is released
free_queue(q);
return 0;
}
int main()
{
int image[] = {
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0
};
int offsets[] = {-10, -9, -8, -1, 1, 8, 9, 10};
int n_offsets = 8;
int start_index = 0;
int seed_value = 3;
int image_size = 81;
int image_dim=9;
flood_fill(image, offsets, n_offsets, start_index, seed_value, image_size,image_dim);
for(int i=0;i<9;i++){
for(int j=0;j<9;j++){
printf("%d ",image[i*9+j]);
}
printf("\n");
}
}
NOTE: I used int instead of size_t, change it later.