How to draw 2D Gaussian blob on an OpenCV image?

Question

There are various available examples with a formula for a 2D Gaussian Blob and drawing it via Pyplot, for example:
How to generate 2D gaussian with Python?
and
How to plot a 2d gaussian with different sigma?

I'm attempting to change this over to OpenCV (in Python).

Some requirements are:

-ability to specify different height and width for the blob, i.e. ability to make the blob an ellipse (not always a circle)
-ability to specify the center point of the blob in the original image
-the value at the exact center of the blob should be 255, and the values should work their way down to 0 towards the edge of the blob
-rotation is not necessary

The final image (depending on settings of course) should look something like this:

In the context of CenterNet (which is my use case for this) the result (image with a Gaussian blob on it) is called a "Heatmap" so that's the term I'm going to use in code for the image.

Here is what I have so far:

import numpy as np
import cv2

def main():
    # suppress numpy printing in scientific notation
    np.set_printoptions(suppress=True)

    hm_width = 1600
    hm_height = 1000

    # create blank heatmap (OpenCV image)
    heatmap = np.zeros((hm_height, hm_width), dtype=np.uint8)

    blob_height = 100
    blob_width = 300
    blob_center_x = 1000
    blob_center_y = 400

    # Create a 2D Gaussian blob
    x, y = np.meshgrid(np.linspace(0, 1, blob_width), np.linspace(0, 1, blob_height))

    print('\n' + 'x: ')
    print(x.dtype)
    print(x.shape)
    print('min = ' + str(np.min(x)) + ' (s/b 0.0)')
    print('max = ' + str(np.max(x)) + ' (s/b 1.0)')
    print(x)
    print('\n' + 'y: ')
    print(y.dtype)
    print(y.shape)
    print('min = ' + str(np.min(y)) + ' (s/b 0.0)')
    print('max = ' + str(np.max(y)) + ' (s/b 1.0)')
    print(y)

    # gaussian_blob = 1.0 / (2.0 * np.pi * blob_width * blob_height) * np.exp(-((x - blob_center_x)**2.0 / (2. * blob_width**2.0) + (y - blob_center_y)**2.0 / (2. * blob_height**2.0)))

    gaussian_x_term = np.power(x - blob_center_x, 2.0) / np.power(blob_width, 2.0)
    gaussian_y_term = np.power(y - blob_center_y, 2.0) / np.power(blob_height, 2.0)
    gaussian_blob = np.exp(-1.0 * (gaussian_x_term + gaussian_y_term))

    print('\n' + 'gaussian_blob before: ')
    print(gaussian_blob.dtype)
    print(gaussian_blob.shape)
    print('min = ' + str(np.min(gaussian_blob)) + ' (s/b 0.0)')
    print('max = ' + str(np.max(gaussian_blob)) + ' (s/b 1.0)')
    print(gaussian_blob)

    # scale up the gaussian blob from the 0.0 to 1.0 range to the 0 to 255 range
    gaussian_blob = gaussian_blob * 255.0
    gaussian_blob = np.clip(gaussian_blob, a_min=0.0, a_max=255.0)
    gaussian_blob = np.rint(gaussian_blob)
    gaussian_blob = np.clip(gaussian_blob, a_min=0, a_max=255)
    gaussian_blob = gaussian_blob.astype(np.uint8)

    print('\n' + 'gaussian_blob after: ')
    print(gaussian_blob.dtype)
    print(gaussian_blob.shape)
    print('min = ' + str(np.min(gaussian_blob)) + ' (s/b 0)')
    print('max = ' + str(np.max(gaussian_blob)) + ' (s/b 255)')
    print(gaussian_blob)

    # show the blob via OpenCV
    cv2.imshow('gaussian blob', gaussian_blob)
    
    # add the gaussian blob image to the heatmap
    blob_left_edge_loc = round(blob_center_x - (0.5 * blob_width))
    blob_right_edge_loc = round(blob_center_x + (0.5 * blob_width))

    blob_top_edge_loc = round(blob_center_y - (0.5 * blob_height))
    blob_bottom_edge_loc = round(blob_center_y + (0.5 * blob_height))

    heatmap[blob_top_edge_loc:blob_bottom_edge_loc, blob_left_edge_loc:blob_right_edge_loc] = gaussian_blob

    # show the heatmap
    cv2.imshow('heatmap', heatmap)

    cv2.waitKey()
