The most simple Difference of Gaussian image pyramid with OpenCL

06 Aug 2023 - tsp
Last update 06 Aug 2023
Reading time 22 mins

The following blog article summarizes a quick implementation of a Difference of Gaussian (DoG) pyramid as approximation of a Laplacian of a Gaussian pyramid. It does not utilize optimizations so performance (around 10ms for a simple 6 octave pyramid) is rather low. Of course all timings taken here are dependent on the CPU and GPU that has been used. To generate the graphs and timing shown in this blog article an AMD Radeon R9 360 has been used with a Xeon E3 host CPU on FreeBSD 13.

Theory

Image Pyramids are often used during machine vision. They usually are used to create multi-scale image representations by repeated application of downsampling. To reduce resolution of a picture one usually performs subsampling after applying (Gaussian) blur. The blurring step is necessary to prevent total loss of high frequency information (for example think about an binary image where every second pixel in every second row is set to white and every other to black where you downsample by a factor of two - then you only would get a single image. When applying blur it would be gray afterwards). This still forms kind of a low pass since one looses high frequency information inside the images - and it forms the different scale spaces. Depending on application one can also device to not perform blurring but applying the mean to build the downsampled version - how one does this depends entirely on the application (calculating a mean for 9 pixels for example would only require only 8 additions and one multiplication, applying an 3x3 Gaussian blur would require 9 multiplications and 8 additions so way more operations).

Inside the scale spaces one often performs bandpass filtering or searching for maxima and minima. This is usually done by calculating the derivative of the image. This can be approximated by applying a Gaussian blur multiple times and calculating differences between the resulting images:

[ f(\sigma_1; \sigma_2) = I * \frac{1}{\sigma_1 \sqrt{2 \pi}} e^{-\frac{x^2 - y^2}{2 \sigma_1^2}} - I * \frac{1}{\sigma_2 \sqrt{2 \pi}} e^{-\frac{x^2 - y^2}{2 \sigma_2^2}} ]

Even though one performs a convolution with the image this yields

[ f(\sigma_1; \sigma_2) = I * \left( \frac{1}{\sigma_1 \sqrt{2 \pi}} e^{-\frac{x^2 - y^2}{2 \sigma_1^2}} - \frac{1}{\sigma_2 \sqrt{2 \pi}} e^{-\frac{x^2 - y^2}{2 \sigma_2^2}} \right) ]

By building those differences from can approximate the Laplacian of a Gaussian well enough - this is an approach thatā€™s often called Loweā€™s Pyramid scheme:

Since one then only wants to perform scale space maxima one can show that itā€™s sufficient to calculate $3$ difference levels inside each octave - this can be verified by checking how many stable keypoints one can find when extending the number of scale spaces which has been done extensively in literature.

Implementing image and difference of Gaussian (DoG) pyramids in OpenCL

It has been shown previously how to implement convolution more efficient on the GPU than on the CPU. The convolution kernel had been pretty simple. Now the next step will be to implement the whole subsampling and blurring for the scale spaces as well as the difference of Gaussian pyramids on the GPU.

Whatā€™s necessary to build a DoG Pyramid?

OpenCL initialization

The first step is to initialize the OpenCL framework. Since Iā€™m going to test the kernels in Python this code uses pyopencl. The very simple testbench simply fetches a list of all platforms and .their devices and then selects a predefined one (usedPlatform and usedDevice). Since this code is only thought as a testbench for the OpenCL code there is no external configuration option. A detailed description can be found in a previous blog article

The initialization creates or initializes:

import pyopencl as cl
import numpy as np
import matplotlib.pyplot as plt

from PIL import Image

from time import time

platforms = cl.get_platforms()
for idxPlatform, platform in enumerate(platforms):
    devices = platform.get_devices()
    print(f"Platform {platform.name}:")
    for idxDevice, dev in enumerate(devices):
        print(f"   {idxDevice}: {dev.name}")

usedPlatform = 0
usedDevice = 0

clPlatform = platforms[usedPlatform]
clDevice = clPlatform.get_devices()[usedDevice]
clContext = cl.Context([ clDevice ])
clQueue = cl.CommandQueue(clContext, clDevice)

Defining the grayscale kernel

The first step for most image processing algorithms is grayscaling. This of course reduces the amount of data and especially the complexity - but also reduces the dimensionality and amount of information in your data. The grayscale kernel will accept an RGBA image with 4 bytes per pixel, one byte per channel - the alpha channel will be ignored. It then grayscales with the factors specified in ITU-R BT.709 standard thatā€™s also used in HDTV to compute the luma component. It then clamps the output value to the 0 to 255 range and writes it into a simple output array. For a quick test itā€™s assumed all images are RGBA images - in a real world scenario one will have to support a variety of input formats (at least RGB, RGBA and YUV frames as well as accepting already grayscaled images).

