     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) {
+__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.x == 0) {
+        // left
+        s_src[(s_row - 2) * s_width + (s_col - 2)] = safe_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)] = safe_global_mem_read(g_src, g_row - 2, g_col - 1, g_width, g_height);
+        s_src[(s_row - 2) * s_width + s_col] = safe_global_mem_read(g_src, g_row - 2, g_col, g_width, g_height);
+
+        s_src[(s_row - 1) * s_width + (s_col - 2)] = safe_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)] = safe_global_mem_read(g_src, g_row - 1, g_col - 1, g_width, g_height);
+        s_src[(s_row - 1) * s_width + s_col] = safe_global_mem_read(g_src, g_row - 1, g_col, g_width, g_height);
+
+        s_src[s_row * s_width + (s_col - 2)] = safe_global_mem_read(g_src, g_row, g_col - 2, g_width, g_height);
+        s_src[s_row * s_width + (s_col - 1)] = safe_global_mem_read(g_src, g_row, g_col - 1, g_width, g_height);
+        s_src[s_row * s_width + s_col] = safe_global_mem_read(g_src, g_row, g_col, g_width, g_height);
+
+        s_src[(s_row + 1) * s_width + (s_col - 2)] = safe_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)] = safe_global_mem_read(g_src, g_row + 1, g_col - 1, g_width, g_height);
+        s_src[(s_row + 1) * s_width + s_col] = safe_global_mem_read(g_src, g_row + 1, g_col, g_width, g_height);
+    } else if (threadIdx.x == (blockDim.x - 1)) {
+        // right
+        s_src[(s_row - 2) * s_width + s_col] = safe_global_mem_read(g_src, g_row - 2, g_col, g_width, g_height);
+        s_src[(s_row - 2) * s_width + (s_col + 1)] = safe_global_mem_read(g_src, g_row - 2, g_col + 1, g_width, g_height);
+
+        s_src[(s_row - 1) * s_width + s_col] = safe_global_mem_read(g_src, g_row - 1, g_col, g_width, g_height);
+        s_src[(s_row - 1) * s_width + (s_col + 1)] = safe_global_mem_read(g_src, g_row - 1, g_col + 1, g_width, g_height);
+
+        s_src[s_row * s_width + s_col] = safe_global_mem_read(g_src, g_row, g_col, g_width, g_height);
+        s_src[s_row * s_width + (s_col + 1)] = safe_global_mem_read(g_src, g_row, g_col + 1, g_width, g_height);
+
+        s_src[(s_row + 1) * s_width + s_col] = safe_global_mem_read(g_src, g_row + 1, g_col, g_width, g_height);
+        s_src[(s_row + 1) * s_width + (s_col + 1)] = safe_global_mem_read(g_src, g_row + 1, g_col + 1, g_width, g_height);
+    } else {
+        // center
+        s_src[(s_row - 2) * s_width + s_col] = safe_global_mem_read(g_src, g_row - 2, g_col, g_width, g_height);
+
+        s_src[(s_row - 1) * s_width + s_col] = safe_global_mem_read(g_src, g_row - 1, g_col, g_width, g_height);
+
+        s_src[s_row * s_width + s_col] = safe_global_mem_read(g_src, g_row, g_col, g_width, g_height);
+
+        s_src[(s_row + 1) * s_width + s_col] = safe_global_mem_read(g_src, g_row + 1, 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_size) {
+    int tid = blockIdx.x * blockDim.x + threadIdx.x;
+
+    while (tid < g_size) {
+        g_out[tid] = (g_in_1[tid] == g_in_2[tid]);
+        tid += (gridDim.x * blockDim.x);
+    }
+}
+
+__device__ uint8_t safe_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) {
+__device__ void safe_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

    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.
