Distort using ftheta equation

How do I apply ftheta distortion equation to an image without using cv2.fisheye?

I am trying to do this to get a more solid grasp of the foundations. But the code produces output like below. What am I doing wrong?

Input:

Output:

import cv2
import numpy as np

def ideal_f_theta_distortion(image):
    h, w = image.shape[:2]
    cx, cy = w // 2, h // 2  # Image center
    f = min(cx, cy)  # Focal length in pixels (approximate)

    # Create mapping arrays
    map_x = np.zeros((h, w), dtype=np.float32)
    map_y = np.zeros((h, w), dtype=np.float32)

    for y in range(h):
        for x in range(w):
            dx, dy = (x - cx), (cy - y) # doing cy - y because my center is top-left
            theta = np.arctan2(dy, dx)  # Angle in radians
            r = np.sqrt(dx**2 + dy**2)  # Original radius

            # Apply f-theta transformation
            r_new = f * theta  # Ideal f-theta projection

            # Convert back to Cartesian
            x_new = int(cx + r_new * np.cos(theta))
            y_new = int(cy - r_new * np.sin(theta)) # my center is top-left

            # Ensure valid mapping
            if 0 <= x_new < w and 0 <= y_new < h:
                map_x[y, x] = x_new
                map_y[y, x] = y_new
            else:
                map_x[y, x] = cx
                map_y[y, x] = cy

    # Apply remap
    distorted = cv2.remap(image, map_x, map_y, interpolation=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT)

    return distorted

# Simulate deformation field
N = 500
sh = (N, N)

# Test image
img = np.zeros(sh)
img[::10, :] = 255
img[:, ::10] = 255

corrected_image = ideal_f_theta_distortion(img)
Image.fromarray(np.array(corrected_image, dtype="uint8"))

This is the code that finally works.

thetashould be the angle between optic axis and the point on the image, so tan(theta) will be r/f assuming ris the distance of the 2D point (on the image) from the center of the image. Looks like I got the core concept wrong in the original post

img = cv2.imread(<impath>)[:,:,::-1]
H, W, _ = img.shape
cX, cY = W//2, H//2 #7, 5 
f = min(cX, cY)
mapX = np.zeros((H, W))
mapY = np.zeros((H, W))
for x in range(W):
    for y in range(H):
        dx, dy = (x - cX), (y - cY)
        r = np.sqrt(dx**2 + dy**2)
        phi = np.arctan2(dy, dx)
        theta = np.arctan(r/f)
        rnew = f*theta
        xnew = cX + rnew*np.cos(phi)
        ynew = cY + rnew*np.sin(phi)
        mapX[H - 1 - int(ynew), int(xnew)] = H - 1 - y
        mapY[H - 1 - int(ynew), int(xnew)] = x
distorted = cv2.remap(np.array(img, "uint8"), np.array(mapY, "float32"), np.array(mapX, "float32"), interpolation=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT)