Coordinate system in digital image processing
Different operating systems, standards, and graphical softwares implement distinct (but quite related) conventions for representing, storing, manipulating, or displaying two-dimensional data on computers. The distinction is about the orientation of the axes.
Basically we have 4 main object types :

$f(x,y)$ , x-y continuous, is a function of mathematical relevance, such as an analog image.
$f[n,m]$, n-m integers, is a sequence obtained by sampling $f(x,y)$, or  synth-generated.
$F(i,j)~$, i-j integers, is a matrix to represent the data in its rows-i and columns-j.
$f[i,j]~~$ , i-j integers, is an array to store the data used in a computer program.

The arguments $x,y,n,m,i,j$ can also be replaced with $x_1,x_2,n_1,n_2,i_1,i_2$, the subscripts indicating the dimension order of the variables, such as the first or the second dimension.
A sequence $f[n_1,n_2]$ can represent a 2D digital image, or a mathematical function defined over the coordinate system of axes $n_1,n_2$. Its samples are placed in the correct elements of the associated matrix $f(i,j)$ for the operation to produce the desired effect at the output.
MATLAB fundamentally uses MATRIX based dimension ordering. The first sample A(1,1) is at the top-left corner, first dimenison $i$ is along vertical-down the rows, and the second dimension, $j$, is along the columns (horizontal-right).
This is equivalent to placing a coordinate system of $n_1 ~, ~n_2$ with its origin (0,0) at the element A(1,1); the first axis $n_1$ points along the first dimension $i$ (vertical down); and the second axis $n_2$ points along the second dimension $j$ (horizontal right) of the matrix A(i,j).
This mapping is obtained when you rotate the conventional coordinate system ($n_1$ horizontal-right, and $n_2$ vertical-up, origin at bottom-left) by 90 degrees clockwise with respect to its origin (0,0) and placing the origin at the element A(1,1) of the matrix.
In a parallel fashion, the 2D-DFT (discrete Fourier transform) $F[k_1,k_2]$ of the sequence $f[n_1,n_2]$ is computed by the statement:  F = fft2( f ) which places the first dimension, $k_1$, along the rows $i$, and the second dimension, $k_2$, along the columns $j$ of the output matrix $F(i,j)$ to represent the DFT sequence $F[k_1,k_2]$. Which is also alligned with the first and second dimenions of the input matrix $f(i,j)$ that represents the seqeunce $f[n_1,n_2]$.
Below is an oldskool discussion of a few mapping modes, and functions in MATLAB related with data orientation. As long as mappings are used consistently, they all yield the same results when interpreted correctly.
% SEQUENCES, BMP IMAGES, MATRICES and MATLAB FUNCTIONS :
% ------------------------------------------------------
%
% All data processing and display in Matlab is done via MATRICES A(i,j).
% But, theory of image processing is based on SEQUENCES f[n1,n2], F[k1,k2]. 
% When processing images, orientation of axes become relevant across
% functions such as CONV2(), FFT2(), IMSHOW(), STEM3(), SURF().
% And a MAPPING convention from f[n1,n2] into A(i,j) should be used.
%
%
% A matrix A(i,j) is indexed by vert rows i, and horz columns j.
% A sequence f[n1,n2] is indexed by horz-right n1, and vert-up n2.
% Which is the most typical, and natural, orientations for n1 and n2.
%
% 
% The mapping convention depends on the functions being called :
% 
%           1- PROCESS functions: conv2(), fft2(), filter2()
%           2- DISPLAY functions: imshow(), surf(), stem3()
%
%
% We consider following mapping modes between the samples of a sequence 
% f[n1,n2] and elements of a matrix A(i,j) as follows:
%
%
% MM-0 : 90d CW ROTATED MAPPING :
% ----------------------------------------------------------------------
%  f[0,0] --&gt; A(1,1), and &quot;n1&quot; grows DOWN from the top-row of A.
%
%    n1 = i-1   , n2 = j-1     ===&gt;          A(i,j) = f[i-1, j-1]
%    i  = n1+1  , j  = n2+1    ===&gt;        f[n1,n2] = A(n1+1,n2+1)
% 
%  ---o---1--------2-----&gt; j (n2)
%   1 | f[0,0]  f[1,0]
%   2 | f[1,0]  f[1,1]    Ex mapping of f[n1,n2] into 3x2 matrix A(i,j)
%   3 | f[2,0]  f[2,1]
%     |
%   i v
%  n1 
%
%
%
% MM-1 : FLIP-DOWN MAPPING :
% ----------------------------------------------------------------------
% f[0,0] ---&gt; A(1,1) and &quot;n2&quot; grows DOWN from the top-row of A.
%
%    n1 = j-1   , n2 = i-1     ===&gt;          A(i,j) = f[j-1, i-1]
%    i  = n2+1  , j  = n1+1    ===&gt;        f[n1,n2] = A(n2+1,n1+1)
% 
%  ---o---1--------2-----3---&gt; j (n1)
%   1 | f[0,0]  f[1,0]  f[2,0]
%   2 | f[0,1]  f[1,1]  f[2,1] Ex mapping of f[n1,n2] into 3x2 matrix A(i,j)
%     |
%   i v
%  n2  
%
%
% MM-2 : BMP MAPPING :
% ----------------------------------------------------------------------
% f[0,0] is stored in A(N2,1) and &quot;n2&quot; grows UP from the bottom-row of A
%
%    n1 = j-1   , n2 = N2-i     ===&gt;          A(i,j) = f[j-1, N2-i]
%    i  = N2-n2 , j  = n1+1     ===&gt;        f[n1,n2] = A(N2-n2,n1+1)
% 
%  n2 ^
%   2 | A(1,1)  A(1,2)         An example mapping into 3x2 matrix A
%   1 | A(2,1)  A(2,2)
%   0 | A(3,1)  A(3,2)
%  ---o---0--------1-----&gt; n1 (j)
%
%
%
% Axis orientations of functions CONV(),FILTER(),FFT2():
% ----------------------------------------------------------------------
% They assume 90d ROTATED / or FLIP DOWN index mapping by default. 
%
%      A(1,1)---------&gt; n2,k2,j (horizontal - column variable)
%           |
%           |  A(i,j)
%           |
%           v n1,k1,i (vertical - row variable)
%
%                  
% NOTE: Strictly speaking, MATLAB does not care which mapping was used to
% generate the matrix A(i,j). Rather it treats &quot;i&quot; as the first dimension
% and &quot;j&quot; as the second dimension, and outputs acordingly.
%
%  
% Axis orientation of DISPLAY functions STEM3(),  SURF()            :
% -------------------------------------------------------------------
% It produces the plots according to NATURAL X-Y orientation in which
% n1 point horizontal right, and n2 point verticcal up and (0,0) being
% at the bottom. But this requires that the MATRIX A(i,j) was filled in    
% according to FLIP-DOWN mapping mode.
%
%                          A(N2,N1)
%                          f22
%                       f12   f21 
%     A(N2,1) (n2-i) f02  f11  f20 (n1-j) A(1,N1)
%                      f01  f10      
%            n2-axis      f00       n1-axis
%                        A(1,1) 
%
%
% Image display function IMSOW(A(i,j)) assumes a BMP mapping:
% -------------------------------- --------------------------
% IMSHOW displays the matrix contents in its row-column order.
% Therefore if you want to display a sequence f[n1,n2], or F[k1,k2] using
% the IMSHOW() function, then in order to have the correct orientation
% according to natural x-y coordinates, you should use BMP mapping mode  
% on the matrix which is to be displayed.
%

SUMMARY
If IMSHOW() iwill be used to display images, or their FFT results, then correct alignment with $n_1,n_2$ and $k_1,k_2$ requires BMP based mapping to be used to fill in the associated matrix. If STEM3() or SURF() will be used to get a 3D display of the sequences, or FFT results, then flipp-down mapping produces correct orientation. Other functions CONV2(), FFT2() etc., work equally well with either Rotated or Flip-Down mapping modes.