The ITU-R BT.709 formula assigns the luma component as:

$Y = 0.2126 * R + 0.7152 * G + 0.0722 * B$

This tries to represent the sensitivity of the human eye to different colors in some way. Since all pixels are independent this is a very simple kernel:

clkernel_Grayscale_Tpl = '''//CL//

#ifndef GRAYSCALE_FACTOR_R
    #define GRAYSCALE_FACTOR_R 0.2126
#endif
#ifndef GRAYSCALE_FACTOR_G
    #define GRAYSCALE_FACTOR_G 0.7152
#endif
#ifndef GRAYSCALE_FACTOR_B
    #define GRAYSCALE_FACTOR_B 0.0722
#endif

__kernel void grayscale_RGBA2Pyrabuffer(
    __global unsigned char* lpRGBAIn,
    __global unsigned char* lpPyraBufferOut
) {
    int y = get_global_id(0);
    int x = get_global_id(1);

    float gray_value = ((float)(lpRGBAIn[(y + x * IMAGEHEIGHT) * 4 + 0])) * GRAYSCALE_FACTOR_R + ((float)(lpRGBAIn[(y + x * IMAGEHEIGHT) * 4 + 1])) * GRAYSCALE_FACTOR_G + ((float)(lpRGBAIn[(y + x * IMAGEHEIGHT) * 4 + 2])) * GRAYSCALE_FACTOR_B;

    if (gray_value > 255) {
        gray_value = 255.0;
    }
    if (gray_value < 1) {
        gray_value = 0;
    }

    lpPyraBufferOut[(y + x * 2 * IMAGEHEIGHT)] = (unsigned char)gray_value;
}
'''

#
# clkernel_Grayscale = '''
# #define IMAGEHEIGHT {0}
# #define IMAGEWIDTH {1}
# '''.format(image_input.shape[1], image_input.shape[0])  + clkernel_Grayscale_Tpl
#
# clProg_Grayscale = cl.Program(clContext, clkernel_Grayscale).build()

Defining the (special) convolution kernel

This kernel is used to perform convolution with a Gaussian in horizontal and vertical direction only and performs simple repeating of border pixels - different modes will require a different implementation, for example using preprocessor definitions. It reads from the same buffer as it writes into and only accepts a source and target index in the horizontal and vertical direction (which is used to apply the convolution only to a window inside a buffer to allow all images to be stored side by side - in the same buffer object) as well as the selector for the horizontal or vertical convolution. The simple implementation of a convolution is described in a previous blog article.

In addition a short Python function is supplied that initializes the buffer containing the Gaussian kernel constants. Those are calculated on the host instead of being calculated on the GPU on the fly.

def gaussiankernel(width, sigma):
    if width % 2 != 1:
        raise ValueError("Require odd kernel size")

    r = np.empty((width,))

    for k in range(int(width/2) + 1):
        r[int(width/2) + k] = 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(- (k*k) / (2 * sigma*sigma))
        r[int(width/2) - k] = 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(- (k*k) / (2 * sigma*sigma))

    # Apply norm ...
    r = r / np.sum(r)

    return r

#imagekernel_GaussianBlurGrayscalePyramid = gaussiankernel(5, np.sqrt(2))
imagekernel_GaussianBlurGrayscalePyramid = gaussiankernel(7, np.sqrt(2))

