Problems in implementing adaptive thresholding using CUDA

If you run the code below, the applied image will look like this: What do I need to fix to make it work correctly?

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace cv;
using namespace std;

// CUDA kernel function
__global__ void adaptiveThresholdCUDA(const uchar* image, uchar* binary, const int* iimage, int nl, int nc, int halfSize, int blockSize, int threshold) {
    // Calculate global indices
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    // Check boundary conditions
    if (j >= nl - halfSize - 1 || i >= nc - halfSize - 1)
        return;

    // Get row addresses
    const uchar* data = image + j * nc;
    uchar* binaryRow = binary + j * nc;
    const int* idata1 = iimage + (j - halfSize) * nc;
    const int* idata2 = iimage + (j + halfSize + 1) * nc;

    // Calculate sum
    int sum = (idata2[i + halfSize + 1] - idata2[i - halfSize] - idata1[i + halfSize + 1] + idata1[i - halfSize]) / (blockSize * blockSize);

    // Apply adaptive threshold
    if (data[i] < (sum - threshold))
        binaryRow[i] = 0;
    else
        binaryRow[i] = 255;
}

int main()
{
    Mat image = imread("image/test.jpg", 0);
    if (!image.data) return 0;

    resize(image, image, Size(), 1.0, 1.0);
    namedWindow("Original Image");
    imshow("Original Image", image);

    /*          Function for Adaptive Thresholding          */
    int blockSize = 16; // Neighborhood size
    int threshold = 10; // Pixel comparison threshold

    Mat binary = image.clone();

    int nl = binary.rows; // Number of lines
    int nc = binary.cols; // Total number of elements per line

    Mat iimage;
> integral(image, iimage, CV_32S);

    // Allocate memory on device (GPU)
    uchar* imageDevice;
    uchar* binaryDevice;
    int* iimageDevice;
    cudaMalloc((void**)&imageDevice, sizeof(uchar) * nl * nc);
    cudaMalloc((void**)&binaryDevice, sizeof(uchar) * nl * nc);
    cudaMalloc((void**)&iimageDevice, sizeof(int) * nl * nc);

    // Copy input data from host to device
    cudaMemcpy(imageDevice, image.data, sizeof(uchar) * nl * nc, cudaMemcpyHostToDevice);
    cudaMemcpy(binaryDevice, binary.data, sizeof(uchar) * nl * nc, cudaMemcpyHostToDevice);
    cudaMemcpy(iimageDevice, iimage.data, sizeof(int) * nl * nc, cudaMemcpyHostToDevice);

    // Define grid and block dimensions
    dim3 blockDim(16, 16);  // You may need to adjust the block dimensions
    dim3 gridDim((nc + blockDim.x - 1) / blockDim.x, (nl + blockDim.y - 1) / blockDim.y);

    // Launch the CUDA kernel
    adaptiveThresholdCUDA << <gridDim, blockDim >> > (imageDevice, binaryDevice, iimageDevice, nl, nc, blockSize / 2, blockSize, threshold);

}

When changing the size of blockDim and applying it, if it is set too large or too small, it may not be applied.

The issue is that iimage has one extra row and one extra column compared to image.

According to the documentation of cv::integral:

>src - input image as W×H, 8-bit or floating-point (32f or 64f).
sum - integral image as (W+1)×(H+1) , 32-bit integer or floating-point (32f or 64f).

The image size of the posted sample is 332 columns by 302 rows.
The image size iimage is 333 columns by 303 rows.

We have to copy iimage to device with one more row and column:

int iimage_nl = iimage.rows; // = nl + 1
int iimage_nc = iimage.cols; // = nc + 1
cudaMalloc((void**)&iimageDevice, sizeof(int) * iimage_nl * iimage_nc);
cudaMemcpy(iimageDevice, iimage.data, sizeof(int) * iimage_nl * iimage_nc, cudaMemcpyHostToDevice);

We have to fix the kernel according to the actual size of iimage:

__global__ void adaptiveThresholdCUDA(...) {
...
int iimage_nc = nc + 1; //Number of columns in iimage equals (nc+1)
const int* idata1 = iimage + (j - halfSize) * iimage_nc;
const int* idata2 = iimage + (j + halfSize + 1) * iimage_nc;
...

Updated code sample:

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <iostream>
//#include <opencv2/core/core.hpp>
//#include <opencv2/highgui/highgui.hpp>
//#include <opencv2/imgproc/imgproc.hpp>
#include "opencv2/opencv.hpp"
using namespace cv;
using namespace std;
// CUDA kernel function
__global__ void adaptiveThresholdCUDA(const uchar* image, uchar* binary, const int* iimage, int nl, int nc, int halfSize, int blockSize, int threshold) {
// Calculate global indices
int j = blockIdx.y * blockDim.y + threadIdx.y;
int i = blockIdx.x * blockDim.x + threadIdx.x;
// Check boundary conditions
if (j >= nl - halfSize - 1 || i >= nc - halfSize - 1)
return;
// Get row addresses
int iimage_nc = nc + 1; //Number of columns in iimage equals (nc+1)
const uchar* data = image + j * nc;
uchar* binaryRow = binary + j * nc;
const int* idata1 = iimage + (j - halfSize) * iimage_nc;
const int* idata2 = iimage + (j + halfSize + 1) * iimage_nc;
// Calculate sum
int sum = (idata2[i + halfSize + 1] - idata2[i - halfSize] - idata1[i + halfSize + 1] + idata1[i - halfSize]) / (blockSize * blockSize);
// Apply adaptive threshold
if (data[i] < (sum - threshold))
binaryRow[i] = 0;
else
binaryRow[i] = 255;
}
int main()
{
//Mat image = imread("image/test.jpg", 0);
Mat image = imread("test.png", 0);
if (!image.data) return 0;
resize(image, image, Size(), 1.0, 1.0);
namedWindow("Original Image");
imshow("Original Image", image);
/*          Function for Adaptive Thresholding          */
int blockSize = 16; // Neighborhood size
int threshold = 10; // Pixel comparison threshold
Mat binary = image.clone();
//nl = 332, nc = 302
int nl = binary.rows; // Number of lines
int nc = binary.cols; // Total number of elements per line
Mat iimage;
integral(image, iimage, CV_32S);
//https://docs.opencv.org/3.4/d7/d1b/group__imgproc__misc.html#gadeaf38d7701d7ad371278d663c50c77d
//src	input image as W×H, 8-bit or floating-point (32f or 64f).
//sum	integral image as (W+1)×(H+1) , 32-bit integer or floating-point (32f or 64f).
//iimage_nl = 333, iimage_nc = 303
int iimage_nl = iimage.rows; // = nl + 1
int iimage_nc = iimage.cols; // = nc + 1
// Allocate memory on device (GPU)
uchar* imageDevice;
uchar* binaryDevice;
int* iimageDevice;
cudaMalloc((void**)&imageDevice, sizeof(uchar) * nl * nc);
cudaMalloc((void**)&binaryDevice, sizeof(uchar) * nl * nc);
cudaMalloc((void**)&iimageDevice, sizeof(int) * iimage_nl * iimage_nc);
// Copy input data from host to device
cudaMemcpy(imageDevice, image.data, sizeof(uchar) * nl * nc, cudaMemcpyHostToDevice);
cudaMemcpy(binaryDevice, binary.data, sizeof(uchar) * nl * nc, cudaMemcpyHostToDevice);
cudaMemcpy(iimageDevice, iimage.data, sizeof(int) * iimage_nl * iimage_nc, cudaMemcpyHostToDevice);
// Define grid and block dimensions
dim3 blockDim(16, 16);  // You may need to adjust the block dimensions
dim3 gridDim((nc + blockDim.x - 1) / blockDim.x, (nl + blockDim.y - 1) / blockDim.y);
// Launch the CUDA kernel
adaptiveThresholdCUDA<<<gridDim, blockDim>>>(imageDevice, binaryDevice, iimageDevice, nl, nc, blockSize / 2, blockSize, threshold);
Mat output = image.clone();
cudaMemcpy(output.data, binaryDevice, sizeof(uchar) * nl * nc, cudaMemcpyDeviceToHost); //Copy from device to host
imshow("Output Image", output);
waitKey();
destroyAllWindows();
imwrite("output.png", output);  //Save the output
}

test.png (input image):

output.png (output image):

Note that the above solution just demonstrates a possible solution (the intention is to show why the output image is "slanted").