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

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?

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()