clkernel_GaussianBlurGrayscalePyramid = '''//CL//
__kernel void convolveGrayscalePyramid(
    __global unsigned char* lpBuffer,
    __global unsigned char* lpTempBuffer,
    __constant float* lpKernel,

    int sourceX, int sourceY,
    int destX, destY
) {
    int x = get_global_id(1);
    int y = get_global_id(0);

    float v = 0.0;

    /* Read the pixel values and apply Gaussian weights */

    for(int i = -KERNELHALFWIDTH; i <= KERNELHALFWIDTH; i=i+1) {
        #if defined(CONVOLVE_HORIZONTAL)
            int realx = x + i;
            if(realx < 0) {
                realx = 0;
            } else if(realx >= IMAGEWIDTH) {
                realx = IMAGEWIDTH-1;
            }

            float kv =  lpKernel[i + KERNELHALFWIDTH];

            realx = realx + sourceX * IMAGEWIDTH;
            int realy = y + sourceY * IMAGEHEIGHT;

            v = v + lpBuffer[realy + realx * 2 * IMAGEHEIGHT] * kv;
        #else
            int realy = y + i;
            if(realy < 0) {
                realy = 0;
            } else if(realy >= IMAGEHEIGHT) {
                realy = IMAGEHEIGHT-1;
            }

            float kv =  lpKernel[i + KERNELHALFWIDTH];

            v = v + lpTempBuffer[realy + x * IMAGEHEIGHT] * kv;
        #endif
    }

    if(v > 255) {
        v = 255;
    } else if(v < 1) {
        v = 0;
    }

    #if defined(CONVOLVE_HORIZONTAL)
        lpTempBuffer[y + (x * IMAGEHEIGHT)] = (unsigned char)v;
    #else
        lpBuffer[(y + (destY * IMAGEHEIGHT)) + ((x + (destX * IMAGEWIDTH)) * 2 * IMAGEHEIGHT)] = (unsigned char)v;
    #endif
}
'''

# clkernel = '''
# #define KERNELHALFWIDTH {0}
# #define IMAGEHEIGHT {1}
# #define IMAGEWIDTH {2}
# #define CONVOLVE_HORIZONTAL 1
# '''.format(int(len(imagekernel)/2), image_input.shape[1], image_input.shape[0])  + clkernel

# clProgH = cl.Program(clContext, clkernel).build()

# clkernel = '''
# #define KERNELHALFWIDTH {0}
# #define IMAGEHEIGHT {1}
# #define IMAGEWIDTH {2}
# '''.format(int(len(imagekernel)/2), image_input.shape[1], image_input.shape[0])  + clkernel

# clProgV = cl.Program(clContext, clkernel).build()

Definition of the subsampling kernel

To decrease the resolution of an image it will be first blurred using a Gaussian (since this will also be done for the scale spaces anyways) - then one can perform simple subsampling. The subsampling kernel reduces the resolution of an image by a factor of two in each direction. It always reads from the third image in the source buffer and writes into the first image in the target buffer. It will simply get launched for half the number of input pixels - again every pixel is independent of each other so the kernel is very simple:

clkernel_SubsampleGrayscalePyramid = '''//CL//
__kernel void subsampleGrayscalePyramid(
    __global unsigned char* lpPyramidLevelIn,
    __global unsigned char* lpPyramidLevelOut
) {
    int x = get_global_id(1);
    int y = get_global_id(0);

    // We source from the third image ...
    float v = lpPyramidLevelIn[(y*2) + IMAGEHEIGHT + ((x*2) * 2 * IMAGEHEIGHT)];

    // And write into the first in the output ... not using IMAGEHEIGHT >> 2 since requiring * 2 again
    lpPyramidLevelOut[y + (x * (IMAGEHEIGHT))] = v;
}
'''

# clkernel = '''
# #define IMAGEHEIGHT {0}
# #define IMAGEWIDTH {1}
# '''.format(image_input.shape[1], image_input.shape[0])  + clkernel

# clProgSubsample = cl.Program(clContext, clkernel).build()

Definition of the difference kernel

The last step for a Difference of Gaussian (DoG) pyramid is the difference step. It will subtract all scale images in a single octave at once and write three different images at the same time. This might not be the wisest implementation due to caching issues - itā€™s just one of the most simple ones and enough for a first try. The kernel accepts the octaves buffer as well as the single difference map output buffer. In addition the source width and height is specified together with the target offset to locate the position of the fourth image of the previous octaves difference images.

# Definitions:
#  IMAGEHEIGHT_DIFFPYRAMID
#  IMAGEWIDTH_DIFFPYRAMID

#define IMAGEHEIGHT_DIFFPYRAMID {0}
#define IMAGEWIDTH_DIFFPYRAMID {1}

clkernel_DifferencePyramidDifference = '''//CL//
__kernel void grayscalePyramidDifference(
    __global unsigned char* lpPyramidLevelIn,
    __global unsigned char* lpDifferenceMapOut,

    int targetOffset_X, int targetOffset_Y,
    int sourceWidth, int sourceHeight
) {
    int x = get_global_id(1);
    int y = get_global_id(0);

    float v1 = lpPyramidLevelIn[y + (x * sourceHeight)];
    float v2 = lpPyramidLevelIn[y + (sourceHeight / 2) + (x * sourceHeight)];
    float v3 = lpPyramidLevelIn[y + ((x + sourceWidth / 2) * sourceHeight)];
    float v4 = lpPyramidLevelIn[y + sourceHeight / 2 + ((x + sourceWidth / 2) * sourceHeight)];

    lpDifferenceMapOut[(y + targetOffset_Y) +  (x + targetOffset_X) * IMAGEHEIGHT_DIFFPYRAMID] = v3 - v1;
    lpDifferenceMapOut[(y + targetOffset_Y + (sourceHeight >> 1)) +  (x + targetOffset_X) * IMAGEHEIGHT_DIFFPYRAMID] = v2 - v3;
    lpDifferenceMapOut[(y + targetOffset_Y) +  (x + targetOffset_X + (sourceWidth >> 1)) * IMAGEHEIGHT_DIFFPYRAMID] = v4 - v2;
}
'''

