hough.py

import numpy
import pyopencl as cl
import pyopencl.array as array

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

program = cl.Program(ctx, """
#define BLOCK_LEN 256

__kernel void HoughCirclesKernel(
    __global int* A,
    __global int* imgData,
    int _width,
    int _height,
    int r
)
{
    __local int imgBuff[BLOCK_LEN];

    int localThreadIndex = get_local_id(0); //threadIdx.x
    int globalThreadIndex = get_local_id(0) + get_group_id(0) * BLOCK_LEN; //threadIdx.x + blockIdx.x * Block_Len
    int width = _width; int height = _height;
    int radius = r;

    A[globalThreadIndex] = 0;
    barrier(CLK_GLOBAL_MEM_FENCE);

    if(globalThreadIndex < width*height)
    {
        imgBuff[localThreadIndex] = imgData[globalThreadIndex];
        barrier(CLK_LOCAL_MEM_FENCE);

        if(imgBuff[localThreadIndex] > 0)
        {
            float s1, c1;
            for(int i = 0; i<180; i++)
            {
                s1 = sincos(i, &c1);
                int centerX = globalThreadIndex % width + radius * c1;
                int centerY = ((globalThreadIndex - centerX) / height) + radius * s1;

                if(centerX < width && centerY < height)
                    atomic_inc(A + centerX + centerY * width);
            }
        }
    }
    barrier(CLK_GLOBAL_MEM_FENCE);
}
""").build()

height = 40
width = 40
block_len = 256

img = numpy.zeros((height, width), numpy.int32)

import numpy
angles = numpy.linspace(0, numpy.pi*2, 200)
x, y = 10, 20
r = 6
for angle in angles:
    img[y + int(r * numpy.sin(angle)), x + int(r * numpy.cos(angle))] = 1

img_dev = array.to_device(queue, img)
a_dev = array.empty_like(img_dev)

program.HoughCirclesKernel(
    queue,
    ((height * width // block_len + 1) * block_len,), (block_len,),
    a_dev.data, img_dev.data, numpy.int32(width), numpy.int32(height), numpy.int32(r))

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig = plt.figure()
s = fig.add_subplot(1, 2, 1)
s.imshow(img, interpolation='none')

s = fig.add_subplot(1, 2, 2)
s.imshow(a_dev.get(), interpolation='none')
fig.savefig('t.png')