The fundamental result in sampling theory states that if a signal $f(t)$ contains no frequencies higher than $\omega$ cycles per second, then $f(t)$ is completely determined by its values $f(\frac{k}{2\omega})$ at a discrete sequence of sample points with spacing $\frac{1}{2\omega}$ and can be reconstructed from these values by the formula$\dagger$
$$\label{1}\tag{1} f(t)=\sum_{k=-\infty}^{\infty}f(\frac{k}{2\omega})\frac{\sin(\pi(2\omega t-k))}{\pi(2\omega t -k)}$$
Is something like $(\ref{1})$ still valid if instead of the signal $f(t)\colon t \in \mathbb{R} \to f(t) \in \mathbb{R}$ we have signal defined in higher dimension? For example the signal
$f(x,y)\colon x,y \in \mathbb{R} \to f(x,y) \in \mathbb{R}$
or
$f(x,y,z)\colon x,y,z \in \mathbb{R} \to f(x,y,z) \in \mathbb{R}$
$\dagger$, Literally taken from the incipit of: CIAURRI, Óscar; VARONA, Juan. A Whittaker-Shannon-Kotel’nikov sampling theorem related to the Dunkl transform. Proceedings of the American Mathematical Society, 2007, 135.9: 2939-2947.
After a bit of search I have found the answer:
THEOREM 1.6$\dagger$
Let $f(t_1 , t_2 ,\ldots , t_n )$ be a function of $n$ real variables, whose $n$-dimensional Fourier integral, $g$, exists and is identically zero outside an $n$-dimensional rectangle and is symmetrical about the origin, that is,
$g(y_1 , y_2 , \ldots , y_n )=0$, $\mathopen|y_k\mathclose|\gt\mathopen|\omega_k\mathclose|$, $k=1,2,\ldots,n$
Then
$$f(t_1 , t_2 ,\ldots , t_n )=f(t)=\sum_{m_1=-\infty}^{\infty}\ldots \sum_{m_n=-\infty}^{\infty}f(\frac{\pi m_1}{\omega_1},\ldots,\frac{\pi m_n}{\omega_n})\times \frac{\sin(\omega_1 t_1-m_1\pi)}{\omega_1 t_1-m_1\pi}\times\ldots\times\frac{\sin(\omega_n t_n-m_n\pi)}{\omega_n t_n-m_n\pi}$$
$\dagger$ POULARIKAS, Alexander D. Transforms and applications handbook. CRC press, 2010.