Quick tests

The following is just a quick run of all kernels after each other to check they really work as expected before assembling a whole pyramid calculation function

Quick test of the grayscale kernel

# Load image

image_input = np.asarray(Image.open("imagetest01.jpg").convert('RGBA'), dtype = np.uint8)

h, w, c = image_input.shape[0], image_input.shape[1], image_input.shape[2]
print(f"Source image has dimensions {h} x {w} x {c}")

# Compile our kernel
clkernel_Grayscale = '''
#define IMAGEHEIGHT {0}
#define IMAGEWIDTH {1}
'''.format(image_input.shape[1], image_input.shape[0])  + clkernel_Grayscale_Tpl

clProg_Grayscale = cl.Program(clContext, clkernel_Grayscale).build()

# Image input buffer

clBufferIn = cl.Buffer(clContext, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf = image_input)

# Would require / 4 for the ARGB -> Grayscale and * 4 for the 4 images so the size stays the same ...
npBuffer_Scales = [
    np.empty(clBufferIn.get_info(cl.mem_info.SIZE), dtype=np.uint8),
    np.empty(int(clBufferIn.get_info(cl.mem_info.SIZE) / 4), dtype=np.uint8),
]
clBuffer_Scales = [
    cl.Buffer(clContext, cl.mem_flags.READ_WRITE | cl.mem_flags.COPY_HOST_PTR, hostbuf = npBuffer_Scales[0]),
    cl.Buffer(clContext, cl.mem_flags.READ_WRITE | cl.mem_flags.COPY_HOST_PTR, hostbuf = npBuffer_Scales[1])
]

krnlGray = clProg_Grayscale.grayscale_RGBA2Pyrabuffer

krnlGray(clQueue, (w, h,), None , clBufferIn, clBuffer_Scales[0])
cl.enqueue_copy(clQueue, npBuffer_Scales[0], clBuffer_Scales[0]) #, origin=(0,0), region=(w,h))

plt.imshow(npBuffer_Scales[0].reshape((h*2, w*2)))

# Release buffers

#for buf in clBuffer_Scales:
#    buf.release()

Simple test of the grayscale kernel

Quick test of Gaussian blur kernel

clkernel = '''
#define KERNELHALFWIDTH {0}
#define IMAGEHEIGHT {1}
#define IMAGEWIDTH {2}
#define CONVOLVE_HORIZONTAL 1
'''.format(int(len(imagekernel_GaussianBlurGrayscalePyramid)/2), image_input.shape[1], image_input.shape[0])  + clkernel_GaussianBlurGrayscalePyramid

clProgH = cl.Program(clContext, clkernel).build()

clkernel = '''
#define KERNELHALFWIDTH {0}
#define IMAGEHEIGHT {1}
#define IMAGEWIDTH {2}
'''.format(int(len(imagekernel_GaussianBlurGrayscalePyramid)/2), image_input.shape[1], image_input.shape[0])  + clkernel_GaussianBlurGrayscalePyramid

clProgV = cl.Program(clContext, clkernel).build()

krnlH = clProgH.convolveGrayscalePyramid
krnlV = clProgV.convolveGrayscalePyramid

npKernel = np.asarray(imagekernel_GaussianBlurGrayscalePyramid, dtype = np.float32)
clBufferKernel = cl.Buffer(clContext, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf = npKernel)

npTempBuffer = np.empty((h, w), dtype = np.uint8)
clTempBuffer = cl.Buffer(clContext, cl.mem_flags.READ_WRITE | cl.mem_flags.COPY_HOST_PTR, hostbuf = npTempBuffer)

krnlH(clQueue, (w, h,), None , clBuffer_Scales[0], clTempBuffer, clBufferKernel, np.int32(0), np.int32(0), np.int32(0), np.int32(1))
krnlV(clQueue, (w, h,), None , clBuffer_Scales[0], clTempBuffer, clBufferKernel, np.int32(0), np.int32(0), np.int32(0), np.int32(1))

