Pinch/bulge distortion using Python OpenCV

I want to apply a pinch/bulge filter on an image using Python OpenCV. The result should be some kind of this example:

https://pixijs.io/pixi-filters/tools/screenshots/dist/bulge-pinch.gif

I've read the following stackoverflow post that should be the correct formula for the filter: Formulas for Barrel/Pincushion distortion

But I'm struggling to implement this in Python OpenCV.

I've read about maps to apply filter on an image: Distortion effect using OpenCv-python

As for my understanding, the code could look something like this:

import numpy as np
import cv2 as cv

f_img = 'example.jpg'
im_cv = cv.imread(f_img)

# grab the dimensions of the image
(h, w, _) = im_cv.shape

# set up the x and y maps as float32
flex_x = np.zeros((h, w), np.float32)
flex_y = np.zeros((h, w), np.float32)

# create map with the barrel pincushion distortion formula
for y in range(h):
    for x in range(w):
        flex_x[y, x] = APPLY FORMULA TO X
        flex_y[y, x] = APPLY FORMULA TO Y

# do the remap  this is where the magic happens
dst = cv.remap(im_cv, flex_x, flex_y, cv.INTER_LINEAR)

cv.imshow('src', im_cv)
cv.imshow('dst', dst)

cv.waitKey(0)
cv.destroyAllWindows()

Is this the correct way to achieve the distortion presented in the example image? Any help regarding useful ressources or preferably examples are much appreciated.

After familiarizing myself with the ImageMagick source code, I've found a way to apply the formula for distortion. With the help of the OpenCV remap function, this is a way to distort an image:

import numpy as np
import cv2 as cv

f_img = 'example.jpg'
im_cv = cv.imread(f_img)

# grab the dimensions of the image
(h, w, _) = im_cv.shape

# set up the x and y maps as float32
flex_x = np.zeros((h, w), np.float32)
flex_y = np.zeros((h, w), np.float32)

# create map with the barrel pincushion distortion formula
for y in range(h):
    delta_y = scale_y * (y - center_y)
    for x in range(w):
        # determine if pixel is within an ellipse
        delta_x = scale_x * (x - center_x)
        distance = delta_x * delta_x + delta_y * delta_y
        if distance >= (radius * radius):
            flex_x[y, x] = x
            flex_y[y, x] = y
        else:
            factor = 1.0
            if distance > 0.0:
                factor = math.pow(math.sin(math.pi * math.sqrt(distance) / radius / 2), -amount)
            flex_x[y, x] = factor * delta_x / scale_x + center_x
            flex_y[y, x] = factor * delta_y / scale_y + center_y

# do the remap  this is where the magic happens
dst = cv.remap(im_cv, flex_x, flex_y, cv.INTER_LINEAR)

cv.imshow('src', im_cv)
cv.imshow('dst', dst)

cv.waitKey(0)
cv.destroyAllWindows()

This has the same effect as using the convert -implode function from ImageMagick.