The difference between using fft2 in matlab for odd and even matrices

I want to perform Fourier transform on a two-dimensional center-symmetric disk in matlab. The specific code is as follows.

N = 20;

matrix = zeros(N, N);

diameter = 20; 
radius = diameter / 2;

center_x = (N + 1) / 2;
center_y = (N + 1) / 2;

[X, Y] = meshgrid(1:N, 1:N);

distances = sqrt((X - center_x).^2 + (Y - center_y).^2);

matrix(distances <= radius) = 1;

figure;
imshow(matrix);

figure;
matrix_ft = fftshift(fft2(ifftshift(matrix))); 
imshow(matrix_ft);

But for N taking odd and even numbers, different results will be obtained. when N=20: when N=21:

From a theoretical point of view, the result should be a first-order Bessel function of the first kind—a series of ring structures. It should be centrally symmetrical.

When N is an odd number, after Fourier transformation, the imaginary parts of the complex elements in the matrix_ft are all 0, and a series of centrally symmetric ring structures are obtained, which is the correct result. When N is equal to an even number, after Fourier transformation, the imaginary part of the complex elements in the matrix_ft will not be 0. Although it also maintains central symmetry (the center of symmetry is offset from the center of the matrix), it does not present a perfect ring structure. The outer ring structure is obviously broken, which is obviously different from the correct result.

I know that with fft, an odd number of samples introduces a DC component. Even-numbered sampling not only enters the DC component, but also introduces the Nyquist frequency.My current thought is, can I get a correct result under the condition that N is an even number? And can it be similar to the matrix at this time and maintain a symmetrical structure(The center of symmetry is at (N/2+0.5=10.5, N/2+0.5=10.5))? When N=20, matrix is:

The center of symmetry has to be on the pixel floor(N/2) + 1 (which ifftshift shifts to the pixel at 1, in MATLAB indexing which starts at 1). If it’s in between pixels, like how you construct the even-sized circle, then it is a shifted disk, which causes the phase in the frequency domain to be modulated. Basically you’re shifting some of the weight of the FFT from the real component (which you display) to the imaginary component, at different rates for different frequencies.