krnlH(clQueue, (w, h,), None , clBuffer_Scales[0], clTempBuffer, clBufferKernel, np.int32(0), np.int32(1), np.int32(1), np.int32(0))
krnlV(clQueue, (w, h,), None , clBuffer_Scales[0], clTempBuffer, clBufferKernel, np.int32(0), np.int32(1), np.int32(1), np.int32(0))

krnlH(clQueue, (w, h,), None , clBuffer_Scales[0], clTempBuffer, clBufferKernel, np.int32(1), np.int32(0), np.int32(1), np.int32(1))
krnlV(clQueue, (w, h,), None , clBuffer_Scales[0], clTempBuffer, clBufferKernel, np.int32(1), np.int32(0), np.int32(1), np.int32(1))

cl.enqueue_copy(clQueue, npBuffer_Scales[0], clBuffer_Scales[0]) #, origin=(0,0), region=(w,h))

plt.imshow(npBuffer_Scales[0].reshape((h*2, w*2)))

Simple test of the Gaussian blur kernel

Quick test of subsampling

clkernel = '''
#define IMAGEHEIGHT {0}
#define IMAGEWIDTH {1}
'''.format(image_input.shape[1], image_input.shape[0])  + clkernel_SubsampleGrayscalePyramid

clProgSubsample = cl.Program(clContext, clkernel).build()

krnlSubsample = clProgSubsample.subsampleGrayscalePyramid

krnlSubsample(clQueue, (int(w / 2), int(h / 2),), None, clBuffer_Scales[0], clBuffer_Scales[1])

cl.enqueue_copy(clQueue, npBuffer_Scales[1], clBuffer_Scales[1])

plt.imshow(npBuffer_Scales[1].reshape((h, w)))

Simple test of the subsampling kernel

#clBufferKernel.release()
for buf in clBuffer_Scales:
    buf.release()

Building a whole Pyramid framework

Since the kernels seem to work - now letā€™s build the whole framework for a pyramid and do some performance tests of this very simple implementation. It will just load the input image using Pillow, convert it into a native array using numpy and transfer it onto the GPU. Note that the data transfer to and from the GPU will not be counted in the performance evaluation. This is done since itā€™s assumed that the images are already available on the GPU when executing such an task and that there will be more processing being done on the GPU before transferring the images back.

The test program works very simple:

Then comes the execution phase. For each repeated iteration - which are done to get a more clean performance measurement:

The execution phase is times multiple times to get a rough feeling about the performance of the OpenCL implementation

nScales = 6
showDebugPlots = True

# Load image

#image_input = np.asarray(Image.open("imagetest01.jpg").convert('RGBA'), dtype = np.uint8)
image_input = np.asarray(Image.open("butterfly_external_testimage.jpg").convert('RGBA'), dtype = np.uint8)
h, w, c = image_input.shape[0], image_input.shape[1], image_input.shape[2]
print(f"Source image has dimensions {h} x {w} x {c}")

# Compile our kernel
clkernel_Grayscale = '''
#define IMAGEHEIGHT {0}
#define IMAGEWIDTH {1}
'''.format(image_input.shape[1], image_input.shape[0])  + clkernel_Grayscale_Tpl

clProg_Grayscale = cl.Program(clContext, clkernel_Grayscale).build()
krnlGray = clProg_Grayscale.grayscale_RGBA2Pyrabuffer

# Image input buffer

clBufferIn = cl.Buffer(clContext, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf = image_input)

# Pyramid output buffer

npPyramidBuffer = np.empty((image_input.shape[1] * 2, image_input.shape[0] * 2), dtype = np.uint8)
clPyramidBuffer = cl.Buffer(clContext, cl.mem_flags.WRITE_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf = npPyramidBuffer)

# Difference Kernel

clkernel_Difference = '''
#define IMAGEHEIGHT_DIFFPYRAMID {0}
#define IMAGEWIDTH_DIFFPYRAMID {1}
'''.format(int(image_input.shape[1] * 2), int(image_input.shape[0] * 2)) + clkernel_DifferencePyramidDifference
clProg_Difference = cl.Program(clContext, clkernel_Difference).build()
krnlDifference = clProg_Difference.grayscalePyramidDifference

# Would require / 4 for the ARGB -> Grayscale and * 4 for the 4 images so the size stays the same ...
curSize = clBufferIn.get_info(cl.mem_info.SIZE)
shapeX = image_input.shape[1]
shapeY = image_input.shape[0]

