 #define PAD_BOTTOM 1
 #define PAD_RIGHT 1
 
-void and_reduction(uint8_t* g_data, int g_width, int g_height, dim3 grid_dim, dim3 block_dim) {
-    int shared_mem_size = block_dim.x * block_dim.y * sizeof(uint8_t);
+void and_reduction(uint8_t* g_data, int g_size, dim3 grid_dim, dim3 block_dim) {
+    int shared_mem_size = block_dim.x * sizeof(uint8_t);
 
     // iterative reductions of g_data
     do {
-        and_reduction<<<grid_dim, block_dim, shared_mem_size>>>(g_data, g_width, g_height);
+        and_reduction<<<grid_dim, block_dim, shared_mem_size>>>(g_data, g_size);
         gpuErrchk(cudaPeekAtLastError());
         gpuErrchk(cudaDeviceSynchronize());
 
-        g_width = grid_dim.x;
-        g_height = grid_dim.y;
-        grid_dim.x = ceil(grid_dim.x / ((double) block_dim.x));
-        grid_dim.y = ceil(grid_dim.y / ((double) block_dim.y));
-    } while ((g_width * g_height) != 1);
-}
-
-// Adapted for 2D arrays from Nvidia cuda SDK samples
-__global__ void and_reduction(uint8_t* g_data, int g_width, int g_height) {
+        g_size = ceil(g_size / ((double) block_dim.x));
+        grid_dim.x = (g_size <= block_dim.x) ? 1 : grid_dim.x;
+    } while (g_size != 1);
+}
+
+__global__ void and_reduction(uint8_t* g_data, int g_size) {
     // shared memory for tile
     extern __shared__ uint8_t s_data[];
 
-    int g_row = blockIdx.y * blockDim.y + threadIdx.y;
-    int g_col = blockIdx.x * blockDim.x + threadIdx.x;
-
-    int tid = threadIdx.y * blockDim.x + threadIdx.x;
-
-    // Load equality values into shared memory tile. We use 1 as the default
-    // value, as it is an AND reduction
-    s_data[tid] = is_outside_image(g_row, g_col, g_width, g_height) ? 1 : global_mem_read(g_data, g_row, g_col, g_width, g_height);
-    __syncthreads();
-
-    // do reduction in shared memory
-    uint8_t write_data = block_and_reduce(s_data);
-
-    // write result for this block to global memory
-    if (tid == 0) {
-        global_mem_write(g_data, blockIdx.y, blockIdx.x, gridDim.x, gridDim.y, write_data);
+    int blockReductionIndex = blockIdx.x;
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+
+    int num_blocks_for_reduction = ceil(g_size / ((double) blockDim.x));
+    int num_iterations_for_reduction = ceil(num_blocks_for_reduction / ((double) gridDim.x));
+
+    for (int iteration = 0; iteration < num_iterations_for_reduction; iteration++) {
+        // Load equality values into shared memory tile. We use 1 as the default
+        // value, as it is an AND reduction
+        s_data[threadIdx.x] = (i < g_size) ? g_data[i] : 1;
+        __syncthreads();
+
+        // do reduction in shared memory
+        block_and_reduce(s_data);
+
+        // write result for this block to global memory
+        if (threadIdx.x == 0) {
+            g_data[blockReductionIndex] = s_data[0];
+        }
+
+        blockReductionIndex += gridDim.x;
+        i += (gridDim.x * blockDim.x);
     }
 }
 

    uint8_t are_identical_bitmaps = 0;
    int iterations = 0;
    do {
        int skeletonize_pass_shared_mem_size = (block_dim.x + PAD_LEFT + PAD_RIGHT) * (block_dim.y + PAD_TOP + PAD_BOTTOM) * sizeof(uint8_t);
        skeletonize_pass<<<grid_dim, block_dim, skeletonize_pass_shared_mem_size>>>(g_src_data, g_dst_data, (*src_bitmap)->width, (*src_bitmap)->height);
        gpuErrchk(cudaPeekAtLastError());
        gpuErrchk(cudaDeviceSynchronize());

        pixel_equality<<<grid_dim, block_dim>>>(g_src_data, g_dst_data, g_equ_data, (*src_bitmap)->width, (*src_bitmap)->height);
        gpuErrchk(cudaPeekAtLastError());
        gpuErrchk(cudaDeviceSynchronize());

        and_reduction(g_equ_data, (*src_bitmap)->width, (*src_bitmap)->height, grid_dim, block_dim);

        // bring reduced bitmap equality information back from device
        gpuErrchk(cudaMemcpy(&are_identical_bitmaps, g_equ_data, 1 * sizeof(uint8_t), cudaMemcpyDeviceToHost));

    // 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.
