 }
 
 __device__ uint8_t block_and_reduce(uint8_t* s_data) {
-    int tid = threadIdx.y * blockDim.x + threadIdx.x;
-
-    for (int s = ((blockDim.x * blockDim.y) / 2); s > 0; s >>= 1) {
-        if (tid < s) {
-            s_data[tid] &= s_data[tid + s];
+    for (int s = (blockDim.x / 2); s > 0; s >>= 1) {
+        if (threadIdx.x < s) {
+            s_data[threadIdx.x] &= s_data[threadIdx.x + s];
         }
         __syncthreads();
     }

    return s_data[0];
}

__device__ uint8_t global_mem_read(uint8_t* g_data, int g_row, int g_col, int g_width, int g_height) {
    return is_outside_image(g_row, g_col, g_width, g_height) ? BINARY_WHITE : g_data[g_row * g_width + g_col];
}

__device__ void global_mem_write(uint8_t* g_data, int g_row, int g_col, int g_width, int g_height, uint8_t write_data) {
    if (!is_outside_image(g_row, g_col, g_width, g_height)) {
        g_data[g_row * g_width + g_col] = write_data;
    }
}

__device__ uint8_t is_outside_image(int g_row, int g_col, int g_width, int g_height) {
    return (g_row < 0) | (g_row > (g_height - 1)) | (g_col < 0) | (g_col > (g_width - 1));
}

__device__ void load_s_src(uint8_t* g_src, int g_row, int g_col, int g_width, int g_height, uint8_t* s_src, int s_row, int s_col, int s_width) {
    if ((threadIdx.y == 0) & (threadIdx.x == 0)) {
        // top-left corner
        s_src[(s_row - 2) * s_width + (s_col - 2)] = global_mem_read(g_src, g_row - 2, g_col - 2, g_width, g_height);
        s_src[(s_row - 2) * s_width + (s_col - 1)] = global_mem_read(g_src, g_row - 2, g_col - 1, g_width, g_height);
        s_src[(s_row - 2) * s_width + s_col] = global_mem_read(g_src, g_row - 2, g_col, g_width, g_height);

        s_src[(s_row - 1) * s_width + (s_col - 2)] = global_mem_read(g_src, g_row - 1, g_col - 2, g_width, g_height);
        s_src[(s_row - 1) * s_width + (s_col - 1)] = global_mem_read(g_src, g_row - 1, g_col - 1, g_width, g_height);
        s_src[(s_row - 1) * s_width + s_col] = global_mem_read(g_src, g_row - 1, g_col, g_width, g_height);

        s_src[s_row * s_width + (s_col - 2)] = global_mem_read(g_src, g_row, g_col - 2, g_width, g_height);
        s_src[s_row * s_width + (s_col - 1)] = global_mem_read(g_src, g_row, g_col - 1, g_width, g_height);
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
    } else if ((threadIdx.y == (blockDim.y - 1)) & (threadIdx.x == 0)) {
        // bottom-left corner
        s_src[s_row * s_width + (s_col - 2)] = global_mem_read(g_src, g_row, g_col - 2, g_width, g_height);
        s_src[s_row * s_width + (s_col - 1)] = global_mem_read(g_src, g_row, g_col - 1, g_width, g_height);
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);

        s_src[(s_row + 1) * s_width + (s_col - 2)] = global_mem_read(g_src, g_row + 1, g_col - 2, g_width, g_height);
        s_src[(s_row + 1) * s_width + (s_col - 1)] = global_mem_read(g_src, g_row + 1, g_col - 1, g_width, g_height);
        s_src[(s_row + 1) * s_width + s_col] = global_mem_read(g_src, g_row + 1, g_col, g_width, g_height);
    } else if ((threadIdx.y == (blockDim.y - 1)) & (threadIdx.x == (blockDim.x - 1))) {
        // bottom-right corner
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
        s_src[s_row * s_width + (s_col + 1)] = global_mem_read(g_src, g_row, g_col + 1, g_width, g_height);

        s_src[(s_row + 1) * s_width + s_col] = global_mem_read(g_src, g_row + 1, g_col, g_width, g_height);
        s_src[(s_row + 1) * s_width + (s_col + 1)] = global_mem_read(g_src, g_row + 1, g_col + 1, g_width, g_height);
    } else if ((threadIdx.y == 0) & (threadIdx.x == (blockDim.x - 1))) {
        // top-right corner
        s_src[(s_row - 2) * s_width + s_col] = global_mem_read(g_src, g_row - 2, g_col, g_width, g_height);
        s_src[(s_row - 2) * s_width + (s_col + 1)] = global_mem_read(g_src, g_row - 2, g_col + 1, g_width, g_height);

        s_src[(s_row - 1) * s_width + s_col] = global_mem_read(g_src, g_row - 1, g_col, g_width, g_height);
        s_src[(s_row - 1) * s_width + (s_col + 1)] = global_mem_read(g_src, g_row - 1, g_col + 1, g_width, g_height);

        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
        s_src[s_row * s_width + (s_col + 1)] = global_mem_read(g_src, g_row, g_col + 1, g_width, g_height);
    } else if (threadIdx.y == 0) {
        // PAD_TOP top rows
        s_src[(s_row - 2) * s_width + s_col] = global_mem_read(g_src, g_row - 2, g_col, g_width, g_height);
        s_src[(s_row - 1) * s_width + s_col] = global_mem_read(g_src, g_row - 1, g_col, g_width, g_height);
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
    } else if (threadIdx.x == 0) {
        // PAD_LEFT left rows
        s_src[s_row * s_width + (s_col - 2)] = global_mem_read(g_src, g_row, g_col - 2, g_width, g_height);
        s_src[s_row * s_width + (s_col - 1)] = global_mem_read(g_src, g_row, g_col - 1, g_width, g_height);
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
    } else if (threadIdx.y == (blockDim.y - 1)) {
        // PAD_BOTTOM bottom rows
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
        s_src[(s_row + 1) * s_width + s_col] = global_mem_read(g_src, g_row + 1, g_col, g_width, g_height);
    } else if (threadIdx.x == (blockDim.x - 1)) {
        // PAD_RIGHT right rows
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
        s_src[s_row * s_width + (s_col + 1)] = global_mem_read(g_src, g_row, g_col + 1, g_width, g_height);
    } else {
        // center pixels
        s_src[s_row * s_width + s_col] = global_mem_read(g_src, g_row, g_col, g_width, g_height);
    }

    __syncthreads();
}

__device__ uint8_t P2_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[(s_row - 1) * s_width + s_col];
}