clBuffer_Scales = []
npBuffer_Scales = []
clProgH = []
clProgV = []
clKrnlH = []
clKrnlV = []
clProgSubsample = []
clKrnlSubsample = []
npTempBuffer = []
clTempBuffer = []

for iScale in range(nScales):
    print(f"Scale {iScale}: {shapeX} x {shapeY} -> {curSize}")

    npBuffer_Scales.append(np.empty(curSize, dtype=np.uint8))
    clBuffer_Scales.append(cl.Buffer(clContext, cl.mem_flags.READ_WRITE | cl.mem_flags.COPY_HOST_PTR, hostbuf = npBuffer_Scales[-1]))

    # Blurring kernels
    clkernel = '''
        #define KERNELHALFWIDTH {0}
        #define IMAGEHEIGHT {1}
        #define IMAGEWIDTH {2}
        #define CONVOLVE_HORIZONTAL 1
    '''.format(int(len(imagekernel_GaussianBlurGrayscalePyramid)/2), shapeX, shapeY)  + clkernel_GaussianBlurGrayscalePyramid
    clProgH.append(cl.Program(clContext, clkernel).build())

    clkernel = '''
        #define KERNELHALFWIDTH {0}
        #define IMAGEHEIGHT {1}
        #define IMAGEWIDTH {2}
    '''.format(int(len(imagekernel_GaussianBlurGrayscalePyramid)/2), shapeX, shapeY)  + clkernel_GaussianBlurGrayscalePyramid
    clProgV.append(cl.Program(clContext, clkernel).build())

    clKrnlH.append(clProgH[-1].convolveGrayscalePyramid)
    clKrnlV.append(clProgV[-1].convolveGrayscalePyramid)

    # Temporary buffers
    npTempBuffer.append(np.empty((shapeY, shapeX), dtype = np.uint8))
    clTempBuffer.append(cl.Buffer(clContext, cl.mem_flags.READ_WRITE | cl.mem_flags.COPY_HOST_PTR, hostbuf = npTempBuffer[-1]))

    # Subsampling kernels

    clkernel = '''
        #define IMAGEHEIGHT {0}
        #define IMAGEWIDTH {1}
    '''.format(shapeX, shapeY)  + clkernel_SubsampleGrayscalePyramid
    clProgSubsample.append(cl.Program(clContext, clkernel).build())
    clKrnlSubsample.append(clProgSubsample[-1].subsampleGrayscalePyramid)

    if shapeX % 2 == 0:
        shapeX = shapeX + 1
    if shapeY % 2 == 0:
        shapeY = shapeY + 1
    shapeX = int(shapeX / 2)
    shapeY = int(shapeY / 2)
    # curSize = int(curSize / 4)
    curSize = (shapeX * shapeY) * 4


# Now execute the kernels on a single input image ...

durations = []
showDebugPlots = False

targetOffset_X = 0
targetOffset_Y = 0

for iIteration in range(500):
    print(f"{iIteration} ... ", end = "")
    stime = time()

    # Grayscale kernel is only required once. Also grayscaling only runs once ...
    krnlGray(clQueue, (w, h,), None , clBufferIn, clBuffer_Scales[0])
    if showDebugPlots:
        cl.enqueue_copy(clQueue, npBuffer_Scales[0], clBuffer_Scales[0]) #, origin=(0,0), region=(w,h))
        plt.imshow(npBuffer_Scales[0].reshape((h*2, w*2)))
        plt.show()

    for iScale in range(nScales):
        clKrnlH[iScale](clQueue, (w, h,), None , clBuffer_Scales[iScale], clTempBuffer[iScale], clBufferKernel, np.int32(0), np.int32(0), np.int32(0), np.int32(1))
        clKrnlV[iScale](clQueue, (w, h,), None , clBuffer_Scales[iScale], clTempBuffer[iScale], clBufferKernel, np.int32(0), np.int32(0), np.int32(0), np.int32(1))

        clKrnlH[iScale](clQueue, (w, h,), None , clBuffer_Scales[iScale], clTempBuffer[iScale], clBufferKernel, np.int32(0), np.int32(1), np.int32(1), np.int32(0))
        clKrnlV[iScale](clQueue, (w, h,), None , clBuffer_Scales[iScale], clTempBuffer[iScale], clBufferKernel, np.int32(0), np.int32(1), np.int32(1), np.int32(0))

        clKrnlH[iScale](clQueue, (w, h,), None , clBuffer_Scales[iScale], clTempBuffer[iScale], clBufferKernel, np.int32(1), np.int32(0), np.int32(1), np.int32(1))
        clKrnlV[iScale](clQueue, (w, h,), None , clBuffer_Scales[iScale], clTempBuffer[iScale], clBufferKernel, np.int32(1), np.int32(0), np.int32(1), np.int32(1))

        cl.enqueue_copy(clQueue, npBuffer_Scales[iScale], clBuffer_Scales[iScale]) #, origin=(0,0), region=(w,h))

        if showDebugPlots:
            plt.imshow(npBuffer_Scales[iScale].reshape((h*2, w*2)))
            plt.show()

        if iScale != (nScales - 1):
            clKrnlSubsample[iScale](clQueue, (int(w / 2), int(h / 2),), None, clBuffer_Scales[iScale], clBuffer_Scales[iScale + 1])
        #    if showDebugPlots:
        #        cl.enqueue_copy(clQueue, npBuffer_Scales[iScale + 1], clBuffer_Scales[iScale + 1])
        #        plt.imshow(npBuffer_Scales[iScale + 1].reshape((h, w)))

        # Perform difference calculation
        krnlDifference(clQueue, (w, h,), None, clBuffer_Scales[iScale], clPyramidBuffer, np.int32(targetOffset_X), np.int32(targetOffset_Y), np.int32(2*h), np.int32(2*w))

        targetOffset_X = targetOffset_X + h
        targetOffset_Y = targetOffset_Y + w

        w = int(w / 2)
        h = int(h / 2)

    etime = time()
    durations.append((etime - stime) * 1000)