# end function

if __name__ == '__main__':
    main()

Currently both images come out almost blank, and based on the output:

x: 
float64
(100, 300)
min = 0.0 (s/b 0.0)
max = 1.0 (s/b 1.0)
[[0.         0.00334448 0.00668896 ... 0.99331104 0.99665552 1.        ]
 [0.         0.00334448 0.00668896 ... 0.99331104 0.99665552 1.        ]
 [0.         0.00334448 0.00668896 ... 0.99331104 0.99665552 1.        ]
 ...
 [0.         0.00334448 0.00668896 ... 0.99331104 0.99665552 1.        ]
 [0.         0.00334448 0.00668896 ... 0.99331104 0.99665552 1.        ]
 [0.         0.00334448 0.00668896 ... 0.99331104 0.99665552 1.        ]]

y: 
float64
(100, 300)
min = 0.0 (s/b 0.0)
max = 1.0 (s/b 1.0)
[[0.         0.         0.         ... 0.         0.         0.        ]
 [0.01010101 0.01010101 0.01010101 ... 0.01010101 0.01010101 0.01010101]
 [0.02020202 0.02020202 0.02020202 ... 0.02020202 0.02020202 0.02020202]
 ...
 [0.97979798 0.97979798 0.97979798 ... 0.97979798 0.97979798 0.97979798]
 [0.98989899 0.98989899 0.98989899 ... 0.98989899 0.98989899 0.98989899]
 [1.         1.         1.         ... 1.         1.         1.        ]]

gaussian_blob before: 
float64
(100, 300)
min = 6.880118208869318e-12 (s/b 0.0)
max = 7.240508138966562e-12 (s/b 1.0)
[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]]

gaussian_blob after: 
uint8
(100, 300)
min = 0 (s/b 0)
max = 0 (s/b 255)
[[0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 ...
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]]

it seems I'm not calculating the Gaussian blob quite right, but I'm not sure how to resolve this. Suggestions?

gaussian blob on _black_? or on some non-trivial background? — Christoph Rackwitz, Jul 19 '23 at 17:41
This is where it goes wrong: `x, y = np.meshgrid(np.linspace(0, 1, blob_width), np.linspace(0, 1, blob_height))`. The blob width or height has nothing to do with this, and the coordinate space you create goes from 0 to 1, your centroid is way far outside this space. — Cris Luengo, Jul 19 '23 at 17:51
You say "I'm attempting to change this over to OpenCV", but you're using only NumPy functions. What does OpenCV have to do with it? — Cris Luengo, Jul 19 '23 at 17:59

Marat · Accepted Answer · 2023-07-19T18:38:42.957

The first issue is that the meshgrid generated incorrectly - it uses 0..1 range while the rest of the code assumes it holds coordinates from the big picture (0..hm_width/height range). Even if we operate in normalized space, 0..1 would only capture a small portion of the curve.

Solution: I'd keep the normalized space, as it will simplify the following code a bit. I'd extend the range to -3..3 standard deviations, though; feel free to adjust.

x, y = np.meshgrid(
    np.linspace(-3, 3, blob_width), 
    np.linspace(-3, 3, blob_height), 
)

The second issue: This is a 2D distribution, you cannot compute the result on x and y separately (in a general case, as pointed out by Chris - here, you can). Normally, the formula would be

sigma = ...  # covariance matrix
pos = np.array([x - blob_center_x, y - blob_center_y])
# I don't want to figure out proper axes here, but I hope you get the idea
gaussian_blob = const * np.exp(-0.5 * np.tensordot(np.tensordot(pos, np.inv(sigma), axes=...), pos, axes=...).sum(axis=...))