__device__ uint8_t P3_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[(s_row - 1) * s_width + (s_col - 1)];
}

__device__ uint8_t P4_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[s_row * s_width + (s_col - 1)];
}

__device__ uint8_t P5_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[(s_row + 1) * s_width + (s_col - 1)];
}

__device__ uint8_t P6_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[(s_row + 1) * s_width + s_col];
}

__device__ uint8_t P7_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[(s_row + 1) * s_width + (s_col + 1)];
}

__device__ uint8_t P8_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[s_row * s_width + (s_col + 1)];
}

__device__ uint8_t P9_f(uint8_t* s_data, int s_row, int s_col, int s_width) {
    return s_data[(s_row - 1) * s_width + (s_col + 1)];
}

__global__ void pixel_equality(uint8_t* g_in_1, uint8_t* g_in_2, uint8_t* g_out, int g_width, int g_height) {
    int g_row = blockIdx.y * blockDim.y + threadIdx.y;
    int g_col = blockIdx.x * blockDim.x + threadIdx.x;

    uint8_t write_data = (global_mem_read(g_in_1, g_row, g_col, g_width, g_height) == global_mem_read(g_in_2, g_row, g_col, g_width, g_height));
    global_mem_write(g_out, g_row, g_col, g_width, g_height, write_data);
}

// Performs an image skeletonization algorithm on the input Bitmap, and stores

    // shared memory for tile
    extern __shared__ uint8_t s_src[];

    int total_size = g_width * g_height;
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    while (tid < total_size) {
        int g_row = tid / g_width;
        int g_col = tid % g_width;

        int s_row = PAD_TOP;
        int s_col = threadIdx.x + PAD_LEFT;
        int s_width = blockDim.x + PAD_LEFT + PAD_RIGHT;

        // load data into shared memory
        load_s_src(g_src, g_row, g_col, g_width, g_height, s_src, s_row, s_col, s_width);

        uint8_t NZ = black_neighbors_around(s_src, s_row, s_col, s_width);
        uint8_t TR_P1 = wb_transitions_around(s_src, s_row, s_col, s_width);
        uint8_t TR_P2 = wb_transitions_around(s_src, s_row - 1, s_col, s_width);
        uint8_t TR_P4 = wb_transitions_around(s_src, s_row, s_col - 1, s_width);
        uint8_t P2 = P2_f(s_src, s_row, s_col, s_width);
        uint8_t P4 = P4_f(s_src, s_row, s_col, s_width);
        uint8_t P6 = P6_f(s_src, s_row, s_col, s_width);
        uint8_t P8 = P8_f(s_src, s_row, s_col, s_width);

        uint8_t thinning_cond_1 = ((2 <= NZ) & (NZ <= 6));
        uint8_t thinning_cond_2 = (TR_P1 == 1);
        uint8_t thinning_cond_3 = (((P2 & P4 & P8) == 0) | (TR_P2 != 1));
        uint8_t thinning_cond_4 = (((P2 & P4 & P6) == 0) | (TR_P4 != 1));
        uint8_t thinning_cond_ok = thinning_cond_1 & thinning_cond_2 & thinning_cond_3 & thinning_cond_4;

        uint8_t write_data = BINARY_WHITE + ((1 - thinning_cond_ok) * s_src[s_row * s_width + s_col]);
        safe_global_mem_write(g_dst, g_row, g_col, g_width, g_height, write_data);

        tid += (gridDim.x * blockDim.x);
    }
}

// Computes the number of white to black transitions around a pixel.