# Fetch LAST Pyramid into host buffer ...
cl.enqueue_copy(clQueue, npPyramidBuffer, clPyramidBuffer)
fig, ax = plt.subplots(figsize = (6.4 * 4, 4.8 * 4))
npPyramidBuffer = (((npPyramidBuffer - np.min(npPyramidBuffer)) / np.max(npPyramidBuffer)) * 255.0).astype(np.uint8)
ax.imshow(npPyramidBuffer.reshape((image_input.shape[0] * 2, image_input.shape[1] * 2)))
#plt.savefig("imagepyramids2_initialpyramid.png")
plt.show()

fig, ax = plt.subplots()
ax.set_xlabel("Iteration")
ax.set_ylabel("Duration [ms]")
ax.plot(durations)
ax.grid()
ax.set_title(f"DoG pyramid, {nScales} scales (average {np.mean(durations)} ms)")
#plt.savefig("imagepyramids2_initialpyramid_timing.png")
#plt.savefig("imagepyramids2_initialpyramid_timing.svg")
plt.show()

Result of the grayscale difference of Gaussian pyramid test

Timing of the kernel runs

In the end one has to clean up allocated resources again:

clBufferKernel.release()
clPyramidBuffer.release()
clBufferIn.release()

for i in range(nScales):
    clBuffer_Scales[i].release()
    clBuffer_Scales[i] = None
    clProgH[i] = None
    clProgV[i] = None
    clKrnlH[i] = None
    clKrnlV[i] = None
    clProgSubsample[i] = None
    clKrnlSubsample[i] = None
    clTempBuffer[i].release()
    clTempBuffer[i] = None

Conclusion

This blog article has demonstrated that itā€™s pretty simple to calculate Lowes image pyramid schema on the GPU using OpenCL. Itā€™s by far not the most performant (I would argue maybe one of the worst) implementations but it can provide a base for optimizing from a working implementation by modifying one kernel after each other. This blog article shows a testbed in which one can play around with image pyramids as well as the optimization of kernels in a very simple way. All kernels have been written in a way to be as readable as possible.

References

This article is tagged: Programming, Math, Basics, Python, Data Mining, Computer Vision, Machine learning, Gaussian

Related articles

Some simple image filters by convolution in Python and OpenCL

This short article contains a summary about image filters and how one can realize them. It will contain a short reminder of how convolution works, what separable filters are and some example of image filters and their (most inefficient but direct) implementation in Python. In the end it will also show how one might naively implement convolution in OpenCL to increase performance in comparison to the naive Python implementations shown. This article provides one with a small playground to quickly try out separable image filter kernels and execute them on the GPU using pyopencl from a Python testing environment.

A quick recap on Householder reflectors

This article is a quick summary on Householder reflectors in real values vector spaces and their application to zero rows or columns in matrices.

Simple ESP8266 blink example with NONOS SDK and custom Makefile

Example on how one can use the ESP8266 NONOS SDK to build applications using a custom Makefile (i.e. not using the nested Makefiles of espressifs SDK) to build one of the most simple microcontroller projects.

