Image Processing

import pyopencl as cl
import numpy as np
import cv2
import time

# Helper function to load OpenCL kernel from file
def load_kernel(filename):
    with open(filename, 'r') as f:
        return f.read()

# Initialize OpenCL context and queue
def init_opencl():
    platform = cl.get_platforms()[0]
    device = platform.get_devices()[0]
    context = cl.Context([device])
    queue = cl.CommandQueue(context)
    return context, queue

# Grayscale conversion kernel
grayscale_kernel = """
__kernel void grayscale(__global const uchar *input, __global uchar *output, const int width, const int height) {
    int x = get_global_id(0);
    int y = get_global_id(1);
    int idx = y * width + x;
    uchar r = input[3 * idx];
    uchar g = input[3 * idx + 1];
    uchar b = input[3 * idx + 2];
    output[idx] = (uchar)(0.299f * r + 0.587f * g + 0.114f * b);
}
"""

# Box blur kernel
box_blur_kernel = """
__kernel void box_blur(__global const uchar *input, __global uchar *output, const int width, const int height) {
    int x = get_global_id(0);
    int y = get_global_id(1);
    int idx = y * width + x;
    float sum = 0.0f;
    int count = 0;
    for (int i = -1; i <= 1; i++) {
        for (int j = -1; j <= 1; j++) {
            int xi = clamp(x + i, 0, width - 1);
            int yj = clamp(y + j, 0, height - 1);
            int kidx = yj * width + xi;
            sum += input[kidx];
            count++;
        }
    }
    output[idx] = (uchar)(sum / count);
}
"""

# Sobel edge detection kernel
sobel_kernel = """
__kernel void sobel(__global const uchar *input, __global uchar *output, const int width, const int height) {
    int x = get_global_id(0);
    int y = get_global_id(1);
    int idx = y * width + x;
    float Gx = 0.0f;
    float Gy = 0.0f;
    int sobel_x[3][3] = {
        {-1, 0, 1},
        {-2, 0, 2},
        {-1, 0, 1}
    };
    int sobel_y[3][3] = {
        {-1, -2, -1},
        {0, 0, 0},
        {1, 2, 1}
    };
    for (int i = -1; i <= 1; i++) {
        for (int j = -1; j <= 1; j++) {
            int xi = clamp(x + i, 0, width - 1);
            int yj = clamp(y + j, 0, height - 1);
            int kidx = yj * width + xi;
            Gx += input[kidx] * sobel_x[i + 1][j + 1];
            Gy += input[kidx] * sobel_y[i + 1][j + 1];
        }
    }
    output[idx] = (uchar)clamp(sqrt(Gx * Gx + Gy * Gy), 0.0f, 255.0f);
}
"""

# Function to run OpenCL kernel
def run_kernel(context, queue, kernel_code, input_image, output_shape):
    program = cl.Program(context, kernel_code).build()
    mf = cl.mem_flags
    input_buf = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=input_image)
    output_buf = cl.Buffer(context, mf.WRITE_ONLY, output_shape.nbytes)
    kernel = program.grayscale if 'grayscale' in kernel_code else program.box_blur if 'box_blur' in kernel_code else program.sobel
    kernel.set_args(input_buf, output_buf, np.int32(input_image.shape[1]), np.int32(input_image.shape[0]))
    cl.enqueue_nd_range_kernel(queue, kernel, input_image.shape[:2], None)
    output_image = np.empty_like(output_shape)
    cl.enqueue_copy(queue, output_image, output_buf)
    return output_image

# Create synthetic images
def create_gradient_image(width, height):
    """Create a grayscale gradient image."""
    gradient = np.linspace(0, 255, width, dtype=np.uint8)
    gradient_image = np.tile(gradient, (height, 1))
    return gradient_image

def create_checkerboard_image(width, height, num_squares):
    """Create a checkerboard pattern image."""
    checkerboard = np.zeros((height, width), dtype=np.uint8)
    square_size = width // num_squares
    for y in range(0, height, square_size * 2):
        for x in range(0, width, square_size * 2):
            checkerboard[y:y+square_size, x:x+square_size] = 255
            checkerboard[y+square_size:y+square_size*2, x+square_size:x+square_size*2] = 255
    return checkerboard

def create_noise_image(width, height):
    """Create a random noise image."""
    noise_image = np.random.randint(0, 256, (height, width), dtype=np.uint8)
    return noise_image

def create_shapes_image(width, height):
    """Create an image with basic geometric shapes."""
    shapes_image = np.zeros((height, width), dtype=np.uint8)
    cv2.circle(shapes_image, (width // 4, height // 4), 50, 255, -1)
    cv2.rectangle(shapes_image, (width // 2, height // 2), (width // 2 + 100, height // 2 + 100), 255, -1)
    cv2.line(shapes_image, (0, height), (width, 0), 255, 5)
    return shapes_image

# Initialize OpenCL
context, queue = init_opencl()

# Create and save synthetic images
gradient_image = create_gradient_image(512, 512)
cv2.imwrite('gradient_image.jpg', gradient_image)

checkerboard_image = create_checkerboard_image(512, 512, 8)
cv2.imwrite('checkerboard_image.jpg', checkerboard_image)

noise_image = create_noise_image(512, 512)
cv2.imwrite('noise_image.jpg', noise_image)

shapes_image = create_shapes_image(512, 512)
cv2.imwrite('shapes_image.jpg', shapes_image)

# Grayscale conversion (using gradient image as an example)
gray_image = run_kernel(context, queue, grayscale_kernel, cv2.cvtColor(gradient_image, cv2.COLOR_GRAY2BGR), np.empty((512, 512), dtype=np.uint8))
cv2.imwrite('gray_image.jpg', gray_image)

# Box blur (using gradient image as an example)
blur_image = run_kernel(context, queue, box_blur_kernel, gray_image, np.empty((512, 512), dtype=np.uint8))
cv2.imwrite('blur_image.jpg', blur_image)

# Sobel edge detection (using gradient image as an example)
edge_image = run_kernel(context, queue, sobel_kernel, gray_image, np.empty((512, 512), dtype=np.uint8))
cv2.imwrite('edge_image.jpg', edge_image)

# Display results
cv2.imshow('Original', gradient_image)
cv2.imshow('Grayscale', gray_image)
cv2.imshow('Box Blur', blur_image)
cv2.imshow('Sobel Edge Detection', edge_image)
cv2.waitKey(0)
cv2.destroyAllWindows()

# Performance testing (using gradient image as an example)
start_time = time.time()
for _ in range(100):
    blur_image = run_kernel(context, queue, box_blur_kernel, gray_image, np.empty((512, 512), dtype=np.uint8))
end_time = time.time()
print(f"Box blur average execution time: {(end_time - start_time) / 100:.4f} seconds")

start_time = time.time()
for _ in range(100):
    edge_image = run_kernel(context, queue, sobel_kernel, gray_image, np.empty((512, 512), dtype=np.uint8))
end_time = time.time()
print(f"Sobel edge detection average execution time: {(end_time - start_time) / 100:.4f} seconds")