But since we're operating in the normalized space now, the mean is zero, and sigma is [[1, 0], [0, 1]], that would reduce to just:

# we dont care about const as it will be normalized later
gaussian_blob = np.exp(-0.5 * (x**2 + y**2))

Finally, taking care of the ignored constant at scaling:

gaussian_blob =  (255.0 * (gaussian_blob - gaussian_blob.min()) / gaussian_blob.max()).astype(np.uint8)

UPD

As pointed out by Chris, you can save a milisecond or so by doing it this way:

X = np.linspace(-3, 3, blob_width)[None, :]
Y = np.linspace(-3, 3, blob_height)[:, None]
gaussian_blob = np.exp(-0.5*X**2) * np.exp(-0.5*Y**2)
gaussian_blob =  255.0 * (gaussian_blob - gaussian_blob.min()) / gaussian_blob.max()

The Gaussian is separable. You actually *can* compute two 1D Gaussians and multiply them together. — Cris Luengo, Jul 19 '23 at 17:53
@CrisLuengo only when your covariance matrix is that trivial, and... what's the advantage? UPD: updated the answer — Marat, Jul 19 '23 at 17:55
OP doesn't define a covariance matrix, only two sigmas and no angle. So the Gaussian is axis-aligned. The advantage is a gigantic reduction in computational cost. — Cris Luengo, Jul 19 '23 at 17:57

cdahms · Answer 2 · 2023-07-20T03:05:31.243

Here is a copy/paste runnable exmaple, including code to make the Gaussian(s) (courtesy of @Marat) and code to add Gaussian(s) to an OpenCV image, including cases where the Gaussians partially overlap. I'm not going to make this as the solution so @Marat gets credit.

There is a slight imperfection in the code to add a Gaussian to the heatmap, in that with close overlapping Gaussians there is a slight seam in the resulting image, see screenshots below. I will try to resolve this at a later time.

# heatmap_drawing_test.py

import numpy as np
import cv2
from termcolor import colored

# ref:
# https://stackoverflow.com/questions/76723027/how-to-draw-2d-gaussian-blob-on-an-opencv-image/76724003#76724003

HEATMAP_WIDTH = 1600
HEATMAP_HEIGHT = 1000

def main():
    # suppress numpy printing in scientific notation
    np.set_printoptions(suppress=True)

    # create blank heatmap (OpenCV image)
    heatmap = np.zeros((HEATMAP_HEIGHT, HEATMAP_WIDTH), dtype=np.uint8)

    blob_1_center_x = 1000
    blob_1_center_y = 400
    blob_1_width = 451
    blob_1_height = 201

    gaussian_blob_1 = make_gaussian_blob(blob_1_width, blob_1_height)

    blob_2_center_x = 1100
    blob_2_center_y = 500
    blob_2_width = 451
    blob_2_height = 201

    gaussian_blob_2 = make_gaussian_blob(blob_2_width, blob_2_height)

    # show the blob via OpenCV
    cv2.imshow('gaussian blob 1', gaussian_blob_1)
    cv2.imshow('gaussian blob 2', gaussian_blob_2)

    heatmap = add_gaussian_blob_to_heatmap(gaussian_blob_1, blob_1_center_x, blob_1_center_y, heatmap)
    heatmap = add_gaussian_blob_to_heatmap(gaussian_blob_2, blob_2_center_x, blob_2_center_y, heatmap)

    print('\n' + 'final heatmap: ')
    print(heatmap.dtype)
    print(heatmap.shape)
    print('min = ' + str(np.min(heatmap)) + ' (s/b 0)')
    print('max = ' + str(np.max(heatmap)) + ' (s/b 255)')
    print(heatmap)

    # show the heatmap
    cv2.imshow('heatmap', heatmap)

    cv2.waitKey()
# end function

