Surface Brightness Partly Black After Convolution

Question

I am convolving the surface brightness of galaxies (one unsheared and one sheared) with a gaussian kernel to simulate the effects of atmospheric and instrumental aberrations.

However, after performing the convolution, there is a black strip on the outer rim of the resulting graph, and this should not happen.

The size of the original array is 64 by 64, the size of the kernel is 5 by 5, the padding is 2, and the stride is 1.

According to the size formula: [(W−K+2P)/S]+1, the output array should have the same dimensions as the original image, so I don't think I got the padding wrong.
However, I can't think of any other reason this phenomenon should happen.

Here's the code:

def generateGalaxy(size, r_disk, f_brightness_ratio, r_bulge):
    x = np.arange(-size/2, size/2)
    y = np.arange(-size/2, size/2)
    X, Y = np.meshgrid(x, y)
    r = (X**2 + Y**2)**0.5
    values = np.exp(-r/r_disk) + f_brightness_ratio*np.exp(-(r/r_bulge)**0.25)
    return values, X, Y

def convolve2D(image, kernel, padding = 2, strides =1):
    
    kernel = np.flipud (np.fliplr(kernel))
    xKernShape = kernel.shape[0]
    yKernShape = kernel.shape[1]
    xImgShape = image.shape[0]
    yImgShape = image.shape[1]
    
    xOutput = int(((xImgShape-xKernShape + 2* padding)/strides)+1)
    yOutput = int(((yImgShape-yKernShape + 2* padding)/strides)+1)
    
    output = np.zeros((xOutput, yOutput))

    if padding != 0:
        imagePadded = np.zeros((image.shape[0] + padding *2, image.shape[1]+padding*2))
        imagePadded[int(padding):int(-1*padding), int(padding):int(-1*padding)] = image
    else:
        imagePadded = image
        
    for y in range(image.shape[1]):
        # y value goes from 0 to image.shape[1] from top to down
        if y > image.shape[1] - yKernShape:
            break
        if y % strides == 0:
            for x in range(image.shape[0]):
                if x > image.shape[0] - xKernShape:
                    break
                try:
                    if x % strides == 0:
                        output [x,y] = (kernel * imagePadded [x:x+xKernShape, y:y+yKernShape]).sum()
                except:
                    break         
    return output

def gkern(sidelength, sigma):

    """
    creates gaussian kernel with side length `l` and a sigma of `sig`
    """
    ax = np.linspace(-(sidelength - 1) / 2., (sidelength - 1) / 2., sidelength)
    gauss = np.exp(-0.5 * np.square(ax) / np.square(sigma))
    kernel = np.outer(gauss, gauss)
    return kernel / np.sum (kernel)

I discovered that if I change the padding to 0, then the black rim will be gone. However, the output size will be reduced. Is there anyway to remove the strip while keep the output size the same as the input size? Also, why does this black strip exist in the first place?