Developing native Python extensions in C - capturing from Video4Linux devices

Since I was not really satisfied with various tutorials out there this is my personal summary on how to develop Python extensions in ANSI C. A simple library that captures frames from Video4Linux devices on various Unices without external dependencies will serve as an example.

The MLX90640 on the RaspberryPi with FreeBSD

This article explores the integration of the MLX90640 far infrared thermal sensor with a RaspberryPi using FreeBSD. It delves into the technicalities of interfacing the sensor, the intricacies of its data processing and compensation algorithms, and demonstrates these concepts with Python code examples.

A first quick look on different sentence embedding methods - playing with word and sentence embeddings

A quick overview on how one can generate word and sentence embeddings using BERT, mpnet and OpenAIs sentence embedding API in Python. The second part also includes a mini semantic search engine (with linear index that only serves demonstration purpose). This article does not offer and in depth explanation but only a quick overview glance on the APIs to play around.

Quick Recap: How to estimate the prediction error for least squares fits

When fitting data with least squares, it's easy to focus solely on getting the best-fit curve - but how confident can one really be in the results? This short article walks through the essential steps to quantify uncertainty: estimating residual noise, calculating parameter variances via the Jacobian and covariance matrix, and computing prediction error at each point. With practical Python code and clear math, it reminds of the tools to go beyond just fitting and toward understanding the reliability of ones model. And if the assumptions of Gaussian noise or local linearity donā€™t hold? A hands-on introduction to bootstrapping shows how to resample your data to build empirical confidence intervals - no derivatives needed. This article provides recipies to generate meaningful error bars and get insight on the reliability of ones fit results.

Manually compiling ESP8266 applications

A short description on what the ESP8266-NONOS-SDK really does (step by step) to build the firmware images that get uploaded to the microcontroller via esptool.

Also on this blog

The Enigma of Consciousness: Are Humans Uniquely Special in the Realm of Mind and Experience?

Are humans metaphysically unique, or is consciousness a matter of organization and dynamics that could, in principle, emerge in other substrates? This essay tackles the long-standing debate by weighing arguments for human exceptionalism against counterpoints from philosophy of mind, neuroscience, and AI. It explores whether qualities like subjective experience, qualia, and the hard problem demand a biological essence, or whether feedback loops, global broadcast, embodiment, intrinsic drives, and learning dynamics suffice to explain why minds feel the way they do. The conclusion points toward a continuum rather than an exclusive divide: human minds feel special, but their vividness and unity can be linked to patterns of information processing, not an untouchable spark. Machines do not yet meet these conditions fully, but no principled barrier rules it out. What makes minds special is excellence of organizationā€”loops, memory, drives, and adaptationā€”not the matter they are built from.

How to handle LoRA messages from thethingsnetwork from JavaEE with TTNClientJ2

This article is a (really) short tutorial on how one can handle messages from LoRA nodes that are routed via thethingsnetwork inside an self-contained JavaEE application.

Harnessing the Power of GPTs - or why GPTs are not better search engines

Generative Pre-trained Transformers (GPTs) like OpenAIā€™s ChatGPT have transformed the landscape of artificial intelligence by offering the ability to generate human-like text, adapt to diverse tasks, and understand context and patterns across various domains. This blog post provides a non-technical, broad overview of GPTs, emphasizing their versatility in applications such as brainstorming, complex problem-solving, data analysis, and more. Despite their impressive capabilities, GPTs are often misunderstood as simple search engines or seen as tools that merely reproduce information from their training data, leading to misconceptions about their potential and legality. The post aims to demystify GPTs by highlighting their unique strengths and differentiating them from traditional tools. By showcasing practical applications and explaining what GPTs can and cannot do, the article seeks to clarify these misunderstandings and demonstrate how GPTs can be powerful partners in fostering innovation, creativity, and effective problem-solving. This understanding will empower readers to leverage GPT technology more effectively in various aspects of their personal and professional lives.

Adding Grafana annotations using HTTP API

Sometimes it's convenient to add annotations to Grafana dashboards and panels using Grafanas flexible HTTP API instead of writing them into a database and fetching them from there using Annotation queries. Since I lacked a complete example and it took me some minutes to find out how this works this blog post summarizes the process.

Data protection policy

Dipl.-Ing. Thomas Spielauer, Wien (webcomplains389t48957@tspi.at)

This webpage is also available via TOR at http://rh6v563nt2dnxd5h2vhhqkudmyvjaevgiv77c62xflas52d5omtkxuid.onion/

Valid HTML 4.01 Strict Powered by FreeBSD IPv6 support