def make_gaussian_blob(blob_width, blob_height):
    assert blob_height % 2 == 1 and blob_width % 2 == 1, \
        colored('\n\n' + 'in make_gaussian_blob, blob_height and blob_width must be odd numbers !!' + '\n', color='red', attrs=['bold'])

    # Create a 2D Gaussian blob
    # +-2.5 was derived from experimentation
    x, y = np.meshgrid(np.linspace(-2.5, 2.5, blob_width), np.linspace(-2.5, 2.5, blob_height))

    print('\n' + 'x: ')
    print(x.dtype)
    print(x.shape)
    print('min = ' + str(np.min(x)) + ' (s/b 0.0)')
    print('max = ' + str(np.max(x)) + ' (s/b 1.0)')
    print(x)

    print('\n' + 'y: ')
    print(y.dtype)
    print(y.shape)
    print('min = ' + str(np.min(y)) + ' (s/b 0.0)')
    print('max = ' + str(np.max(y)) + ' (s/b 1.0)')
    print(y)

    gaussian_blob = np.exp(-0.5 * (x**2 + y**2))

    print('\n' + 'gaussian_blob before: ')
    print(gaussian_blob.dtype)
    print(gaussian_blob.shape)
    print('min = ' + str(np.min(gaussian_blob)) + ' (s/b 0.0)')
    print('max = ' + str(np.max(gaussian_blob)) + ' (s/b 1.0)')
    print(gaussian_blob)

    # scale up the gaussian blob from the 0.0 to 1.0 range to the 0 to 255 range
    gaussian_blob = gaussian_blob * 255.0
    gaussian_blob = np.clip(gaussian_blob, a_min=0.0, a_max=255.0)
    gaussian_blob = np.rint(gaussian_blob).astype(np.uint8)

    print('\n' + 'gaussian_blob after: ')
    print(gaussian_blob.dtype)
    print(gaussian_blob.shape)
    print('min = ' + str(np.min(gaussian_blob)) + ' (s/b 0)')
    print('max = ' + str(np.max(gaussian_blob)) + ' (s/b 255)')
    print(gaussian_blob)

    return gaussian_blob
# end function

def add_gaussian_blob_to_heatmap(gaussian_blob, blob_center_x, blob_center_y, heatmap):
    # ToDo: this function is not perfect, there is a slight seam when two blobs overlap each other, eventually should resolve this

    blob_height, blob_width = gaussian_blob.shape[0:2]
    blob_left_edge_loc = round(blob_center_x - ((blob_width - 1) * 0.5))
    blob_right_edge_loc = round(blob_center_x + ((blob_width - 1) * 0.5))

    print('\n' + 'blob_left_edge_loc = ' + str(blob_left_edge_loc))
    print('\n' + 'blob_right_edge_loc = ' + str(blob_right_edge_loc))

    blob_top_edge_loc = round(blob_center_y - ((blob_height - 1) * 0.5))
    blob_bottom_edge_loc = round(blob_center_y + ((blob_height - 1) * 0.5))

    print('\n' + 'blob_top_edge_loc = ' + str(blob_top_edge_loc))
    print('\n' + 'blob_bottom_edge_loc = ' + str(blob_bottom_edge_loc))

    heatmap = heatmap.astype(np.uint16)
    gaussian_blob = gaussian_blob.astype(np.uint16)

    heatmap[blob_top_edge_loc:blob_bottom_edge_loc+1, blob_left_edge_loc:blob_right_edge_loc+1] += gaussian_blob

    heatmap = np.where(heatmap > 255, 255, heatmap).astype(np.uint8)

    return heatmap
# end function

if __name__ == '__main__':
    main()

result:

The seam will go away if you create them larger than 2.5*sigma. You are cutting off the Gaussian too early, creating a sharp edge. Also, all that clipping to uint8 range, and casting back and fourth between 8 and 16 bits… I would just create the image in a floating-point array, and quantize and clip only once at the very end. — Cris Luengo, Jul 20 '23 at 01:45

How to draw 2D Gaussian blob on an OpenCV image?

2 Answers2